Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commitaa92b65

Browse files
Disk margin calculations (#1146)
* Initial version of disk margin calculation and example/test script* Comment updates: update margins.py header, clarify import exception handler comment, fix typo in skew description of disk_margins docstring* More work in progress on disk margin calculation, adding new prototype function to plot allowable gain/phase variations.* Add disk_margin_plot to subroutine list in comment header in margins.py* Follow-on toba15789, add disk_margin_plot to list of functions within the margins module* More work in progress on disk_margin_plot. Corrected a typo/bug in thecalculation of 'f', the bounding complex curve. Seems to look correct for balanced(skew = 0) case, still verifying the skewed equivalent.* Further progress/debugging on disk margin calculation + plot utility* Clean up docstring/code for disk_margin_plot* Clean up docstring/code for disk_margin_plot* Remove debugging statements, update comments, add unit tests.* Minor change to fix logic to find minimum across DGM, DPM numpy vectors* Rename disk margin example, since unit tests are now written in control/tests/margin_test.py* Remove unneeded dependencies from margins.py, used for debugging* Minor updates to docstrings* Undod92fb20* Minor tweaks to plots in example script for readability* Fix typo in disk_margin_plot.* Fix mag2db import hack/workaround and trim down disk_margin docstring.* Add input handling to disk_margin, clean up column width/comments* Move disk_margin_plot out of the library into the example script* Recommended changes from the linter* Follow-on to5f34a7b* Add disk_margins to function list* Whittle down the docstring from disk_margins* Put more comments in the disk margin example, add example to documentation* Fixing docstrings* Corrected expected values for 'no-slycot' condition in newly-added unit tests* Attempt#2 at397efab, based on linter recommendation* Address@murrayrm review comments.* Update formatting per PEP8/@murrayrm review comments. Add additional reference on disk/ellipse-based margin calculations.* Follow-on toe8897f6: remove now-unnecessary import of importlib* Update formatting per@murrayrm review comments* Remove temporarily-added string from docstring* Minor tweak to docstring to fit the word 'function' back into the description of skew = 0.0
1 parent7b06774 commitaa92b65

File tree

5 files changed

+814
-4
lines changed

5 files changed

+814
-4
lines changed

‎control/margins.py

Lines changed: 141 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -11,12 +11,17 @@
1111
importnumpyasnp
1212
importscipyassp
1313

14-
from .importfrdata,freqplot,xferfcn
14+
from .importfrdata,freqplot,xferfcn,statesp
1515
from .exceptionimportControlMIMONotImplemented
1616
from .iosysimportissiso
17+
from .ctrlutilimportmag2db
18+
try:
19+
fromslycotimportab13md
20+
exceptImportError:
21+
ab13md=None
1722

18-
__all__= ['stability_margins','phase_crossover_frequencies','margin']
19-
23+
__all__= ['stability_margins','phase_crossover_frequencies','margin',
24+
'disk_margins']
2025

2126
# private helper functions
2227
def_poly_iw(sys):
@@ -168,6 +173,7 @@ def fun(wdt):
168173

169174
returnz,w
170175

176+
171177
def_likely_numerical_inaccuracy(sys):
172178
# crude, conservative check for if
173179
# num(z)*num(1/z) << den(z)*den(1/z) for DT systems
@@ -517,3 +523,135 @@ def margin(*args):
517523
%len(args))
518524

