From d9736a4db3d37f16f0ec4fedbcfa8cf782d54748 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 13 Feb 2025 08:30:24 +0100 Subject: [PATCH 1/5] handle phase wrapping in margin --- lib/ControlSystemsBase/src/analysis.jl | 10 ++++++++-- lib/ControlSystemsBase/src/plotting.jl | 8 ++++---- lib/ControlSystemsBase/test/test_analysis.jl | 19 +++++++++++++++++++ 3 files changed, 31 insertions(+), 6 deletions(-) diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index 568303188..975659097 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/lib/ControlSystemsBase/src/analysis.jl @@ -496,12 +496,18 @@ function sisomargin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMa wgm, = _allPhaseCrossings(w, phase) gm = similar(wgm) for i = eachindex(wgm) - gm[i] = 1 ./ abs(freqresp(sys,[wgm[i]])[1][1]) + gm[i] = 1 ./ abs(freqresp(sys,wgm[i])[1]) end wpm, fi = _allGainCrossings(w, mag) pm = similar(wpm) for i = eachindex(wpm) - pm[i] = mod(rad2deg(angle(freqresp(sys,[wpm[i]])[1][1])),360)-180 + # We have to access the actual phase value from the `phase` array to get unwrapped phase. This value is not fully accurate since it is computed at a grid point, so we compute the more accurate phase at the interpolated frequency. This accurate value is not unwrapped, so we add an integer multiple of 360 to get the closest unwrapped phase. + φ_nom = rad2deg(angle(freqresp(sys,wpm[i])[1])) + φ_rounded = phase[clamp(round(Int, fi[i]), 1, length(phase))] # fi is interpolated, so we round to the closest integer + φ_int = φ_nom - 360 * round( (φ_nom - φ_rounded) / 360 ) + + # Now compute phase margin relative to -180: + pm[i] = 180 + φ_int end if !allMargins #Only output the smallest margins gm, idx = findmin([gm;Inf]) diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl index 6a18ddfd0..9d9dbd414 100644 --- a/lib/ControlSystemsBase/src/plotting.jl +++ b/lib/ControlSystemsBase/src/plotting.jl @@ -771,10 +771,10 @@ marginplot mag = 1 ./ gm oneLine = 1 end - titles[j,i,1,1] *= "["*join([Printf.@sprintf("%2.2f",v) for v in gm],", ")*"] " - titles[j,i,1,2] *= "["*join([Printf.@sprintf("%2.2f",v) for v in wgm],", ")*"] " - titles[j,i,2,1] *= "["*join([Printf.@sprintf("%2.2f",v) for v in pm],", ")*"] " - titles[j,i,2,2] *= "["*join([Printf.@sprintf("%2.2f",v) for v in wpm],", ")*"] " + titles[j,i,1,1] *= "["*join([Printf.@sprintf("%3.2g",v) for v in gm],", ")*"] " + titles[j,i,1,2] *= "["*join([Printf.@sprintf("%3.2g",v) for v in wgm],", ")*"] " + titles[j,i,2,1] *= "["*join([Printf.@sprintf("%3.2g",v) for v in pm],", ")*"] " + titles[j,i,2,2] *= "["*join([Printf.@sprintf("%3.2g",v) for v in wpm],", ")*"] " subplot := min(s2i((plotphase ? (2i-1) : i),j), prod(plotattributes[:layout])) if si == length(systems) diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index f21a1dbe1..aecf508c1 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/lib/ControlSystemsBase/test/test_analysis.jl @@ -250,6 +250,25 @@ P = tf(1,[5, 10.25, 6.25, 1]) w_180, gm, w_c, pm = margin(50P) @test pm[] ≈ -35.1 rtol=1e-2 +## Tricky case from https://www.reddit.com/r/ControlTheory/comments/1inhxsz/understanding_stability_in_highorder/ +s = tf("s") +kpu = -10.593216768722073; kiu = -0.00063; t = 1000; tau = 180; a = 1/8.3738067325406132e-5; +kpd = 15.92190277847431; kid = 0.000790960718241793; +kpo = -10.39321676872207317; kio = -0.00063; +kpb = kpd; kib = kid; + +C1 = (kpu + kiu/s)*(1/(t*s + 1)) +C2 = (kpu + kiu/s)*(1/(t*s + 1)) +C3 = (kpo + kio/s)*(1/(t*s + 1)) +Cb = (kpb + kib/s)*(1/(t*s + 1)) +OL = (ss(Cb)*ss(C1)*ss(C2)*ss(C3)*exp(-3*tau*s))/((C1 - a*s)*(C2 - a*s)*(C3 - a*s)); + +wgm, gm, ωϕₘ, ϕₘ = margin(OL; full=true, allMargins=true) +@test ϕₘ[][] ≈ -320 rtol=1e-2 +for wgm in wgm[] + @test mod(rad2deg(angle(freqresp(OL, wgm)[])), 360)-180 ≈ 0 atol=1e-1 +end + # RGA a = 10 P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) From f4d497d0c47a37873cfa0ef6697bd3d3cefdcebc Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 13 Feb 2025 08:30:46 +0100 Subject: [PATCH 2/5] handle method ambiguity ss/tf --- lib/ControlSystemsBase/src/connections.jl | 4 ++++ lib/ControlSystemsBase/test/test_connections.jl | 2 ++ 2 files changed, 6 insertions(+) diff --git a/lib/ControlSystemsBase/src/connections.jl b/lib/ControlSystemsBase/src/connections.jl index b2478a5bf..3434e2fae 100644 --- a/lib/ControlSystemsBase/src/connections.jl +++ b/lib/ControlSystemsBase/src/connections.jl @@ -215,6 +215,10 @@ function /(sys1::Union{StateSpace,AbstractStateSpace}, sys2::LTISystem) sys1new, sys2new = promote(sys1, 1/sys2) return sys1new*sys2new end +function /(sys1::Union{StateSpace,AbstractStateSpace}, sys2::TransferFunction) # This method is handling ambiguity between method above and one with explicit TF as second argument, hit by ss(1)/tf(1) + sys1new, sys2new = promote(sys1, 1/sys2) + return sys1new*sys2new +end @static if VERSION >= v"1.8.0-beta1" blockdiag(anything...) = cat(anything..., dims=Val((1,2))) diff --git a/lib/ControlSystemsBase/test/test_connections.jl b/lib/ControlSystemsBase/test/test_connections.jl index 7a34969df..d944386d9 100644 --- a/lib/ControlSystemsBase/test/test_connections.jl +++ b/lib/ControlSystemsBase/test/test_connections.jl @@ -438,4 +438,6 @@ Pr = input_resolvent(P) @test Pr.C == I @test iszero(Pr.D) +@test ss(1) / tf(1) == ss(1) # Test no method ambiguity + end \ No newline at end of file From 280ed4585edd986fa8512e58fee8daba15522956 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 13 Feb 2025 08:31:01 +0100 Subject: [PATCH 3/5] tidy up latex math --- lib/ControlSystemsBase/src/pid_design.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/ControlSystemsBase/src/pid_design.jl b/lib/ControlSystemsBase/src/pid_design.jl index af10e2226..1e223a73e 100644 --- a/lib/ControlSystemsBase/src/pid_design.jl +++ b/lib/ControlSystemsBase/src/pid_design.jl @@ -150,8 +150,8 @@ r ┌─────┐ ┌─────┐ │ │ │ ``` The `form` can be chosen as one of the following (determines how the arguments `param_p, param_i, param_d` are interpreted) -* `:standard` - ``K_p*(br-y + (r-y)/(T_i s) + T_d s (cr-y)/(T_f s + 1))`` -* `:parallel` - ``K_p*(br-y) + K_i (r-y)/s + K_d s (cr-y)/(Tf s + 1)`` +* `:standard` - ``K_p(br-y + (r-y)/(T_i s) + T_d s (cr-y)/(T_f s + 1))`` +* `:parallel` - ``K_p(br-y) + K_i (r-y)/s + K_d s (cr-y)/(Tf s + 1)`` - `b` is a set-point weighting for the proportional term - `c` is a set-point weighting for the derivative term, this defaults to 0. From 74d7e0d43bd4498a252188627e106230f13f24be Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 13 Feb 2025 09:49:10 +0100 Subject: [PATCH 4/5] adjust for integrator excess in marginplot --- lib/ControlSystemsBase/src/plotting.jl | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl index 9d9dbd414..76e3222f0 100644 --- a/lib/ControlSystemsBase/src/plotting.jl +++ b/lib/ControlSystemsBase/src/plotting.jl @@ -293,6 +293,9 @@ end else sbal = s end + if plotphase && adjust_phase_start && isrational(sbal) + intexcess = integrator_excess(sbal) + end mag, phase = bode(sbal, w; unwrap=false) if _PlotScale == "dB" # Set by setPlotScale(str) globally mag = 20*log10.(mag) @@ -330,7 +333,6 @@ end plotphase || continue if adjust_phase_start == true && isrational(sbal) - intexcess = integrator_excess(sbal) if intexcess != 0 # Snap phase so that it starts at -90*intexcess nineties = round(Int, phasedata[1] / 90) @@ -725,11 +727,12 @@ Plot all the amplitude and phase margins of the system(s) `sys`. - A frequency vector `w` can be optionally provided. - `balance`: Call [`balance_statespace`](@ref) on the system before plotting. +- `adjust_phase_start`: If true, the phase will be adjusted so that it starts at -90*intexcess degrees, where `intexcess` is the integrator excess of the system. `kwargs` is sent as argument to RecipesBase.plot. """ marginplot -@recipe function marginplot(p::Marginplot; plotphase=true, hz=false, balance=true) +@recipe function marginplot(p::Marginplot; plotphase=true, hz=false, balance=true, adjust_phase_start=true) systems, w = _processfreqplot(Val{:bode}(), p.args...) ny, nu = size(systems[1]) s2i(i,j) = LinearIndices((nu,(plotphase ? 2 : 1)*ny))[j,i] @@ -746,6 +749,11 @@ marginplot s = balance_statespace(s)[1] end bmag, bphase = bode(s, w) + + if plotphase && adjust_phase_start && isrational(s) + intexcess = integrator_excess(s) + end + for j=1:nu for i=1:ny wgm, gm, wpm, pm, fullPhase = sisomargin(s[i,j],w, full=true, allMargins=true) @@ -801,6 +809,15 @@ marginplot [wgm wgm]', [ones(length(mag)) mag]' end plotphase || continue + + phasedata = bphase[i, j, :] + if plotphase && adjust_phase_start && isrational(s) + if intexcess != 0 + # Snap phase so that it starts at -90*intexcess + nineties = round(Int, phasedata[1] / 90) + phasedata .+= ((90*(-intexcess-nineties)) ÷ 360) * 360 + end + end # Phase margins subplot := s2i(2i,j) @@ -810,7 +827,7 @@ marginplot @series begin primary := true seriestype := :bodephase - w, bphase[i, j, :] + w, phasedata end @series begin color --> :gray From a3d9b57c085eb9102e1443720692459c43be960f Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 13 Feb 2025 09:49:42 +0100 Subject: [PATCH 5/5] better hover labels in nyquistplot --- lib/ControlSystemsBase/src/plotting.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl index 76e3222f0..c3f0ddd65 100644 --- a/lib/ControlSystemsBase/src/plotting.jl +++ b/lib/ControlSystemsBase/src/plotting.jl @@ -431,7 +431,7 @@ nyquistplot if lab !== nothing label --> lab end - hover --> [hz ? Printf.@sprintf("f = %.3f", w/2π) : Printf.@sprintf("ω = %.3f", w) for w in w] + hover --> [hz ? Printf.@sprintf("f = %.3g", w/2π) : Printf.@sprintf("ω = %.3g", w) for w in w] (redata, imdata) end