I know this seems like a dumb question because you should know by your system size, but I can't seem seem to get my lsim function to work with an array of inputs.
My code seems to work fine with one input.
Here is my code:
The part where I define my inputs is found at:
P = zeros(timesize(1),16);P(:,7) = -10;I feel like my number of inputs should be 8 or 16 but niether seems to work.
% Lec 12 39 min
clear; clctime = (0:.01:22.5)';% -------------------------------------------------------
%{
% Project Demo (1) Requirements
l = 3; % length in inchesEI1 = 5*10^6; % EI lbf x in^2EI2 = 2*EI1; % EIEI3 = EI1; % EIEI4 = EI2; % EIk1 = -200; % spring constant 1 lbf per inchk2 = k1; % spring constant 2m = 268.3; % lbf per inchpo = 0; % distributed loadd1 = 2; % distance to Pd2 = 1; % distance to spring on Element 3d3 = 1.5; % distance to spring on Element 4P = 0*(time+1)./(time+1);P(1) = 500;P(2) = 500;end%}
% -------------------------------------------------------% -------------------------------------------------------%{% Project Demo (2) Requirements
l = 3; % length in inchesEI1 = 5*10^6; % EI lbf x in^2EI2 = 2*EI1; % EIEI3 = EI1; % EIEI4 = EI2; % EIk1 = -200; % spring constant 1 lbf per inchk2 = k1; % spring constant 2m = 268.3; % lbf per inchpo = 0; % distributed loadd1 = 2; % distance to Pd2 = 1; % distance to spring on Element 3d3 = 1.5; % distance to spring on Element 4P = -500*sin(2*pi*10*time); % point load lbfP = P';%}% -------------------------------------------------------% -------------------------------------------------------%{% Project Demo (3) Requirements
l = 3; % length in inchesEI1 = 5*10^6; % EI lbf x in^2EI2 = 2*EI1; % EIEI3 = EI1; % EIEI4 = EI2; % EIk1 = -200; % spring constant 1 lbf per inchk2 = k1; % spring constant 2m = 268.3; % lbf per inchpo = 100; % distributed loadd1 = 2; % distance to Pd2 = 1; % distance to spring on Element 3d3 = 1.5; % distance to spring on Element 4P = zeros(time,1); % point load lbf%}% -------------------------------------------------------% -------------------------------------------------------% Project Code Validation Requirements
l = 3; % length in inches
EI1 = 5*10^6; % EI lbf x in^2
EI2 = EI1; % EI
EI3 = EI1; % EIEI4 = EI1; % EIk1 = 0; % spring constant 1 lbf per inch
k2 = k1; % spring constant 2
m = 268.3; % lbf per inch
po = 0; % distributed load
d1 = 3; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = -10; % point load lbf
P = P';% -------------------------------------------------------% -------------------------------------------------------% Unit conversion
l = l/12;EI1 = EI1/144;EI2 = EI2/144;EI3 = EI3/144;EI4 = EI4/144;k1 = k1*12;k2 = k2*12;m = m*12;po = po*12;d1 = d1/12;d2 = d2/12;d3 = d3/12;P = P;% -------------------------------------------------------% -------------------------------------------------------% Phi function handles
phi1 = @(zeta) 1-3*zeta^2+2*zeta^3;phi2 = @(zeta) l*zeta-2*l*zeta^2+l*zeta^3;phi3 = @(zeta) 3*zeta^2-2*zeta^3;phi4 = @(zeta) -l*zeta^2+l*zeta^3;Phi = @(zeta)[... phi1(zeta)*phi1(zeta) phi1(zeta)*phi2(zeta) phi1(zeta)*phi3(zeta) phi1(zeta)*phi4(zeta);... phi2(zeta)*phi1(zeta) phi2(zeta)*phi2(zeta) phi2(zeta)*phi3(zeta) phi2(zeta)*phi4(zeta);... phi3(zeta)*phi1(zeta) phi3(zeta)*phi2(zeta) phi3(zeta)*phi3(zeta) phi3(zeta)*phi4(zeta);... phi4(zeta)*phi1(zeta) phi4(zeta)*phi2(zeta) phi4(zeta)*phi3(zeta) phi4(zeta)*phi4(zeta) ];%PPhi = Phi(zeta)
%}% This function needs to be moved to where it is needed
% -------------------------------------------------------% -------------------------------------------------------% Displacement function handles
Qpo = @(po,l) [po*l/2 po*l^2/12 po*l/2 -po*l^2/12]';QP = @(zeta,P) [-P*phi1(zeta) -P*phi2(zeta) -P*phi3(zeta) -P*phi4(zeta)]';Qspring3 = @(k1,zeta) k1*Phi(zeta);Qspring4 = @(k2,zeta) k2*Phi(zeta);% -------------------------------------------------------% -------------------------------------------------------% Calculation simplification
a=zeros(4,6);aa=zeros(4);aaa=zeros(4,2);I4=eye(4);A1=[I4 a];A2=[aaa I4 aa];A3=[aa I4 aaa];A4=[a I4];AA11 = zeros(10);AA12 = eye(10);% -------------------------------------------------------% -------------------------------------------------------% load transformations
% p
% P
% -------------------------------------------------------% -------------------------------------------------------% Mass matrix, Stiffness matrix
M1 = (m*l/420)*[... 156 22*l 54 -13*l; 22*l 4*l^2 13*l -3*l^2 54 13*l 156 -22*l -13*l -3*l^2 -22*l 4*l^2];M2 = M1; % note, this can be different
M3 = M1; % note, this can be differentM4 = M1; % note, this can be differentK1 = (EI1/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];K2 = (EI2/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];K3 = (EI3/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];K4 = (EI4/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];% -------------------------------------------------------% calculations
K3 = K3 + Qspring3(k1,d2/l); % Khat
K4 = K4 + Qspring4(k1,d3/l); % KhatQpo4 = Qpo(po,l); % distributed loadQP4 = QP(d1/l,P); % Point load
Q(1) = 0;Q(2) = 0;Q(3) = 0;Q(4) = 0;Q(5) = 0;Q(6) = 0;Q(7) = QP4(1)+Qpo4(1);Q(8) = QP4(2)+Qpo4(2);Q(9) = QP4(3)+Qpo4(3);Q(10) = QP4(4)+Qpo4(4);Q = Q';% -------------------------------------------------------% -------------------------------------------------------% global calculations
M = A1'*M1*A1+A2'*M2*A2+A3'*M3*A3+A4'*M4*A4;K = A1'*K1*A1+A2'*K2*A2+A3'*K3*A3+A4'*K4*A4;% -------------------------------------------------------% -------------------------------------------------------% Value Saves
Ktrue = K;Mtrue = M;Qtrue = Q;% -------------------------------------------------------% -------------------------------------------------------% Finds F1,M1
NN = 8;%C = zeros(NN); % no damping
C = 200*eye(NN); % Damping, better to find steady state
A11 = zeros(NN,NN);A12 = eye(NN);M(1,:) = [];M(1,:) = [];M(:,1) = [];M(:,1) = [];K(1,:) = [];K(1,:) = [];K(:,1) = [];K(:,1) = [];Q(1) = [];Q(1) = [];A21 = -inv(M)*K;A22 = -inv(M)*C;AA = [A11 A12; A21 A22];BB = [zeros(NN,1); M\Q];%CC = [zeros(1,2*NN)];
%CC(7) = 1; % decides which value to pull
% For CC, first NN are position, 2nd NN are acceleration
CC = [eye(NN), zeros(NN)];DD = [0];% -------------------------------------------------------% -------------------------------------------------------% Short Simulation
time = (0:.01:22.5)';%P = 0*(time+1)./(time+1);
%P(1) = 500;
%P(2) = 500;
timesize = size(time);P = zeros(timesize(1),16);P(:,7) = -10;%P(2) = 500;%P = zeros(timesize(1),1);
SYS = ss(AA,BB,CC,DD);[YY, TT] = lsim(SYS, P, time);% -------------------------------------------------------plot(TT,YY(:,7));xlabel('Time(sec)');ylabel('Displacement');title('Stepped Response(Underdamped)');% -------------------------------------------------------% Resets True Values
K = Ktrue;M = Mtrue;Q = Qtrue;% -------------------------------------------------------% -------------------------------------------------------% Sets real Q values
%Q(1) = dot(K(1,:),x1(1000,7:16));
%Q(2) = dot(K(2,:),x1(1000,7:16));
% -------------------------------------------------------% -------------------------------------------------------% After F1, M1 are Found
NN = 10;C = zeros(NN);A11 = zeros(NN,NN);A12 = eye(NN);A21 = -inv(M)*K;A22 = -inv(M)*C;AA = [A11 A12; A21 A22];BB = [zeros(NN,1); M\Q];CC = [zeros(1,2*NN)]; % decides which value to pull
% For CC, first NN are position, 2nd NN are acceleration% -------------------------------------------------------
Best Answer