some fixes to uncertainty

docs/src/uncertainty.md

@@ -259,7 +259,7 @@ Ps = P*W
 ```
 `Ps` is now represented as a upper linear fractional transform (upper LFT).
 
-We can draw samples from this uncertainty representation like so
+We can draw samples from this uncertainty representation (sampling of $\Delta$ and closing the loop `starprod(Δ, Ps)`) like so
 ```@example satellite
 Psamples = rand(Ps, 100)
 sigmaplot(Psamples, w)

src/RobustAndOptimalControl.jl

@@ -53,7 +53,7 @@ export glover_mcfarlane, hanus
 include("glover_mcfarlane.jl")
 
 
-export diskmargin, Diskmargin, Disk, sim_diskmargin, loop_diskmargin, structured_singular_value, broken_feedback
+export diskmargin, Diskmargin, Disk, sim_diskmargin, loop_diskmargin, structured_singular_value, broken_feedback, robstab
 include("diskmargin.jl")
 include("mimo_diskmargin.jl")
 

src/diskmargin.jl

@@ -105,7 +105,7 @@ function Disk(γmin, γmax, c, r)
  ϕm = 90.0
  else
  ϕm = (1 + γmin*γmax) / (γmin + γmax)
- ϕm = ϕm >= 1 ? Inf : rad2deg(acos(ϕm))
+ ϕm = ϕm >= 1 ? 0 : ϕm <= -1 ? Inf : rad2deg(acos(ϕm))
  end
  Disk(γmin, γmax, c, r, ϕm)
 end
@@ -161,8 +161,8 @@ plot(w, dms)
 """
 function diskmargin(L::LTISystem, σ::Real=0; kwargs...)
  issiso(L) || return sim_diskmargin(L, σ)
- S̄ = 1/(1 + L) + (σ-1)/2
- n,ω0 = hinfnorm(S̄; kwargs...)
+ M = 1/(1 + L) + (σ-1)/2
+ n,ω0 = hinfnorm(M; kwargs...)
  diskmargin(L, σ, ω0)
 end
 
@@ -175,9 +175,9 @@ diskmargin(L::LTISystem, σ::Real, ω::AbstractArray) = map(w->diskmargin(L, σ,
 
 function diskmargin(L::LTISystem, σ::Real, ω0::Real)
  issiso(L) || return sim_diskmargin(L, σ, [ω0])[]
- S̄ = 1/(1 + L) + (σ-1)/2
+ M = 1/(1 + L) + (σ-1)/2
  freq = isdiscrete(L) ? cis(ω0*L.Ts) : complex(0, ω0)
- Sω = evalfr(S̄, freq)[]
+ Sω = evalfr(M, freq)[]
  αmax = 1/abs(Sω)
  δ0 = inv(Sω)
  dp = Disk(; α = αmax, σ)
@@ -243,25 +243,26 @@ end
  end
 end
 
-@recipe function plot(dm::AbstractVector{<:Diskmargin})
+@recipe function plot(dm::AbstractVector{<:Diskmargin}; lower=true, phase=true)
  w = [dm.ω0 for dm in dm]
- layout --> (2, 1)
+ layout --> (phase ? 2 : 1, 1)
  link --> :x
  @series begin
  subplot --> 1
  title --> "Gain margin"
- label --> ["Lower" "Upper"]
+ label --> (lower ? ["Lower" "Upper"] : "Upper")
  # xguide --> "Frequency"
  xscale --> :log10
  yscale --> :log10
  gma = getfield.(dm, :γmax)
  gmi = getfield.(dm, :γmin)
  # ylims --> (0, min(10, maximum(gma)))
- data = [gmi gma]
+ data = (lower ? [gmi gma] : gma)
  data = max.(0, data)
  replace!(data, 0 => -Inf)
  w, data
  end
+ phase || return
  @series begin
  subplot --> 2
  title --> "Phase margin"
@@ -278,11 +279,11 @@ end
 end
 
 
-@recipe function plot(dm::AbstractMatrix{Diskmargin})
+@recipe function plot(dm::AbstractMatrix{Diskmargin}; lower=true, phase=true)
  w = getfield.(dm[:,1], :ω0)
  ny = size(dm, 2)
  length(w) == size(dm, 1) || throw(ArgumentError("Frequency vector and diskmargin vector must have the same lengths."))
- layout --> (2, 1)
+ layout --> (phase ? 2 : 1, 1)
  link --> :x
  neg2inf(x) = x <= 0 ? Inf : x
  for i = 1:ny
@@ -291,16 +292,17 @@ end
  subplot --> 1
  title --> "Gain margin"
  # label --> permutedims([["Lower $i" for i = 1:ny]; ["Upper $i" for i = 1:ny]])
- label --> ["Upper"*is "Lower"*is]
+ label --> (lower ? ["Upper"*is "Lower"*is] : "Upper"*is)
  # xguide --> "Frequency"
  xscale --> :log10
  yscale --> :log10
  gma = getfield.(dm[:, i], :γmax) .|> neg2inf
  gmi = getfield.(dm[:, i], :γmin) .|> neg2inf
  replace!(gmi, 0 => -Inf)
  # ylims --> (0, min(10, maximum(gma)))
- w, [gmi gma]
+ w, (lower ? [gmi gma] : gma)
  end
+ phase || continue
  @series begin
  subplot --> 2
  title --> "Phase margin"

src/mimo_diskmargin.jl

@@ -14,7 +14,7 @@ end
 any0det(D::Matrix{<:Complex{<:Interval}}) = 0 ∈ det(D)
 
 """
