How do you know how many inputs your lsim function has?

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
% -------------------------------------------------------

Answer 1

[Mark Sherstan]

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.

[Jesse Crotts]

Thanks Mark,

This helped me to figure out how to debug to get to be able to put in multiple inputs.

What determines the number of inputs to the lsim function is the number of column in the BB matrix of sys = ss(AA,BB,CC,DD).

So before I had a column vector for BB. The fix was to diagonalize this matrix to give me a 16x8 matrix instead of a 16x1 vector.