From f7089f0e2a87f66860c32548d5bbda2884a4b4a2 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 6 Aug 2025 13:46:37 +0200 Subject: [PATCH] add more robustness options in causal simplification --- Project.toml | 6 +++--- src/ode_system.jl | 31 +++++++++++-------------------- test/test_ODESystem.jl | 26 ++++++++++++++++++++------ 3 files changed, 34 insertions(+), 29 deletions(-) diff --git a/Project.toml b/Project.toml index b317ca4..ba376e7 100644 --- a/Project.toml +++ b/Project.toml @@ -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" @@ -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"] diff --git a/src/ode_system.jl b/src/ode_system.jl index f9aeae8..8689d8a 100644 --- a/src/ode_system.jl +++ b/src/ode_system.jl @@ -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)' @@ -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") diff --git a/test/test_ODESystem.jl b/test/test_ODESystem.jl index b3613fa..a4d8594 100644 --- a/test/test_ODESystem.jl +++ b/test/test_ODESystem.jl @@ -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 @@ -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))) \ No newline at end of file +# bodeplot([G, G2, G3, G4, G5], exp10.(LinRange(-2, 2, 200))) \ No newline at end of file