10000 Merge pull request #566 from sawyerbfuller/frd-margin · python-control/python-control@b41574f · GitHub
[go: up one dir, main page]
Skip to content

Commit b41574f

Browse files
bnavigatorbnavigator
authored
Merge pull request #566 from sawyerbfuller/frd-margin
new method keyword in stability_margins and auto-fallback back to FRD ..
2 parents 258c9c2 + 72cb23f commit b41574f

File tree

3 files changed

+73
-11
lines changed

3 files changed

+73
-11
lines changed

‎control/freqplot.py

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@
8888

8989
def bode_plot(syslist, omega=None,
9090
plot=True, omega_limits=None, omega_num=None,
91-
margins=None, *args, **kwargs):
91+
margins=None, method='best', *args, **kwargs):
9292
"""Bode plot for a system
9393
9494
Plots a Bode plot for the system over a (optional) frequency range.
@@ -117,6 +117,7 @@ def bode_plot(syslist, omega=None,
117117
config.defaults['freqplot.number_of_samples'].
118118
margins : bool
119119
If True, plot gain and phase margin.
120+
method : method to use in computing margins (see :func:`stability_margins`)
120121
*args : :func:`matplotlib.pyplot.plot` positional properties, optional
121122
Additional arguments for `matplotlib` plots (color, linestyle, etc)
122123
**kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
@@ -373,7 +374,7 @@ def bode_plot(syslist, omega=None,
373374
# Show the phase and gain margins in the plot
374375
if margins:
375376
# Compute stability margins for the system
376-
margin = stability_margins(sys)
377+
margin = stability_margins(sys, method=method)
377378
gm, pm, Wcg, Wcp = (margin[i] for i in (0, 1, 3, 4))
378379

379380
# Figure out sign of the phase at the first gain crossing

‎control/margins.py

Lines changed: 48 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -51,11 +51,13 @@
5151
"""
5252

5353
import math
54+
from warnings import warn
5455
import numpy as np
5556
import scipy as sp
5657
from . import xferfcn
5758
from .lti import issiso, evalfr
5859
from . import frdata
60+
from . import freqplot
5961
from .exception import ControlMIMONotImplemented
6062

6163
__all__ = ['stability_margins', 'phase_crossover_frequencies', 'margin']
@@ -206,6 +208,17 @@ def fun(wdt):
206208

207209
return z, w
208210

211+
def _likely_numerical_inaccuracy(sys):
212+
# crude, conservative check for if
213+
# num(z)*num(1/z) << den(z)*den(1/z) for DT systems
214+
num, den, num_inv_zp, den_inv_zq, p_q, dt = _poly_z_invz(sys)
215+
p1 = np.polymul(num, num_inv_zp)
216+
p2 = np.polymul(den, den_inv_zq)
217+
if p_q < 0:
218+
# * z**(-p_q)
219+
x = [1] + [0] * (-p_q)
220+
p1 = np.polymul(p1, x)
221+
return np.linalg.norm(p1) < 1e-3 * np.linalg.norm(p2)
209222

210223
# Took the framework for the old function by
211224
# Sawyer B. Fuller <minster@uw.edu>, removed a lot of the innards
@@ -237,25 +250,34 @@ def fun(wdt):
237250
# systems
238251

239252

240-
def stability_margins(sysdata, returnall=False, epsw=0.0):
253+
def stability_margins(sysdata, returnall=False, epsw=0.0, method='best'):
241254
"""Calculate stability margins and associated crossover frequencies.
242255
243256
Parameters
244257
----------
245-
sysdata: LTI system or (mag, phase, omega) sequence
258+
sysdata : LTI system or (mag, phase, omega) sequence
246259
sys : LTI system
247260
Linear SISO system representing the loop transfer function
248261
mag, phase, omega : sequence of array_like
249262
Arrays of magnitudes (absolute values, not dB), phases (degrees),
250263
and corresponding frequencies. Crossover frequencies returned are
251264
in the same units as those in `omega` (e.g., rad/sec or Hz).
252-
returnall: bool, optional
265+
returnall : bool, optional
253266
If true, return all margins found. If False (default), return only the
254267
minimum stability margins. For frequency data or FRD systems, only
255268
margins in the given frequency region can be found and returned.
256-
epsw: float, optional
269+
epsw : float, optional
257270
Frequencies below this value (default 0.0) are considered static gain,
258271
and not returned as margin.
272+
method : string, optional
273+
Method to use (default is 'best'):
274+
'poly': use polynomial method if passed a :class:`LTI` system.
275+
'frd': calculate crossover frequencies using numerical interpolation
276+
of a :class:`FrequencyResponseData` representation of the system if
277+
passed a :class:`LTI` system.
278+
'best': use the 'poly' method if possible, reverting to 'frd' if it is
279+
detected that numerical inaccuracy is likey to arise in the 'poly'
280+
method for for discrete-time systems.
259281
260282
Returns
261283
-------
@@ -279,6 +301,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
279301
determined by the frequency of maximum sensitivity (given by the
280302
magnitude of 1/(1+L)).
281303
"""
304+
# TODO: FRD method for cont-time systems doesn't work
282305
try:
283306
if isinstance(sysdata, frdata.FRD):
284307
sys = frdata.FRD(sysdata, smooth=True)
@@ -300,6 +323,27 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
300323
raise ControlMIMONotImplemented(
301324
"Can only do margins for SISO system")
302325

