interface DescriptorSystems

Author: baggepinnen

Project.toml

@@ -6,6 +6,7 @@ version = "0.2.3"
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 ControlSystems = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
+DescriptorSystems = "a81e2ce2-54d1-11eb-2c75-db236b00f339"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,7 +22,8 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 ComponentArrays = "0.9, 0.10, 0.11"
-ControlSystems = "0.11.6"
+ControlSystems = "0.11.9"
+DescriptorSystems = "1.2"
 Distributions = "0.25"
 IntervalArithmetic = "0.20"
 MatrixEquations = "2"

src/RobustAndOptimalControl.jl

@@ -14,7 +14,11 @@ import IntervalArithmetic
 import IntervalArithmetic: Interval
 
 
-import MatrixPencils, MatrixEquations
+import MatrixPencils, MatrixEquations, DescriptorSystems
+using DescriptorSystems: dss
+
+export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap
+include("descriptor.jl")
 
 export ExtendedStateSpace, system_mapping, performance_mapping, noise_mapping, ssdata_e, partition, ss
 include("ExtendedStateSpace.jl")

src/descriptor.jl

@@ -0,0 +1,84 @@
+
+"""
+    DescriptorSystems.dss(sys::AbstractStateSpace)
+
+Convert `sys` to a descriptor statespace system from [DescriptorSystems.jl](https://andreasvarga.github.io/DescriptorSystems.jl/dev/index.html)
+"""
+function DescriptorSystems.dss(sys::AbstractStateSpace)
+    A,B,C,D = ssdata(sys)
+    DescriptorSystems.dss(A, B, C, D; Ts = isdiscrete(sys) ? sys.Ts : 0)
+end
+
+function ControlSystems.ss(G::DescriptorSystems.DescriptorStateSpace)
+    try
+        G,r = DescriptorSystems.gss2ss(G)
+        if r < size(G.A, 1)
+            G = DescriptorSystems.dss2ss(G, fast=false)
+        end
+    catch
+        G,_ = DescriptorSystems.dss2ss(G, simple_infeigs = false, fast=false)
+    end
+    G.Ts >= 0 ? ss(G.A, G.B, G.C, G.D) : ss(G.A, G.B, G.C, G.D, G.Ts)
+end
+
+"""
+    n, ω = hinfnorm2(sys::LTISystem; kwargs...)
+
+A numerically robust version of `hinfnorm` using DescriptorSystems.jl
+
+For keyword arguments, see the docstring of DescriptorSystems.ghinfnorm, reproduced below
+$(@doc(DescriptorSystems.ghinfnorm))
+"""
+function hinfnorm2(sys::LTISystem; kwargs...)
+    DescriptorSystems.ghinfnorm(dss(ss(sys)); kwargs...)
+end
+
+"""
+    n = h2norm(sys::LTISystem; kwargs...)
+
+A numerically robust version of `norm` using DescriptorSystems.jl
+
+For keyword arguments, see the docstring of DescriptorSystems.gh2norm, reproduced below
+$(@doc(DescriptorSystems.gh2norm))
+"""
+function h2norm(sys::LTISystem; kwargs...)
+    DescriptorSystems.gh2norm(dss(ss(sys)); kwargs...)
+end
+
+"""
+    n, hsv = hankelnorm(sys::LTISystem; kwargs...)
+
+Compute the hankelnorm and the hankel singular values
+
+For keyword arguments, see the docstring of DescriptorSystems.ghanorm, reproduced below
+$(@doc(DescriptorSystems.ghanorm))
+"""
+function hankelnorm(sys::LTISystem; kwargs...)
+    DescriptorSystems.ghanorm(dss(ss(sys)); kwargs...)
+end
+
+"""
+    nugap(sys0::LTISystem, sys1::LTISystem; kwargs...)
+
+Compute the ν-gap metric between two systems.
+
+For keyword arguments, see the docstring of DescriptorSystems.gnugap, reproduced below
+$(@doc(DescriptorSystems.gnugap))
+"""
+function nugap(sys0::LTISystem, sys1::LTISystem; kwargs...)
+    DescriptorSystems.gnugap(dss(ss(sys0)), dss(ss(sys1)); kwargs...)
+end
+
+const νgap = nugap
+
+
+##
+
+function Base.:\(G1::AbstractStateSpace, G2::AbstractStateSpace)
+    g1,g2 = dss(G1), dss(G2)
+    try
+        ss(g1\g2)
+    catch
+        error("Failed to solve G1\\G2, the product might not be proper. To keep the product on descriptor form, use dss(G1)\\dss(G2)")
+    end
+end

test/test_descriptor.jl

