From dcc7f999cb906e9233ff1b313723f23d78592539 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 26 Jan 2022 08:41:15 +0100 Subject: [PATCH 1/3] improve rlocusplot closes #611 --- src/pid_design.jl | 47 ++++++++++++++++++++++++++++++----------------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/src/pid_design.jl b/src/pid_design.jl index 64be02c3d..216431ae0 100644 --- a/src/pid_design.jl +++ b/src/pid_design.jl @@ -84,17 +84,27 @@ end -function getpoles(G, K) # If OrdinaryDiffEq is installed, we override getpoles with an adaptive method - P = numpoly(G)[1] - Q = denpoly(G)[1] - f = (y,_,k) -> sort(ComplexF64.(Polynomials.roots(k[1]*P+Q)), by=imag) - prob = OrdinaryDiffEq.ODEProblem(f,f(0.,0.,0.),(0.,K[end])) +function getpoles(G, K) + P = numpoly(G)[1] + Q = denpoly(G)[1] + T = float(eltype(K)) + ϵ = eps(T) + f = function (y,_,k) + if k == 0 && length(P) > length(Q) + # More zeros than poles, make sure the vector of roots is of correct length when k = 0 + # When this happens, there are fewer poles for k = 0, these poles can be seen as beeing located somewhere at Inf + # We get around the problem by not allowing k = 0 for non-proper systems. + k = ϵ + end + sort(ComplexF64.(Polynomials.roots(k[1]*P+Q)), by=imag) + end + prob = OrdinaryDiffEq.ODEProblem(f,f(0.,0.,0),(0,K[end])) integrator = OrdinaryDiffEq.init(prob,OrdinaryDiffEq.Tsit5(),reltol=1e-8,abstol=1e-8) ts = Vector{Float64}() poleout = Vector{Vector{ComplexF64}}() for i in integrator - push!(poleout,integrator.k[1]) - push!(ts,integrator.t[1]) + push!(poleout,integrator.k[1]) + push!(ts,integrator.t[1]) end poleout = hcat(poleout...)' poleout, ts @@ -111,16 +121,21 @@ range(1e-6,stop=500,length=10000) is used. If `OrdinaryDiffEq.jl` is installed and loaded by the user (`using OrdinaryDiffEq`), `rlocusplot` will use an adaptive step-size algorithm to select values of `K`. A scalar `Kmax` can then be given as second argument. """ -rlocus -@recipe function rlocus(p::Rlocusplot; K=500) +rlocusplot +@recipe function rlocusplot(p::Rlocusplot; K=500) P = p.args[1] K = K isa Number ? range(1e-6,stop=K,length=10000) : K Z = tzeros(P) - poles, K = getpoles(P,K) - redata = real.(poles) - imdata = imag.(poles) - ylim = (max(-50,minimum(imdata)), min(50,maximum(imdata))) - xlim = (max(-50,minimum(redata)), min(50,maximum(redata))) + roots, K = getpoles(P,K) + redata = real.(roots) + imdata = imag.(roots) + + openloop_features = [poles(P); tzeros(P)] + + reth = 1.2 * maximum(abs, real(openloop_features)) # A nice threshold for limits + imth = 1.2 * maximum(abs, imag(openloop_features)) # A nice threshold for limits + ylims --> (max(-imth,minimum(imdata)-1), min(imth,maximum(imdata)+1)) + xlims --> (max(-reth,minimum(redata)-1), min(reth,maximum(redata)+1)) framestyle --> :zerolines title --> "Root locus" xguide --> "Re(roots)" @@ -128,9 +143,7 @@ rlocus form(k, p) = Printf.@sprintf("%.4f", k) * " pole=" * Printf.@sprintf("%.3f%+.3fim", real(p), imag(p)) @series begin legend --> false - ylims --> ylim - xlims --> xlim - hover := "K=" .* form.(K,poles) + hover := "K=" .* form.(K,roots) label := "" redata, imdata end From 0e713dace6a483273eadde2c2307e919882fbd38 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 26 Jan 2022 16:53:36 +0100 Subject: [PATCH 2/3] tweak default limits --- src/pid_design.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/pid_design.jl b/src/pid_design.jl index 216431ae0..81a819ed3 100644 --- a/src/pid_design.jl +++ b/src/pid_design.jl @@ -129,13 +129,9 @@ rlocusplot roots, K = getpoles(P,K) redata = real.(roots) imdata = imag.(roots) - - openloop_features = [poles(P); tzeros(P)] - reth = 1.2 * maximum(abs, real(openloop_features)) # A nice threshold for limits - imth = 1.2 * maximum(abs, imag(openloop_features)) # A nice threshold for limits - ylims --> (max(-imth,minimum(imdata)-1), min(imth,maximum(imdata)+1)) - xlims --> (max(-reth,minimum(redata)-1), min(reth,maximum(redata)+1)) + ylims --> (max(-50,minimum(imdata) - 1), min(50,maximum(imdata) + 1)) + xlims --> (max(-50,minimum(redata) - 1), clamp(maximum(redata) + 1, 1, 50)) framestyle --> :zerolines title --> "Root locus" xguide --> "Re(roots)" From d059caad0560c038b61965b8c876a6e9b90ad50c Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 26 Jan 2022 16:55:02 +0100 Subject: [PATCH 3/3] convert ss to tf in rlocusplot --- src/pid_design.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pid_design.jl b/src/pid_design.jl index 81a819ed3..aa8516f7d 100644 --- a/src/pid_design.jl +++ b/src/pid_design.jl @@ -83,8 +83,8 @@ end @deprecate rlocus(args...;kwargs...) rlocusplot(args...;kwargs...) - function getpoles(G, K) + G isa TransferFunction || (G = tf(G)) P = numpoly(G)[1] Q = denpoly(G)[1] T = float(eltype(K))