How do you know how many inputs your lsim function has?
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Mark Sherstan
2018-12-10
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 System Identification Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!