Disk margin calculations #1146

josiahdelange · 2025-04-24T12:06:11Z

Am relatively new to this toolbox, have been using it here and there and found #726.

This is an initial prototype Python implementation of disk margins, built on top of python-control and slycot. The code should work both for SISO and MIMO systems, the latter of which requires slycot (for SLICOT's AB13MD) to compute the upper bound of $\mu$ at each discrete frequency. The function disk_margins computes disk margins (and corresponding disk-based gain/phase margins), optionally returning the whole frequency-dependent vectors for further plotting.

It's been verified against the example SISO loop transfer functions in the published paper and the "spinning satellite" MIMO example from the MathWorks documentation. I also confirmed SISO disk margins match the relevant output of existing function stability_margins corresponding to the Nyquist plot's distance to -1. All seem to match within a few significant digits, so the general behavior seems correct. Might be good to get another set of eyes to double check/test further.

I tried to base as much as possible (e.g. style conventions) on existing code in control/margins.py. Example usage (see examples/disk_margin.py):

import os, sys, math
import numpy as np
import control

import math
import matplotlib
import matplotlib.pyplot as plt
from warnings import warn

import numpy as np
import scipy as sp

# Frequencies of interest
omega = np.logspace(-1, 3, 1001)

# Laplace variable
s = control.tf('s')

# MIMO loop transfer gain for the spinning satellite example
# https://www.mathworks.com/help/robust/ug/mimo-stability-margins-for-spinning-satellite.html
P = control.ss([[0, 10],[-10, 0]], np.eye(2), [[1, 10], [-10, 1]], [[0, 0],[0, 0]]) # plant
C = control.ss([],[],[], [[1, -2], [0, 1]]) # controller
L = P*C # output loop gain

print(f"------------- Sensitivity function (S) -------------")
DM, GM, PM = control.disk_margins(L, omega, skew = 1.0, returnall = True) # S-based (S)
print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
print(f"GM = {GM[np.argmin(DM)]} dB")
print(f"PM = {PM[np.argmin(DM)]} deg")
print(f"min(GM) = {min(GM)} dB")
print(f"min(PM) = {min(PM)} deg\n")

plt.figure(3)
plt.subplot(3,3,1)
plt.semilogx(omega, DM, label='$\\alpha$')
plt.legend()
plt.title('Disk Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 2])

plt.figure(3)
plt.subplot(3,3,4)
plt.semilogx(omega, GM, label='$\\gamma_{m}$')
plt.ylabel('Gain Margin (dB)')
plt.legend()
plt.title('Gain-Only Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 40])

plt.figure(3)
plt.subplot(3,3,7)
plt.semilogx(omega, PM, label='$\\phi_{m}$')
plt.ylabel('Phase Margin (deg)')
plt.legend()
plt.title('Phase-Only Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 90])

print(f"------------- Complementary sensitivity function (T) -------------")
DM, GM, PM = control.disk_margins(L, omega, skew = -1.0, returnall = True) # T-based (T)
print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
print(f"GM = {GM[np.argmin(DM)]} dB")
print(f"PM = {PM[np.argmin(DM)]} deg")
print(f"min(GM) = {min(GM)} dB")
print(f"min(PM) = {min(PM)} deg\n")

plt.figure(3)
plt.subplot(3,3,2)
plt.semilogx(omega, DM, label='$\\alpha$')
plt.legend()
plt.title('Disk Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 2])

plt.figure(3)
plt.subplot(3,3,5)
plt.semilogx(omega, GM, label='$\\gamma_{m}$')
plt.ylabel('Gain Margin (dB)')
plt.legend()
plt.title('Gain-Only Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 40])

plt.figure(3)
plt.subplot(3,3,8)
plt.semilogx(omega, PM, label='$\\phi_{m}$')
plt.ylabel('Phase Margin (deg)')
plt.legend()
plt.title('Phase-Only Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 90])

print(f"------------- Balanced sensitivity function (S - T) -------------")
DM, GM, PM = control.disk_margins(L, omega, skew = 0.0, returnall = True) # balanced (S - T)
print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
print(f"GM = {GM[np.argmin(DM)]} dB")
print(f"PM = {PM[np.argmin(DM)]} deg")
print(f"min(GM) = {min(GM)} dB")
print(f"min(PM) = {min(PM)} deg\n")

plt.figure(3)
plt.subplot(3,3,3)
plt.semilogx(omega, DM, label='$\\alpha$')
plt.legend()
plt.title('Disk Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 2])

plt.figure(3)
plt.subplot(3,3,6)
plt.semilogx(omega, GM, label='$\\gamma_{m}$')
plt.ylabel('Gain Margin (dB)')
plt.legend()
plt.title('Gain-Only Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 40])

plt.figure(3)
plt.subplot(3,3,9)
plt.semilogx(omega, PM, label='$\\phi_{m}$')
plt.ylabel('Phase Margin (deg)')
plt.legend()
plt.title('Phase-Only Margin')
plt.grid()
plt.xlim([omega[0], omega[-1]])
plt.ylim([0, 90])

Output:

------------- Sensitivity function (S) -------------
min(DM) = 0.3697398698410271 (omega = 0.1)
GM = 2.732761944304522 dB
PM = 21.30709883385843 deg
min(GM) = 2.732761944304522 dB
min(PM) = 21.30709883385843 deg

------------- Complementary sensitivity function (T) -------------
min(DM) = 0.3683035923739884 (omega = 0.1)
GM = 2.7236493433653735 dB
PM = 21.223368586091564 deg
min(GM) = 2.7236493433653735 dB
min(PM) = 21.223368586091564 deg

------------- Balanced sensitivity function (S - T) -------------
min(DM) = 0.3769872636944355 (omega = 0.1)
GM = 3.3140985364257256 dB
PM = 21.349285542574695 deg
min(GM) = 3.3140985364257256 dB
min(PM) = 21.349285542574695 deg

The example script also shows how to plot the allowable/stable region of gain and phase variations which will not destabilize the loop, in a local function plot_allowable_region, e.g.

.
.
.
DM_plot = []
DM_plot.append(control.disk_margins(L, omega, skew = -1.0)[0]) # T-based (T)
DM_plot.append(control.disk_margins(L, omega, skew = 0.0)[0]) # balanced (S - T)
DM_plot.append(control.disk_margins(L, omega, skew = 1.0)[0]) # S-based (S)
plt.figure(30)
plot_allowable_region(DM_plot, skew = [-1.0, 0.0, 1.0])

…andler comment, fix typo in skew description of disk_margins docstring

…e function to plot allowable gain/phase variations.

…n the margins module

calculation of 'f', the bounding complex curve. Seems to look correct for balanced (skew = 0) case, still verifying the skewed equivalent.

…on-control into jdelange/disk-margins

…ol/tests/margin_test.py

slivingston · 2025-04-24T20:57:12Z

control/margins.py

+try:
+    from . import mag2db
+except ImportError:
+    # Likely due the following circular import issue:
+    #
+    # ImportError: cannot import name 'mag2db' from partially initialized module
+    # 'control' (most likely due to a circular import) (control/__init__.py)
+    #
+    def mag2db(mag):
+        return 20*np.log10(mag)


Try to import from the module where this is defined:

from .ctrlutil import mag2db

slivingston · 2025-04-24T21:02:18Z

control/margins.py

+    >>
+    >> omega = np.logspace(-1, 3, 1001)
+    >>
+    >> P = control.ss([[0, 10],[-10, 0]], np.eye(2), [[1, 10], [-10, 1]], [[0, 0],[0, 0]])


This example seems too big for a docstring, and anyway, it already is given in examples/disk_margins.py. Is there a smaller example that can be given in this docstring?

Agree. Yes - will remove all the plotting.

…ation

josiahdelange · 2025-04-26T16:03:21Z

I fixed the docstrings and implemented the linter's recommendations so that automated workflows finally pass. One thing I previously hadn't appreciated was the complexity of the control plotting code, so to make this PR simpler I'm only proposing to add a library function to calculate the disk-based margins. All other plotting is left to user code (or future PRs), e.g. the local function plot_allowable_region in examples/disk_margins.py.

…it tests

coveralls · 2025-04-28T02:17:19Z

coverage: 94.732% (-0.01%) from 94.746%
when pulling fe79760 on josiahdelange:jdelange/disk-margins
into 632391c on python-control:main.

murrayrm

Some additional suggestions attached, mainly focused on coding style (to be consistent with the standard style in the package).

murrayrm · 2025-04-28T02:07:53Z

control/tests/margin_test.py

+            assert_allclose([DMo], [0.3754], atol = 0.1) # disk margin of 0.3754
+            assert_allclose([DGMo], [3.3], atol = 0.1) # disk-based gain margin of 3.3 dB 
+            assert_allclose([DPMo], [21.26], atol = 0.1) # disk-based phase margin of 21.26 deg
+
+            # Balanced (S - T) disk-based stability margins at plant input
+            DMi, DGMi, DPMi = disk_margins(Li, omega, skew = 0.0)
+            assert_allclose([DMi], [0.3754], atol = 0.1) # disk margin of 0.3754
+            assert_allclose([DGMi], [3.3], atol = 0.1) # disk-based gain margin of 3.3 dB 
+            assert_allclose([DPMi], [21.26], atol = 0.1) # disk-based phase margin of 21.26 deg


This code never runs if slycot is not loaded => can you just remove it?

Yes - will fix.

murrayrm · 2025-04-28T02:09:49Z

control/tests/margin_test.py

+    Lo = P*K # loop transfer function, broken at plant output
+    Li = K*P # loop transfer function, broken at plant input
+
+    if importlib.util.find_spec('slycot') == None:


You can just use ct.slycot_check() here and avoid using the importlib call.

murrayrm · 2025-04-28T02:12:04Z

control/margins.py

@@ -517,3 +521,130 @@ def margin(*args):
                         % len(args))

    return margin[0], margin[1], margin[3], margin[4]
+
+def disk_margins(L, omega, skew = 0.0, returnall = False):


Remove spaces around = to be consistent with python-control (PEP8) coding conventions.

murrayrm · 2025-04-28T02:12:53Z

control/margins.py

+        skew = 0 uses the "balanced" sensitivity function 0.5*(S - T)
+        skew = 1 uses the sensitivity function S
+        skew = -1 uses the complementary sensitivity function T


Which of these is the default value?

Will update. Default of 0.0 is consistent with the MATLAB function diskmargin.

murrayrm · 2025-04-28T02:13:17Z

control/margins.py

+        skew = 1 uses the sensitivity function S
+        skew = -1 uses the complementary sensitivity function T
+    returnall : bool, optional
+        If true, return frequency-dependent margins. If False (default),


true → True

murrayrm · 2025-04-28T02:18:57Z

control/margins.py

+    DM = np.zeros(omega.shape, np.float64)
+    DGM = np.zeros(omega.shape, np.float64)
+    DPM = np.zeros(omega.shape, np.float64)
+    for ii in range(0,len(omega)):


Missing space after comma (PEP8).

murrayrm · 2025-04-28T02:21:53Z

control/margins.py

+    DPM = np.zeros(omega.shape, np.float64)
+    for ii in range(0,len(omega)):
+        # Disk margin (a.k.a. "alpha") vs. frequency
+        if L.issiso() and (ab13md == None):


Parentheses around second expression are not needed.

murrayrm · 2025-04-28T02:22:24Z

control/margins.py

+                np.array(num_loops*[2]))[0]
+
+        # Disk-based gain margin (dB) and phase margin (deg)
+        with np.errstate(divide = 'ignore', invalid = 'ignore'):


Remove excess spaces around keyword assignment operator.

murrayrm · 2025-04-28T02:22:55Z

control/margins.py

+
+    if returnall:
+        # Frequency-dependent disk margin, gain margin and phase margin
+        return (DM, DGM, DPM)


Parentheses not needed here.

murrayrm · 2025-04-28T02:31:59Z

control/margins.py

+        return ((not DM.shape[0] and float('inf')) or np.amin(DM),
+            (not gmidx != -1 and float('inf')) or DGM[gmidx][0],
+            (not DPM.shape[0] and float('inf')) or DPM[pmidx][0])


I found this very hard to parse. I wonder whether using if-else syntax might be simpler:

return ( float('inf') if DM.shape[0] == 0 else np.amin(DM), float('inf') if gmidx == -1 else DGM[gmidx][0], float('inf') if DPM.shape[0] == 0 else DPM[pmidx][0])

Agree - thanks for the snippet, that is much cleaner. I tried to emulate stability_margins, but had whittled it down somewhat. Will update.

…reference on disk/ellipse-based margin calculations.

josiahdelange · 2025-04-29T01:08:04Z

@murrayrm (and @slivingston) thanks for taking time to review. I think everything's addressed - but let me know of anything else!

josiahdelange added 22 commits January 11, 2025 01:12

Initial version of disk margin calculation and example/test script

f6dada2

Merge remote-tracking branch 'origin' into jdelange/disk-margins

bdb82a5

Comment updates: update margins.py header, clarify import exception h…

e46f824

…andler comment, fix typo in skew description of disk_margins docstring

More work in progress on disk margin calculation, adding new prototyp…

1e3af88

…e function to plot allowable gain/phase variations.

Add disk_margin_plot to subroutine list in comment header in margins.py

ba15789

Follow-on to ba15789, add disk_margin_plot to list of functions withi…

e47ae02

…n the margins module

Merge branch 'python-control:main' into jdelange/disk-margins

20686fe

More work in progress on disk_margin_plot. Corrected a typo/bug in the

f84221d

calculation of 'f', the bounding complex curve. Seems to look correct for balanced (skew = 0) case, still verifying the skewed equivalent.

Merge branch 'jdelange/disk-margins' of github.com:josiahdelange/pyth…

edf8040

…on-control into jdelange/disk-margins

Merge remote-tracking branch 'origin/main' into jdelange/disk-margins

9d55419

Further progress/debugging on disk margin calculation + plot utility

2cf1545

Merge remote-tracking branch 'origin/main' into jdelange/disk-margins

bbf37f0

Clean up docstring/code for disk_margin_plot

c3efe75

Clean up docstring/code for disk_margin_plot

63c8523

Remove debugging statements, update comments, add unit tests.

cffc3e5

Minor change to fix logic to find minimum across DGM, DPM numpy vectors

91517f9

Rename disk margin example, since unit tests are now written in contr…

86329e0

…ol/tests/margin_test.py

Remove unneeded dependencies from margins.py, used for debugging

d92fb20

Minor updates to docstrings

b2a2edc

Undo d92fb20

1f0ee52

Minor tweaks to plots in example script for readability

ba41e8c

Fix typo in disk_margin_plot.

14eb315

slivingston reviewed Apr 24, 2025

View reviewed changes

josiahdelange added 3 commits April 25, 2025 06:57

Fix mag2db import hack/workaround and trim down disk_margin docstring.

0bebc1d

Add input handling to disk_margin, clean up column width/comments

87714bd

Move disk_margin_plot out of the library into the example script

c17910f

josiahdelange marked this pull request as draft April 26, 2025 14:29

josiahdelange added 2 commits April 26, 2025 10:32

Recommended changes from the linter

5f34a7b

Follow-on to 5f34a7b

f0e2d74

josiahdelange added 4 commits April 26, 2025 11:09

Add disk_margins to function list

a5fcb91

Whittle down the docstring from disk_margins

077d538

Put more comments in the disk margin example, add example to document…

8f0c037

…ation

Fixing docstrings

ce80819

josiahdelange marked this pull request as ready for review April 26, 2025 16:03

josiahdelange added 2 commits April 26, 2025 14:49

Corrected expected values for 'no-slycot' condition in newly-added un…

397efab

…it tests

Attempt python-control#2 at 397efab, based on linter recommendation

cc02712

murrayrm reviewed Apr 28, 2025

View reviewed changes

josiahdelange added 3 commits April 28, 2025 20:04

Address @murrayrm review comments.

e8897f6

Update formatting per PEP8/@murrayrm review comments. Add additional …

579f24b

…reference on disk/ellipse-based margin calculations.

Follow-on to e8897f6: remove now-unnecessary import of importlib

fe79760

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Disk margin calculations #1146

Disk margin calculations #1146

josiahdelange commented Apr 24, 2025 •

edited

Loading

slivingston Apr 24, 2025

slivingston Apr 24, 2025

josiahdelange Apr 25, 2025

josiahdelange commented Apr 26, 2025

coveralls commented Apr 28, 2025 •

edited

Loading

murrayrm left a comment

murrayrm Apr 28, 2025

josiahdelange Apr 28, 2025

murrayrm Apr 28, 2025

murrayrm Apr 28, 2025

josiahdelange Apr 28, 2025

murrayrm Apr 28, 2025

josiahdelange Apr 28, 2025

murrayrm Apr 28, 2025

josiahdelange Apr 28, 2025

murrayrm Apr 28, 2025

murrayrm Apr 28, 2025

murrayrm Apr 28, 2025

murrayrm Apr 28, 2025

murrayrm Apr 28, 2025

josiahdelange Apr 28, 2025

josiahdelange commented Apr 29, 2025

Disk margin calculations #1146

Are you sure you want to change the base?

Disk margin calculations #1146

Conversation

josiahdelange commented Apr 24, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

josiahdelange commented Apr 26, 2025

coveralls commented Apr 28, 2025 • edited Loading

murrayrm left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

josiahdelange commented Apr 29, 2025

josiahdelange commented Apr 24, 2025 •

edited

Loading

coveralls commented Apr 28, 2025 •

edited

Loading