minreal works -- but gives "wrong" result

[B-LIE]

Today, minreal works on Windows 11/ Julia v1.9. THANKS!

The result from minreal is, however, "wrong" -- or perhaps extremely sensitive to tolerances...

Here is my system eigenvalues -- I insert copyable matrices for the system below:

• From physics, there should be 4 states (independent differential equations) and 2 algebraic equations.
• If I apply the minreal function to the system with default tolerances, the system is reduced to a 3 state ODE.
• If I specify the relative tolerance rtol, the system minimal realization has 3 states all the way to rtol=1e-16.
• If I specify relative tolerance at rtol=1e-17, the minimal realization gives 1 state
• If I specify relative tolerance rtol at 1e-18 and smaller, the minimal realization gives 6 states

A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]

B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]

C = [0.0 0.0 0.0 1.0 0.0 0.0]

D = [0.0]

linsys = ss(A,B,C,D)

[baggepinnen]

I would encourage you to open this issue in MatrixPencils.jl where minreal is implemented.

[B-LIE]

OK -- I can do that. (In MatrixPencils, the tolerance atol1 is probably the same as tolerance atol in ControlSystems??)

[B-LIE]

I have opened the issue in MatrixPencils.jl.

[baggepinnen]

You have two pole-zero cancellations in the origin and one in -7.27070, so it makes sense that minreal removes three state variables

julia> poles(linsys)
6-element Vector{ComplexF64}:
     -7.270709095526688 + 0.0im
     -6.683187769600364 + 0.0im
    -0.5184873504279723 + 0.9245137064074253im
    -0.5184873504279723 - 0.9245137064074253im
 1.7969420165981442e-17 + 0.0im
  4.638847626343683e-16 + 0.0im

julia> tzeros(linsys)
3-element Vector{Float64}:
  3.0660800180583606e-16
 -7.270709095526687
 -1.4243805976205752e-16

minreal does not really care about the physics, only the input-output mapping of the linear operator and this mapping is subject to the pole-zero cancellations

This is our implementation of minreal, i.e., it just wraps MatrixPencils.lsminreal

function minreal(sys::T, tol=nothing; fast=false, atol=0.0, kwargs...) where T <: AbstractStateSpace
    A,B,C,D = ssdata(sys)
    if tol !== nothing
        atol == 0 || atol == tol || error("Both positional argument `tol` and keyword argument `atol` were set but were not equal. `tol` is provided for backwards compat and can not be set to another value than `atol`.")
        atol = tol
    end
    Ar, Br, Cr = MatrixPencils.lsminreal(A,B,C; atol, fast, kwargs...)
    basetype(T)(Ar,Br,Cr,D, ntuple(i->getfield(sys, i+4), fieldcount(T)-4)...)
end

[B-LIE]

Hm... OK -- but the bode plots before and after minimal realization are very different.

[baggepinnen]

I think that's because your system is very poorly balanced, if you perform balancing, you notice that the reduced-order system has the same frequency response

using ControlSystemsBase, Plots
A = [-6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9 0.0; 0.0 -6.537773175952662 0.0 0.0 0.0 -9.892378564622923e-9; 2.0163803998106024e8 2.0163803998106024e8 -0.006223894167415392 -1.551620418759878e8 0.002358202548321148 0.002358202548321148; 0.0 0.0 5.063545034365582e-9 -0.4479539754649166 0.0 0.0; -2.824060629317756e8 2.0198389074625736e8 -0.006234569427701143 -1.5542817673286995e8 -0.7305736722226711 0.0023622473513548576; 2.0198389074625736e8 -2.824060629317756e8 -0.006234569427701143 -1.5542817673286995e8 0.0023622473513548576 -0.7305736722226711]

B = [0.004019511633336128; 0.004019511633336128; 0.0; 0.0; 297809.51426114445; 297809.51426114445]

C = [0.0 0.0 0.0 1.0 0.0 0.0]

D = [0.0]

linsys = ss(A,B,C,D)
bsys, T = balance_statespace(linsys)
rsys = minreal(bsys)

w = exp10.(LinRange(-2, 2, 200))

bodeplot(linsys, w)
bodeplot!(bsys, w)
bodeplot!(rsys, w, l=(:dash, ))

[baggepinnen]

Maybe we should perform balancing by default in the plotting functions?

[B-LIE]

If I do not do balancing, the original system has 2 real poles, 2 complex conjugate poles, and 2 stray poles.

The original system has 1 real zero and 2 stray zeros.

A correct minimal realization should cancel one of the real poles with the real zero, and lead to a minimal realization with 1 real pole and 2 complex conjugate poles.

However, minreal leads to 3 real poles.

[baggepinnen]

minreal on the balanced model, like in the example above, does the right thing

julia> poles(rsys)
3-element Vector{ComplexF64}:
  -6.683187769600371 + 0.0im
 -0.5184873504279726 + 0.9245137064074256im
 -0.5184873504279726 - 0.9245137064074256im

in general, we do not perform automatic balancing everywhere since it's in some situations an expensive operation, and also changes the state representation, instead, we require the user to explicitly opt into the balancing when needed.

[B-LIE]

I get an error message when I try to use balance_statespace... my set-up doesn't find this method

[baggepinnen]

I can't do much without more information, the function is exported so it should work if you have loaded CS or CSBase

[B-LIE]

OK. Could it be that I generated the linear system using ControlSystemsMTK and that balance_statespace doesn't recognize that system, but recognizes a system generated from "pure" matrices?

I'll check.

[baggepinnen]

The problem is that you have not loaded ControlSystemsBase. ControlSystemsMTK does not contain this functionality

[B-LIE]

