diff --git a/lib/ControlSystemsBase/src/analysis.jl b/lib/ControlSystemsBase/src/analysis.jl index d455da77c..bbf35e762 100644 --- a/lib/ControlSystemsBase/src/analysis.jl +++ b/lib/ControlSystemsBase/src/analysis.jl @@ -494,12 +494,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/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/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. diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl index 6a18ddfd0..c3f0ddd65 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) @@ -429,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 @@ -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) @@ -771,10 +779,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) @@ -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 diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index bf5e3a4a9..a225900a8 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/lib/ControlSystemsBase/test/test_analysis.jl @@ -278,6 +278,7 @@ nintbig = ControlSystemsBase.count_integrators(big(1.0)*pade(OL, 2)) @test ControlSystemsBase.count_integrators(pade(OL, 4)) == nintbig @test ControlSystemsBase.count_integrators(pade(OL, 10)) == nintbig + # RGA a = 10 P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0) 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