8000 add warning for nyquist_plot() when indent_radius is too large · python-control/python-control@18e997c · GitHub
[go: up one dir, main page]
Skip to content

Commit 18e997c

Browse files
murrayrmmurrayrm
committed
add warning for nyquist_plot() when indent_radius is too large
1 parent 2102181 commit 18e997c

File tree

2 files changed

+68
-4
lines changed

2 files changed

+68
-4
lines changed

‎control/freqplot.py

Lines changed: 41 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -723,35 +723,72 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
723723
and indent_direction != 'none':
724724
if sys.isctime():
725725
splane_poles = sys.poles()
726+
splane_cl_poles = sys.feedback().poles()
726727
else:
727728
# map z-plane poles to s-plane, ignoring any at the origin
728729
# because we don't need to indent for them
729730
zplane_poles = sys.poles()
730731
zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)]
731732
splane_poles = np.log(zplane_poles)/sys.dt
732733

734+
zplane_cl_poles = sys.feedback().poles()
735+
zplane_cl_poles = zplane_cl_poles[
736+
~np.isclose(abs(zplane_poles), 0.)]
737+
splane_cl_poles = np.log(zplane_cl_poles)/sys.dt
738+
739+
#
740+
# Check to make sure indent radius is small enough
741+
#
742+
# If there is a closed loop pole that is near the imaginary access
743+
# at a point that is near an open loop pole, it is possible that
744+
# indentation might skip or create an extraneous encirclement.
745+
# We check for that situation here and generate a warning if that
746+
# could happen.
747+
#
748+
for p_cl in splane_cl_poles:
749+
# See if any closed loop poles are near the imaginary axis
750+
if abs(p_cl.real) <= indent_radius:
751+
# See if any open loop poles are close to closed loop poles
752+
p_ol = splane_poles[
753+
(np.abs(splane_poles - p_cl)).argmin()]
754+
755+
if abs(p_ol - p_cl) <= indent_radius:
756+
warnings.warn(
757+
"indented contour may miss closed loop pole; "
758+
8000 "consider reducing indent_radius to be less than "
759+
f"{abs(p_ol - p_cl):5.2g}", stacklevel=2)
760+
761+
# See if we should add some frequency points near the origin
733762
if splane_contour[1].imag > indent_radius \
734763
and np.any(np.isclose(abs(splane_poles), 0)) \
735764
and not omega_range_given:
736765
# add some points for quarter circle around poles at origin
766+
# (these will get indented left or right below)
737767
splane_contour = np.concatenate(
738768
(1j * np.linspace(0., indent_radius, 50),
739769
splane_contour[1:]))
770+
740771
for i, s in enumerate(splane_contour):
741772
# Find the nearest pole
742773
p = splane_poles[(np.abs(splane_poles - s)).argmin()]
774+
743775
# See if we need to indent around it
744776
if abs(s - p) < indent_radius:
777+
# Figure out how much to offset (simple trigonometry)
778+
offset = np.sqrt(indent_radius ** 2 - (s-p).imag ** 2) \
779+
-(s-p).real
780+
781+
# Figure out which way to offset the contour point
745782
if p.real < 0 or (np.isclose(p.real, 0) \
746783
and indent_direction == 'right'):
747784
# Indent to the right
748-
splane_contour[i] += \
749-
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
785+
splane_contour[i] += offset
786+
750787
elif p.real > 0 or (np.isclose(p.real, 0) \
751788
and indent_direction == 'left'):
752789
# Indent to the left
753-
splane_contour[i] -= \
754-
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
790+
splane_contour[i] -= offset
791+
755792
else:
756793
ValueError("unknown value for indent_direction")
757794

‎control/tests/nyquist_test.py

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,33 @@ def test_nyquist_basic():
4141
N_sys = ct.nyquist_plot(sys)
4242
assert _Z(sys) == N_sys + _P(sys)
4343

44+
# Previously identified bug
45+
#
46+
# This example has an open loop pole at -0.06 and a closed loop pole at
47+
# 0.06, so if you use an indent_radius of larger than 0.12, then the
48+
# encirclements computed by nyquist_plot() will not properly predict
49+
# stability. A new warning messages was added to catch this case.
50+
#
51+
A = np.array([
52+
[-3.56355873, -1.22980795, -1.5626527 , -0.4626829 , -0.16741484],
53+
[-8.52361371, -3.60331459, -3.71574266, -0.43839201, 0.41893656],
54+
[-2.50458726, -0.72361335, -1.77795489, -0.4038419 , 0.52451147],
55+
[-0.281183 , 0.23391825, 0.19096003, -0.9771515 , 0.66975606],
56+
[-3.04982852, -1.1091943 , -1.40027242, -0.1974623 , -0.78930791]])
57+
B = np.array([[-0.], [-1.42827213], [ 0.76806551], [-1.07987454], [0.]])
58+
C = np.array([[-0., 0.35557249, 0.35941791, -0., -1.42320969]])
59+
D = np.array([[0]])
60+
sys = ct.ss(A, B, C, D)
61+
62+
# With a small indent_radius, all should be fine
63+
N_sys = ct.nyquist_plot(sys, indent_radius=0.001)
64+
assert _Z(sys) == N_sys + _P(sys)
65+
66+
# With a larger indent_radius, we get a warning message + wrong answer
67+
with pytest.warns(UserWarning, match="contour may miss closed loop pole"):
68+
N_sys = ct.nyquist_plot(sys, indent_radius=0.2)
69+
assert _Z(sys) != N_sys + _P(sys)
70+
4471
# Unstable system
4572
sys = ct.tf([10], [1, 2, 2, 1])
4673
N_sys = ct.nyquist_plot(sys)

0 commit comments

Comments
 (0)
0