diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index eb34e599..127881c2 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,7 +10,7 @@ jobs: fail-fast: false matrix: version: - - '1.7' + - '1.6' os: - ubuntu-latest arch: diff --git a/Project.toml b/Project.toml index 86c8fbf2..83923b27 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/uncertainty.md b/docs/src/uncertainty.md index 85b8b2c5..4cf036a5 100644 --- a/docs/src/uncertainty.md +++ b/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 diff --git a/src/uncertainty_interface.jl b/src/uncertainty_interface.jl index a397eb87..ae24c842 100644 --- a/src/uncertainty_interface.jl +++ b/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) diff --git a/test/test_diskmargin.jl b/test/test_diskmargin.jl index 4414853c..66b1fa0b 100644 --- a/test/test_diskmargin.jl +++ b/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)) diff --git a/test/test_uncertainty.jl b/test/test_uncertainty.jl index 48b583ff..da2b869e 100644 --- a/test/test_uncertainty.jl +++ b/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])