8000 update nyquist_plot to accomodate discrete-time transfer functions wi… · python-control/python-control@d88c67c · GitHub
[go: up one dir, main page]
Skip to content

Commit d88c67c

Browse files
sawyerbfullersawyerbfuller
committed
update nyquist_plot to accomodate discrete-time transfer functions with poles at the origin and unity
1 parent 999189c commit d88c67c

File tree

2 files changed

+56
-39
lines changed

2 files changed

+56
-39
lines changed

‎control/freqplot.py

Lines changed: 38 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -819,29 +819,29 @@ def _parse_linestyle(style_name, allow_false=False):
819819
splane_contour = 1j * omega_sys
820820

821821
# Bend the contour around any poles on/near the imaginary axis
822-
# TODO: smarter indent radius that depends on dcgain of system
823-
# and timebase of discrete system.
824822
if isinstance(sys, (StateSpace, TransferFunction)) \
825823
and indent_direction != 'none':
826824
if sys.isctime():
827825
splane_poles = sys.poles()
828826
splane_cl_poles = sys.feedback().poles()
829827
else:
830-
# map z-plane poles to s-plane, ignoring any at the origin
831-
# because we don't need to indent for them
828+
# map z-plane poles to s-plane. We ignore any at the origin
829+
# to avoid numerical warnings because we know we
830+
# don't need to indent for them
832831
zplane_poles = sys.poles()
833832
zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)]
834833
splane_poles = np.log(zplane_poles) / sys.dt
835834

836835
zplane_cl_poles = sys.feedback().poles()
836+
# eliminate z-plane poles at the origin to avoid warnings
837837
zplane_cl_poles = zplane_cl_poles[
838-
~np.isclose(abs(zplane_poles), 0.)]
838+
~np.isclose(abs(zplane_cl_poles), 0.)]
839839
splane_cl_poles = np.log(zplane_cl_poles) / sys.dt
840840

841841
#
842842
# Check to make sure indent radius is small enough
843843
#
844-
# If there is a closed loop pole that is near the imaginary access
844+
# If there is a closed loop pole that is near the imaginary axis
845845
# at a point that is near an open loop pole, it is possible that
846846
# indentation might skip or create an extraneous encirclement.
847847
# We check for that situation here and generate a warning if that
@@ -851,15 +851,16 @@ def _parse_linestyle(style_name, allow_false=False):
851851
# See if any closed loop poles are near the imaginary axis
852852
if abs(p_cl.real) <= indent_radius:
853853
# See if any open loop poles are close to closed loop poles
854-
p_ol = splane_poles[
855-
(np.abs(splane_poles - p_cl)).argmin()]
854+
if len(splane_poles) > 0:
855+
p_ol = splane_poles[
856+
(np.abs(splane_poles - p_cl)).argmin()]
856857

857-
if abs(p_ol - p_cl) <= indent_radius and \
858-
warn_encirclements:
859-
warnings.warn(
860-
"indented contour may miss closed loop pole; "
861-
"consider reducing indent_radius to be less than "
862-
f"{abs(p_ol - p_cl):5.2g}", stacklevel=2)
858+
if abs(p_ol - p_cl) <= indent_radius and \
859+
warn_encirclements:
860+
warnings.warn(
861+
"indented contour may miss closed loop pole; "
862+
"consider reducing indent_radius to below "
863+
f"{abs(p_ol - p_cl):5.2g}", stacklevel=2)
863864

864865
#
865866
# See if we should add some frequency points near imaginary poles
@@ -897,29 +898,30 @@ def _parse_linestyle(style_name, allow_false=False):
897898
splane_contour[last_point:]))
898899

899900
# Indent points that are too close to a pole
900-
for i, s in enumerate(splane_contour):
901-
# Find the nearest pole
902-
p = splane_poles[(np.abs(splane_poles - s)).argmin()]
903-
904-
# See if we need to indent around it
905-
if abs(s - p) < indent_radius:
906-
# Figure out how much to offset (simple trigonometry)
907-
offset = np.sqrt(indent_radius ** 2 - (s - p).imag ** 2) \
908-
- (s - p).real
909-
910-
# Figure out which way to offset the contour point
911-
if p.real < 0 or (p.real == 0 and
912-
indent_direction == 'right'):
913-
# Indent to the right
914-
splane_contour[i] += offset
915-
916-
elif p.real > 0 or (p.real == 0 and
917-
indent_direction == 'left'):
918-
# Indent to the left
919-
splane_contour[i] -= offset
901+
if len(splane_poles) > 0: # accomodate no splane poles if dtime sys
902+
for i, s in enumerate(splane_contour):
903+
# Find the nearest pole
904+
p = splane_poles[(np.abs(splane_poles - s)).argmin()]
905+
906+
# See if we need to indent around it
907+
if abs(s - p) < indent_radius:
908+
# Figure out how much to offset (simple trigonometry)
909+
offset = np.sqrt(indent_radius ** 2 - (s - p).imag ** 2) \
910+
- (s - p).real
911+
912+
# Figure out which way to offset the contour point
913+
if p.real < 0 or (p.real == 0 and
914+
indent_direction == 'right'):
915+
# Indent to the right
916+
splane_contour[i] += offset
917+
918+
elif p.real > 0 or (p.real == 0 and
919+
indent_direction == 'left'):
920+
# Indent to the left
921+
splane_contour[i] -= offset
920922

921-
else:
922-
raise ValueError("unknown value for indent_direction")
923+
else:
924+
raise ValueError("unknown value for indent_direction")
923925

924926
# change contour to z-plane if necessary
925927
if sys.isctime():

‎control/tests/nyquist_test.py

Lines changed: 18 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -370,8 +370,23 @@ def test_nyquist_legacy():
370370
def test_discrete_nyquist():
371371
# Make sure we can handle discrete time systems with negative poles
372372
sys = ct.tf(1, [1, -0.1], dt=1) * ct.tf(1, [1, 0.1], dt=1)
373-
ct.nyquist_plot(sys)
374-
373+
ct.nyquist_plot(sys, plot=False)
374+
375+
# system with a pole at the origin
376+
sys = ct.zpk([1,], [.3, 0], 1, dt=True)
377+
ct.nyquist_plot(sys, plot=False)
378+
sys = ct.zpk([1,], [0], 1, dt=True)
379+
ct.nyquist_plot(sys, plot=False)
380+
381+
# only a pole at the origin
382+
sys = ct.zpk([], [0], 2, dt=True)
383+
ct.nyquist_plot(sys, plot=False)
384+
385+
# pole at zero (pure delay)
386+
sys = ct.zpk([], [1], 1, dt=True)
387+
ct.nyquist_plot(sys, plot=False)
388+
389+
375390
if __name__ == "__main__":
376391
#
377392
# Interactive mode: generate plots for manual viewing
@@ -427,5 +442,5 @@ def test_discrete_nyquist():
427442
np.array2string(sys.poles(), precision=2, separator=','))
428443
count = ct.nyquist_plot(sys)
429444

430-
445+
431446

0 commit comments

Comments
 (0)
0