`minreal` works -- but gives "wrong" result

B-LIE commented 1 year ago

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 commented 1 year ago

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

B-LIE commented 1 year ago

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

B-LIE commented 1 year ago

I have opened the issue in MatrixPencils.jl.

baggepinnen commented 1 year ago

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 commented 1 year ago

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

baggepinnen commented 1 year ago

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 commented 1 year ago

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

B-LIE commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

baggepinnen commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

B-LIE commented 1 year ago

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 commented 1 year ago

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

B-LIE commented 1 year ago

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 commented 1 year ago

I've fixed so that balance_statespace works for named systems in #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 commented 1 year ago

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 commented 1 year ago

The PR I linked above has not been merged yet 😉

837

baggepinnen commented 1 year ago

A new release will be out in a few minutes

B-LIE commented 1 year ago

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:

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).
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.
Now suddenly bodeplot doesn't work on "named ss" structure, but it works on the standard ss structure.
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 commented 1 year ago

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 commented 1 year ago

Named systems will work after the merge of

https://github.com/JuliaControl/RobustAndOptimalControl.jl/pull/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, ))

JuliaControl / ControlSystems.jl

`minreal` works -- but gives "wrong" result #835

837