add gainphaseplot

Author: baggepinnen

docs/src/index.md

@@ -253,6 +253,12 @@ dms = diskmargin(L, 0, w)
 plot(dms)
 ```
 
+We can also visualize the plane of complex perturbations that are allowed to be simultaneously applied to the system:
+```@example diskmargin
+gainphaseplot(L)
+```
+The green regions are stable perturbations while red regions are unstable. The diskmargin is the largest disk that can be placed entirely inside the green area. The center of the disk is determined by the skew of the diskmargin. The classical gain margin is the length of the green area along the x-axis starting at the point 1, while the classical phase margin is the length of the green area along the unit circle starting a the point 1.
+
 
 ## Closed-loop analysis
 - [`output_sensitivity`](@ref)

src/diskmargin.jl

@@ -257,59 +257,64 @@ end
     end
 end
 
-@recipe function plotdiskmargin(dm::AbstractVector{<:Diskmargin}; lower=true, phase=true)
+@recipe function plotdiskmargin(dm::AbstractVector{<:Diskmargin}; lower=true, phase=true, gain=true)
+    gain || phase || error("Must plot either gain or phase")
     w = [dm.ω0 for dm in dm]
-    layout --> (phase ? 2 : 1, 1)
+    layout --> (phase + gain, 1)
     link --> :x
     gma = getfield.(dm, :γmax)
     gmi = getfield.(dm, :γmin)
     # ylims --> (0, min(10, maximum(gma)))
     data = (lower ? [gmi gma] : gma)
     data = max.(0, data)
     replace!(data, 0 => -Inf)
-    @series begin
-        subplot --> 1
-        title --> "Gain margin"
-        label --> (lower ? ["Lower" "Upper"] : "Upper")
-        # xguide --> "Frequency"
-        xscale --> :log10
-        yscale --> :log10
-        w, data
-    end
-    replace!(data, -Inf => Inf)
-    m,i = findmin(data[:, end])
-    @series begin
-        subplot --> 1
-        primary := true
-        seriestype := :scatter
-        seriescolor := :red
-        label := string(round(m, sigdigits=3))
-        [w[i]], [m]
-    end
-    @series begin
-        primary := false
-        seriestype := :hline
-        linecolor := :black
-        [1]
-    end
-    phase || return
-    ϕm = getfield.(dm, :ϕm)
-    m,i = findmin(ϕm)
-    @series begin
-        subplot --> 2
-        primary := true
-        seriestype := :scatter
-        seriescolor := :red
-        label := string(round(m, sigdigits=3))
-        [w[i]], [m]
+    @show gain, phase
+    if gain
+        @series begin
+            # subplot --> 1
+            title --> "Gain margin"
+            label --> (lower ? ["Lower" "Upper"] : "Upper")
+            # xguide --> "Frequency"
+            xscale --> :log10
+            yscale --> :log10
+            w, data
+        end
+        replace!(data, -Inf => Inf)
+        m,i = findmin(data[:, end])
+        @series begin
+            # subplot --> 1
+            primary := true
+            seriestype := :scatter
+            seriescolor := :red
+            label := string(round(m, sigdigits=3))
+            [w[i]], [m]
+        end
+        @series begin
+            primary := false
+            seriestype := :hline
+            linecolor := :black
+            [1]
+        end
     end
-    @series begin
-        subplot --> 2
-        title --> "Phase margin (deg)"
-        xguide --> "Frequency"
-        xscale --> :log10
-        label --> ""
-        w, ϕm
+    if phase
+        ϕm = getfield.(dm, :ϕm)
+        m,i = findmin(ϕm)
+        @series begin
+            subplot --> (phase + gain)
+            primary := true
+            seriestype := :scatter
+            seriescolor := :red
+            label := string(round(m, sigdigits=3))
+            [w[i]], [m]
+        end
+        @series begin
+            subplot --> (phase + gain)
+            title --> "Phase margin (deg)"
+            xguide --> "Frequency"
+            xscale --> :log10
+            label --> ""
+            w, ϕm
+        end
     end
 end
 
@@ -354,6 +359,64 @@ end
     end
 end
 
+@userplot GainPhasePlot
+
+@recipe function f(gp::GainPhasePlot)
+    if length(gp.args) == 3
+        P, re, im = gp.args
+    elseif length(gp.args) == 1
+        P = gp.args[1]
+        re = LinRange(-1.5, 10, 300)
+        im = LinRange(-3, 3, 300)
+    else
+        error("gainphaseplot takes either 1 (sys) or 3 arguments (sys, re, im).")
+    end
+
+    res = map(Iterators.product(re,im)) do (r,i)
+        isstable(feedback(P*complex(r,i)))
+    end
+    colorbar --> false
+    aspect_ratio --> 1
+    @series begin
+        color --> [:red, :green]
+        seriesalpha --> 0.5
+        seriestype := :heatmap
+        re, im, res'
+    end
+    θ = LinRange(0, 2π, 100)
+    S = sin.(θ)
+    C = cos.(θ)
+    framestyle := :zerolines
+    @series begin
+        primary := false
+        linestyle := :dash
+        seriescolor := :black
+        C, S
+    end
+    r = [0, 1.5]
+    c = [:red, :orange, :green]
+    for (i, θ) = enumerate([30, 45, 60])
+        points = r*cis(deg2rad(θ))
+        @series begin
+            linestyle := :dash
+            linewidth --> 2
+            xguide --> "Re"
+            yguide --> "Im"
+            seriescolor --> c[i]
+            label --> string(θ)*"°"
+            real(points), imag(points)
+        end
+    end
+end
+
+"""
+    gainphaseplot(P)
+    gainphaseplot(P, re, im)
+
+Plot complex perturbantions to the plant `P` and indicate whether or not the closed-loop system is stable. The diskmargin is the largest disk that can be fit inside the green region that only contains stable variations.
+"""
+gainphaseplot
+
 
 """
     passivity_index(P; kwargs...)

test/test_diskmargin.jl

@@ -14,6 +14,7 @@ plot(dm)
 nyquistplot(L)
 plot!(dm, nyquist=true)
 plot!(Disk(dm), nyquist=true)
+gainphaseplot(L)
 
 @test dm.alpha == dm.α
 @test length(dm.gainmargin) == 2