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

Commit618ea45

Browse files
authored
Merge pull request#620 from bnavigator/freqplot-omega
Refine automatic contour determination in Nyquist plot
2 parents608093e +5cf0358 commit618ea45

File tree

2 files changed

+78
-82
lines changed

2 files changed

+78
-82
lines changed

‎control/freqplot.py

Lines changed: 33 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -628,7 +628,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
628628
629629
2. If a continuous-time system contains poles on or near the imaginary
630630
axis, a small indentation will be used to avoid the pole. The radius
631-
of the indentation is given by `indent_radius` and it is takenthe the
631+
of the indentation is given by `indent_radius` and it is takento the
632632
right of stable poles and the left of unstable poles. If a pole is
633633
exactly on the imaginary axis, the `indent_direction` parameter can be
634634
used to set the direction of indentation. Setting `indent_direction`
@@ -683,26 +683,11 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
683683
ifnothasattr(syslist,'__iter__'):
684684
syslist= (syslist,)
685685

686-
# Decide whether to go above Nyquist frequency
687-
omega_range_given=TrueifomegaisnotNoneelseFalse
688-
689-
# Figure out the frequency limits
690-
ifomegaisNone:
691-
ifomega_limitsisNone:
692-
# Select a default range if none is provided
693-
omega=_default_frequency_range(
694-
syslist,number_of_samples=omega_num)
695-
696-
# Replace first point with the origin
697-
omega[0]=0
698-
else:
699-
omega_range_given=True
700-
omega_limits=np.asarray(omega_limits)
701-
iflen(omega_limits)!=2:
702-
raiseValueError("len(omega_limits) must be 2")
703-
omega=np.logspace(np.log10(omega_limits[0]),
704-
np.log10(omega_limits[1]),num=omega_num,
705-
endpoint=True)
686+
omega,omega_range_given=_determine_omega_vector(
687+
syslist,omega,omega_limits,omega_num)
688+
ifnotomega_range_given:
689+
# Start contour at zero frequency
690+
omega[0]=0.
706691

707692
# Go through each system and keep track of the results
708693
counts,contours= [], []
@@ -734,9 +719,15 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
734719
contour=1j*omega_sys
735720

736721
# Bend the contour around any poles on/near the imaginary axis
737-
ifisinstance(sys, (StateSpace,TransferFunction))and\
738-
sys.isctime()andindent_direction!='none':
722+
ifisinstance(sys, (StateSpace,TransferFunction)) \
723+
andsys.isctime()andindent_direction!='none':
739724
poles=sys.pole()
725+
ifcontour[1].imag>indent_radius \
726+
and0.inpolesandnotomega_range_given:
727+
# add some points for quarter circle around poles at origin
728+
contour=np.concatenate(
729+
(1j*np.linspace(0.,indent_radius,50),
730+
contour[1:]))
740731
fori,sinenumerate(contour):
741732
# Find the nearest pole
742733
p=poles[(np.abs(poles-s)).argmin()]
@@ -1242,17 +1233,15 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num):
12421233
omega_in or through omega_limits. False if both omega_in
12431234
and omega_limits are None.
12441235
"""
1245-
1246-
# Decide whether to go above Nyquist frequency
1247-
omega_range_given=Trueifomega_inisnotNoneelseFalse
1236+
omega_range_given=True
12481237

12491238
ifomega_inisNone:
12501239
ifomega_limitsisNone:
1240+
omega_range_given=False
12511241
# Select a default range if none is provided
12521242
omega_out=_default_frequency_range(syslist,
12531243
number_of_samples=omega_num)
12541244
else:
1255-
omega_range_given=True
12561245
omega_limits=np.asarray(omega_limits)
12571246
iflen(omega_limits)!=2:
12581247
raiseValueError("len(omega_limits) must be 2")
@@ -1267,11 +1256,12 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num):
12671256
# Compute reasonable defaults for axes
12681257
def_default_frequency_range(syslist,Hz=None,number_of_samples=None,
12691258
feature_periphery_decades=None):
1270-
"""Compute a reasonable default frequency range for frequency
1271-
domain plots.
1259+
"""Compute a default frequency range for frequency domain plots.
12721260
1273-
Finds a reasonable default frequency range by examining the features
1274-
(poles and zeros) of the systems in syslist.
1261+
This code looks at the poles and zeros of all of the systems that
1262+
we are plotting and sets the frequency range to be one decade above
1263+
and below the min and max feature frequencies, rounded to the nearest
1264+
integer. If no features are found, it returns logspace(-1, 1)
12751265
12761266
Parameters
12771267
----------
@@ -1302,12 +1292,6 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
13021292
>>> omega = _default_frequency_range(sys)
13031293
13041294
"""
1305-
# This code looks at the poles and zeros of all of the systems that
1306-
# we are plotting and sets the frequency range to be one decade above
1307-
# and below the min and max feature frequencies, rounded to the nearest
1308-
# integer. It excludes poles and zeros at the origin. If no features
1309-
# are found, it turns logspace(-1, 1)
1310-
13111295
# Set default values for options
13121296
number_of_samples=config._get_param(
13131297
'freqplot','number_of_samples',number_of_samples)
@@ -1329,30 +1313,31 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
13291313
features_=np.concatenate((np.abs(sys.pole()),
13301314
np.abs(sys.zero())))
13311315
# Get rid of poles and zeros at the origin
1332-
features_=features_[features_!=0.0]
1333-
features=np.concatenate((features,features_))
1316+
toreplace=features_==0.0
1317+
ifnp.any(toreplace):
1318+
features_=features_[~toreplace]
13341319
elifsys.isdtime(strict=True):
13351320
fn=math.pi*1./sys.dt
13361321
# TODO: What distance to the Nyquist frequency is appropriate?
13371322
freq_interesting.append(fn*0.9)
13381323

