add more robustness options in causal simplification

Author: baggepinnen

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"]

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")

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)))
+# bodeplot([G, G2, G3, G4, G5], exp10.(LinRange(-2, 2, 200)))