From dffd6167e22fe0cdd6e583d52f7171cb02c97ab9 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 24 Jun 2026 08:47:25 +0200 Subject: [PATCH] diskmargin: derive MIMO frequency-grid bounds from the system by default MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The MIMO diskmargin searched a fixed grid l=1e-3..u=1e3 rad/s. When the worst-case margin lay near a bound, the reported worst-case frequency pinned to the grid edge (e.g. exactly 1000.0 rad/s), corrupting derived quantities such as the delay margin (phase_margin/ω_crit). Default l and u to `nothing`; when unset, derive them from the poles/zeros of L via ControlSystemsBase._bounds_and_features(L, Val{:nyquist}()), the same heuristic used for Bode/Nyquist plots. Explicit l/u still override. Co-Authored-By: Claude Opus 4.8 (1M context) --- Project.toml | 2 +- src/diskmargin.jl | 3 ++- src/mimo_diskmargin.jl | 15 +++++++++++---- test/test_diskmargin.jl | 5 ++++- 4 files changed, 18 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 29388928..9b4a38aa 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.4.50" +version = "0.4.51" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/diskmargin.jl b/src/diskmargin.jl index b5eec87d..0d40d710 100644 --- a/src/diskmargin.jl +++ b/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 diff --git a/src/mimo_diskmargin.jl b/src/mimo_diskmargin.jl index db1a4af1..2f344af1 100644 --- a/src/mimo_diskmargin.jl +++ b/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 diff --git a/test/test_diskmargin.jl b/test/test_diskmargin.jl index 6286d0f3..bf538c17 100644 --- a/test/test_diskmargin.jl +++ b/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