Skip to content

Ease cutoff between poly and frd method for stability margins #575

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.

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Mar 19, 2021
Merged

Ease cutoff between poly and frd method for stability margins #575

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4381a7e
merge the second example from #523 into test_stability_margins_discrete
bnavigator Mar 12, 2021
52d8fc3
improve error message display, ease poly to frd check
bnavigator Mar 17, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions control/margins.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,7 @@ def _likely_numerical_inaccuracy(sys):
# * z**(-p_q)
x = [1] + [0] * (-p_q)
p1 = np.polymul(p1, x)
return np.linalg.norm(p1) < 1e-3 * np.linalg.norm(p2)
return np.linalg.norm(p1) < 1e-4 * np.linalg.norm(p2)

# Took the framework for the old function by
# Sawyer B. Fuller <minster@uw.edu>, removed a lot of the innards
Expand Down Expand Up @@ -337,7 +337,8 @@ def stability_margins(sysdata, returnall=False, epsw=0.0, method='best'):
if isinstance(sys, xferfcn.TransferFunction) and not sys.isctime():
if _likely_numerical_inaccuracy(sys):
warn("stability_margins: Falling back to 'frd' method "
"because of chance of numerical inaccuracy in 'poly' method.")
"because of chance of numerical inaccuracy in 'poly' method.",
stacklevel=2)
omega_sys = freqplot._default_frequency_range(sys)
omega_sys = omega_sys[omega_sys < np.pi / sys.dt]
sys = frdata.FRD(sys, omega_sys, smooth=True)
Expand Down
54 changes: 27 additions & 27 deletions control/tests/margin_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -336,41 +336,41 @@ def test_zmore_stability_margins(tsys_zmore):
@pytest.mark.parametrize(
'cnum, cden, dt,'
'ref,'
'rtol',
'rtol, poly_is_inaccurate',
[( # gh-465
[2], [1, 3, 2, 0], 1e-2,
[2.9558, 32.390, 0.43584, 1.4037, 0.74951, 0.97079],
2e-3), # the gradient of the function reduces numerical precision
[ 2.955761, 32.398492, 0.429535, 1.403725, 0.749367, 0.923898],
1e-5, True),
( # 2/(s+1)**3
[2], [1, 3, 3, 1], .1,
[3.4927, 65.4212, 0.5763, 1.6283, 0.76625, 1.2019],
1e-4),
( # gh-523
1e-4, True),
( # gh-523 a
[1.1 * 4 * np.pi**2], [1, 2 * 0.2 * 2 * np.pi, 4 * np.pi**2], .05,
[2.3842, 18.161, 0.26953, 11.712, 8.7478, 9.1504],
1e-4),
1e-4, False),
( # gh-523 b
# H1 = w1**2 / (z**2 + 2*zt*w1 * z + w1**2)
# H2 = w2**2 / (z**2 + 2*zt*w2 * z + w2**2)
# H = H1 * H2
# w1 = 1, w2 = 100, zt = 0.5
[5e4], [1., 101., 10101., 10100., 10000.], 1e-3,
[18.8766, 26.3564, 0.406841, 9.76358, 2.32933, 2.55986],
1e-5, True),
])
def test_stability_margins_discrete(cnum, cden, dt, ref, rtol):
@pytest.mark.filterwarnings("error")
def test_stability_margins_discrete(cnum, cden, dt,
ref,
rtol, poly_is_inaccurate):
"""Test stability_margins with discrete TF input"""
tf = TransferFunction(cnum, cden).sample(dt)
out = stability_margins(tf, method='poly')
if poly_is_inaccurate:
with pytest.warns(UserWarning, match="numerical inaccuracy in 'poly'"):
out = stability_margins(tf)
# cover the explicit frd branch and make sure it yields the same
# results as the fallback mechanism
out_frd = stability_margins(tf, method='frd')
assert_allclose(out, out_frd)
else:
out = stability_margins(tf)
assert_allclose(out, ref, rtol=rtol)


def test_stability_margins_methods():
# the following system gives slightly inaccurate result for DT systems
# because of numerical issues
omegan = 1
zeta = 0.5
resonance = TransferFunction(omegan**2, [1, 2*zeta*omegan, omegan**2])
omegan2 = 100
resonance2 = TransferFunction(omegan2**2, [1, 2*zeta*omegan2, omegan2**2])
sys = 5 * resonance * resonance2
sysd = sys.sample(0.001, 'zoh')
"""Test stability_margins() function with different methods"""
out = stability_margins(sysd, method='best')
# confirm getting reasonable results using FRD method
assert_allclose(
(18.876634740386308, 26.356358386241055, 0.40684127995261044,
9.763585494645046, 2.3293357226374805, 2.55985695034263),
stability_margins(sysd, method='frd'), rtol=1e-5)