test on v1.6 (#31)

* test on v1.6

* 1.6 compat

Author: baggepinnen

.github/workflows/CI.yml

@@ -10,7 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.7'
+          - '1.6'
         os:
           - ubuntu-latest
         arch:

Project.toml

@@ -1,7 +1,7 @@
 name = "RobustAndOptimalControl"
 uuid = "21fd56a4-db03-40ee-82ee-a87907bee541"
 authors = ["Fredrik Bagge Carlson", "Marcus Greiff"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
@@ -21,7 +21,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 ComponentArrays = "0.9, 0.10"
-ControlSystems = "0.11.2"
+ControlSystems = "0.11.6"
 Distributions = "0.25"
 IntervalArithmetic = "0.20"
 MatrixEquations = "2"
@@ -30,7 +30,7 @@ MonteCarloMeasurements = "1.0"
 Optim = "1.5"
 RecipesBase = "1"
 UnPack = "1.0"
-julia = "1.7"
+julia = "1.6"
 
 [extras]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/src/uncertainty.md

@@ -235,13 +235,26 @@ P_{32} & P_{33}\\
 ```
 into $P_{22}$ for the purposes of uncertainty analysis (use `ss` to convert it to a standard statespace object), and later use [`partition`](@ref) to recover the internal block structure. 
 
-Given an [`UncertainSS`](@ref) $P$, we can close the loop around $\Delta$ by calling `starprod(Δ, P)` or `lft(P, Δ, :u)` (note the different order of the arguments), and given an [`ExtendedStateSpace`](@ref), we can close the loop around `K` by calling `starprod(P, K)` or `lft(P, K)` (using positive feedback). This works even if `P` is a regular statespace object, in which case the convention is that the inputs and outputs are ordered as in the block diagrams above. The number of signals that will be connected by [`lft`](@ref) is determined by the input-output arity of $K$ and $\Delta$ respectively.
+Given an [`UncertainSS`](@ref) $P$, we can close the loop around $\Delta$ by calling `lft(P, Δ, :u)`, and given an [`ExtendedStateSpace`](@ref), we can close the loop around `K` by calling `starprod(P, K)` or `lft(P, K)` (using positive feedback). This works even if `P` is a regular statespace object, in which case the convention is that the inputs and outputs are ordered as in the block diagrams above. The number of signals that will be connected by [`lft`](@ref) is determined by the input-output arity of $K$ and $\Delta$ respectively.
 
 We have the following methods for `lft` (in addition to the standard ones in ControlSystems.jl)
 - `lft(G::UncertainSS, K::LTISystem)` forms the lower LFT closing the loop around $K$.
 - `lft(G::UncertainSS, Δ::AbstractArray=G.Δ)` forms the upper LFT closing the loop around $\Delta$.
 - `lft(G::ExtendedStateSpace, K)` forms the lower LFT closing the loop around $K$.
 
+### Robust stability and performance
+To check robust stability of the system in the last block diagram (with or without $z$ and $w$), we can use the functions [`structured_singular_value`](@ref), [`robstab`](@ref) and [`diskmargin`](@ref).
+
+Currently, [`structured_singular_value`](@ref) is rather limited and supports diagonal complex blocks only. If $\Delta$ is a single full complex block, `opnorm(P.M) < 1` is the condition for stability.
+
+Robust performance can be verified by introducing an additional fictitious "performance perturbation" $\Delta_p$ which is a full complex block, around which we close the loop from $z$ to $w$ and check the [`structured_singular_value`](@ref) with the augmented perturbation block
+```math
+\Delta_a = \begin{bmatrix}
+\Delta & 0\\
+0      & \Delta_p\\
+\end{bmatrix}
+```
+
 
 
 ### Examples

src/uncertainty_interface.jl

@@ -148,10 +148,18 @@ function Base.promote_rule(::Type{U}, ::Type{L}) where {U <: UncertainSS, L <: L
     UncertainSS
 end
 
+function Base.promote_rule(::Type{U}, ::Type{L}) where {U <: UncertainSS, L <: AbstractArray}
+    UncertainSS
+end
+
 function Base.convert(::Type{U}, d::δ) where {U <: UncertainSS}
     uss(d)
 end
 
+function Base.convert(::Type{U}, d::AbstractArray) where {U <: UncertainSS}
+    uss(d)
+end
+
 function Base.convert(::Type{U}, s::LTISystem) where {U <: UncertainSS}
     s isa UncertainSS && return s
     sys = ss(s)
@@ -165,8 +173,8 @@ Base.:+(n::δ, sys::UncertainSS{TE}) where TE <: ControlSystems.TimeEvolution =
 Base.:+(sys::UncertainSS{TE}, n::Number) where TE <: ControlSystems.TimeEvolution = +(sys, uss(n))
 Base.:+(n::Number, sys::UncertainSS{TE}) where TE <: ControlSystems.TimeEvolution = +(uss(n), sys)
 
-Base.:*(G::LTISystem, d::UncertainElement) = uss(G) * uss(d)
-Base.:*(d::UncertainElement, G::LTISystem) = uss(d) * uss(G)
+Base.:*(G::LTISystem, d::UncertainElement) = uss(ss(G)) * uss(d)
+Base.:*(d::UncertainElement, G::LTISystem) = uss(d) * uss(ss(G))
 
 """
     uss(D11, D12, D21, D22, Δ, Ts = nothing)
@@ -207,11 +215,16 @@ uss(n::Number, Ts=nothing) = uss(zeros(0,0), zeros(0,1), zeros(1,0), 1, [], Ts)
 """
     uss(D::AbstractArray, Δ, Ts = nothing)
 
-If only a single `D` matrix is provided, it's treated as `D11`.
+If only a single `D` matrix is provided, it's treated as `D11` if Δ is given, and as `D22` if no Δ is provided.
 """
-function uss(D::AbstractArray, Δ, Ts=nothing)
-    length(Δ) == size(D,1) || throw(DimensionMismatch("length(Δ) != size(D,1)"))
-    uss(D, zeros(size(D,1),0), zeros(0,size(D,2)), zeros(0,0), Δ, Ts)
+function uss(D::AbstractArray, Δ=[], Ts=nothing)
+    if length(Δ) == size(D,1)
+        uss(D, zeros(size(D,1),0), zeros(0,size(D,2)), zeros(0,0), Δ, Ts)
+    elseif isempty(Δ)
+        uss(zeros(0,0), zeros(0,size(D,2)), zeros(size(D,1),0), D, Δ, Ts)
+    else
+        throw(DimensionMismatch("length(Δ) != size(D,1)"))
+    end
 end
 
 uss(s::UncertainSS) = s
@@ -474,7 +487,7 @@ function ss2particles(G::Vector{<:AbstractStateSpace})
     B = reduce(hcat, vec.(getproperty.(G, :B))) |> pdp
     C = reduce(hcat, vec.(getproperty.(G, :C))) |> pdp
     D = reduce(hcat, vec.(getproperty.(G, :D))) |> pdp
-    (; nx,ny,nu) = G[1]
+    @unpack nx,ny,nu = G[1]
     A = reshape(A, nx, nx)
     B = reshape(B, nx, nu)
     C = reshape(C, ny, nx)

test/test_diskmargin.jl

@@ -55,7 +55,7 @@ plot(dms)
 ## Loop at a time
 a = 10
 P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
-K = ss(1.0I(2))
+K = ss(I(2))
 
 Li = K*P
 Lo = P*K
@@ -86,7 +86,7 @@ w = 2π .* exp10.(LinRange(-2, 2, 300))
 dm = loop_diskmargin(L3, 0, 4.05)
 @test dm[1].α ≈ 0.794418036911981 rtol=1e-3
 
-dm = diskmargin(L3, ss(1.0I(3), L3.Ts), 0, w)
+dm = diskmargin(L3, ss(I(3), L3.Ts), 0, w)
 plot(dm.simultaneous)
 plot!(dm.simultaneous_input)
 plot!(dm.simultaneous_output)
@@ -98,7 +98,7 @@ plot!(dm.output)
 
 a = 10
 P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
-K = ss(1.0I(2))
+K = ss(I(2))
 
 w = 2π .* exp10.(LinRange(-2, 2, 300))
 # @time bisect_a(P, K, w)
@@ -146,8 +146,8 @@ L3 = let
 end
 
 
-a = bisect_a(L3, ss(1.0I(3), L3.Ts), w; tol=2e-3)
-au = bisect_a(L3, ss(1.0I(3), L3.Ts), w; tol=2e-3, upper=true, N=256)
+a = bisect_a(L3, ss(I(3), L3.Ts), w; tol=2e-3)
+au = bisect_a(L3, ss(I(3), L3.Ts), w; tol=2e-3, upper=true, N=256)
 plot(w, a, xscale=:log10, xlabel="Frequency", ylims=(0,Inf))
 plot!(w, au, xscale=:log10, xlabel="Frequency", ylims=(0,Inf))
 

test/test_uncertainty.jl

@@ -117,13 +117,13 @@ end
 
 
 ##
-@test uss(ss(I(2))).M.D == I
+@test uss(I(2)).M.D == I
 @test sum(δss(2, 2).D) == 4
 @test sum(δss(4, 4).D) == 8
 
 ##
 H = δss(2, 3)
-W = ss(tf([1, .1],[.1, 1]))*I(2)
+W = tf([1, .1],[.1, 1]) .* I(2)
 WH = W*H
 @test WH.D == [zeros(3,2) I(3); 10I(2) zeros(2,3)]
 
@@ -155,14 +155,14 @@ temp = (W*d)
 @test temp.nu == temp.ny == 1
 @test temp.nz == temp.nw == 1
 
-temp = (ss(I(1)) + W*d)
+temp = I(1) + W*d
 @test temp.nu == temp.ny == 1
 @test temp.nz == temp.nw == 1
 
 @test length(d.Δ) == 1
 
 
-Pd = Pn*(ss(I(1)) + W*d)
+Pd = Pn*(I(1) + W*d)
 
 
 usyss = sminreal(system_mapping(Pd))
@@ -182,7 +182,7 @@ end
 
 Pn = ssrand(3,4,5)
 @test δss(4,4, bound=0.2).D == [0I(4) sqrt(0.2)*I(4); sqrt(0.2)*I(4) 0I(4)]
-mu = ss(I(4)) + δss(4,4, bound=0.2)
+mu = I(4) + δss(4,4, bound=0.2)
 
 
 @test mu.D == [0I(4) sqrt(0.2)*I(4); sqrt(0.2)*I(4) I(4)]
@@ -199,8 +199,8 @@ Pn2 = system_mapping(P)
 ## this time mimo real
 delta = uss([δr(), δr()])
 a = 1
-P = ss([0 a; -a -1], I(2), [1 a; 0 1], 0)* (ss(1.0*I(2)) + delta)
-K = ss(1.0I(2))
+P = ss([0 a; -a -1], I(2), [1 a; 0 1], 0)* (I(2) + delta)
+K = ss(I(2))
 
 G = lft(P, -K)
 hn = norm(G, Inf)
@@ -260,7 +260,7 @@ blocks, M = RobustAndOptimalControl.blocksort(P)
 
 w = exp10.(LinRange(-2, 2, 500))
 delta = δss(1,1)
-P = ss(tf(1,[1, .2, 1])) * (1+0.2*delta)
+P = (tf(1,[1, .2, 1])) * (1+0.2*delta)
 s = tf("s")
 K = ss(1 + 2/s + 0.9s/(0.1s+1))
 Gcl = lft(P, -K)
@@ -278,7 +278,7 @@ mu = structured_singular_value(Gcl, w)
 ## same as above but with scalar instead of 1×1 system
 w = exp10.(LinRange(-2, 2, 500))
 delta = δc()
-P = ss(tf(1,[1, .2, 1])) * (1+0.2*delta)
+P = (tf(1,[1, .2, 1])) * (1+0.2*delta)
 s = tf("s")
 K = ss(1 + 2/s + 0.9s/(0.1s+1))
 Gcl = lft(P, -K)
@@ -295,7 +295,7 @@ mu = structured_singular_value(Gcl, w)
 
 ## this time mimo complex
 delta = uss([δc(), δc()])
-P = ss([0 1; 0 0], I(2), [1 0], 0) * (ss(1.0I(2)) + delta)
+P = ss([0 1; 0 0], I(2), [1 0], 0) * (I(2) + delta)
 # diagonal input uncertainty
 
 K = ss([1;1])