diff --git a/Project.toml b/Project.toml index 924b04d6d..64c8d689b 100644 --- a/Project.toml +++ b/Project.toml @@ -16,10 +16,10 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4" MatrixPencils = "48965c70-4690-11ea-1f13-43a2532b2fa8" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] @@ -33,14 +33,15 @@ MacroTools = "0.5" MatrixEquations = "1, 2.1" MatrixPencils = "1.6" OrdinaryDiffEq = "5.2, 6.0" -Plots = "0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 1.0" Polynomials = "1.1.10, 2.0" +RecipesBase = "1" julia = "1.6" [extras] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Documenter", "GR"] +test = ["Test", "Documenter", "GR", "Plots"] diff --git a/README.md b/README.md index d4cec0e98..08ba48ca4 100644 --- a/README.md +++ b/README.md @@ -103,6 +103,7 @@ CLs = TransferFunction[kp*P/(1 + kp*P) for kp = [1, 5, 15]]; # Plot the step response of the controllers # Any keyword arguments supported in Plots.jl can be supplied +using Plots plot(step.(CLs, 5), label=["Kp = 1" "Kp = 5" "Kp = 15"]) ``` diff --git a/docs/src/examples/example.md b/docs/src/examples/example.md index 2b66e77a1..46892d222 100644 --- a/docs/src/examples/example.md +++ b/docs/src/examples/example.md @@ -13,6 +13,7 @@ end # LQR design ```jldoctest; output = false using LinearAlgebra # For identity matrix I +using Plots Ts = 0.1 A = [1 Ts; 0 1] B = [0 1]' # To handle bug TODO @@ -26,7 +27,7 @@ u(x,t) = -L*x .+ 1.5(t>=2.5)# Form control law (u is a function of t and x), a t =0:Ts:5 x0 = [1,0] y, t, x, uout = lsim(sys,u,t,x0=x0) -Plots.plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]") +plot(t,x', lab=["Position" "Velocity"], xlabel="Time [s]") save_docs_plot("lqrplot.svg"); # hide @@ -53,13 +54,14 @@ we notice that the sensitivity function is a bit too high around frequencies ω function `loopshapingPI` and tell it that we want 60 degrees phase margin at this frequency. The resulting gang of four is plotted for both the constructed controller and for unit feedback. ```jldoctest PIDDESIGN; output = false +using Plots ωp = 0.8 kp,ki,C = loopshapingPI(P,ωp,phasemargin=60) p1 = gangoffourplot(P, [tf(1), C]); p2 = nyquistplot([P, P*C], ylims=(-1,1), xlims=(-1.5,1.5)); -Plots.plot(p1,p2, layout=(2,1), size=(800,800)) +plot(p1,p2, layout=(2,1), size=(800,800)) # save_docs_plot("pidgofplot2.svg") # hide # save_docs_plot("pidnyquistplot.svg"); # hide save_docs_plot("pidgofnyquistplot.svg") # hide @@ -71,13 +73,14 @@ save_docs_plot("pidgofnyquistplot.svg") # hide We could also cosider a situation where we want to create a closed-loop system with the bandwidth ω = 2 rad/s, in which case we would write something like ```jldoctest PIDDESIGN; output = false +using Plots ωp = 2 kp,ki,C60 = loopshapingPI(P,ωp,rl=1,phasemargin=60, doplot=true) p1 = gangoffourplot(P, [tf(1), C60]); p2 = nyquistplot([P, P*C60], ylims=(-2,2), xlims=(-3,3)); -Plots.plot(p1,p2, layout=(2,1), size=(800,800)) +plot(p1,p2, layout=(2,1), size=(800,800)) # gangoffourplot(P, [tf(1), C60]) # hide # save_docs_plot("pidgofplot3.svg") # hide diff --git a/docs/src/man/introduction.md b/docs/src/man/introduction.md index ed57bac97..f52876b20 100644 --- a/docs/src/man/introduction.md +++ b/docs/src/man/introduction.md @@ -49,12 +49,13 @@ TransferFunction{Continuous, ControlSystems.SisoRational{Float64}} Continuous-time transfer function model ``` ## Plotting -Plotting requires some extra care. The ControlSystems package is using `Plots.jl` ([link](https://github.com/tbreloff/Plots.jl)) as interface to generate all the plots. This means that the user is able to freely choose back-end. The plots in this manual are generated using `GR`. If you have several back-ends for plotting then you can select the one you want to use with the corresponding `Plots` call (for `GR` this is `Plots.gr()`, some alternatives are `pyplot(), plotly(), pgfplots()`). A simple example where we generate a plot and save it to a file is +Plotting requires some extra care. The ControlSystems package is using `RecipesBase.jl` ([link](https://github.com/JuliaPlots/RecipesBase.jl)) as interface to generate all the plots. This means that it is up to the user to choose a plotting library that supports `RecipesBase.jl`, a suggestion would be `Plots.jl` with which the user is also able to freely choose a back-end. The plots in this manual are generated using `Plots.jl` with the `GR` backend. If you have several back-ends for plotting then you can select the one you want to use with the corresponding `Plots` call (for `GR` this is `Plots.gr()`, some alternatives are `pyplot(), plotly(), pgfplots()`). A simple example where we generate a plot and save it to a file is ```jldoctest; output=false +using Plots fig = bodeplot(tf(1,[1,2,1])) -Plots.savefig(fig, "myfile.svg") +savefig(fig, "myfile.svg") save_docs_plot(fig, "intro_bode.svg") # hide diff --git a/example/autodiff.jl b/example/autodiff.jl index ba47343c9..1bd22d4a1 100644 --- a/example/autodiff.jl +++ b/example/autodiff.jl @@ -1,4 +1,4 @@ -using ControlSystems, OrdinaryDiffEq, NLopt, BlackBoxOptim, ForwardDiff +using ControlSystems, OrdinaryDiffEq, NLopt, BlackBoxOptim, ForwardDiff, Plots p0 = [0.2,0.8,1] # Initial guess K(kp,ki,kd) = pid(kp=kp, ki=ki, kd=kd) K(p) = K(p...) diff --git a/example/dc_motor_lqg_design.jl b/example/dc_motor_lqg_design.jl index 19d5430bf..53e72e7d2 100644 --- a/example/dc_motor_lqg_design.jl +++ b/example/dc_motor_lqg_design.jl @@ -1,4 +1,4 @@ -using ControlSystems +using ControlSystems, Plots """ Example for designing an LQG speed controller for an electrical DC motor. """ @@ -29,7 +29,7 @@ function motor(Ke, Kt, L, R, J, b=1e-3) end p60 = motor(Ke, Kt, L, Rel, J) -f1 = stepplot(p60, 1) +f1 = plot(step(p60, 1)) f2 = bodeplot(p60) # LQR control @@ -60,5 +60,5 @@ S = 1-T # 1000 logarithmically spaced values from -3 to 3 f3 = bodeplot([Gcl, S, T], exp10.(range(-3, stop=3, length=1000))) -f4 = stepplot(Gcl, 1, label="Closed loop system using LQG") -Plots.plot(f1, f2, f3, f4, layout=(2,2), size=(800, 600)) +f4 = plot(step(Gcl, 1), label="Closed loop system using LQG") +plot(f1, f2, f3, f4, layout=(2,2), size=(800, 600)) diff --git a/example/delayed_lti_timeresp.jl b/example/delayed_lti_timeresp.jl index 0753c3480..22a55ae4d 100644 --- a/example/delayed_lti_timeresp.jl +++ b/example/delayed_lti_timeresp.jl @@ -19,9 +19,9 @@ function u0(out,t) return end -@time y, t, x = lsim(sys, u0, t) +@time res = lsim(sys, u0, t) -plot(t, y') +plot(res) s = tf("s") P = delay(2.6)*ss((s+3.0)/(s^2+0.3*s+1)) @@ -30,8 +30,8 @@ P*C T = feedback(P*C,1.0) t = 0:0.1:70 -y, t, x = lsim(T, t -> (t<0 ? 0 : 1 ), t) -plot(t, y, c = :blue) +res = lsim(T, t -> (t<0 ? 0 : 1 ), t) +plot(res, c = :blue) w = 10 .^ (-2:0.01:2) marginplot(P*C, w) @@ -41,16 +41,16 @@ notch = ss(tf([1, 0.2, 1],[1, .8, 1])); C = ss(0.05 * (1 + 1/s)); Tnotch = feedback(P*C*notch, 1.0) -stepplot(Tnotch) +plot(step(Tnotch)) y, t, x = step(C, method=:zoh) y2, t2, x2 = step(Tnotch) -stepplot(Tnotch) +plot(step(Tnotch)) -stepplot(Tnotch, 40, 0.1) +plot(step(Tnotch, 40)) -stepplot(T, 100) +plot(step(T, 100)) G = delay(5)/(s+1) T = feedback(G, 0.5) @@ -58,7 +58,6 @@ w = 10 .^ (-2:0.01:3) bodeplot(T, w, plotphase=false) # Test conversion, promotion -delay(1,Int64) + 3.5 G = 1 + 0.5 * delay(3) w = 10 .^(-2:0.001:2) @@ -67,14 +66,13 @@ bodeplot(G, w, plotphase=false) G = delay(1) * ((0.8*s^2+s+2)/(s^2+s)) T = feedback(G,1) # Not possible with direct term -stepplot(T) +plot(step(T)) bodeplot(T) G = 1/(s+1) + delay(4) T = feedback(1,G) # Not possible to lsim with direct term -stepplot(T) +plot(step(T)) bodeplot(T) -s = tf("s") diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 39c9977d5..3b09588b1 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -99,11 +99,10 @@ export LTISystem, # QUESTION: are these used? LaTeXStrings, Requires, IterTools -using Plots, LaTeXStrings, LinearAlgebra +using RecipesBase, LaTeXStrings, LinearAlgebra import Polynomials import Polynomials: Polynomial, coeffs using OrdinaryDiffEq -export Plots import Base: +, -, *, /, (==), (!=), isapprox, convert, promote_op import Base: getproperty, getindex import Base: exp # for exp(-s) diff --git a/src/pid_design.jl b/src/pid_design.jl index a7cf6ce06..64be02c3d 100644 --- a/src/pid_design.jl +++ b/src/pid_design.jl @@ -225,7 +225,7 @@ See also `Leadlink, leadlinkat` function leadlinkcurve(start=1) N = range(start, stop=10, length=50) dph = 180/pi*map(Ni->atan(sqrt(Ni))-atan(1/sqrt(Ni)), N) - Plots.plot(N,dph, xlabel="N", ylabel="Phase advance [deg]") + RecipesBase.plot(N,dph, xlabel="N", ylabel="Phase advance [deg]") end @@ -253,7 +253,7 @@ function stabregionPID(P, ω = _default_freq_vector(P,Val{:bode}()); kd=0, doplo phi = angle.(Pv) kp = -cos.(phi)./r ki = kd.*ω.^2 .- ω.*sin.(phi)./r - Plots.plot(kp,ki,linewidth = 1.5, xlabel=L"k_p", ylabel=L"k_i", title="Stability region of P, k_d = $(round(kd, digits=4))"), kp, ki + RecipesBase.plot(kp,ki,linewidth = 1.5, xlabel=L"k_p", ylabel=L"k_i", title="Stability region of P, k_d = $(round(kd, digits=4))"), kp, ki end @@ -263,7 +263,7 @@ function stabregionPID(P::Function, ω = exp10.(range(-3, stop=1, length=50)); k phi = angle.(Pv) kp = -cos.(phi)./r ki = kd.*ω.^2 .- ω.*sin.(phi)./r - Plots.plot(kp,ki,linewidth = 1.5, xlabel=L"k_p", ylabel=L"k_i", title="Stability region of P, k_d = $(round(kd, digits=4))"), kp, ki + RecipesBase.plot(kp,ki,linewidth = 1.5, xlabel=L"k_p", ylabel=L"k_i", title="Stability region of P, k_d = $(round(kd, digits=4))"), kp, ki end diff --git a/src/plotting.jl b/src/plotting.jl index 565692fcb..a19fa733f 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -68,11 +68,7 @@ function getPhaseTicks(x, minmax) ## this helps identifying at the edges. major = [(min-0.5);min:max;(max+0.5)].*90 end - if Plots.backend() != Plots.GRBackend() - majorText = [latexstring("\$ $(round(Int64,i))\$") for i = major] - else - majorText = ["$(round(Int64,i))" for i = major] - end + majorText = ["$(round(Int64,i))" for i = major] return major, majorText @@ -85,19 +81,12 @@ function getLogTicks(x, minmax) min = minx <= 0 ? minimum(x) : ceil(log10(minx)) max = floor(log10(maxx)) major = exp10.(min:max) - if Plots.backend() ∉ [Plots.GRBackend(), Plots.PlotlyBackend()] - majorText = [latexstring("\$10^{$(round(Int64,i))}\$") for i = min:max] - else - majorText = ["10^{$(round(Int64,i))}" for i = min:max] - end + + majorText = ["10^{$(round(Int64,i))}" for i = min:max] + if max - min < major_minor_limit minor = [j*exp10(i) for i = (min-1):(max+1) for j = 2:9] - if Plots.backend() ∉ [Plots.GRBackend(), Plots.PlotlyBackend()] - minorText = [latexstring("\$$j\\cdot10^{$(round(Int64,i))}\$") for i = (min-1):(max+1) for j = 2:9] - else - minorText = ["$j*10^{$(round(Int64,i))}" for i = (min-1):(max+1) for j = 2:9] - end - + minorText = ["$j*10^{$(round(Int64,i))}" for i = (min-1):(max+1) for j = 2:9] ind = findall(minx .<= minor .<= maxx) minor = minor[ind] minorText = minorText[ind] @@ -189,7 +178,7 @@ optionally provided. To change the Magnitude scale see `setPlotScale(str)` If `hz=true`, the plot x-axis will be displayed in Hertz, the input frequency vector is still treated as rad/s. -`kwargs` is sent as argument to Plots.plot. +`kwargs` is sent as argument to RecipesBase.plot. """ bodeplot @@ -258,37 +247,26 @@ end w = x magdata = y seriestype := :path - primary := false - @series begin - grid --> true - yscale --> :log10 - xscale --> :log10 - yguide --> "Magnitude" - xticks --> getLogTicks(w, getlims(:xlims, plotattributes, w)) - yticks --> getLogTicks(magdata, getlims(:ylims, plotattributes, magdata)) - x := w; y := magdata - () - end - x := [] - y := [] + primary --> false + grid --> true + yscale --> :log10 + xscale --> :log10 + yguide --> "Magnitude" + x := w + y := magdata () end @recipe function f(::Type{Val{:bodephase}}, x, y, z) w = x phasedata = y seriestype := :path - primary := false - @series begin - grid --> true - xscale --> :log10 - yguide --> "Phase (deg)" - xguide --> "Frequency (rad/s)" - xticks --> getLogTicks(w, getlims(:xlims, plotattributes, w)) - x := w; y := phasedata - () - end - x := [] - y := [] + primary --> false + grid --> true + xscale --> :log10 + yguide --> "Phase (deg)" + xguide --> "Frequency (rad/s)" + x := w + y := phasedata () end @@ -402,7 +380,7 @@ fontsize = 10 `val` ∈ [0,1] determines the brightness of the gain lines Additional keyword arguments are sent to the function plotting the systems and can be -used to specify colors, line styles etc. using regular Plots.jl syntax +used to specify colors, line styles etc. using regular RecipesBase.jl syntax This function is based on code subject to the two-clause BSD licence Copyright 2011 Will Robertson @@ -468,7 +446,7 @@ nicholsplot offset = (l+1) TextX = Niϕ(k,210) .+offset TextY = Ni_Ga(k,210) - annotations := (TextX,TextY,Plots.text("$(string(k)) dB",fontsize)) + annotations := (TextX,TextY,("$(string(k)) dB")) end ϕVals .+ 360(l+1),GVals end @@ -513,7 +491,7 @@ nicholsplot end end TextX - annotations := (TextX,TextY,Plots.text("$(string(k))°",fontsize)) + annotations := (TextX,TextY,("$(string(k))°")) title --> "Nichols chart" grid --> false @@ -581,6 +559,7 @@ sigmaplot end end +@userplot Marginplot """ fig = marginplot(sys::LTISystem [,w::AbstractVector]; kwargs...) marginplot(sys::Vector{LTISystem}, w::AbstractVector; kwargs...) @@ -588,19 +567,20 @@ end Plot all the amplitude and phase margins of the system(s) `sys`. A frequency vector `w` can be optionally provided. -`kwargs` is sent as argument to Plots.plot. +`kwargs` is sent as argument to RecipesBase.plot. """ -function marginplot(systems::Union{AbstractVector{T},T}, args...; kwargs...) where T<:LTISystem - systems, w = _processfreqplot(Val{:bode}(), systems, args...) +@recipe function marginplot(p::Marginplot) + systems, w = _processfreqplot(Val{:bode}(), p.args...) ny, nu = size(systems[1]) - fig = bodeplot(systems, w; kwargs...) s2i(i,j) = LinearIndices((ny,2nu))[j,i] titles = Array{AbstractString}(undef, nu,ny,2,2) titles[:,:,1,1] .= "Gm: " titles[:,:,2,1] .= "Pm: " titles[:,:,1,2] .= "Wgm: " titles[:,:,2,2] .= "Wpm: " + layout --> (2ny, nu) for (si, s) in enumerate(systems) + bmag, bphase = bode(s, w) for j=1:nu for i=1:ny wgm, gm, wpm, pm, fullPhase = sisomargin(s[i,j],w, full=true, allMargins=true) @@ -625,32 +605,65 @@ function marginplot(systems::Union{AbstractVector{T},T}, args...; kwargs...) whe mag = 1 ./ gm oneLine = 1 end - for k=1:length(wgm) - #Plot gain margins - Plots.plot!(fig, [wgm[k];wgm[k]], [1;mag[k]]; lab="", subplot=s2i(2i-1,j), group=si) - end - #Plot gain line at 1 - Plots.plot!(fig, [w[1],w[end]], [oneLine,oneLine], l=:dash, c=:gray, lab="", subplot=s2i(2i-1,j)) 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],", ")*"] " - for k=1:length(wpm) - #Plot the phase margins - Plots.plot!(fig, [wpm[k];wpm[k]],[fullPhase[k];fullPhase[k]-pm[k]]; lab="", subplot=s2i(2i,j)) - #Plot the line at 360*k - Plots.plot!(fig, [w[1],w[end]],(fullPhase[k]-pm[k])*ones(2); l=:dash, c=:gray, lab="", subplot=s2i(2i,j)) - end 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],", ")*"] " + + @series begin + primary := true + subplot --> s2i(2i-1,j) + seriestype := :bodemag + w, bmag[:, i, j] + end + + primary --> false + #Plot gain margins + @series begin + subplot --> s2i(2i-1,j) + primary --> false + color --> :gray + linestyle --> :dash + [w[1],w[end]], [oneLine,oneLine] + end + @series begin + subplot --> s2i(2i-1,j) + title --> titles[j,i,1,1]*" "*titles[j,i,1,2] + [wgm wgm]', [ones(length(mag)) mag]' + end + + # Phase margins + @series begin + primary := true + subplot --> s2i(2i,j) + seriestype := :bodephase + w, bphase[:, i, j] + end + @series begin + subplot --> s2i(2i,j) + primary --> false + color --> :gray + linestyle --> :dash + [w[1],w[end]], [oneLine,oneLine] + end + @series begin + primary --> false + subplot --> s2i(2i,j) + [wpm wpm]', [fullPhase fullPhase-pm]' + end + @series begin + title --> titles[j,i,2,1]*" "*titles[j,i,2,2] + subplot --> s2i(2i,j) + primary --> false + color --> :gray + linestyle --> :dash + [w[1] w[end]]', ((fullPhase .- pm) .* ones(1, 2))' + end + end end end - for j = 1:nu - for i = 1:ny - Plots.title!(fig, titles[j,i,1,1]*" "*titles[j,i,1,2], subplot=s2i(2i-1,j)) - Plots.title!(fig, titles[j,i,2,1]*" "*titles[j,i,2,2], subplot=s2i(2i,j)) - end - end - return fig + end # HELPERS: @@ -722,7 +735,7 @@ pzmap!(sys::LTISystem; kwargs...) = pzmap!([sys]; kwargs...) fig = gangoffourplot(P::LTISystem, C::LTISystem; minimal=true, plotphase=false, kwargs...) gangoffourplot(P::Union{Vector, LTISystem}, C::Vector; minimal=true, plotphase=false, kwargs...) -Gang-of-Four plot. `kwargs` is sent as argument to Plots.plot. +Gang-of-Four plot. `kwargs` is sent as argument to RecipesBase.plot. """ function gangoffourplot(P::Union{Vector, LTISystem}, C::Vector, args...; minimal=true, plotphase=false, kwargs...) if P isa LTISystem # Don't broadcast over scalar (with size?) @@ -730,12 +743,12 @@ function gangoffourplot(P::Union{Vector, LTISystem}, C::Vector, args...; minimal end sys = gangoffour.(P,C; minimal=minimal) fig = bodeplot([[sys[i][1] sys[i][2]; sys[i][3] sys[i][4]] for i = 1:length(C)], args..., plotphase=plotphase; kwargs...) - hline!([1 1 1 1], l=(:black, :dash), primary=false) + RecipesBase.plot!(fig, [x-> _PlotScale == "dB" ? 0 : 1 for _ in 1:4], l=(:black, :dash), primary=false) titles = fill("", 1, plotphase ? 8 : 4) # Empty titles on phase titleIdx = plotphase ? [1,2,5,6] : [1,2,3,4] titles[titleIdx] = ["S = 1/(1+PC)", "P/(1+PC)", "C/(1+PC)", "T = PC/(1+PC)"] - Plots.plot!(fig, title = titles) + RecipesBase.plot!(fig, title = titles) return fig end diff --git a/src/simulators.jl b/src/simulators.jl index 069c2fd0d..ae950c4c1 100644 --- a/src/simulators.jl +++ b/src/simulators.jl @@ -22,7 +22,7 @@ end Used to simulate continuous-time systems. See function `?solve` for additional info. # Usage: ``` -using OrdinaryDiffEq +using OrdinaryDiffEq, Plots dt = 0.1 tfinal = 20 t = 0:dt:tfinal diff --git a/src/synthesis.jl b/src/synthesis.jl index 36b37fd93..261a2d839 100644 --- a/src/synthesis.jl +++ b/src/synthesis.jl @@ -19,6 +19,7 @@ See also `LQG` Usage example: ```julia using LinearAlgebra # For identity matrix I +using Plots A = [0 1; 0 0] B = [0;1] C = [1 0] @@ -86,6 +87,7 @@ The `args...; kwargs...` are sent to the Riccati solver, allowing specification Usage example: ```julia using LinearAlgebra # For identity matrix I +using Plots Ts = 0.1 A = [1 Ts; 0 1] B = [0;1]