How do you know how many inputs your lsim function has?
2 views (last 30 days)
Show older comments
Jesse Crotts
on 10 Dec 2018
Commented: Jesse Crotts
on 10 Dec 2018
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; clc
time = (0:.01:22.5)';
% -------------------------------------------------------
%{
% Project Demo (1) Requirements
l = 3; % length in inches
EI1 = 5*10^6; % EI lbf x in^2
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = -200; % spring constant 1 lbf per inch
k2 = k1; % spring constant 2
m = 268.3; % lbf per inch
po = 0; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = 0*(time+1)./(time+1);
P(1) = 500;
P(2) = 500;
end
%}
% -------------------------------------------------------
% -------------------------------------------------------
%{
% Project Demo (2) Requirements
l = 3; % length in inches
EI1 = 5*10^6; % EI lbf x in^2
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = -200; % spring constant 1 lbf per inch
k2 = k1; % spring constant 2
m = 268.3; % lbf per inch
po = 0; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = -500*sin(2*pi*10*time); % point load lbf
P = P';
%}
% -------------------------------------------------------
% -------------------------------------------------------
%{
% Project Demo (3) Requirements
l = 3; % length in inches
EI1 = 5*10^6; % EI lbf x in^2
EI2 = 2*EI1; % EI
EI3 = EI1; % EI
EI4 = EI2; % EI
k1 = -200; % spring constant 1 lbf per inch
k2 = k1; % spring constant 2
m = 268.3; % lbf per inch
po = 100; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = 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; % EI
EI4 = EI1; % EI
k1 = 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 different
M4 = M1; % note, this can be different
K1 = (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); % Khat
Qpo4 = Qpo(po,l); % distributed load
QP4 = 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
% -------------------------------------------------------
0 Comments
Accepted Answer
Mark Sherstan
on 10 Dec 2018
Run:
size(SYS)
and you get:
>> State-space model with 8 outputs, 1 inputs, and 16 states.
According to lsim documentation: The input u is an array having as many rows as time samples (length(t)) and as many columns as system inputs. As you your time vector is 2251 long your input u (or P in your case) should be a 2551x1 matrix.
More Answers (0)
See Also
Categories
Find more on General Applications in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!