MATLAB: Strange behaviour of step response (unstable system)

controlControl System ToolboxMATLABsimulinktheoryunstanble step response

Hi!
The system defned in the following code should be stable (verified also in simulink). The step response looks correct in the beginning 25sec) but becomes unstable for simulation up to 50sec. Whats happening? Nummerical problems?
s = tf('s');
ks=2;
a=1;
b=2;
 
G=ks/s/(s-a)/(s+b);
 
kp=0.4;
Tv=2*a*b/(b-a)
 
Gr=kp*(1+Tv*s)
 
Gcl=Gr*G/(1+Gr*G);
t=0:0.0001:50;
figure;
step(Gcl,t);

Best Answer

When you do transfer function math like this, matlab doesn't attempt to deal with poles and zeros that should cancel. Then you can run into problems when the pole/zero cancellations that should be exact but aren't because of rounding errors. You can see from your expression for Gcl that Gcl should have a third order denominator. But Gcl as computed is:
>> Gcl
Gcl =
 
     3.2 s^4 + 4 s^3 - 5.6 s^2 - 1.6 s
  ---------------------------------------
  s^6 + 2 s^5 + 0.2 s^4 - 1.6 s^2 - 1.6 s
 
Continuous-time transfer function.
Now look at the poles and zeros:
>> [p,z]=pzmap((Gcl))
p =
   0.0000 + 0.0000i
  -2.0000 + 0.0000i
   1.0000 + 0.0000i
  -0.1107 + 1.0076i
  -0.1107 - 1.0076i
  -0.7785 + 0.0000i
z =
         0
   -2.0000
    1.0000
   -0.2500
Note that the first three poles seem to be equivalent to the first three zeros. Mathematically, they should be. But are they numerically:
>> [z(1:3)-p(1:3)].'
ans =
   1.0e-15 *
         0         0   -0.8882
So we see that third pole at s=1 doesn't perfectly cancel with the zero. As a result, your step response will eventually go unstable. If you're confident that the those three poles and zeros should cancel based on your problem, you can do:
>> minreal(Gcl)
ans =
 
        3.2 s + 0.8
  -----------------------
  s^3 + s^2 + 1.2 s + 0.8
 
Continuous-time transfer function.
Which is probably the answer you're looking for. Or, you could have gotten that answer directly using:
>> feedback(Gr*G,1)
ans =
 
        3.2 s + 0.8
  -----------------------
  s^3 + s^2 + 1.2 s + 0.8
 
Continuous-time transfer function.
Related Question