- bisect_a(P, K, w; W = (:), Z = (:), au0 = 3.0, tol = 0.001, N = 32, upper = true, δ = δc)
+ bisect_a(P, K, w; W = (:), Z = (:), au0 = 3.0, tol = 0.001, N = 32, upper = true, δ = δc, σ=-1)
 
 EXPERIMENTAL AND SUBJECT TO BUGS, BREAKAGE AND REMOVAL
 
@@ -58,7 +58,7 @@ function bisect_a(args...; au0 = 3.0, tol=1e-3, kwargs...)
 end
 
 """
- M,D = get_M(P, K, w; W = (:), Z = (:), N = 32, upper = false, δ = δc)
+ M,D = get_M(P, K, w; W = (:), Z = (:), N = 32, upper = false, δ = δc, σ=-1)
 
 Return the frequency response of `M` in the `M-Δ` formulation that arises when individual, complex/real perturbations are introduced on inputs `W` and outputs `Z` (defaults to all).
 
@@ -72,22 +72,49 @@ Return the frequency response of `M` in the `M-Δ` formulation that arises when
 - `upper`: Indicate whether an upper or lower bound is to be computed
 - `δ = δc` for complex perturbations and `δ = δr` for real perturbations.
 """
-function get_M(P, K, w; W = (:), Z = (:), N = 32, upper=false, δ = δc)
+function get_M(P::LTISystem, K::LTISystem, w::AbstractVector; W = (:), Z = (:), N = 32, upper=false, δ = δc, σ = -1)
  Z1 = W2 = Z == (:) ? (1:P.ny) : Z
  W1 = Z2 = W == (:) ? (1:P.nu) : W
- ny,nu = length(Z2), length(W2)
+ ny,nu = length(Z1), length(W1)
  D = Δ(ny+nu, δ)
  if upper
  D = rand(D, N)
  else
  D = Diagonal([Interval(d) for d in diag(D)])
  end
