What is Input for lsim?

10 views (last 30 days)
Jesse Crotts
Jesse Crotts on 10 Dec 2018
Commented: Jesse Crotts on 10 Dec 2018
When using the lsim function, the second variable to input is the "input" to the system. This seems very straightforward for a single dimensional system. However, when you have a multi-degree of freedom system. What is the input, inputing too?
1) For example, if I had two masses between three springs, then would the input be toward the first or second mass or to the container that they are attached to?
2) Is the input to the system, a force or a displacement?
3) Why do you need an input variable? In a code like below, you define your inputs in the BB and DD matrices.
AA = [A11 A12; A21 A22];
BB = [zeros(NN,1); M\Q];
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 = 0*(time+1)./(time+1);
P(1:timesize) = 500;
%P(2) = 500;
%P = zeros(timesize(1),1);
SYS = ss(AA,BB,CC,DD);
[YY, TT] = lsim(SYS, P, time);
% -------------------------------------------------------
4) Should I be defining my inputs in the BB and DD matrices or the "P" input or both?
  2 Comments
Star Strider
Star Strider on 10 Dec 2018
I can’t run your code.
The documentation states:
  • The input u is an array having as many rows as time samples (length(t)) and as many columns as system inputs.
Jesse Crotts
Jesse Crotts on 10 Dec 2018
Here is the full length of my code:
% 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);
DD = zeros(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)
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 = 0*(time+1)./(time+1);
P(1:timesize) = -10;
%P(2) = 500;
%P = zeros(timesize(1),1);
SYS = ss(AA,BB,CC,DD);
[YY, TT] = lsim(SYS, P, time);
% -------------------------------------------------------
plot(TT,YY(:,2));
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
% -------------------------------------------------------

Sign in to comment.

Answers (0)

Categories

Find more on Robust Control Toolbox 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!