State Space Model with PID controller unstable in Julia compared to Matlab

I just realized that we can actually perform the multiplication by s in ControlSystems without using DescriptorSystems.jl. Multiplying by s is the same as dividing by 1/s, and division of two statespace systems can be performed even if the denominator is strictly proper (inverse is non-proper) as long as the quotient is proper:

julia> using ControlSystemsBase

julia> P = DemoSystems.double_mass_model(); # An example system with poles in the origin

julia> C = ss(1/tf('s')) # The inverse of a pure derivative term
StateSpace{Continuous, Float64}
A = 
 0.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

julia> P/C # This is now P*s
StateSpace{Continuous, Float64}
A = 
  48.499999999999474   37.123106012292105     79.21505538721604
 -70.00357133746986   -51.99999999999945    -111.35618527949073
   0.0                  1.1180339887498683    -0.5000000000000259
B = 
  0.7071067811865231
 -1.5000000000000109
  1.1180339887498945
C = 
 0.0  0.0  0.8944271909999159
D = 
 0.0

Continuous-time state-space model

The fraction \dfrac{P}{C} above is now equal to P(s)s will all computations perfomed in statespace form.

I played a bit with the code in Matlab and I can only confirm @baggepinnen 's explanation. The PID controller in this case is formed solely by the D term, for which the transfer function is just 278.7369 s and no state equations can be formulated to model it. Instead, differential-algebraic (also called descriptor) equations can be formulated. And that is what Matlab does automatically:

>> K = 278.7369 * tf([1 0],1)

K =
 
  278.7 s
 
Continuous-time transfer function.
Model Properties
>> ss(K)

ans =
 
  A = 
       x1  x2
   x1   1   0
   x2   0   1
 
  B = 
        u1
   x1    0
   x2  -16
 
  C = 
          x1     x2
   y1  17.42      0
 
  D = 
       u1
   y1   0
 
  E = 
       x1  x2
   x1   0   1
   x2   0   0
 
Continuous-time state-space model.
Model Properties

Matlab then proceeds with the descriptor description (dss instead of the more familiar ss class) when forming the open- and whater closed-loop transfer functions. All these are then modelled by the five matrices A, B, C, D, and E. It is then incorrect to ignore the E matrix. For example when analysing the stability of the closed-loop system, checking the eigenvalues of the A matrix leads to incorrect conclusions.

In this particular case, since the original plant is strictly proper (more poles than zeros), it seems useful to form the open-loop transfer function in a way that avoids the need to form a standalone model of the derivative term first, as @baggepinnen suggests.

When doing this in Matlab, I first reduced the model of the plant a bit. The original model is rather redundant – most poles can be either exactly or approximately cancelled by zeros (I was just wondering, wasn’t the model obtained by some ARX-like identification procedure?).

>> G = ss(A,B,C,D); % SISO, order 163
G = minreal(G);  %       order 146
19 states removed.
>> Gred = reduce(G,3);
>> zpk(Gred)

ans =
 
  -0.0049294 (s+8.172) (s+33.99)
  ------------------------------
  (s+2.323) (s-2.047) (s+13.72)
 
Continuous-time zero/pole/gain model.
Model Properties

It appears that this simple model captures the main aspects of the original model quite well. It is strictly proper, it does contain a single unstable pole, and the Bode plots reasonably agree too. I found it more convenient to play with this simple model.

Now, I obtained the open-loop transfer function by manipulating the transfer functions instead of state-space models. Perhaps this could be done more efficiently (transforming the model description from the state-space format to the transfer function format and back is not particularly recommendable), anyway:

>> L = tf(Gred)*K 

L =
 
  -1.374 s^3 - 57.93 s^2 - 381.6 s
  ---------------------------------
  s^3 + 13.99 s^2 - 0.9693 s - 65.2
 
Continuous-time transfer function.
Model Properties
>> zpk(L)

ans =
 
  -1.374 s (s+33.99) (s+8.172)
  -----------------------------
  (s+13.72) (s+2.323) (s-2.047)
 
Continuous-time zero/pole/gain model.
Model Properties

Obviously, the open-loop transfer function is proper and can be conveniently described by a state-space model. But what we are ultimately after is that one particular closed-loop transfer function

>> U = feedback(tf(Gred),K)

U =
 
     0.004929 s^2 + 0.2078 s + 1.369
  --------------------------------------
  0.374 s^3 + 43.94 s^2 + 382.6 s + 65.2
 
Continuous-time transfer function.
Model Properties
>> zpk(U)

ans =
 
  0.01318 (s+33.99) (s+8.172)
  ----------------------------
  (s+108) (s+9.282) (s+0.1739)
 
Continuous-time zero/pole/gain model.
Model Properties

All the poles are obviously stable. This can also be seen upon drawing a root locus:

For sufficiently large gain (just by trying it appears the treshold is around 0.75) the close-loop pole that lies on the locus emanating from the only unstable open-loop pole “jumps” over the infinity to the left half plane.

1 Like

That explains why my calculations with DescriptorSystems.jl would not yield a stable closed-loop system, I created the derivative term s instead of 278.7369*s.

So it seems one must indeed have a pure derivative for the unstable pole to jump over the infinity here, when using the filtered derivative, I would only ever get it to move “far to the right” but never actually jump over.

@lstagner What does this derivative term correspond to in reality? An ideal derivative term isn’t really possible to implement in practice, and if you indeed require it to be ideal, you might not get a stable system if you try to implement this closed loop.


Some additional notes:

3 Likes