326+
if method == 'frd':
327+
# convert to FRD if we got a transfer function
328+
if isinstance(sys, xferfcn.TransferFunction):
329+
omega_sys = freqplot._default_frequency_range(sys)
330+
if sys.isctime():
331+
sys = frdata.FRD(sys, omega_sys)
332+
else:
333+
omega_sys = omega_sys[omega_sys < np.pi / sys.dt]
334+
sys = frdata.FRD(sys, omega_sys, smooth=True)
335+
elif method == 'best':
336+
# convert to FRD if anticipated numerical issues
337+
if isinstance(sys, xferfcn.TransferFunction) and not sys.isctime():
338+
if _likely_numerical_inaccuracy(sys):
339+
warn("stability_margins: Falling back to 'frd' method "
340+
"because of chance of numerical inaccuracy in 'poly' method.")
341+
omega_sys = freqplot._default_frequency_range(sys)
342+
omega_sys = omega_sys[omega_sys < np.pi / sys.dt]
343+
sys = frdata.FRD(sys, omega_sys, smooth=True)
344+
elif method != 'poly':
345+
raise ValueError("method " + method + " unknown")
346+
303347
if isinstance(sys, xferfcn.TransferFunction):
304348
if sys.isctime():
305349
num_iw, den_iw = _poly_iw(sys)
@@ -366,7 +410,6 @@ def _dstab(w):
366410
w_180 = np.array(
367411
[sp.optimize.brentq(_arg, sys.omega[i], sys.omega[i+1])
368412
for i in widx])
369-
# TODO: replace by evalfr(sys, 1J*w) or sys(1J*w), (needs gh-449)
370413
w180_resp = sys(1j * w_180)
371414

372415
# Find all crossings, note that this depends on omega having

‎control/tests/margin_test.py

Lines changed: 22 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -104,7 +104,6 @@ def test_margin_sys(tsys):
104104
out = margin(sys)
105105
assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
106106

107-
108107
def test_margin_3input(tsys):
109108
sys, refout, refoutall = tsys
110109
"""Test margin() function with mag, phase, omega input"""
@@ -339,11 +338,11 @@ def test_zmore_stability_margins(tsys_zmore):
339338
'ref,'
340339
'rtol',
341340
[( # gh-465
342-
[2], [1, 3, 2, 0], 1e-2,
341+
[2], [1, 3, 2, 0], 1e-2,
343342
[2.9558, 32.390, 0.43584, 1.4037, 0.74951, 0.97079],
344343
2e-3), # the gradient of the function reduces numerical precision
345344
( # 2/(s+1)**3
346-
[2], [1, 3, 3, 1], .1,
345+
[2], [1, 3, 3, 1], .1,
347346
[3.4927, 65.4212, 0.5763, 1.6283, 0.76625, 1.2019],
348347
1e-4),
349348
( # gh-523
@@ -354,5 +353,24 @@ def test_zmore_stability_margins(tsys_zmore):
354353
def test_stability_margins_discrete(cnum, cden, dt, ref, rtol):
355354
"""Test stability_margins with discrete TF input"""
356355
tf = TransferFunction(cnum, cden).sample(dt)
357-
out = stability_margins(tf)
356+
out = stability_margins(tf, method='poly')
358357
assert_allclose(out, ref, rtol=rtol)
358+
359+
360+
def test_stability_margins_methods():
361+
# the following system gives slightly inaccurate result for DT systems
362+
# because of numerical issues
363+
omegan = 1
364+
zeta = 0.5
365+
resonance = TransferFunction(omegan**2, [1, 2*zeta*omegan, omegan**2])
366+
omegan2 = 100
367+
resonance2 = TransferFunction(omegan2**2, [1, 2*zeta*omegan2, omegan2**2])
368+
sys = 5 * resonance * resonance2
369+
sysd = sys.sample(0.001, 'zoh')
370+
"""Test stability_margins() function with different methods"""
371+
out = stability_margins(sysd, method='best')
372+
# confirm getting reasonable results using FRD method
373+
assert_allclose(
374+
(18.876634740386308, 26.356358386241055, 0.40684127995261044,
375+
9.763585494645046, 2.3293357226374805, 2.55985695034263),
376+
stability_margins(sysd, method='frd'), rtol=1e-5)

0 commit comments

Comments
 (0)
0