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

Disk margin calculations#1146

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to ourterms of service andprivacy statement. We’ll occasionally send you account related emails.

Already on GitHub?Sign in to your account

Merged

Conversation

josiahdelange
Copy link
Contributor

@josiahdelangejosiahdelange commentedApr 24, 2025•
edited
Loading

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 ofpython-control andslycot. The code should work both for SISO and MIMO systems, the latter of which requiresslycot (for SLICOT'sAB13MD) to compute the upper bound of$\mu$ at each discrete frequency. The functiondisk_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 thepublished paper and the "spinning satellite" MIMOexample from the MathWorks documentation. I also confirmed SISO disk margins match the relevant output of existing functionstability_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 incontrol/margins.py. Example usage (seeexamples/disk_margin.py):

importos,sys,mathimportnumpyasnpimportcontrolimportmathimportmatplotlibimportmatplotlib.pyplotaspltfromwarningsimportwarnimportnumpyasnpimportscipyassp# Frequencies of interestomega=np.logspace(-1,3,1001)# Laplace variables=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.htmlP=control.ss([[0,10],[-10,0]],np.eye(2), [[1,10], [-10,1]], [[0,0],[0,0]])# plantC=control.ss([],[],[], [[1,-2], [0,1]])# controllerL=P*C# output loop gainprint(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 dBPM = 21.30709883385843 degmin(GM) = 2.732761944304522 dBmin(PM) = 21.30709883385843 deg------------- Complementary sensitivity function (T) -------------min(DM) = 0.3683035923739884 (omega = 0.1)GM = 2.7236493433653735 dBPM = 21.223368586091564 degmin(GM) = 2.7236493433653735 dBmin(PM) = 21.223368586091564 deg------------- Balanced sensitivity function (S - T) -------------min(DM) = 0.3769872636944355 (omega = 0.1)GM = 3.3140985364257256 dBPM = 21.349285542574695 degmin(GM) = 3.3140985364257256 dBmin(PM) = 21.349285542574695 deg

Figure_3

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 functionplot_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])

Figure_30

…andler comment, fix typo in skew description of disk_margins docstring
…e function to plot allowable gain/phase variations.
calculation of 'f', the bounding complex curve.  Seems to look correct for balanced(skew = 0) case, still verifying the skewed equivalent.
@josiahdelangejosiahdelange marked this pull request as draftApril 26, 2025 14:29
@josiahdelange
Copy link
ContributorAuthor

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 functionplot_allowable_region inexamples/disk_margins.py.

@josiahdelangejosiahdelange marked this pull request as ready for reviewApril 26, 2025 16:03
@coveralls
Copy link

coveralls commentedApr 28, 2025•
edited
Loading

Coverage Status

coverage: 94.745% (-0.01%) from 94.759%
when pullingbb06c9e on josiahdelange:jdelange/disk-margins
into34c6d59 on python-control:main.

Copy link
Member

@murrayrmmurrayrm left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others.Learn more.

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

@josiahdelange
Copy link
ContributorAuthor

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

Copy link
Member

@murrayrmmurrayrm left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others.Learn more.

Found a few more style issues.

@josiahdelange
Copy link
ContributorAuthor

@murrayrm thanks and latest formatting issues are now fixed. The one failing check "Slycot from source / build-linux (pull_request)" does not seem to be the disk margin tests.

@murrayrm
Copy link
Member

Agreed. Tracking the failed unit test in#1161. I'll go ahead and merge this since the fail check is unrelated.

@murrayrmmurrayrm merged commitaa92b65 intopython-control:mainJun 25, 2025
23 of 24 checks passed
Sign up for freeto join this conversation on GitHub. Already have an account?Sign in to comment
Reviewers

@murrayrmmurrayrmmurrayrm left review comments

@slivingstonslivingstonslivingston left review comments

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

4 participants
@josiahdelange@coveralls@murrayrm@slivingston
[8]ページ先頭
©2009-2025 Movatter.jp