Merge pull request #915 from JuliaControl/margifix

handle rounding error in Nyquist freq for default freq vec
JuliaControl · Feb 6, 2024 · d5bf7c6 · d5bf7c6
2 parents 01c8e08 + d4638ea
commit d5bf7c6
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 12 deletions.
diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl
@@ -410,18 +410,19 @@ See also [`delaymargin`](@ref) and [`RobustAndOptimalControl.diskmargin`](https:
 function margin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false)
     ny, nu = size(sys)
 
+    T = float(numeric_type(sys))
     if allMargins
-        wgm         = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
-        gm          = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
-        wpm         = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
-        pm          = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
-        fullPhase   = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
+        wgm         = Array{Array{T,1}}(undef, ny,nu)
+        gm          = Array{Array{T,1}}(undef, ny,nu)
+        wpm         = Array{Array{T,1}}(undef, ny,nu)
+        pm          = Array{Array{T,1}}(undef, ny,nu)
+        fullPhase   = Array{Array{T,1}}(undef, ny,nu)
     else
-        wgm         = Array{numeric_type(sys),2}(undef, ny, nu)
-        gm          = Array{numeric_type(sys),2}(undef, ny, nu)
-        wpm         = Array{numeric_type(sys),2}(undef, ny, nu)
-        pm          = Array{numeric_type(sys),2}(undef, ny, nu)
-        fullPhase   = Array{numeric_type(sys),2}(undef, ny, nu)
+        wgm         = Array{T,2}(undef, ny, nu)
+        gm          = Array{T,2}(undef, ny, nu)
+        wpm         = Array{T,2}(undef, ny, nu)
+        pm          = Array{T,2}(undef, ny, nu)
+        fullPhase   = Array{T,2}(undef, ny, nu)
     end
     for j=1:nu
         for i=1:ny

diff --git a/lib/ControlSystemsBase/src/freqresp.jl b/lib/ControlSystemsBase/src/freqresp.jl
@@ -382,7 +382,11 @@ function _default_freq_vector(systems::Vector{<:LTISystem}, plot)
     w2 = maximum(maximum, bounds)
 
     nw = round(Int, max(min_pt_total, min_pt_per_dec*(w2 - w1)))
-    return exp10.(range(w1, stop=w2, length=nw))
+    w = exp10.(range(w1, stop=w2, length=nw))
+    if length(systems) == 1 && isdiscrete(systems[1])
+        w[end] = π/systems[1].Ts # To account for numerical rounding problems from exp(log())
+    end
+    w
 end
 _default_freq_vector(sys::LTISystem, plot) = _default_freq_vector(
         [sys], plot)

diff --git a/lib/ControlSystemsBase/test/test_freqresp.jl b/lib/ControlSystemsBase/test/test_freqresp.jl
@@ -173,6 +173,9 @@ mag, mag, ws2 = bode(sys2)
 @test maximum(ws2) >= 5max(p,z)
 @test minimum(ws2) <= 0.2min(p,z)
 @test length(ws2) > 100
+
+@test margin(tf(1, [1, -1], 0.01)).gm == [2;;]
+
 end
 
 
@@ -204,4 +207,5 @@ end
 
 # f2 = plot(sizes, last.(times1), scale=:log10, lab="Allocations freqresp", m=:o)
 # plot!(sizes, last.(times2), scale=:log10, lab="Allocations freqresp_large", xlabel="Model order", m=:o)
-# plot(f1, f2)
+# plot(f1, f2)
+