diskmargin: auto-compute MIMO frequency-grid bounds by default

Project.toml

@@ -1,7 +1,7 @@
 name = "RobustAndOptimalControl"
 uuid = "21fd56a4-db03-40ee-82ee-a87907bee541"
 authors = ["Fredrik Bagge Carlson", "Marcus Greiff"]
-version = "0.4.50"
+version = "0.4.51"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/diskmargin.jl

@@ -152,6 +152,7 @@ The margins are available as fields of the returned objects, see [`Diskmargin`](
 - `σ`: If little is known about the distribution of gain variations then σ = 0 is a reasonable choice as it allows for a gain increase or decrease by the same relative amount. *The choice σ < 0* is justified if the gain can decrease by a larger factor than it can increase. Similarly, *the choice σ > 0* is justified when the gain can increase by a larger factor than it can decrease. *If σ = −1* then the disk margin condition is αmax = inv(MT). This margin is related to the robust stability condition for models with multiplicative uncertainty of the form P (1 + δ). If σ = +1 then the disk margin condition is αmax = inv(MS)
 - `kwargs`: Are sent to the [`hinfnorm`](@ref) calculation
 - `ω`: If a vector of frequencies is supplied, the frequency-dependent disk margin will be computed, see example below.
+- `l, u`: Lower and upper frequency bounds [rad/s] for the search grid used for MIMO systems (ignored for SISO systems, which are solved analytically). The default `nothing` derives the bounds automatically from the dynamics of `L`; pass explicit values to override. If the reported worst-case frequency lands exactly on `u`, the true worst case may lie outside the grid — extend the range.
 
 # Example: 
 ```
@@ -186,7 +187,7 @@ isapprox(dm.ω0, wMt, rtol=1e-1)
 
 See also [`ncfmargin`](@ref) and [`loop_diskmargin`](@ref).
 """
-function diskmargin(L::LTISystem, σ::Real=0; l=1e-3, u=1e3, kwargs...)
+function diskmargin(L::LTISystem, σ::Real=0; l=nothing, u=nothing, 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; kwargs...)
     M = feedback(1, L) + (σ-1)/2

src/mimo_diskmargin.jl

@@ -286,13 +286,20 @@ function sim_diskmargin(P::LTISystem, C::LTISystem, σ::Real=0, args...; kwargs.
 end
 
 """
-    sim_diskmargin(L, σ::Real = 0, l=1e-3, u=1e3)
+    sim_diskmargin(L, σ::Real = 0, l=nothing, u=nothing)
 
-Return the smallest simultaneous diskmargin over the grid l:u
+Return the smallest simultaneous diskmargin over a log-spaced frequency grid spanning `l` to `u` rad/s.
+If `l` or `u` is `nothing` (the default), the corresponding bound is derived from the dynamics of `L`
+(poles and zeros) using the same heuristic as the frequency-domain plots, with high-frequency headroom.
 See also [`ncfmargin`](@ref).
 """
-function sim_diskmargin(L, σ::Real=0, l::Real=1e-3, u::Real=1e3; kwargs...)
-    # For discrete L the default upper bound 1e3 typically exceeds the Nyquist frequency π/Ts.
+function sim_diskmargin(L, σ::Real=0, l::Union{Real,Nothing}=nothing, u::Union{Real,Nothing}=nothing; kwargs...)
+    if l === nothing || u === nothing
+        b = ControlSystemsBase._bounds_and_features(L, Val{:nyquist}())[1] # log10 frequency bounds [rad/s]
+        l === nothing && (l = exp10(b[1]))
+        u === nothing && (u = exp10(b[2]))
+    end
+    # For discrete L the upper bound may exceed the Nyquist frequency π/Ts.
     # Above Nyquist `freqresp` aliases (cis(ω·Ts) is 2π/Ts-periodic) so margins computed there
     # reflect the response at an aliased frequency rather than the intended one.
     u = isdiscrete(L) ? min(u, π/L.Ts) : u

test/test_diskmargin.jl

@@ -137,7 +137,10 @@ let
     dm_default = sim_diskmargin(Pd, Kd)
     @test dm_default isa Diskmargin
     @test dm_default.ω0 ≤ π/Ts + sqrt(eps())
-    wgrid = exp10.(LinRange(log10(1e-3), log10(π/Ts), 500))
+    # Dense reference grid up to Nyquist; the default grid now derives its lower
+    # bound from the dynamics, so use a fine explicit grid to make the comparison
+    # robust to the exact bounds.
+    wgrid = exp10.(LinRange(log10(1e-5*π/Ts), log10(π/Ts), 4000))
     dms_explicit = sim_diskmargin(Pd, Kd, 0, wgrid)
     αs = [d.α for d in dms_explicit]
     @test dm_default.α ≈ minimum(αs) rtol=1e-3