- M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
+ if isempty(Z1) # No outputs
+ L = (K*P)[W1, W1]
+ elseif isempty(W1) # No inputs
+ L = (P*K)[Z1, Z1]
+ else
+ L = [ss(zeros(P.ny, P.ny), P.timeevol) P;-K ss(zeros(K.ny, K.ny), P.timeevol)]
+ L = L[[Z1; P.ny .+ Z2], [W1; P.nu .+ W2]]
+ end
+ n = L.ny
+ X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
+ M = starprod(X,L)
+ # M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
  M0 = freqresp(M, w)
  M0 = permutedims(M0, (2,3,1))
  M0, D
 end
 
+function get_M(L::LTISystem, w::AbstractVector; N = 32, upper=false, δ = δc, σ = -1)
+ ny,nu = size(L)
+ D = Δ(ny+nu, δ)
+ if upper
+ D = rand(D, N)
+ else
+ D = Diagonal([Interval(d) for d in diag(D)])
+ end
+ n = L.ny
+ X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
+ M = starprod(X,L)
+ # M = feedback(P, K; W2, Z2, Z1, W1) # TODO: this is probably not correct
+ M0 = freqresp(M, w)
+ M0 = permutedims(M0, (2,3,1))
+ M0, D
+end
 
 # function muloss(M, d)
 # D = Diagonal(d)
@@ -110,14 +137,15 @@ function structured_singular_value(M::AbstractArray{T}; tol=1e-4) where T
  ldiv!(Ms2, D, Ms1)
  opnorm(Ms2)
  end
