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

Commitd88c67c

Browse files
committed
update nyquist_plot to accomodate discrete-time transfer functions with poles at the origin and unity
1 parent999189c commitd88c67c

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
ifisinstance(sys, (StateSpace,TransferFunction)) \
825823
andindent_direction!='none':
826824
ifsys.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 imaginaryaccess
844+
# If there is a closed loop pole that is near the imaginaryaxis
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
ifabs(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+
iflen(splane_poles)>0:
855+
p_ol=splane_poles[
856+
(np.abs(splane_poles-p_cl)).argmin()]
856857

857-
ifabs(p_ol-p_cl)<=indent_radiusand \
858-
warn_encirclements:
859-
warnings.warn(
860-
"indented contour may miss closed loop pole; "
861-
"consider reducing indent_radius tobe less than "
862-
f"{abs(p_ol-p_cl):5.2g}",stacklevel=2)
858+
ifabs(p_ol-p_cl)<=indent_radiusand \
859+
warn_encirclements:
860+
warnings.warn(
861+
"indented contour may miss closed loop pole; "
862+
"consider reducing indent_radius tobelow "
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-
fori,sinenumerate(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-
ifabs(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-
ifp.real<0or (p.real==0and
912-
indent_direction=='right'):
913-
# Indent to the right
914-
splane_contour[i]+=offset
915-
916-
elifp.real>0or (p.real==0and
917-
indent_direction=='left'):
918-
# Indent to the left
919-
splane_contour[i]-=offset
901+
iflen(splane_poles)>0:# accomodate no splane poles if dtime sys
902+
fori,sinenumerate(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+
ifabs(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+
ifp.real<0or (p.real==0and
914+
indent_direction=='right'):
915+
# Indent to the right
916+
splane_contour[i]+=offset
917+
918+
elifp.real>0or (p.real==0and
919+
indent_direction=='left'):
920+
# Indent to the left
921+
splane_contour[i]-=offset
920922

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

924926
# change contour to z-plane if necessary
925927
ifsys.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
deftest_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)
[8]ページ先頭
©2009-2025 Movatter.jp