13391324
features_=np.concatenate((sys.pole(),
13401325
sys.zero()))
1341-
# Get rid of poles and zeros
1342-
# *at theorigin and real <= 0 & imag==0: log!
1326+
# Get rid of poles and zeros on the real axis (imag==0)
1327+
# * origin and real < 0
13431328
# * at 1.: would result in omega=0. (logaritmic plot!)
1344-
features_=features_[
1345-
(features_.imag!=0.0)|(features_.real>0.)]
1346-
features_=features_[
1347-
np.bitwise_not((features_.imag==0.0)&
1348-
(np.abs(features_.real-1.0)<1.e-10))]
1329+
toreplace=(features_.imag==0.0)& (
1330+
(features_.real<=0.)|
1331+
(np.abs(features_.real-1.0)<1.e-10))
1332+
ifnp.any(toreplace):
1333+
features_=features_[~toreplace]
13491334
# TODO: improve
1350-
features__=np.abs(np.log(features_)/ (1.j*sys.dt))
1351-
features=np.concatenate((features,features__))
1335+
features_=np.abs(np.log(features_)/ (1.j*sys.dt))
13521336
else:
13531337
# TODO
13541338
raiseNotImplementedError(
13551339
"type of system in not implemented now")
1340+
features=np.concatenate((features,features_))
13561341
exceptNotImplementedError:
13571342
pass
13581343

‎control/tests/nyquist_test.py

Lines changed: 45 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -10,12 +10,10 @@
1010

1111
importpytest
1212
importnumpyasnp
13-
importscipyassp
1413
importmatplotlib.pyplotasplt
1514
importcontrolasct
1615

17-
# In interactive mode, turn on ipython interactive graphics
18-
plt.ion()
16+
pytestmark=pytest.mark.usefixtures("mplcleanup")
1917

2018

2119
# Utility function for counting unstable poles of open loop (P in FBS)
@@ -37,7 +35,6 @@ def _Z(sys):
3735

3836

3937
# Basic tests
40-
@pytest.mark.usefixtures("mplcleanup")
4138
deftest_nyquist_basic():
4239
# Simple Nyquist plot
4340
sys=ct.rss(5,1,1)
@@ -112,7 +109,6 @@ def test_nyquist_basic():
112109

113110

114111
# Some FBS examples, for comparison
115-
@pytest.mark.usefixtures("mplcleanup")
116112
deftest_nyquist_fbs_examples():
117113
s=ct.tf('s')
118114

@@ -154,7 +150,6 @@ def test_nyquist_fbs_examples():
154150
1,2,3,4,# specified number of arrows
155151
[0.1,0.5,0.9],# specify arc lengths
156152
])
157-
@pytest.mark.usefixtures("mplcleanup")
158153
deftest_nyquist_arrows(arrows):
159154
sys=ct.tf([1.4], [1,2,1])*ct.tf(*ct.pade(1,4))
160155
plt.figure();
@@ -163,7 +158,6 @@ def test_nyquist_arrows(arrows):
163158
assert_Z(sys)==count+_P(sys)
164159

165160