- d0 = ones(real(T), size(M, 2))
+ n = size(M, 2)
+ d0 = ones(real(T), n)
  mu = map(axes(M, 3)) do i
  @views M0 = M[:,:,i]
  res = Optim.optimize(
  d->muloss(M0,d),
  d0,
  # i == 1 ? ParticleSwarm() : BFGS(alphaguess = LineSearches.InitialStatic(alpha=0.9), linesearch = LineSearches.HagerZhang()),
- i == 1 ? ParticleSwarm() : NelderMead(), # Initialize using Particle Swarm
+ n == 1 ? BFGS() : i == 1 ? ParticleSwarm() : NelderMead(), # Initialize using Particle Swarm
  Optim.Options(
  store_trace = false,
  show_trace = false,
@@ -139,6 +167,30 @@ function structured_singular_value(M::AbstractArray{T}; tol=1e-4) where T
  M isa AbstractMatrix ? mu[] : mu
 end
 
+"""
+ robstab(M0::UncertainSS, w=exp10.(LinRange(-3, 3, 1500)); kwargs...)
+
+Return the robust stability margin of an uncertain model, defined as the inverse of the structured singular value.
+"""
+robstab(M0::UncertainSS, args...; kwargs...) = 1/norm(structured_singular_value(M0, args...; kwargs...), Inf)
+
+"""
+ structured_singular_value(M0::UncertainSS, [w::AbstractVector]; kwargs...)
+
+- `w`: Frequency vector, if none is provided, the maximum μ over a brid 1e-3 : 1e3 will be returned.
+"""
+function structured_singular_value(M0::UncertainSS, w::AbstractVector; kwargs...)
+ M = freqresp(M0.M, w)
+ M = permutedims(M, (2,3,1))
+ μ = structured_singular_value(M; kwargs...)
+end
+
+function structured_singular_value(M0::UncertainSS; kwargs...)
+ w = exp10.(LinRange(-3, 3, 1500))
+ μ = structured_singular_value(M0, w; kwargs...)
+ maximum(μ)
+end
+
 """
  sim_diskmargin(L, σ::Real, w::AbstractVector)
  sim_diskmargin(L, σ::Real = 0)
@@ -147,11 +199,11 @@ Simultaneuous diskmargin at the outputs of `L`.
 Uses should consider using [`diskmargin`](@ref).
 """
 function sim_diskmargin(L::LTISystem,σ::Real,w::AbstractVector)
- # X = S+(σ-1)/2*I = lft([(1+σ)/2 -1;1 -1], L)
+ # S̄ = S+(σ-1)/2*I = lft([(1+σ)/2 -1;1 -1], L)
  n = L.ny
  X = ss(kron([(1+σ)/2 -1;1 -1], I(n)), L.timeevol)
- S̄ = starprod(X,L)
- M0 = permutedims(freqresp(S̄, w), (2,3,1))
+ M = starprod(X,L)
+ M0 = permutedims(freqresp(M, w), (2,3,1))
  mu = structured_singular_value(M0)
  imu = inv.(structured_singular_value(M0))
  simultaneous = [Diskmargin(imu, σ; ω0 = w, L) for (imu, w) in zip(imu,w)]

src/plotting.jl

@@ -110,14 +110,16 @@ nyquistplot
  θ = range(0, stop=2π, length=100)
  S, C = sin.(θ), cos.(θ)
  L = freqresp(sys, w)
- dets = mapslices(L, dims=(2,3)) do L
- det(I + L)
+ dets = map(axes(L, 1)) do i
+ det(I + L[i,:,:])
  end
  dets = vec(dets)
  redata = real.(dets)
  imdata = imag.(dets)
- ylims --> (clamp(minimum(imdata), -10, -1), clamp(maximum(imdata), 1, 10))
- xlims --> (clamp(minimum(redata), -10, -1), clamp(maximum(redata), 1, 10))
+ if eltype(imdata) <: AbstractFloat
+ ylims --> (clamp(minimum(imdata), -10, -1), clamp(maximum(imdata), 1, 10))
+ xlims --> (clamp(minimum(redata), -10, -1), clamp(maximum(redata), 1, 10))
+ end
  title --> "Mutivariable Nyquist plot"
 
  @series begin

src/uncertainty_interface.jl

@@ -128,8 +128,8 @@ function Base.getproperty(sys::UncertainSS, s::Symbol)
  elseif s === :uinds
  return sys.nw .+ (1:size(sys.B2, 2))
  elseif s ===:M
- # return sminreal(performance_mapping(sys.sys))
- return sys.sys
+ return sminreal(performance_mapping(sys.sys))
+ # return sys.sys
  elseif s ===:delta
  return sys.Δ
  end
@@ -168,19 +168,33 @@ Base.:+(n::Number, sys::UncertainSS{TE}) where TE <: ControlSystems.TimeEvolutio
 Base.:*(G::LTISystem, d::UncertainElement) = uss(G) * uss(d)
 Base.:*(d::UncertainElement, G::LTISystem) = uss(d) * uss(G)
 
+"""
+ uss(D11, D12, D21, D22, Δ, Ts = nothing)
 
+Create an uncertain statespace object with only gin matrices.
+"""
 function uss(D11, D12, D21, D22, Δ, Ts=nothing)
  D11, D12, D21, D22 = vcat.((D11, D12, D21, D22))
  sys = ExtendedStateSpace(zeros(0,0), zeros(0,size(D21,2)), zeros(0,size(D12,2)), zeros(size(D12,1),0), zeros(size(D21,1),0), D11, D12, D21, D22, Ts)
  # B1 B2 C1 C2
  UncertainSS(sys, Δ)
 end
 
+"""
+ uss(d::δ{C, F}, Ts = nothing)
+
+Convert a δ object to an UncertainSS
+"""
 function uss(d::δ{C,F}, Ts = nothing) where {C,F}
  # sys = partition(ss([d.val d.radius; 1 0]), u=2, y=2, w=1, z=1)
  uss(0, 1, d.radius, d.val, [normalize(d)], Ts)
 end
 
+"""
+ uss(d::AbstractVector{<:δ}, Ts = nothing)
+
+Create a diagonal uncertain statespace object with the uncertain elements `d` on the diagonal.
+"""
 function uss(d::AbstractVector{<:δ}, Ts = nothing)
  vals = getfield.(d, :val)
  rads = getfield.(d, :radius)
@@ -190,6 +204,11 @@ end
 
 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`.
+"""
 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)