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

**URL:** <https://discourse.julialang.org/t/state-space-model-with-pid-controller-unstable-in-julia-compared-to-matlab/117688>\
**Category:** Specific Domains\
**Tags:** control, controlsystems\
**Created:** [July 31, 2024, 5:14pm UTC](https://discourse.julialang.org/t/state-space-model-with-pid-controller-unstable-in-julia-compared-to-matlab/117688 "2024-07-31T17:14:13Z")\
**Posts on this page:** 3\
**Page:** 2

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [August 2, 2024, 8:42am UTC](https://discourse.julialang.org/t/state-space-model-with-pid-controller-unstable-in-julia-compared-to-matlab/117688/22 "2024-08-02T08:42:25Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [August 2, 2024, 10:49am UTC](https://discourse.julialang.org/t/state-space-model-with-pid-controller-unstable-in-julia-compared-to-matlab/117688/23 "2024-08-02T10:49:22Z")

</div>

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:

```Matlab
>> 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?).

```Matlab
>> 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.

 ![bode](https://global.discourse-cdn.com/julialang/original/3X/8/b/8b82f5cd6bc9f55c710bb64f944b2658cf3e4b6f.png)

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:

```matlab
>> 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

```matlab
>> 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:

 ![rlocus](https://global.discourse-cdn.com/julialang/original/3X/2/f/2f299aa9ba3cb75e8dcd87d8382480d4dfc7bdb9.png)

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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [August 2, 2024, 11:23am UTC](https://discourse.julialang.org/t/state-space-model-with-pid-controller-unstable-in-julia-compared-to-matlab/117688/24 "2024-08-02T11:23:21Z")

</div>

> [@zdenek\_hurak](#):
>
> 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.

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:

- To reduce the order of unstable systems, we have the function [`RobustAndOptimalControl.baltrunc_unstab`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/api/#RobustAndOptimalControl.baltrunc_unstab)
- A [new PR](https://github.com/JuliaControl/ControlSystems.jl/pull/930) added support multiplication of statespace systems and improper transfer functions when the product is proper.

[Previous page](https://discourse.julialang.org/t/state-space-model-with-pid-controller-unstable-in-julia-compared-to-matlab/117688.md?page=1)