519525
returnmargin[0],margin[1],margin[3],margin[4]
526+
527+
528+
defdisk_margins(L,omega,skew=0.0,returnall=False):
529+
"""Disk-based stability margins of loop transfer function.
530+
531+
Parameters
532+
----------
533+
L : `StateSpace` or `TransferFunction`
534+
Linear SISO or MIMO loop transfer function.
535+
omega : sequence of array_like
536+
1D array of (non-negative) frequencies (rad/s) at which
537+
to evaluate the disk-based stability margins.
538+
skew : float or array_like, optional
539+
skew parameter(s) for disk margin (default = 0.0).
540+
skew = 0.0 "balanced" sensitivity function 0.5*(S - T).
541+
skew = 1.0 sensitivity function S.
542+
skew = -1.0 complementary sensitivity function T.
543+
returnall : bool, optional
544+
If True, return frequency-dependent margins.
545+
If False (default), return worst-case (minimum) margins.
546+
547+
Returns
548+
-------
549+
DM : float or array_like
550+
Disk margin.
551+
DGM : float or array_like
552+
Disk-based gain margin.
553+
DPM : float or array_like
554+
Disk-based phase margin.
555+
556+
Example
557+
--------
558+
>> omega = np.logspace(-1, 3, 1001)
559+
>> P = control.ss([[0, 10], [-10, 0]], np.eye(2), [[1, 10],
560+
[-10, 1]], 0)
561+
>> K = control.ss([], [], [], [[1, -2], [0, 1]])
562+
>> L = P * K
563+
>> DM, DGM, DPM = control.disk_margins(L, omega, skew=0.0)
564+
"""
565+
566+
# First argument must be a system
567+
ifnotisinstance(L, (statesp.StateSpace,xferfcn.TransferFunction)):
568+
raiseValueError(
569+
"Loop gain must be state-space or transfer function object")
570+
571+
# Loop transfer function must be square
572+
ifstatesp.ss(L).B.shape[1]!=statesp.ss(L).C.shape[0]:
573+
raiseValueError("Loop gain must be square (n_inputs = n_outputs)")
574+
575+
# Need slycot if L is MIMO, for mu calculation
576+
ifnotL.issiso()andab13md==None:
577+
raiseControlMIMONotImplemented(
578+
"Need slycot to compute MIMO disk_margins")
579+
580+
# Get dimensions of feedback system
581+
num_loops=statesp.ss(L).C.shape[0]
582+
I=statesp.ss([], [], [],np.eye(num_loops))
583+
584+
# Loop sensitivity function
585+
S=I.feedback(L)
586+
587+
# Compute frequency response of the "balanced" (according
588+
# to the skew parameter "sigma") sensitivity function [1-2]
589+
ST=S+0.5* (skew-1)*I
590+
ST_mag,ST_phase,_=ST.frequency_response(omega)
591+
ST_jw= (ST_mag*np.exp(1j*ST_phase))
592+
ifnotL.issiso():
593+
ST_jw=ST_jw.transpose(2,0,1)
594+
595+
# Frequency-dependent complex disk margin, computed using
596+
# upper bound of the structured singular value, a.k.a. "mu",
597+
# of (S + (skew - I)/2).
598+
DM=np.zeros(omega.shape)
599+
DGM=np.zeros(omega.shape)
600+
DPM=np.zeros(omega.shape)
601+
foriiinrange(0,len(omega)):
602+
# Disk margin (a.k.a. "alpha") vs. frequency
603+
ifL.issiso()andab13md==None:
604+
# For the SISO case, the norm on (S + (skew - I)/2) is
605+
# unstructured, and can be computed as the magnitude
606+
# of the frequency response.
607+
DM[ii]=1.0/ST_mag[ii]
608+
else:
609+
# For the MIMO case, the norm on (S + (skew - I)/2)
610+
# assumes a single complex uncertainty block diagonal
611+
# uncertainty structure. AB13MD provides an upper bound
612+
# on this norm at the given frequency omega[ii].
613+
DM[ii]=1.0/ab13md(ST_jw[ii],np.array(num_loops* [1]),
614+
np.array(num_loops* [2]))[0]
615+
616+
# Disk-based gain margin (dB) and phase margin (deg)
617+
withnp.errstate(divide='ignore',invalid='ignore'):
618+
# Real-axis intercepts with the disk
619+
gamma_min= (1-0.5*DM[ii]* (1-skew)) \
620+
/ (1+0.5*DM[ii]* (1+skew))
621+
gamma_max= (1+0.5*DM[ii]* (1-skew)) \
622+
/ (1-0.5*DM[ii]* (1+skew))
623+
624+
# Gain margin (dB)
625+
DGM[ii]=mag2db(np.minimum(1/gamma_min,gamma_max))
626+
ifnp.isnan(DGM[ii]):
627+
DGM[ii]=float('inf')
628+
629+
# Phase margin (deg)
630+
ifnp.isinf(gamma_max):
631+
DPM[ii]=90.0
632+
else:
633+
DPM[ii]= (1+gamma_min*gamma_max) \
634+
/ (gamma_min+gamma_max)
635+
ifabs(DPM[ii])>=1.0:
636+
DPM[ii]=float('Inf')
637+
else:
638+
DPM[ii]=np.rad2deg(np.arccos(DPM[ii]))
639+
640+
ifreturnall:
641+
# Frequency-dependent disk margin, gain margin and phase margin
642+
returnDM,DGM,DPM
643+
else:
644+
# Worst-case disk margin, gain margin and phase margin
645+
ifDGM.shape[0]andnotnp.isinf(DGM).all():
646+
withnp.errstate(all='ignore'):
647+
gmidx=np.where(DGM==np.min(DGM))
648+
else:
649+
gmidx=-1
650+
651+
ifDPM.shape[0]:
652+
pmidx=np.where(DPM==np.min(DPM))
653+
654+
return (
655+
float('inf')ifDM.shape[0]==0elsenp.amin(DM),
656+
float('inf')ifgmidx==-1elseDGM[gmidx][0],
657+
float('inf')ifDPM.shape[0]==0elseDPM[pmidx][0])

