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