Skip to content

add more robustness options in causal simplification #79

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
baggepinnen merged 1 commit into from
Aug 6, 2025
Merged

add more robustness options in causal simplification #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
name = "ControlSystemsMTK"
uuid = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
authors = ["Fredrik Bagge Carlson"]
version = "2.4.1"

version = "2.5.0"

[deps]
ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
Expand All @@ -28,9 +27,10 @@ julia = "1.6"

[extras]
ControlSystems = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test", "ControlSystems", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqNonlinearSolve"]
test = ["Test", "ControlSystems", "GenericLinearAlgebra", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqNonlinearSolve"]
31 changes: 11 additions & 20 deletions src/ode_system.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,26 +151,28 @@ function get_x_names(lsys, sys; descriptor, simple_infeigs, balance)
end

"""
causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; descriptor=true, simple_infeigs=true, balance=false)
causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; descriptor=true, simple_infeigs=true, balance=false, big=false)

If `descriptor = true`, the function `DescriptorSystems.dss2ss` is used. In this case,
- `balance`: indicates whether to balance the system using `DescriptorSystems.gprescale` before conversion to `StateSpace`. Balancing changes the state realization (through scaling).
- `simple_infeigs`: if set to false, further simplification may be performed in some cases.

If `descriptor = false`, the argument `big = true` performs computations in `BigFloat` precision, useful for poorly scaled systems.
The argument `big = true` performs computations in `BigFloat` precision, useful for poorly scaled systems. This may require the user to install and load `GenericLinearAlgebra` (if you get error `no method matching svd!(::Matrix{BigFloat})`).
"""
function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; balance=false, descriptor=true, simple_infeigs = true, big = false)
fm(x) = convert(Matrix{Float64}, x)
T = big ? BigFloat : Float64
b1 = big ? Base.big(1.0) : 1.0
fm(x) = convert(Matrix{T}, x)
nx = size(sys.A, 1)
ny = size(sys.C, 1)
ndu = length(u2duinds)
nu = size(sys.B, 2) - ndu
u_with_du_inds = first.(u2duinds)
duinds = last.(u2duinds)
B = sys.B[:, 1:nu]
B̄ = sys.B[:, duinds]
D = sys.D[:, 1:nu]
D̄ = sys.D[:, duinds]
B = b1*sys.B[:, 1:nu]
B̄ = b1*sys.B[:, duinds]
D = b1*sys.D[:, 1:nu]
D̄ = b1*sys.D[:, duinds]
iszero(fm(D̄)) || error("Nonzero feedthrough matrix from input derivative not supported")
if descriptor
Iu = u_with_du_inds .== (1:nu)'
Expand All @@ -191,22 +193,11 @@ function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; balance=fa

return ss(RobustAndOptimalControl.DescriptorSystems.dss2ss(dsys; simple_infeigs, fast=false)[1])
else
if big
b1 = Base.big(1.0)
ss(b1*sys.A, b1*B, b1*sys.C, b1*D) + diff_out(ss(b1*sys.A, b1*B̄, b1*sys.C, b1*D̄))
else
b = balance ? s->balance_statespace(sminreal(s))[1] : identity
b(ss(sys.A, B, sys.C, D)) + b(ss(sys.A, B̄, sys.C, D̄))*tf('s')
end
b = balance ? s->balance_statespace(sminreal(s))[1] : identity
b(ss(sys.A, B, sys.C, D)) + b(ss(sys.A, B̄, sys.C, D̄))*tf('s')
end
end

function diff_out(sys)
A,B,C,D = ssdata(sys)
iszero(D) || error("Nonzero feedthrough matrix not supported")
ss(A,B,C*A, C*B, sys.timeevol)
end


for f in [:sensitivity, :comp_sensitivity, :looptransfer]
fnn = Symbol("get_named_$f")
Expand Down
26 changes: 20 additions & 6 deletions test/test_ODESystem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -335,18 +335,25 @@ op2[cart.f] = 0
@test G.ny == length(lin_outputs)

## Test difficult `named_ss` simplification
using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl
using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl, Test, GenericLinearAlgebra
lsys = (A = [0.0 2.778983834717109e8 1.4122312296634873e6 0.0; 0.0 0.0 0.0 0.037848975765016724; 0.0 24.837541148074962 0.12622006230897712 0.0; -0.0 -4.620724819774693 -0.023481719514324866 -0.6841991610512456], B = [-5.042589978197361e8 0.0; -0.0 0.0; -45.068824982602656 -0.0; 8.384511049369085 54.98555939873381], C = [0.0 0.0 0.954929658551372 0.0], D = [0.0 0.0])

# lsys = (A = [-0.0075449237853825925 1.6716817118020731e-6 0.0; 1864.7356343162514 -0.4131578457122937 0.0; 0.011864343330426718 -2.6287085638214332e-6 0.0], B = [0.0 0.0; 0.0 52566.418015009294; 0.0 0.3284546792274811], C = [1.4683007399899215e8 0.0 0.0], D = [-9.157636303058283e7 0.0])

G = ControlSystemsMTK.causal_simplification(lsys, [1=>2])
G = minreal(G, 1e-12)
G2 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=false)
G2 = minreal(G2, 1e-12)
G3 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], big=true)
G3 = minreal(G3, 1e-12)
G4 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=true, balance=true)
G4 = minreal(G4, 1e-12)

w_test = [1e-5, 1, 10, 100]
@test freqresp(G, w_test) ≈ freqresp(G2, w_test) rtol=1e-6
@test freqresp(G, w_test) ≈ freqresp(G3, w_test) rtol=1e-6
@test freqresp(G, w_test) ≈ freqresp(G4, w_test) rtol=1e-6

@test dcgain(G, 1e-5)[] ≈ dcgain(G2, 1e-5)[] rtol=1e-3
@test freqresp(G, 1)[] ≈ freqresp(G2, 1)[]
@test freqresp(G, 10)[] ≈ freqresp(G2, 10)[]

z = 0.462726166562343204837317130554462562

Expand All @@ -366,15 +373,22 @@ w = exp10.(LinRange(-12, 2, 2000))

## artificial fully dense test

lsys = ssrand(2,2,3,proper=true)
lsys = let
tempA = [-0.8101413207021115 0.905054220094988 -0.12282983628692003; 0.9066960393438117 -1.3378400902008518 -0.0610101630607034; -0.07253623899748166 1.4118402411983302 -1.7957106725392422]
tempB = [-0.8984842931160537 1.4012460461970973; -0.6086273804588082 0.9336277670619954; -0.1376448663325406 0.19241847208402496]
tempC = [1.1228782382333735 0.041349950615147096 -1.3629397459682921; 0.10499710831314196 -1.228883432780842 0.044498460069701234]
tempD = [0.0 0.0; 0.0 0.0]
ss(tempA, tempB, tempC, tempD)
end
G = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=false)
G2 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=true)
G3 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=false, big=true)
G4 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=true, simple_infeigs=false, balance=true)
G5 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=true, big=true)

w_test = [1e-5, 1, 100]
@test freqresp(G, w_test) ≈ freqresp(G2, w_test)
@test freqresp(G, w_test) ≈ freqresp(G3, w_test)
@test freqresp(G, w_test) ≈ freqresp(G4, w_test)

# bodeplot([G, G2, G3, G4], exp10.(LinRange(-2, 2, 200)))
# bodeplot([G, G2, G3, G4, G5], exp10.(LinRange(-2, 2, 200)))