‎control/tests/margin_test.py

Lines changed: 105 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,8 @@
1414

1515
fromcontrolimportControlMIMONotImplemented,FrequencyResponseData, \
1616
StateSpace,TransferFunction,margin,phase_crossover_frequencies, \
17-
stability_margins
17+
stability_margins,disk_margins,tf,ss
18+
fromcontrol.exceptionimportslycot_check
1819

1920
s=TransferFunction.s
2021

@@ -372,3 +373,106 @@ def test_stability_margins_discrete(cnum, cden, dt,
372373
else:
373374
out=stability_margins(tf)
374375
assert_allclose(out,ref,rtol=rtol)
376+
377+
deftest_siso_disk_margin():
378+
# Frequencies of interest
379+
omega=np.logspace(-1,2,1001)
380+
381+
# Loop transfer function
382+
L=tf(25, [1,10,10,10])
383+
384+
# Balanced (S - T) disk-based stability margins
385+
DM,DGM,DPM=disk_margins(L,omega,skew=0.0)
386+
assert_allclose([DM], [0.46],atol=0.1)# disk margin of 0.46
387+
assert_allclose([DGM], [4.05],atol=0.1)# disk-based gain margin of 4.05 dB
388+
assert_allclose([DPM], [25.8],atol=0.1)# disk-based phase margin of 25.8 deg
389+
390+
# For SISO systems, the S-based (S) disk margin should match the third output
391+
# of existing library "stability_margins", i.e., minimum distance from the
392+
# Nyquist plot to -1.
393+
_,_,SM=stability_margins(L)[:3]
394+
DM=disk_margins(L,omega,skew=1.0)[0]
395+
assert_allclose([DM], [SM],atol=0.01)
396+
397+
deftest_mimo_disk_margin():
398+
# Frequencies of interest
399+
omega=np.logspace(-1,3,1001)
400+
401+
# Loop transfer gain
402+
P=ss([[0,10], [-10,0]],np.eye(2), [[1,10], [-10,1]],0)# plant
403+
K=ss([], [], [], [[1,-2], [0,1]])# controller
404+
Lo=P*K# loop transfer function, broken at plant output
405+
Li=K*P# loop transfer function, broken at plant input
406+
407+
ifslycot_check():
408+
# Balanced (S - T) disk-based stability margins at plant output
409+
DMo,DGMo,DPMo=disk_margins(Lo,omega,skew=0.0)
410+
assert_allclose([DMo], [0.3754],atol=0.1)# disk margin of 0.3754
411+
assert_allclose([DGMo], [3.3],atol=0.1)# disk-based gain margin of 3.3 dB
412+
assert_allclose([DPMo], [21.26],atol=0.1)# disk-based phase margin of 21.26 deg
413+
414+
# Balanced (S - T) disk-based stability margins at plant input
415+
DMi,DGMi,DPMi=disk_margins(Li,omega,skew=0.0)
416+
assert_allclose([DMi], [0.3754],atol=0.1)# disk margin of 0.3754
417+
assert_allclose([DGMi], [3.3],atol=0.1)# disk-based gain margin of 3.3 dB
418+
assert_allclose([DPMi], [21.26],atol=0.1)# disk-based phase margin of 21.26 deg
419+
else:
420+
# Slycot not installed. Should throw exception.
421+
withpytest.raises(ControlMIMONotImplemented,\
422+
match="Need slycot to compute MIMO disk_margins"):
423+
DMo,DGMo,DPMo=disk_margins(Lo,omega,skew=0.0)
424+
425+
deftest_siso_disk_margin_return_all():
426+
# Frequencies of interest
427+
omega=np.logspace(-1,2,1001)
428+
429+
# Loop transfer function
430+
L=tf(25, [1,10,10,10])
431+
432+
# Balanced (S - T) disk-based stability margins
433+
DM,DGM,DPM=disk_margins(L,omega,skew=0.0,returnall=True)
434+
assert_allclose([omega[np.argmin(DM)]], [1.94],\
435+
atol=0.01)# sensitivity peak at 1.94 rad/s
436+
assert_allclose([min(DM)], [0.46],atol=0.1)# disk margin of 0.46
437+
assert_allclose([DGM[np.argmin(DM)]], [4.05],\
438+
atol=0.1)# disk-based gain margin of 4.05 dB
439+
assert_allclose([DPM[np.argmin(DM)]], [25.8],\
440+
atol=0.1)# disk-based phase margin of 25.8 deg
441+
442+
deftest_mimo_disk_margin_return_all():
443+
# Frequencies of interest
444+
omega=np.logspace(-1,3,1001)
445+
446+
# Loop transfer gain
447+
P=ss([[0,10], [-10,0]],np.eye(2),\
448+
[[1,10], [-10,1]],0)# plant
449+
K=ss([], [], [], [[1,-2], [0,1]])# controller
450+
Lo=P*K# loop transfer function, broken at plant output
451+
Li=K*P# loop transfer function, broken at plant input
452+
453+
ifslycot_check():
454+
# Balanced (S - T) disk-based stability margins at plant output
455+
DMo,DGMo,DPMo=disk_margins(Lo,omega,skew=0.0,returnall=True)
456+
assert_allclose([omega[np.argmin(DMo)]], [omega[0]],\
457+
atol=0.01)# sensitivity peak at 0 rad/s (or smallest provided)
458+
assert_allclose([min(DMo)], [0.3754],atol=0.1)# disk margin of 0.3754
459+
assert_allclose([DGMo[np.argmin(DMo)]], [3.3],\
460+
atol=0.1)# disk-based gain margin of 3.3 dB
461+
assert_allclose([DPMo[np.argmin(DMo)]], [21.26],\
462+
atol=0.1)# disk-based phase margin of 21.26 deg
463+
464+
# Balanced (S - T) disk-based stability margins at plant input
465+
DMi,DGMi,DPMi=disk_margins(Li,omega,skew=0.0,returnall=True)
466+
assert_allclose([omega[np.argmin(DMi)]], [omega[0]],\
467+
atol=0.01)# sensitivity peak at 0 rad/s (or smallest provided)
468+
assert_allclose([min(DMi)], [0.3754],\
469+
atol=0.1)# disk margin of 0.3754
470+
assert_allclose([DGMi[np.argmin(DMi)]], [3.3],\
471+
atol=0.1)# disk-based gain margin of 3.3 dB
472+
assert_allclose([DPMi[np.argmin(DMi)]], [21.26],\
473+
atol=0.1)# disk-based phase margin of 21.26 deg
474+
else:
475+
# Slycot not installed. Should throw exception.
476+
withpytest.raises(ControlMIMONotImplemented,\
477+
match="Need slycot to compute MIMO disk_margins"):
478+
DMo,DGMo,DPMo=disk_margins(Lo,omega,skew=0.0,returnall=True)

‎doc/examples/disk_margins.rst

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
Disk margin example
2+
------------------------------------------
3+
4+
This example demonstrates the use of the `disk_margins` routine
5+
to compute robust stability margins for a feedback system, i.e.,
6+
variation in gain and phase one or more loops. The SISO examples
7+
are drawn from the published paper and the MIMO example is the
8+
"spinning satellite" example from the MathWorks documentation.
9+
10+
Code
11+
....
12+
..literalinclude::disk_margins.py
13+
:language: python
14+
:linenos:
15+
16+
Notes
17+
.....
18+
1. The environment variable `PYCONTROL_TEST_EXAMPLES` is used for
19+
testing to turn off plotting of the outputs.

‎doc/functions.rst

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,7 @@ Frequency domain analysis:
142142

143143
bandwidth
144144
dcgain
145+
disk_margins
145146
linfnorm
146147
margin
147148
stability_margins

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp