Lsim gives NaN for square input

Hi,

I want to simulate 6 state, 6 output, 1input linear system. But lsim gives NaN. Here is my code. What can I do to fix this problem? Thanks in advice.

m=3255;

ms =2956;

a =1.459;

b= 1.895;

h =1.234;

hsr= 0.781;

Cf=-120000;

Cr= -120000;

Tr =1.615;

Ixx= 1830;

Iyy= 6488;

Izz= 7913;

Ixz= 500;

D= 4500;

K= 145330;

f= 0.7;

r= 0.23;

g =9.81;

u=13.4;

M=[m 0 ms*hsr 0 0 0;

ms*hsr -Ixz Ixx+ms*hsr*hsr 0 0 0;

0 Izz 0 -Ixz 0 0;

0 0 0 1 0 0;

0 0 0 0 1 0;

0 0 0 0 0 1];

N=[0 m*u 0 0 -1 -1;

0 ms*hsr*u D K-ms*hsr*g 0 0;

0 0 0 0 -a -b;

0 0 -1 0 0 0;

-Cf/f -a*Cf/f 0 0 u/f 0;

-Cr/r b*Cr/r 0 0 0 u/r];

F=[0 0 0 0 Cf*u/-f 0].';

A=-M*N;

B=-M*F;

C=diag(ones(6,1));

sys=ss(A,B,C,0);

k=0:0.1:10;

i=1;

for t=0:0.1:10

Sf(i)=0;

if t>=1

if t<=3

Sf(i)=0.1;

end

end

i=i+1;

end

y=lsim(sys,Sf,k);

Answer 1

[Paul]

Look at the eigenvalues of the plant matrix:

>> eig(sys)

ans =

    -6.846331220081878e+07
    -5.237359879589575e+07
     5.902674140623139e+04
    -5.905163677506245e+04
    -2.946129710311468e+01
    -2.726234629459069e+01

To simulate this model, the step size must be orders of magnitude smaller than the value of 0.1 specified in k. But that won't really help because the system is very unstable. Are the equations correct? If the system were stable, there are model reduction techniques that could be used to eliminate the fast dynamics in the plant and approximate its response to the pulse defined in Sf.

[tokyo]

First of all, thank you for answer. I reduced the step size of k and it gives me a graphics like this. This is one of the response graphic and the others are similar to this one. I could not figure it out. There is no value after 1 sec. But the simulation should be 10 seconds. Also I know this is unstable system. Can Lsim not simulate unstable system?

[Paul]

Yes, lsim can simulate an unstable system. However, for just about any input the output of an ustable system eventually becomes too large to represent. But, to show that lsim can simulate an ustable system and what happens if you let it go too long, let's consider a simple example. First, I'm going to define the input Sf as a function of time, because we'll need that a lot:

>> Sf = @(t) 0.1*(t>= 1 & t <= 3);

Now, look at what happens using lsim to simulate a simple, unstable system, in response to Sf:

>> t=0:.001:10;
>> lsim(tf(1,conv([1 1],[1 -1])),Sf(t),t);

As you'll see, you get what looks like exponential growth, which is expected for an unstable system with a real pole in the RHP. But exponentials eventually grow quickly. Try it again with

>> t=0:.001:100;

If you extend the time vector out long enough, eventually the output will be become inf, i.e., too large to represent in floating point. But, over a reasonable time scale, lsim can simulate an unstable system just fine, where reasonable is relative to the growth rate due to the unstable pole(s).

Now, let's look at your system and focus only on the first output:

>> zpk(sys(1,1))
ans =
 
          -1.5215e10 (s-7.28e07) (s+3.053e07) (s+65.03) (s+27.26)
  -----------------------------------------------------------------------
  (s+6.846e07) (s+5.237e07) (s+5.905e04) (s-5.903e04) (s+29.46) (s+27.26)

We see a RHP pole at s = 5.903e4 ! That exponential gets excited by the step in Sf(t) at t = 1 and it only takes a fraction of a second before exp(5.903e4*(t - 1)) becomes way too large, which is what the plot shows (assuming it was generated with an extremely small dt in the time vector).