margin() does not compute all margins even when allMargins=true

[hurak]

The function margin() does not compute all the margins even when the setting allMargins=true. Below I insert a code aimed at designing an LQG controller. The open loop has one unstable pole, hence the Nyquist close has to encircle the point -1 just once, and the particular character of the Nyquist curve implies that there are both GM_plus and GM_minus margins. But the function computes only the former.

using ControlSystems: ss, lqr, margin, marginplot, bodeplot, bodeplot!, nyquistplot, nyquistplot!, kalman, poles
using LinearAlgebra
using MatrixEquations: arec
using Plots
using Printf

# Plant model
A  = [0.0 1.0; 0.5 0.0]
Bu = [0.0; 0.1]
Bw = [0.0; 0.1]
C  = [1.0 0.0]

# Plant state-space model
G = ss(A, Bu, C, 0.0)

# Spectral densities
Sw = 1.0    # Disturbance (aka process noise) spectral density
Sv = 1.0    # Measurement noise spectral density

# LQR gain
Q = Diagonal([1.0, 0.0])
R = 16.0
K = lqr(A, Bu, Q, R)

# Kalman filter gain
L = kalman(A, C, Bw*Sw*Bw', Sv)   # Note that the function assumes Bw = I.

# Build the output-feedback controller and the open-loop transfer function L(S) = K(s)·G(s)
Klqg  = ss(A - L*C - Bu*K, L, K, zeros(size(K,1), size(C,1)))  # Observer-controller (output u = -K*x̂)
OL = Klqg*G                                                    # Open-loop transfer function

# The open loop has one unstable pole
poles(OL) 

# Compute the GM and PM margins
wgm, gm, wpm, pm = margin(OL, allMargins=true)

# Margin plotted
ω = logrange(1e-2, 1e3, length=500)
p_lqg = marginplot(OL, ω, title="Stability margins LQG")
display(plot(p_lqg))

# Nyquist plots
# The plant has one unstable open-loop pole, so a stable closed loop
# requires exactly one counter-clockwise encirclement of the critical point −1.
pn = nyquistplot(OL, ω)
θ = range(0, 2π, length=300)
plot!(pn, cos.(θ), sin.(θ), color=:grey, lw=1, alpha=0.5, label="", aspect_ratio=:equal, xlabel="Re", ylabel="Im", title="Nyquist plot")
display(pn)

Here I show the output of the margin function:

julia> wgm, gm, wpm, pm = margin(OL, allMargins=true)
(wgm = [[1.2311038829891778];;], gm = [[2.005125180939021];;], wpm = [[0.3830306271727272];;], pm = [[10.600059869029764];;])

and the nyquistplot:

Indeed, the computation only gives the GMplus, while the plot confirms that there should also be the "attenuation GM". And it should be around 0.8747 (based on Matlab).

[baggepinnen]

I believe this to be because we do not include negative frequencies in the frequency grid used for the margin computations, if I include one negative frequency the other margin is found

julia> wgm, gm, wpm, pm = margin(OL, [-1; ω], allMargins=true)
(wgm = [[-0.04241083165812409, 1.2311117077345854];;], gm = [[0.8749978316484006, 2.0051395673666277];;], wpm = [[-0.2653375947456835, 0.38305494108339083];;], pm = [[351.73329694522033, 370.6004299256333];;])