add functions relating Ms peak to gain and phase margins

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -80,6 +80,9 @@ export  LTISystem,
         output_comp_sensitivity,
         G_PS,
         G_CS,
+        margin_bounds,
+        Ms_from_phase_margin,
+        Ms_from_gain_margin,
         resolvent,
         input_resolvent,
         # Discrete

lib/ControlSystemsBase/src/analysis.jl

@@ -73,7 +73,7 @@ end
 
 Count the number of integrators in the system by finding the difference between the number of poles in the origin and the number of zeros in the origin. If the number of zeros in the origin is greater than the number of poles in the origin, the count is negative.
 
-For discrete-tiem systems, the origin ``s = 0`` is replaced by the point ``z = 1``.
+For discrete-time systems, the origin ``s = 0`` is replaced by the point ``z = 1``.
 """
 function integrator_excess(P::LTISystem)
     p = poles(P)

lib/ControlSystemsBase/src/plotting.jl

@@ -416,7 +416,7 @@ Create a Nyquist plot of the `LTISystem`(s). A frequency vector `w` can be
 optionally provided.
 
 - `unit_circle`: if the unit circle should be displayed. The Nyquist curve crosses the unit circle at the gain crossover frequency.
-- `Ms_circles`: draw circles corresponding to given levels of sensitivity (circles around -1 with  radii `1/Ms`). `Ms_circles` can be supplied as a number or a vector of numbers. A design staying outside such a circle has a phase margin of at least `2asin(1/(2Ms))` rad and a gain margin of at least `Ms/(Ms-1)`.
+- `Ms_circles`: draw circles corresponding to given levels of sensitivity (circles around -1 with  radii `1/Ms`). `Ms_circles` can be supplied as a number or a vector of numbers. A design staying outside such a circle has a phase margin of at least `2asin(1/(2Ms))` rad and a gain margin of at least `Ms/(Ms-1)`. See also [`margin_bounds`](@ref), [`Ms_from_phase_margin`](@ref) and [`Ms_from_gain_margin`](@ref).
 - `Mt_circles`: draw circles corresponding to given levels of complementary sensitivity. `Mt_circles` can be supplied as a number or a vector of numbers.
 - `critical_point`: point on real axis to mark as critical for encirclements
 - If `hz=true`, the hover information will be displayed in Hertz, the input frequency vector is still treated as rad/s.

lib/ControlSystemsBase/src/sensitivity_functions.jl

@@ -27,6 +27,7 @@ The output [sensitivity function](https://en.wikipedia.org/wiki/Sensitivity_(con
 ```math
 ϕ_m ≥ 2\\text{sin}^{-1}(\\frac{1}{2M_S}), g_m ≥ \\frac{M_S}{M_S-1}
 ```
+(see [`margin_bounds`](@ref) for a function that computes these bounds, and [`Ms_from_phase_margin`](@ref) and [`Ms_from_gain_margin`](@ref) for the inverse functions.)
 Generally, bounding ``M_S`` is a better objective than looking at gain and phase margins due to the possibility of combined gain and pahse variations, which may lead to poor robustness despite large gain and pahse margins.
 
 $sensdoc
@@ -35,6 +36,39 @@ function sensitivity(args...)# Sensitivity function
     return output_sensitivity(args...)
 end
 
+
+"""
+    g_m, ϕ_m = margin_bounds(M_S)
+
+Compute the phase margin lower bound ϕ_m (in radians) and gain margin lower bound g_m given a maximum sensitivity peak ``M_S = ||S||_∞``. These bounds are derived from the fact that the inverse of the sensitivity function is the distance from the open-loop Nyquist curve to the critical point -1.
+
+See also [`Ms_from_phase_margin`](@ref) and [`Ms_from_gain_margin`](@ref) for the inverse functions. The circle corresponding to the maximum sensitivity peak ``M_S`` can be plotted in [`nyquistplot`](@ref) by passing the keyword argument `Ms_circles = [Ms]`.
+"""
+function margin_bounds(M_S)
+    M_S < 1 && error("Maximum sensitivity peak M_S must be greater than or equal to 1, got $M_S")
+    ϕ_m = 2 * asin(1 / (2 * M_S))
+    g_m = M_S / (M_S - 1)
+    return (; g_m, ϕ_m, ϕ_m_deg = rad2deg(ϕ_m))
+end
+
+"""
+    Ms_from_phase_margin(ϕ_m)
+
+Compute the maximum sensitivity peak ``M_S = ||S||_∞`` such that if respected, gives a phase margin of at least ϕ_m (in radians).
+
+See also [`Ms_from_gain_margin`](@ref) and [`margin_bounds`](@ref).
+"""
+Ms_from_phase_margin(ϕ_m) = 1 / (2 * sin(ϕ_m / 2))
+
+"""
+    Ms_from_gain_margin(g_m)
+
+Compute the maximum sensitivity peak ``M_S = ||S||_∞`` such that if respected, gives a gain margin of at least g_m.
+
+See also [`Ms_from_phase_margin`](@ref) and [`margin_bounds`](@ref).
+"""
+Ms_from_gain_margin(g_m) = g_m / (g_m - 1)
+
 """
 See [`output_comp_sensitivity`](@ref)
 $sensdoc

lib/ControlSystemsBase/test/test_analysis.jl

@@ -370,6 +370,11 @@ gof = gangoffour(P,C)
 F = ssrand(2,2,2,proper=true)
 @test_nowarn gangofseven(P, C, F)
 
+## Test bound functions
+M_S = 1.3
+g_m, ϕ_m = margin_bounds(M_S)
+@test Ms_from_gain_margin(g_m) ≈ M_S
+@test Ms_from_phase_margin(ϕ_m) ≈ M_S
 
 ## Approximate double integrator
 P = let