why are the output simulation results of the following LTI system different by using 'lsim' command and theorical method calculations?

Dear all,

I would like to obtain the output response of an LTI state space model to an exponential input compatible with system invariant zero frequency and adjusted amplitude to get a zero oputput in all the time (input signal : u(t)=-2*exp(-3*t); initial state vector =x0=[1;-1;0];)

so I wrote the following code by two different ways :

1-by using 'lsim 'command

2- by using output response theoretical formulation.

unfortunatley ,the system output results weren't compatible with each other. what is the reason for this difference?is there any limitation on using "lsim "command for LTI systems?

Thank you for your consideration.

clc

clear all

A=[1 4 0;0 -1 0;0 2 -3]

B=[0;-1;-1]

C=[-1 -1 0];

D=0;

sys1=ss(A,B,C,D);

G1=tf(sys1)

time=0:0.1:10;

for i=1:length(time)

u(i)=-2*exp(-3*time(i));

end

x0=[1;-1;0]; %% initial state vector

% % % % % % % % % % % % % % % %simulation 1

y=lsim(sys1,u',time,x0);

figure(1)

plot(time,y)

xlabel('time')

ylabel('output')

title('output response to exponetial inout(simulation 1)')

grid on

% % % % % % % % % % % % % % % %simulation 2

syms t to

expm(A*t)

y3=(C*expm(A*t)*(x0))+(C*expm(A*t)*(int(expm(-A*to)*B*-2*exp(-3*to),0,t)))

for i=1:length(time)

y3n(i)=subs(y3,t,time(i));

end

y3t=double(y3n)

figure(2)

plot(time,y3t,'r');

xlabel('time')

ylabel('output')

title('output response to exponetial inout(simulation 2)')

grid on

Answer 1

[Paul]

I suspect the problem is that lsim() is ineherently an approximation to the LTI system response. I think that lsim() generates a discrete time approximation, and so its accuracy will be a function of the time step. And of course there will always be numerical inaccuracy. I'm pretty sure, but not certain, that lsim computes the output due to the input and the initial conditions separately, and then adds them together at the end. For this example, it would be subtracting two very large numbers (because the system is unstable) in the hopes of getting zero as the net output.

Let's look at the code:

A=[1 4 0;0 -1 0;0 2 -3];
B=[0;-1;-1];
C=[-1 -1 0];
D=0;
sys1=ss(A,B,C,D);
G1=tf(sys1);
time = 0:0.1:10;
u = -2*exp(-3*time);
x0 = [1;-1;0]; %% initial state vector
% % % % % % % % % % % % % % % %simulation 1
y = lsim(sys1,u,time,x0);
figure
plot(time,y)
xlabel('time')
ylabel('output')
title('output response to exponential inout(simulation 1)')
grid on

Now try it with a much smaller time step

time = 0:0.0001:10;
u = -2*exp(-3*time);
x0 = [1;-1;0]; %% initial state vector
% % % % % % % % % % % % % % % %simulation 1
y = lsim(sys1,u,time,x0);
figure
plot(time,y)
xlabel('time')
ylabel('output')
grid on

Although the shape of the curve is the same, note the scale on the y-axis where the order of magnitude is much lower.

Now the exact answer:

% % % % % % % % % % % % % % % %simulation 2
syms t to
expm(A*t);
y3(t) =(C*expm(A*t)*(x0))+(C*expm(A*t)*(int(expm(-A*to)*B*-2*exp(-3*to),0,t)));
y3(t)

ans = 

which can be simplified to

y3(t) = simplify(y3(t))

y3(t) = 

0

which is the expected result. The plot with the un-simplified version of y3(t) was jsut suffering from numerical inaccuracy.

[sara elmi]

Dear Paul,

I really appreciate that you took the time to answer my question. It was comprehensive and exactly what I needed to know. Thank you so much.