# 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.