Interface improvements + discrete in diskmargin

src/diskmargin.jl

@@ -13,7 +13,7 @@ The notation follows "An Introduction to Disk Margins", Peter Seiler, Andrew Pac
 `ϕm`: is analogous to the classical phase margin.
 `σ`: The skew parameter that was used to calculate the margin
 
-Note, `γmax` and `ϕm` are in smaller than the classical gain and phase margins sicne the classical margins do not consider simultaneous perturbations in gain and phase. 
+Note, `γmax` and `ϕm` are smaller than the classical gain and phase margins since the classical margins do not consider simultaneous perturbations in gain and phase.
 
 The "disk" margin becomes a half plane for `α = 2` and an inverted circle for `α > 2`. In this case, the upper gain margin is infinite. See the paper for more details, in particular figure 6.
 """
@@ -29,7 +29,7 @@ struct Diskmargin
     L
 end
 
-function Diskmargin(α, σ=0; ω0=mising, f0=missing, δ0=missing, L=nothing)
+function Diskmargin(α, σ=0; ω0=missing, f0=missing, δ0=missing, L=nothing)
     d = Disk(; α, σ)
     γmin = d.γmin
     γmax = d.γmax
@@ -145,7 +145,7 @@ The implementation and notation follows
 ["An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet](https://arxiv.org/abs/2003.04771).
 
 
-The margins are aviable as fields of the returned objects, see [`Diskmargin`](@ref).
+The margins are available as fields of the returned objects, see [`Diskmargin`](@ref).
 
 # Arguments:
 - `L`: A loop-transfer function.
@@ -188,7 +188,7 @@ See also [`ncfmargin`](@ref) and [`loop_diskmargin`](@ref).
 """
 function diskmargin(L::LTISystem, σ::Real=0; l=1e-3, u=1e3, kwargs...)
     L isa DelayLtiSystem && @warn "To compute the diskmargin of delay systems, consider approximating the delay with a Pade-approximation (`pade`) or discretize the system (`c2d`)."
-    issiso(L) || return sim_diskmargin(L, σ, l, u)
+    issiso(L) || return sim_diskmargin(L, σ, l, u; kwargs...)
     M = feedback(1, L) + (σ-1)/2
     n,ω0 = hinfnorm2(M; kwargs...)
     diskmargin(L, σ, ω0)
@@ -200,11 +200,11 @@ end
 Calculate the diskmargin at a particular frequency or vector of frequencies. If `ω` is a vector, you get a frequency-dependent diskmargin plot if you plot the returned value.
 See also [`ncfmargin`](@ref).
 """
-diskmargin(L::LTISystem, σ::Real, ω::AbstractArray) = map(w->diskmargin(L, σ, w), ω)
+diskmargin(L::LTISystem, σ::Real, ω::AbstractArray; kwargs...) = map(w->diskmargin(L, σ, w; kwargs...), ω)
 
-function diskmargin(L::LTISystem, σ::Real, ω0::Real)
+function diskmargin(L::LTISystem, σ::Real, ω0::Real; kwargs...)
     L isa DelayLtiSystem && @warn "To compute the diskmargin of delay systems, consider approximating the delay with a Pade-approximation (`pade`) or discretize the system (`c2d`)."
-    issiso(L) || return sim_diskmargin(L, σ, [ω0])[]
+    issiso(L) || return sim_diskmargin(L, σ, [ω0]; kwargs...)[]
     M = feedback(1, L) + (σ-1)/2
     freq = isdiscrete(L) ? cis(ω0*L.Ts) : complex(0, ω0)
     Sω = evalfr(M, freq)[]
@@ -428,7 +428,7 @@ end
     gainphaseplot(P)
     gainphaseplot(P, re, im)
 
-Plot complex perturbantions to the plant `P` and indicate whether or not the closed-loop system is stable. The diskmargin is the largest disk that can be fit inside the green region that only contains stable variations.
+Plot complex perturbations to the plant `P` and indicate whether or not the closed-loop system is stable. The diskmargin is the largest disk that can be fit inside the green region that only contains stable variations.
 """
 gainphaseplot
 

src/mimo_diskmargin.jl

@@ -194,7 +194,7 @@ function structured_singular_value(M::AbstractArray{T}; tol=1e-4, scalings=false
         return scalings ? (res.minimum, d0) : res.minimum
     else
         if scalings
-            Dm = Matrix{Float64}(undef, n, size(M,3))
+            Dm = Matrix{real(T)}(undef, n, size(M,3))
         end
         mu = map(axes(M, 3)) do i
             @views M0 = M[:,:,i]
@@ -256,30 +256,31 @@ end
     sim_diskmargin(L, σ::Real, w::AbstractVector)
     sim_diskmargin(L, σ::Real = 0)
 
-Simultaneuous diskmargin at the outputs of `L`. 
+Simultaneous diskmargin at the outputs of `L`. 
 Users should consider using [`diskmargin`](@ref).
 """
-function sim_diskmargin(L::LTISystem,σ::Real,w::AbstractVector)
+function sim_diskmargin(L::LTISystem, σ::Real, w::AbstractVector; kwargs...)
     # 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)
     M = starprod(X,L)
     M0 = freqresp(M, w)
-    imu = inv.(structured_singular_value(M0))
+    imu = inv.(structured_singular_value(M0; kwargs...))
     [Diskmargin(imu, σ; ω0 = w, L) for (imu, w) in zip(imu,w)]
 end
 
 """
     sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real = 0)
 