That is not the problem -- I loaded ControlSystemsBase, but still balance_statespace does not work on the the linear system generated by ControlSystemsMTK. I'll check a little bit more.

[B-LIE]

OK -- now it seems to work. Don't know what I did wrong "in the heat of the moment"...

[B-LIE]

I get an error message when I try to use balance_statespace... my set-up doesn't find this method

The problem persists, and I found the problem...

• Method balance_statespace works on structures of type StateSpace{ControlSystemsBase.Continuous, Float64}

• Method balance_statespace does not work on structures of type NamedStateSpace{ControlSystemsBase.Continuous, Float64}.

[baggepinnen]

• I've fixed so that balance_statespace works for named systems in make balance_statespace work for more system types #837
• I've also added a section to the docs that talk about numerical balancing https://juliacontrol.github.io/ControlSystems.jl/dev/man/numerical/#State-space-balancing
• bodeplot and other plots now perform the balancing automatically

[B-LIE]

Hm... still doesn't work properly. I have updated the packages, and this time, function balance_statespace does run on a structure of type NamedStateSpace{ControlSystemsBase.Continuous, Float64}.

However, doing balance_statespace on the named statespace structure followed by minreal gives the same poles as if I drop doing balance_statespace, while if I instead create an unnamed statespace structure before doing balance_statespace, things work.

I know you don't have my NamedStateSpace structure, but here is what I get:

# Linearization using ControlSystemsMTK
linsys = named_ss(sys_p, sys_in, sys_out; op=op_0)
# poles
poles(linsys)

6-element Vector{ComplexF64}:
     -7.270709095526688 + 0.0im
     -6.683187769600364 + 0.0im
    -0.5184873504279723 + 0.9245137064074253im
    -0.5184873504279723 - 0.9245137064074253im
 1.7969420165981442e-17 + 0.0im
  4.638847626343683e-16 + 0.0im

# minimal realization (unbalanced) + poles
linsys_m = minreal(linsys)
poles(linsys_m)

3-element Vector{Float64}:
 -0.850985115940591
 -6.546113999987366
 -4.460027611772281

# doing balance_statespace on Named SS
linsys_b = balance_statespace(linsys)[1]
poles(linsys_b)

3-element Vector{Float64}:
 -0.850985115940591
 -6.546113999987366
 -4.460027611772281

# The above shows that doing balance_statespace on
# *Named* statespace does not improve the accuracy

# Instead, build an unnamed SS:
linsys_unm = ss(linsys.A,linsys.B,linsys.C,linsys.D)
poles(linsys_unm)

6-element Vector{ComplexF64}:
     -7.270709095526688 + 0.0im
     -6.683187769600364 + 0.0im
    -0.5184873504279723 + 0.9245137064074253im
    -0.5184873504279723 - 0.9245137064074253im
 1.7969420165981442e-17 + 0.0im
  4.638847626343683e-16 + 0.0im

# Above is the same as original named SS

# Using balance_statespace on unnamed SS + minreal
balance_statespace(linsys_unm)[1] |> minreal |> poles

3-element Vector{ComplexF64}:
  -6.683187769600371 + 0.0im
 -0.5184873504279726 + 0.9245137064074256im
 -0.5184873504279726 - 0.9245137064074256im

# Notice that *this time*, the balancing has an effect

[baggepinnen]

The PR I linked above has not been merged yet 😉

• make balance_statespace work for more system types #837

[baggepinnen]

A new release will be out in a few minutes

[B-LIE]

Hm... I still have problems with "named statespace" structures.

I'm on Julia v1.9.0, with packages:

status `[C:\Users\Bernt\.julia\environments\v1.9\Project.toml](file:///C:/Users/Bernt/.julia/environments/v1.9/Project.toml)`
  [a6e380b2] ControlSystems v1.7.2
  [aaaaaaaa] ControlSystemsBase v1.4.4
  [687d7614] ControlSystemsMTK v0.1.12
  [0c46a032] DifferentialEquations v7.7.0
  [31c24e10] Distributions v0.25.90
  [b964fa9f] LaTeXStrings v1.3.0
  [961ee093] ModelingToolkit v8.55.1
  [91a5bcdd] Plots v1.38.11

Here is what does and doesn't work:

1. minreal now works on both "named ss" (from MTK) and the standard unnamed ss structure. The result seems to be the same, and it suffers from numeric problems, i.e., finds "wrong" approximations (with solely real poles).
2. balance_statespace() doesn't work on a "named ss" structure, but it works on a standard/unnamed ss object and gives the expected 1 real and 2 complex conjugate poles.
3. Now suddenly bodeplot doesn't work on "named ss" structure, but it works on the standard ss structure.
4. However, the "bode" command works on both structures.

Perhaps I should send my code. I realize that it is difficult to debug/check when you don't have the named linsys structure available.

My code is in a Jupyter notebook in VSCode. I guess you don't like to work with such notebooks. It will take me some effort to ,convert the code into a standard .jl file. Is there another way to store the named ss file so that I can send it to you guys?

[B-LIE]

I notice that now, the original "named ss" converted to an unnamed ss (ss(linsys.A,linsys.B,linsys.C,linsys.D)) and the balance_statespace followed by minimal realization seem to give identical bodeplots (which they should have). Hm...

[baggepinnen]

Named systems will work after the merge of

• add method balance_statespace(::NamedStateSpace) RobustAndOptimalControl.jl#91

I now get

for the code posted above, and the same result if I additionally do

using RobustAndOptimalControl
nss = named_ss(linsys)

bsys, T = balance_statespace(nss)
rsys = minreal(bsys)

w = exp10.(LinRange(-2, 2, 200))

bodeplot(nss, w)
bodeplot!(bsys, w)
bodeplot!(rsys, w, l=(:dash, ))