From 559767ee885a5129f9313a2eb044865f4b4e4347 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 15 May 2026 19:42:38 +0200 Subject: [PATCH] =?UTF-8?q?Clamp=20default=20frequency=20grid=20to=20Nyqui?= =?UTF-8?q?st=20for=20discrete=20diskmargin/=CE=BC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `sim_diskmargin(L, σ, l, u)` and `structured_singular_value(M0::LTISystem)` both built their default frequency grid up to 1e3 rad/s. For discrete-time inputs whose Nyquist frequency π/Ts is below 1e3, that grid extended beyond Nyquist and `freqresp` aliased (`cis(ω·Ts)` is 2π/Ts-periodic), so the worst-case margin / μ computed on the default grid reflected an aliased frequency rather than the intended one. The MIMO branch of `diskmargin(L, σ; l, u)` is also fixed because it forwards to `sim_diskmargin`. Clamp the upper bound to `min(u, π/L.Ts)` (or `min(3.0, log10(π/M0.Ts))` in log-frequency) when the system is discrete. Continuous-time behavior is unchanged. Users who deliberately want frequencies past Nyquist can still pass an explicit `w` vector. Adds a regression test in `test/test_diskmargin.jl` using `Ts=1.0` (Nyquist ≈ 3.14 rad/s, well below the default 1e3) that confirms the default-grid result agrees with an explicit Nyquist-respecting grid and that the worst-case frequency is at or below Nyquist. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/mimo_diskmargin.jl | 8 +++++++- test/test_diskmargin.jl | 19 +++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/mimo_diskmargin.jl b/src/mimo_diskmargin.jl index e5be5865..db1a4af1 100644 --- a/src/mimo_diskmargin.jl +++ b/src/mimo_diskmargin.jl @@ -247,7 +247,9 @@ function structured_singular_value(M0::LTISystem, w::AbstractVector; kwargs...) end function structured_singular_value(M0::LTISystem; kwargs...) - w = exp10.(LinRange(-3, 3, 1500)) + # Clamp the upper end of the default grid to the Nyquist frequency for discrete `M0`. + u = isdiscrete(M0) ? min(3.0, log10(π/M0.Ts)) : 3.0 + w = exp10.(LinRange(-3, u, 1500)) μ = structured_singular_value(M0, w; kwargs...) maximum(μ) end @@ -290,6 +292,10 @@ 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; kwargs...) + # For discrete L the default upper bound 1e3 typically exceeds 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 m = sim_diskmargin(L, σ, exp10.(LinRange(log10(l), log10(u), 500)); kwargs...) m = argmin(d->d.α, m) end diff --git a/test/test_diskmargin.jl b/test/test_diskmargin.jl index a1dcb6aa..6286d0f3 100644 --- a/test/test_diskmargin.jl +++ b/test/test_diskmargin.jl @@ -124,6 +124,25 @@ let @test maximum(abs, αs_c .- αs_d) < 0.02 end +# Discrete-time sim_diskmargin: default frequency grid must respect Nyquist. +# Regression: with slow sampling (here Ts=1.0, Nyquist ≈ 3.14 rad/s) the default +# upper bound u=1e3 used to alias under freqresp. The default-grid result must +# agree with an explicit Nyquist-respecting grid. +let + a = 10 + Pc = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) + Ts = 1.0 + Pd = c2d(Pc, Ts) + Kd = ss(1.0I(2), Ts) + 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)) + dms_explicit = sim_diskmargin(Pd, Kd, 0, wgrid) + αs = [d.α for d in dms_explicit] + @test dm_default.α ≈ minimum(αs) rtol=1e-3 +end + ## MIMO a = 10