@@ -0,0 +1,23 @@
+using RobustAndOptimalControl
+G1 = ssrand(2,3,4, proper=true)
+G2 = ssrand(2,3,4, proper=true)
+
+
+@test h2norm(G1-G1) < 1e-6
+@test hinfnorm2(G1-G1)[1] < 1e-6
+@test nugap(G1, G1)[1] < 1e-6
+@test hankelnorm(G1-G1)[1] < 1e-6
+
+
+
+
+@test count(1:100) do _
+    G1 = ssrand(1,1,3, proper=true)
+    G2 = ssrand(1,1,3, proper=true)
+    try
+        νgap(G1\G2, ss(inv(tf(G1))*tf(G2)))[1] < 1e-8
+    catch
+        false
+    end
+end >= 95
+

test/test_reduction.jl

@@ -1,5 +1,6 @@
 using RobustAndOptimalControl: ControlSystems
 using RobustAndOptimalControl, ControlSystems
+hinorm(x) = hinfnorm2(x)[1]
 
 sys = let
     _A = [-2.467030078096955 -1.7579086436230489 0.2642212618287617 -1.8139967744112642 0.9581654870157803 -1.1573282803258818 -2.409554294298775 1.2006959533714947 0.2661208318398341 0.9999177692959037 -0.3558865668962797 0.9416230396873052 -0.04719968633569187 1.9048861053182367 0.04105479772871481; 1.7913628319734 -2.762572481036251 1.1195768845996525 -0.5685580615144388 0.6726614103668934 1.0083399722116977 0.5812702064616616 -0.7829110624242633 -0.5724789585031453 0.2739159931739683 0.2757716078817949 -1.1772430269658822 -0.6381809965413997 -0.5567081736321586 0.3625746395689792; -0.005143392441108892 -1.7086991517875587 -1.717585106651333 1.5009119207685369 1.8738081845828842 -1.0939749683914888 -0.5264010960315267 0.22497823448076093 1.246044372509282 -0.5303445465251139 -0.5548425351067188 0.33884707949495035 0.5101563763162115 0.20539506245514677 0.26057781689650356; 0.604556580615853 1.2138308974005823 -2.059380841591183 -1.191665990395876 1.150735331516496 0.8110103414198591 -0.154085845809975 1.1119072467025846 0.11430616747450531 1.2819431451075372 0.921519933726758 0.4697140300022479 0.005096755657399703 -0.6758568389946026 0.23258764043708394; -0.3279280635059737 -0.3899618604207743 0.12657097470583745 0.8299061544706142 -2.6661102102915484 -0.33139258711550823 -0.012925000635212677 -0.42098301162071433 1.5144416025925511 -0.796229814982033 -0.0657636569172887 1.2193694935643702 -1.700730366738557 -0.480614226127818 0.5646449944635483; 0.2450793035368965 -1.8805339894110746 0.30667982435893143 0.9479179427233335 -1.2301806768500736 -3.136669168734015 -0.2481694826668052 0.380557579742639 -1.416193374954648 0.6987044103237035 -0.8853533466891377 -1.030250851383474 1.4525483702515944 0.36288849878688056 0.385057654662975; -2.0022001785391392 -0.2835451807794339 -0.26028704624052934 0.07217513465675776 0.08414471225474451 0.024542390758292842 -2.6763138550617667 0.011591697680752557 1.7133530919443034 0.17602929220838714 -0.6726695861811227 1.5441809892864893 -0.00015597605070589597 1.301000782058375 0.03551676106571688; -0.16249013289239514 0.13198122357867412 1.3236770560172366 0.3301510019655505 0.45179991512285006 0.07486588236983925 -1.6288494827036406 -1.5392724756838398 -0.6217318765063355 0.2858777881749128 -1.4166215418299406 -0.35825226416886896 1.8697114631330003 -0.3059163232549132 1.934355556146854; -0.3097986537718366 1.9026121572125902 -0.6714075159339921 -0.2504951318869757 -1.0960812903814519 0.551662724714743 -0.5206166054513874 -0.10057595046577751 -4.293450503414131 0.24751394898682735 0.5853194271595069 1.2157066833567851 -0.5437491323644809 0.03154589360470189 -1.3601989491821653; -1.8758967387335015 -0.08090769024414653 1.703446780809511e-5 -0.15226980080851651 -0.4474269369804122 -1.0627549839203736 1.1462306659249069 0.4825852860273927 1.1227319809099574 -2.5566746127781825 0.2744599619933535 1.514041120509779 0.8781499781073471 -0.7005275854271636 0.8540288603061493; 1.3102799390114126 2.0268412603985713 -0.26445369495117343 0.1673277757620624 0.33995168379334895 1.3424350975481887 -0.5262869773852599 1.772445455267548 0.7427203416097469 -0.6894548775913741 -4.648846973081957 0.6160480706048977 0.27693416312465485 -0.5843844597444985 -0.05046410580231136; -0.2667927151895562 -0.8394567537703295 0.7356804516400021 0.021352734365649183 1.6248465860949257 1.0375631106350516 -0.44557649819583367 -1.318671094466728 0.1198337174336245 0.028769936034846715 1.1651225420369067 -2.5802262184585376 -0.000961942820959872 0.3531297512430618 0.8801838685743447; -0.2340108920478918 -0.47020830048023965 0.0905085276240309 -1.159187923078835 -0.32263237957704327 -0.20398562485896923 -1.3882045935474685 0.6799676640930128 0.34600863452034686 -0.8950627077509554 0.37192323444628994 -0.17983996993004395 -3.254198004959321 -0.12498721862068969 -0.7674453121406706; -0.48890022859302995 0.7856693831160491 0.5530727968750053 0.28679439003867857 -0.353796629578632 1.3881939385941933 0.04935762134562992 -0.6880042763129538 1.3414828049193512 -1.6887109143288497 -0.846493928570519 -2.8372692994687916 -0.8703006934200097 -2.0491003542049913 -0.04580815825737827; 0.14776326405901308 -1.0181952707299127 0.17409898161459436 0.5045755253054325 1.2756549793646714 0.0513146326521474 1.9399237561334002 -1.3740316671021375 -1.6331290401813572 -0.4569927755409647 -0.6655757376846796 -0.11273958657543114 0.17478370341836297 0.9495603386771154 -2.069076828262977]