166-
@pytest.mark.usefixtures("mplcleanup")
167161
deftest_nyquist_encirclements():
168162
# Example 14.14: effect of friction in a cart-pendulum system
169163
s=ct.tf('s')
@@ -188,21 +182,34 @@ def test_nyquist_encirclements():
188182
assert_Z(sys)==count+_P(sys)
189183

190184

191-
@pytest.mark.usefixtures("mplcleanup")
192185
deftest_nyquist_indent():
193186
# FBS Figure 10.10
194187
s=ct.tf('s')
195188
sys=3* (s+6)**2/ (s* (s+1)**2)
189+
# poles: [-1, -1, 0]
196190

197191
plt.figure();
198192
count=ct.nyquist_plot(sys)
199193
plt.title("Pole at origin; indent_radius=default")
200194
assert_Z(sys)==count+_P(sys)
201195

196+
# first value of default omega vector was 0.1, replaced by 0. for contour
197+
# indent_radius is larger than 0.1 -> no extra quater circle around origin
198+
count,contour=ct.nyquist_plot(sys,plot=False,indent_radius=.1007,
199+
return_contour=True)
200+
np.testing.assert_allclose(contour[0],.1007+0.j)
201+
# second value of omega_vector is larger than indent_radius: not indented
202+
assertnp.all(contour.real[2:]==0.)
203+
202204
plt.figure();
203-
count=ct.nyquist_plot(sys,indent_radius=0.01)
205+
count,contour=ct.nyquist_plot(sys,indent_radius=0.01,
206+
return_contour=True)
204207
plt.title("Pole at origin; indent_radius=0.01; encirclements = %d"%count)
205208
assert_Z(sys)==count+_P(sys)
209+
# indent radius is smaller than the start of the default omega vector
210+
# check that a quarter circle around the pole at origin has been added.
211+
np.testing.assert_allclose(contour[:50].real**2+contour[:50].imag**2,
212+
0.01**2)
206213

207214
plt.figure();
208215
count=ct.nyquist_plot(sys,indent_direction='left')
@@ -255,34 +262,38 @@ def test_nyquist_exceptions():
255262
ct.nyquist_plot(sys,np.logspace(-2,3))
256263

257264

258-
#
259-
# Interactive mode: generate plots for manual viewing
260-
#
261-
# Running this script in python (or better ipython) will show a collection of
262-
# figures that should all look OK on the screeen.
263-
#
265+
if__name__=="__main__":
266+
#
267+
# Interactive mode: generate plots for manual viewing
268+
#
269+
# Running this script in python (or better ipython) will show a collection of
270+
# figures that should all look OK on the screeen.
271+
#
272+
273+
# In interactive mode, turn on ipython interactive graphics
274+
plt.ion()
264275

265-
# Start by clearing existing figures
266-
plt.close('all')
276+
# Start by clearing existing figures
277+
plt.close('all')
267278

268-
print("Nyquist examples from FBS")
269-
test_nyquist_fbs_examples()
279+
print("Nyquist examples from FBS")
280+
test_nyquist_fbs_examples()
270281

271-
print("Arrow test")
272-
test_nyquist_arrows(None)
273-
test_nyquist_arrows(1)
274-
test_nyquist_arrows(3)
275-
test_nyquist_arrows([0.1,0.5,0.9])
282+
print("Arrow test")
283+
test_nyquist_arrows(None)
284+
test_nyquist_arrows(1)
285+
test_nyquist_arrows(3)
286+
test_nyquist_arrows([0.1,0.5,0.9])
276287

277-
print("Stability checks")
278-
test_nyquist_encirclements()
288+
print("Stability checks")
289+
test_nyquist_encirclements()
279290

280-
print("Indentation checks")
281-
test_nyquist_indent()
291+
print("Indentation checks")
292+
test_nyquist_indent()
282293

283-
print("Unusual Nyquist plot")
284-
sys=ct.tf([1], [1,3,2])*ct.tf([1], [1,0,1])
285-
plt.figure()
286-
plt.title("Poles: %s"%np.array2string(sys.pole(),precision=2,separator=','))
287-
count=ct.nyquist_plot(sys)
288-
assert_Z(sys)==count+_P(sys)
294+
print("Unusual Nyquist plot")
295+
sys=ct.tf([1], [1,3,2])*ct.tf([1], [1,0,1])
296+
plt.figure()
297+
plt.title("Poles: %s"%np.array2string(sys.pole(),precision=2,separator=','))
298+
count=ct.nyquist_plot(sys)
299+
assert_Z(sys)==count+_P(sys)

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp