Disk margin example
This example demonstrates the use of the disk_margins routine
to compute robust stability margins for a feedback system, i.e.,
variation in gain and phase one or more loops. The SISO examples
are drawn from the published paper and the MIMO example is the
“spinning satellite” example from the MathWorks documentation.
Code
1"""disk_margins.py
2
3Demonstrate disk-based stability margin calculations.
4
5References:
6[1] Blight, James D., R. Lane Dailey, and Dagfinn Gangsaas. “Practical
7 Control Law Design for Aircraft Using Multivariable Techniques.”
8 International Journal of Control 59, no. 1 (January 1994): 93-137.
9 https://doi.org/10.1080/00207179408923071.
10
11[2] Seiler, Peter, Andrew Packard, and Pascal Gahinet. “An Introduction
12 to Disk Margins [Lecture Notes].” IEEE Control Systems Magazine 40,
13 no. 5 (October 2020): 78-95.
14
15[3] P. Benner, V. Mehrmann, V. Sima, S. Van Huffel, and A. Varga, "SLICOT
16 - A Subroutine Library in Systems and Control Theory", Applied and
17 Computational Control, Signals, and Circuits (Birkhauser), Vol. 1, Ch.
18 10, pp. 505-546, 1999.
19
20[4] S. Van Huffel, V. Sima, A. Varga, S. Hammarling, and F. Delebecque,
21 "Development of High Performance Numerical Software for Control", IEEE
22 Control Systems Magazine, Vol. 24, Nr. 1, Feb., pp. 60-76, 2004.
23
24[5] Deodhare, G., & Patel, V. (1998, August). A "Modern" Look at Gain
25 and Phase Margins: An H-Infinity/mu Approach. In Guidance, Navigation,
26 and Control Conference and Exhibit (p. 4134).
27"""
28
29import os
30import control
31import matplotlib.pyplot as plt
32import numpy as np
33
34def plot_allowable_region(alpha_max, skew, ax=None):
35 """Plot region of allowable gain/phase variation, given worst-case disk margin.
36
37 Parameters
38 ----------
39 alpha_max : float (scalar or list)
40 worst-case disk margin(s) across all frequencies. May be a scalar or list.
41 skew : float (scalar or list)
42 skew parameter(s) for disk margin calculation.
43 skew=0 uses the "balanced" sensitivity function 0.5*(S - T)
44 skew=1 uses the sensitivity function S
45 skew=-1 uses the complementary sensitivity function T
46 ax : axes to plot bounding curve(s) onto
47
48 Returns
49 -------
50 DM : ndarray
51 1D array of frequency-dependent disk margins. DM is the same
52 size as "omega" parameter.
53 GM : ndarray
54 1D array of frequency-dependent disk-based gain margins, in dB.
55 GM is the same size as "omega" parameter.
56 PM : ndarray
57 1D array of frequency-dependent disk-based phase margins, in deg.
58 PM is the same size as "omega" parameter.
59 """
60
61 # Create axis if needed
62 if ax is None:
63 ax = plt.gca()
64
65 # Allow scalar or vector arguments (to overlay plots)
66 if np.isscalar(alpha_max):
67 alpha_max = np.asarray([alpha_max])
68 else:
69 alpha_max = np.asarray(alpha_max)
70
71 if np.isscalar(skew):
72 skew=np.asarray([skew])
73 else:
74 skew=np.asarray(skew)
75
76 # Add a plot for each (alpha, skew) pair present
77 theta = np.linspace(0, np.pi, 500)
78 legend_list = []
79 for ii in range(0, skew.shape[0]):
80 legend_str = "$\\sigma$ = %.1f, $\\alpha_{max}$ = %.2f" %(\
81 skew[ii], alpha_max[ii])
82 legend_list.append(legend_str)
83
84 # Complex bounding curve of stable gain/phase variations
85 f = (2 + alpha_max[ii] * (1 - skew[ii]) * np.exp(1j * theta))\
86 /(2 - alpha_max[ii] * (1 + skew[ii]) * np.exp(1j * theta))
87
88 # Allowable combined gain/phase variations
89 gamma_dB = control.ctrlutil.mag2db(np.abs(f)) # gain margin (dB)
90 phi_deg = np.rad2deg(np.angle(f)) # phase margin (deg)
91
92 # Plot the allowable combined gain/phase variations
93 out = ax.plot(gamma_dB, phi_deg, alpha=0.25, label='_nolegend_')
94 ax.fill_between(ax.lines[ii].get_xydata()[:,0],\
95 ax.lines[ii].get_xydata()[:,1], alpha=0.25)
96
97 plt.ylabel('Phase Variation (deg)')
98 plt.xlabel('Gain Variation (dB)')
99 plt.title('Range of Gain and Phase Variations')
100 plt.legend(legend_list)
101 plt.grid()
102 plt.tight_layout()
103
104 return out
105
106def test_siso1():
107 #
108 # Disk-based stability margins for example
109 # SISO loop transfer function(s)
110 #
111
112 # Frequencies of interest
113 omega = np.logspace(-1, 2, 1001)
114
115 # Loop transfer gain
116 L = control.tf(25, [1, 10, 10, 10])
117
118 print("------------- Python control built-in (S) -------------")
119 GM_, PM_, SM_ = control.stability_margins(L)[:3] # python-control default (S-based...?)
120 print(f"SM_ = {SM_}")
121 print(f"GM_ = {GM_} dB")
122 print(f"PM_ = {PM_} deg\n")
123
124 print("------------- Sensitivity function (S) -------------")
125 DM, GM, PM = control.disk_margins(L, omega, skew=1.0, returnall=True) # S-based (S)
126 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
127 print(f"GM = {GM[np.argmin(DM)]} dB")
128 print(f"PM = {PM[np.argmin(DM)]} deg")
129 print(f"min(GM) = {min(GM)} dB")
130 print(f"min(PM) = {min(PM)} deg\n")
131
132 plt.figure(1)
133 plt.subplot(3, 3, 1)
134 plt.semilogx(omega, DM, label='$\\alpha$')
135 plt.ylabel('Disk Margin (abs)')
136 plt.legend()
137 plt.title('S-Based Margins')
138 plt.grid()
139 plt.xlim([omega[0], omega[-1]])
140 plt.ylim([0, 2])
141
142 plt.figure(1)
143 plt.subplot(3, 3, 4)
144 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
145 plt.ylabel('Gain Margin (dB)')
146 plt.legend()
147 #plt.title('Gain-Only Margin')
148 plt.grid()
149 plt.xlim([omega[0], omega[-1]])
150 plt.ylim([0, 40])
151
152 plt.figure(1)
153 plt.subplot(3, 3, 7)
154 plt.semilogx(omega, PM, label='$\\phi_{m}$')
155 plt.ylabel('Phase Margin (deg)')
156 plt.legend()
157 #plt.title('Phase-Only Margin')
158 plt.grid()
159 plt.xlim([omega[0], omega[-1]])
160 plt.ylim([0, 90])
161 plt.xlabel('Frequency (rad/s)')
162
163 print("------------- Complementary sensitivity function (T) -------------")
164 DM, GM, PM = control.disk_margins(L, omega, skew=-1.0, returnall=True) # T-based (T)
165 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
166 print(f"GM = {GM[np.argmin(DM)]} dB")
167 print(f"PM = {PM[np.argmin(DM)]} deg")
168 print(f"min(GM) = {min(GM)} dB")
169 print(f"min(PM) = {min(PM)} deg\n")
170
171 plt.figure(1)
172 plt.subplot(3, 3, 2)
173 plt.semilogx(omega, DM, label='$\\alpha$')
174 plt.ylabel('Disk Margin (abs)')
175 plt.legend()
176 plt.title('T_Based Margins')
177 plt.grid()
178 plt.xlim([omega[0], omega[-1]])
179 plt.ylim([0, 2])
180
181 plt.figure(1)
182 plt.subplot(3, 3, 5)
183 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
184 plt.ylabel('Gain Margin (dB)')
185 plt.legend()
186 #plt.title('Gain-Only Margin')
187 plt.grid()
188 plt.xlim([omega[0], omega[-1]])
189 plt.ylim([0, 40])
190
191 plt.figure(1)
192 plt.subplot(3, 3, 8)
193 plt.semilogx(omega, PM, label='$\\phi_{m}$')
194 plt.ylabel('Phase Margin (deg)')
195 plt.legend()
196 #plt.title('Phase-Only Margin')
197 plt.grid()
198 plt.xlim([omega[0], omega[-1]])
199 plt.ylim([0, 90])
200 plt.xlabel('Frequency (rad/s)')
201
202 print("------------- Balanced sensitivity function (S - T) -------------")
203 DM, GM, PM = control.disk_margins(L, omega, skew=0.0, returnall=True) # balanced (S - T)
204 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
205 print(f"GM = {GM[np.argmin(DM)]} dB")
206 print(f"PM = {PM[np.argmin(DM)]} deg")
207 print(f"min(GM) = {min(GM)} dB")
208 print(f"min(PM) = {min(PM)} deg\n")
209
210 plt.figure(1)
211 plt.subplot(3, 3, 3)
212 plt.semilogx(omega, DM, label='$\\alpha$')
213 plt.ylabel('Disk Margin (abs)')
214 plt.legend()
215 plt.title('Balanced Margins')
216 plt.grid()
217 plt.xlim([omega[0], omega[-1]])
218 plt.ylim([0, 2])
219
220 plt.figure(1)
221 plt.subplot(3, 3, 6)
222 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
223 plt.ylabel('Gain Margin (dB)')
224 plt.legend()
225 #plt.title('Gain-Only Margin')
226 plt.grid()
227 plt.xlim([omega[0], omega[-1]])
228 plt.ylim([0, 40])
229
230 plt.figure(1)
231 plt.subplot(3, 3, 9)
232 plt.semilogx(omega, PM, label='$\\phi_{m}$')
233 plt.ylabel('Phase Margin (deg)')
234 plt.legend()
235 #plt.title('Phase-Only Margin')
236 plt.grid()
237 plt.xlim([omega[0], omega[-1]])
238 plt.ylim([0, 90])
239 plt.xlabel('Frequency (rad/s)')
240
241 # Disk margin plot of admissible gain/phase variations for which
242 DM_plot = []
243 DM_plot.append(control.disk_margins(L, omega, skew=-2.0)[0])
244 DM_plot.append(control.disk_margins(L, omega, skew=0.0)[0])
245 DM_plot.append(control.disk_margins(L, omega, skew=2.0)[0])
246 plt.figure(10); plt.clf()
247 plot_allowable_region(DM_plot, skew=[-2.0, 0.0, 2.0])
248
249 return
250
251def test_siso2():
252 #
253 # Disk-based stability margins for example
254 # SISO loop transfer function(s)
255 #
256
257 # Frequencies of interest
258 omega = np.logspace(-1, 2, 1001)
259
260 # Laplace variable
261 s = control.tf('s')
262
263 # Loop transfer gain
264 L = (6.25 * (s + 3) * (s + 5)) / (s * (s + 1)**2 * (s**2 + 0.18 * s + 100))
265
266 print("------------- Python control built-in (S) -------------")
267 GM_, PM_, SM_ = control.stability_margins(L)[:3] # python-control default (S-based...?)
268 print(f"SM_ = {SM_}")
269 print(f"GM_ = {GM_} dB")
270 print(f"PM_ = {PM_} deg\n")
271
272 print("------------- Sensitivity function (S) -------------")
273 DM, GM, PM = control.disk_margins(L, omega, skew=1.0, returnall=True) # S-based (S)
274 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
275 print(f"GM = {GM[np.argmin(DM)]} dB")
276 print(f"PM = {PM[np.argmin(DM)]} deg")
277 print(f"min(GM) = {min(GM)} dB")
278 print(f"min(PM) = {min(PM)} deg\n")
279
280 plt.figure(2)
281 plt.subplot(3, 3, 1)
282 plt.semilogx(omega, DM, label='$\\alpha$')
283 plt.ylabel('Disk Margin (abs)')
284 plt.legend()
285 plt.title('S-Based Margins')
286 plt.grid()
287 plt.xlim([omega[0], omega[-1]])
288 plt.ylim([0, 2])
289
290 plt.figure(2)
291 plt.subplot(3, 3, 4)
292 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
293 plt.ylabel('Gain Margin (dB)')
294 plt.legend()
295 #plt.title('Gain-Only Margin')
296 plt.grid()
297 plt.xlim([omega[0], omega[-1]])
298 plt.ylim([0, 40])
299
300 plt.figure(2)
301 plt.subplot(3, 3, 7)
302 plt.semilogx(omega, PM, label='$\\phi_{m}$')
303 plt.ylabel('Phase Margin (deg)')
304 plt.legend()
305 #plt.title('Phase-Only Margin')
306 plt.grid()
307 plt.xlim([omega[0], omega[-1]])
308 plt.ylim([0, 90])
309 plt.xlabel('Frequency (rad/s)')
310
311 print("------------- Complementary sensitivity function (T) -------------")
312 DM, GM, PM = control.disk_margins(L, omega, skew=-1.0, returnall=True) # T-based (T)
313 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
314 print(f"GM = {GM[np.argmin(DM)]} dB")
315 print(f"PM = {PM[np.argmin(DM)]} deg")
316 print(f"min(GM) = {min(GM)} dB")
317 print(f"min(PM) = {min(PM)} deg\n")
318
319 plt.figure(2)
320 plt.subplot(3, 3, 2)
321 plt.semilogx(omega, DM, label='$\\alpha$')
322 plt.ylabel('Disk Margin (abs)')
323 plt.legend()
324 plt.title('T-Based Margins')
325 plt.grid()
326 plt.xlim([omega[0], omega[-1]])
327 plt.ylim([0, 2])
328
329 plt.figure(2)
330 plt.subplot(3, 3, 5)
331 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
332 plt.ylabel('Gain Margin (dB)')
333 plt.legend()
334 #plt.title('Gain-Only Margin')
335 plt.grid()
336 plt.xlim([omega[0], omega[-1]])
337 plt.ylim([0, 40])
338
339 plt.figure(2)
340 plt.subplot(3, 3, 8)
341 plt.semilogx(omega, PM, label='$\\phi_{m}$')
342 plt.ylabel('Phase Margin (deg)')
343 plt.legend()
344 #plt.title('Phase-Only Margin')
345 plt.grid()
346 plt.xlim([omega[0], omega[-1]])
347 plt.ylim([0, 90])
348 plt.xlabel('Frequency (rad/s)')
349
350 print("------------- Balanced sensitivity function (S - T) -------------")
351 DM, GM, PM = control.disk_margins(L, omega, skew=0.0, returnall=True) # balanced (S - T)
352 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
353 print(f"GM = {GM[np.argmin(DM)]} dB")
354 print(f"PM = {PM[np.argmin(DM)]} deg")
355 print(f"min(GM) = {min(GM)} dB")
356 print(f"min(PM) = {min(PM)} deg\n")
357
358 plt.figure(2)
359 plt.subplot(3, 3, 3)
360 plt.semilogx(omega, DM, label='$\\alpha$')
361 plt.ylabel('Disk Margin (abs)')
362 plt.legend()
363 plt.title('Balanced Margins')
364 plt.grid()
365 plt.xlim([omega[0], omega[-1]])
366 plt.ylim([0, 2])
367
368 plt.figure(2)
369 plt.subplot(3, 3, 6)
370 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
371 plt.ylabel('Gain Margin (dB)')
372 plt.legend()
373 #plt.title('Gain-Only Margin')
374 plt.grid()
375 plt.xlim([omega[0], omega[-1]])
376 plt.ylim([0, 40])
377
378 plt.figure(2)
379 plt.subplot(3, 3, 9)
380 plt.semilogx(omega, PM, label='$\\phi_{m}$')
381 plt.ylabel('Phase Margin (deg)')
382 plt.legend()
383 #plt.title('Phase-Only Margin')
384 plt.grid()
385 plt.xlim([omega[0], omega[-1]])
386 plt.ylim([0, 90])
387 plt.xlabel('Frequency (rad/s)')
388
389 # Disk margin plot of admissible gain/phase variations for which
390 # the feedback loop still remains stable, for each skew parameter
391 DM_plot = []
392 DM_plot.append(control.disk_margins(L, omega, skew=-1.0)[0]) # T-based (T)
393 DM_plot.append(control.disk_margins(L, omega, skew=0.0)[0]) # balanced (S - T)
394 DM_plot.append(control.disk_margins(L, omega, skew=1.0)[0]) # S-based (S)
395 plt.figure(20)
396 plot_allowable_region(DM_plot, skew=[-1.0, 0.0, 1.0])
397
398 return
399
400def test_mimo():
401 #
402 # Disk-based stability margins for example
403 # MIMO loop transfer function(s)
404 #
405
406 # Frequencies of interest
407 omega = np.logspace(-1, 3, 1001)
408
409 # Loop transfer gain
410 P = control.ss([[0, 10],[-10, 0]], np.eye(2), [[1, 10], [-10, 1]], 0) # plant
411 K = control.ss([], [], [], [[1, -2], [0, 1]]) # controller
412 L = P * K # loop gain
413
414 print("------------- Sensitivity function (S) -------------")
415 DM, GM, PM = control.disk_margins(L, omega, skew=1.0, returnall=True) # S-based (S)
416 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
417 print(f"GM = {GM[np.argmin(DM)]} dB")
418 print(f"PM = {PM[np.argmin(DM)]} deg")
419 print(f"min(GM) = {min(GM)} dB")
420 print(f"min(PM) = {min(PM)} deg\n")
421
422 plt.figure(3)
423 plt.subplot(3, 3, 1)
424 plt.semilogx(omega, DM, label='$\\alpha$')
425 plt.ylabel('Disk Margin (abs)')
426 plt.legend()
427 plt.title('S-Based Margins')
428 plt.grid()
429 plt.xlim([omega[0], omega[-1]])
430 plt.ylim([0, 2])
431
432 plt.figure(3)
433 plt.subplot(3, 3, 4)
434 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
435 plt.ylabel('Gain Margin (dB)')
436 plt.legend()
437 #plt.title('Gain-Only Margin')
438 plt.grid()
439 plt.xlim([omega[0], omega[-1]])
440 plt.ylim([0, 40])
441
442 plt.figure(3)
443 plt.subplot(3, 3, 7)
444 plt.semilogx(omega, PM, label='$\\phi_{m}$')
445 plt.ylabel('Phase Margin (deg)')
446 plt.legend()
447 #plt.title('Phase-Only Margin')
448 plt.grid()
449 plt.xlim([omega[0], omega[-1]])
450 plt.ylim([0, 90])
451 plt.xlabel('Frequency (rad/s)')
452
453 print("------------- Complementary sensitivity function (T) -------------")
454 DM, GM, PM = control.disk_margins(L, omega, skew=-1.0, returnall=True) # T-based (T)
455 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
456 print(f"GM = {GM[np.argmin(DM)]} dB")
457 print(f"PM = {PM[np.argmin(DM)]} deg")
458 print(f"min(GM) = {min(GM)} dB")
459 print(f"min(PM) = {min(PM)} deg\n")
460
461 plt.figure(3)
462 plt.subplot(3, 3, 2)
463 plt.semilogx(omega, DM, label='$\\alpha$')
464 plt.ylabel('Disk Margin (abs)')
465 plt.legend()
466 plt.title('T-Based Margins')
467 plt.grid()
468 plt.xlim([omega[0], omega[-1]])
469 plt.ylim([0, 2])
470
471 plt.figure(3)
472 plt.subplot(3, 3, 5)
473 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
474 plt.ylabel('Gain Margin (dB)')
475 plt.legend()
476 #plt.title('Gain-Only Margin')
477 plt.grid()
478 plt.xlim([omega[0], omega[-1]])
479 plt.ylim([0, 40])
480
481 plt.figure(3)
482 plt.subplot(3, 3, 8)
483 plt.semilogx(omega, PM, label='$\\phi_{m}$')
484 plt.ylabel('Phase Margin (deg)')
485 plt.legend()
486 #plt.title('Phase-Only Margin')
487 plt.grid()
488 plt.xlim([omega[0], omega[-1]])
489 plt.ylim([0, 90])
490 plt.xlabel('Frequency (rad/s)')
491
492 print("------------- Balanced sensitivity function (S - T) -------------")
493 DM, GM, PM = control.disk_margins(L, omega, skew=0.0, returnall=True) # balanced (S - T)
494 print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})")
495 print(f"GM = {GM[np.argmin(DM)]} dB")
496 print(f"PM = {PM[np.argmin(DM)]} deg")
497 print(f"min(GM) = {min(GM)} dB")
498 print(f"min(PM) = {min(PM)} deg\n")
499
500 plt.figure(3)
501 plt.subplot(3, 3, 3)
502 plt.semilogx(omega, DM, label='$\\alpha$')
503 plt.ylabel('Disk Margin (abs)')
504 plt.legend()
505 plt.title('Balanced Margins')
506 plt.grid()
507 plt.xlim([omega[0], omega[-1]])
508 plt.ylim([0, 2])
509
510 plt.figure(3)
511 plt.subplot(3, 3, 6)
512 plt.semilogx(omega, GM, label='$\\gamma_{m}$')
513 plt.ylabel('Gain Margin (dB)')
514 plt.legend()
515 #plt.title('Gain-Only Margin')
516 plt.grid()
517 plt.xlim([omega[0], omega[-1]])
518 plt.ylim([0, 40])
519
520 plt.figure(3)
521 plt.subplot(3, 3, 9)
522 plt.semilogx(omega, PM, label='$\\phi_{m}$')
523 plt.ylabel('Phase Margin (deg)')
524 plt.legend()
525 #plt.title('Phase-Only Margin')
526 plt.grid()
527 plt.xlim([omega[0], omega[-1]])
528 plt.ylim([0, 90])
529 plt.xlabel('Frequency (rad/s)')
530
531 # Disk margin plot of admissible gain/phase variations for which
532 # the feedback loop still remains stable, for each skew parameter
533 DM_plot = []
534 DM_plot.append(control.disk_margins(L, omega, skew=-1.0)[0]) # T-based (T)
535 DM_plot.append(control.disk_margins(L, omega, skew=0.0)[0]) # balanced (S - T)
536 DM_plot.append(control.disk_margins(L, omega, skew=1.0)[0]) # S-based (S)
537 plt.figure(30)
538 plot_allowable_region(DM_plot, skew=[-1.0, 0.0, 1.0])
539
540 return
541
542if __name__ == '__main__':
543 #test_siso1()
544 #test_siso2()
545 test_mimo()
546 if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
547 #plt.tight_layout()
548 plt.show()
Notes
1. The environment variable PYCONTROL_TEST_EXAMPLES is used for
testing to turn off plotting of the outputs.