@@ -8,12 +9,10 @@ sys = let
     _D = [0.0]
     ss(_A, _B, _C, _D)
 end
-sysr = frequency_weighted_reduction(sys, 1, 1, 10)
+sysr = frequency_weighted_reduction(sys, 1, 1, 10, residual=false)
 sysr2 = baltrunc(sys, n=10)[1]
 @test sysr.nx == 10
-@test norm(sys-sysr) < 1e-5
-@test norm(sys-sysr, Inf) < 1e-5
-@test norm(sys-sysr, Inf) ≈ norm(sys-sysr2, Inf) rtol=1e-6
+@test hinorm(sys-sysr) < 1e-5
 # bodeplot([sys, sysr])
 
 
@@ -26,36 +25,35 @@ sys = let
     ss(_A, _B, _C, _D)
 end
 sysi = fudge_inv(sys)
-sysr = frequency_weighted_reduction(sys, sysi, 1, 5)
+sysr = frequency_weighted_reduction(sys, sysi, 1, 5, residual=false)
 sysr2 = baltrunc(sys, n=5)[1]
 @test sysr.nx == 5
 @test norm(sys-sysr) < 0.1
 @test norm(sysi*(sys-sysr)) < 0.05
 @test norm(sysi*(sys-sysr)) <= norm(sysi*(sys-sysr2))
-@test_broken norm(minreal(sysi*(sys-sysr)), Inf) <= norm(minreal(sysi*(sys-sysr2)), Inf) # broken due to non-convergence of hinfnorm
-
+@test hinorm(minreal(sysi*(sys-sysr))) <= hinorm(minreal(sysi*(sys-sysr2))) 
 
 
 Wi = tf(1, [1, 1])
-sysr = frequency_weighted_reduction(sys, 1, Wi, 5)
+sysr = frequency_weighted_reduction(sys, 1, Wi, 5, residual=false)
 sysr2 = baltrunc(sys, n=5)[1]
 @test sysr.nx == 5
 @test norm(sys-sysr) < 0.1
 @test norm((sys-sysr)*Wi) < 0.05
 @test norm((sys-sysr)*Wi) <= norm((sys-sysr2)*Wi)
-@test norm((sys-sysr)*Wi, Inf) <= norm((sys-sysr2)*Wi, Inf) # broken due to non-convergence of hinfnorm
+@test hinorm((sys-sysr)*Wi) <= hinorm((sys-sysr2)*Wi) 
 
 
 
 Wo = sysi
 Wi = tf(1, [1, 1])
-sysr = frequency_weighted_reduction(sys, Wo, Wi, 5)
+sysr = frequency_weighted_reduction(sys, Wo, Wi, 5, residual=false)
 sysr2 = baltrunc(sys, n=5)[1]
 @test sysr.nx == 5
 @test_broken norm(sys-sysr) < 0.1
 @test norm(Wo*(sys-sysr)*Wi) < 0.05
 # @test_broken norm(Wo*(sys-sysr)*Wi) <= norm(Wo*(sys-sysr2)*Wi)
-@test_broken norm(Wo*(sys-sysr)*Wi, Inf) <= norm(Wo*(sys-sysr2)*Wi, Inf) 
+@test_broken hinorm(Wo*(sys-sysr)*Wi) <= hinorm(Wo*(sys-sysr2)*Wi) 
 
 
 # residual
@@ -65,7 +63,7 @@ sysr2 = baltrunc(sys, n=3, residual=true)[1]
 @test sysr.nx == 3
 @test norm(sys-sysr, Inf) < 3
 @test norm(sysi*(sys-sysr), Inf) < 0.4
-@test norm(minreal(sysi*(sys-sysr)), Inf) <= norm(minreal(sysi*(sys-sysr2)), Inf) # broken due to non-convergence of hinfnorm
+@test hinorm(minreal(sysi*(sys-sysr))) <= hinorm(minreal(sysi*(sys-sysr2))) 
 @test dcgain(sys)[] ≈ dcgain(sysr)[] rtol=1e-10 # test the residual property