-Simultaneuous diskmargin at both outputs and inputs of `P`.
+Simultaneous diskmargin at both outputs and inputs of `P`.
 Ref: "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet
 https://arxiv.org/abs/2003.04771
 See also [`ncfmargin`](@ref).
 """
-function sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real=0, args...)
-    L = [ss(zeros(P.ny, P.ny)) P;-C ss(zeros(C.ny, C.ny))]
-    sim_diskmargin(L,σ, args...)
+function sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real=0, args...; kwargs...)
+    te = P.timeevol
+    L = [ss(zeros(P.ny, P.ny), te) P;-C ss(zeros(C.ny, C.ny), te)]
+    sim_diskmargin(L, σ, args...; kwargs...)
 end
 
 """
@@ -288,16 +289,16 @@ end
 Return the smallest simultaneous diskmargin over the grid l:u
 See also [`ncfmargin`](@ref).
 """
-function sim_diskmargin(L, σ::Real=0, l::Real=1e-3, u::Real=1e3)
-    m = sim_diskmargin(L, σ, exp10.(LinRange(log10(l), log10(u), 500)))
+function sim_diskmargin(L, σ::Real=0, l::Real=1e-3, u::Real=1e3; kwargs...)
+    m = sim_diskmargin(L, σ, exp10.(LinRange(log10(l), log10(u), 500)); kwargs...)
     m = argmin(d->d.α, m)
 end
 
 
 """
     diskmargin(P::LTISystem, C::LTISystem, σ, w::AbstractVector, args...; kwargs...)
 
-Simultaneuous diskmargin at outputs, inputs and input/output simultaneously of `P`. 
+Simultaneous diskmargin at outputs, inputs and input/output simultaneously of `P`. 
 Returns a named tuple with the fields `input, output, simultaneous_input, simultaneous_output, simultaneous` where `input` and `output` represent loop-at-a-time margins, `simultaneous_input` is the margin for simultaneous perturbations on all inputs and `simultaneous` is the margin for perturbations on all inputs and outputs simultaneously.
 
 Note: simultaneous margins are more conservative than single-loop margins and are likely to be much lower than the single-loop margins. Indeed, with several simultaneous perturbations, it's in general easier to make the system unstable. It's not uncommon for a simultaneous margin involving two signals to be on the order of half the size of the single-loop margins.
@@ -355,20 +356,18 @@ end
 Closes all loops in square MIMO system `L` except for loops `i`.
 Forms L1 in fig 14. of ["An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet](https://arxiv.org/abs/2003.04771)
 """
-function broken_feedback(L::LTISystem, i)
+function broken_feedback(L::LTISystem, i::Integer)
     ny, nu = size(L)
     ny == nu || throw(ArgumentError("Only square loop-transfer functions supported"))
     connection_inds = setdiff(1:ny, i)
-    i isa AbstractVector || (i = [i])
     open_L = feedback(
         L,
         ss(I(ny-1), L.timeevol), # open one loop
         U1 = connection_inds,
         Y1 = connection_inds,
-        Z1 = i,
-        W1 = i,
+        Z1 = [i],
+        W1 = [i],
     )
-    @assert issiso(open_L)
     open_L
 end
 

test/test_diskmargin.jl

@@ -99,6 +99,31 @@ plot!(dm.simultaneous_output)
 plot(dm.input)
 plot!(dm.output)
 
+# Discrete-time sim_diskmargin(P, C, ...) — covers timeevol fix.
+# Compute the simultaneous I/O margin for a continuous-time (P, C) pair and
+# for a fast-sampled discretization of the same pair; the margins should be
+# close.
+let
+    a = 10
+    Pc = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
+    Kc = ss(1.0I(2))
+    dm_c = sim_diskmargin(Pc, Kc)
+    Ts = 1e-3
+    Pd = c2d(Pc, Ts)
+    Kd = ss(1.0I(2), Ts)
+    dm_d = sim_diskmargin(Pd, Kd)
+    @test dm_d isa Diskmargin
+    @test dm_d.α ≈ dm_c.α rtol=1e-2
+
+    wgrid = exp10.(LinRange(-1, 2, 50))
+    dms_c = sim_diskmargin(Pc, Kc, 0, wgrid)
+    dms_d = sim_diskmargin(Pd, Kd, 0, wgrid)
+    @test dms_d isa AbstractVector{<:Diskmargin}
+    αs_c = [d.α for d in dms_c]
+    αs_d = [d.α for d in dms_d]
+    @test maximum(abs, αs_c .- αs_d) < 0.02
+end
+
 ## MIMO
 
 a = 10