10000 Merge pull request #688 from bnavigator/omega_default_round · python-control/python-control@df6b352 · GitHub
[go: up one dir, main page]
Skip to content

Commit df6b352

Browse files
bnavigatorbnavigator
authored
Merge pull request #688 from bnavigator/omega_default_round
Round to nearest integer decade for default omega vector
2 parents a111b03 + 021372f commit df6b352

File tree

3 files changed

+46
-28
lines changed

3 files changed

+46
-28
lines changed

‎control/freqplot.py

Lines changed: 19 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -209,7 +209,7 @@ def bode_plot(syslist, omega=None,
209209
syslist = (syslist,)
210210

211211
omega, omega_range_given = _determine_omega_vector(
212-
syslist, omega, omega_limits, omega_num)
212+
syslist, omega, omega_limits, omega_num, Hz=Hz)
213213

214214
if plot:
215215
# Set up the axes with labels so that multiple calls to
@@ -965,7 +965,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
965965
# Select a default range if none is provided
966966
# TODO: This needs to be made more intelligent
967967
if omega is None:
968-
omega = _default_frequency_range((P, C, S))
968+
omega = _default_frequency_range((P, C, S), Hz=Hz)
969969

970970
# Set up the axes with labels so that multiple calls to
971971
# gangof4_plot will superimpose the data. See details in bode_plot.
@@ -1115,7 +1115,7 @@ def singular_values_plot(syslist, omega=None,
11151115
syslist = (syslist,)
11161116

11171117
omega, omega_range_given = _determine_omega_vector(
1118-
syslist, omega, omega_limits, omega_num)
1118+
syslist, omega, omega_limits, omega_num, Hz=Hz)
11191119

11201120
omega = np.atleast_1d(omega)
11211121

@@ -1210,7 +1210,8 @@ def singular_values_plot(syslist, omega=None,
12101210

12111211

12121212
# Determine the frequency range to be used
1213-
def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num):
1213+
def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num,
1214+
Hz=None):
12141215
"""Determine the frequency range for a frequency-domain plot
12151216
according to a standard logic.
12161217
@@ -1236,6 +1237,10 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num):
12361237
omega_num : int
12371238
Number of points to be used for the frequency
12381239
range (if the frequency range is not user-specified)
1240+
Hz : bool, optional
1241+
If True, the limits (first and last value) of the frequencies
1242+
are set to full decades in Hz so it fits plotting with logarithmic
1243+
scale in Hz otherwise in rad/s. Omega is always returned in rad/sec.
12391244
12401245
Returns
12411246
-------
@@ -1253,7 +1258,8 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num):
12531258
omega_range_given = False
12541259
# Select a default range if none is provided
12551260
omega_out = _default_frequency_range(syslist,
1256-
number_of_samples=omega_num)
1261+
number_of_samples=omega_num,
1262+
Hz=Hz)
12571263
else:
12581264
omega_limits = np.asarray(omega_limits)
12591265
if len(omega_limits) != 2:
@@ -1280,7 +1286,7 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
12801286
----------
12811287
syslist : list of LTI
12821288
List of linear input/output systems (single system is OK)
1283-
Hz : bool
1289+
Hz : bool, optional
12841290
If True, the limits (first and last value) of the frequencies
12851291
are set to full decades in Hz so it fits plotting with logarithmic
12861292
scale in Hz otherwise in rad/s. Omega is always returned in rad/sec.
@@ -1326,7 +1332,7 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
13261332
features_ = np.concatenate((np.abs(sys.pole()),
13271333
np.abs(sys.zero())))
13281334
# Get rid of poles and zeros at the origin
1329-
toreplace = features_ == 0.0
1335+
toreplace = np.isclose(features_, 0.0)
13301336
if np.any(toreplace):
13311337
features_ = features_[~toreplace]
13321338
elif sys.isdtime(strict=True):
@@ -1339,7 +1345,7 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
13391345
# Get rid of poles and zeros on the real axis (imag==0)
13401346
# * origin and real < 0
13411347
# * at 1.: would result in omega=0. (logaritmic plot!)
1342-
toreplace = (features_.imag == 0.0) & (
1348+
toreplace = np.isclose(features_.imag, 0.0) & (
13431349
(features_.real <= 0.) |
13441350
(np.abs(features_.real - 1.0) < 1.e-10))
13451351
if np.any(toreplace):
@@ -1360,15 +1366,13 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
13601366

13611367
if Hz:
13621368
features /= 2. * math.pi
1363-
features = np.log10(features)
1364-
lsp_min = np.floor(np.min(features) - feature_periphery_decades)
1365-
lsp_max = np.ceil(np.max(features) + feature_periphery_decades)
1369+
features = np.log10(features)
1370+
lsp_min = np.rint(np.min(features) - feature_periphery_decades)
1371+
lsp_max = np.rint(np.max(features) + feature_periphery_decades)
1372+
if Hz:
13661373
lsp_min += np.log10(2. * math.pi)
13671374
lsp_max += np.log10(2. * math.pi)
1368-
else:
1369-
features = np.log10(features)
1370-
lsp_min = np.floor(np.min(features) - feature_periphery_decades)
1371-
lsp_max = np.ceil(np.max(features) + feature_periphery_decades)
1375+
13721376
if freq_interesting:
13731377
lsp_min = min(lsp_min, np.log10(min(freq_interesting)))
13741378
lsp_max = max(lsp_max, np.log10(max(freq_interesting)))

‎control/tests/nyquist_test.py

Lines changed: 25 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -182,42 +182,56 @@ def test_nyquist_encirclements():
182182
assert _Z(sys) == count + _P(sys)
183183

184184

185-
def test_nyquist_indent():
185+
@pytest.fixture
186+
def indentsys():
186187
# FBS Figure 10.10
187-
s = ct.tf('s')
188-
sys = 3 * (s+6)**2 / (s * (s+1)**2)
189188
# poles: [-1, -1, 0]
189+
s = ct.tf('s')
190+
return 3 * (s+6)**2 / (s * (s+1)**2)
190191

192+
193+
def test_nyquist_indent_default(indentsys):
191194
plt.figure();
192-
count = ct.nyquist_plot(sys)
195+
count = ct.nyquist_plot(indentsys)
193196
plt.title("Pole at origin; indent_radius=default")
194-
assert _Z(sys) == count + _P(sys)
197+
assert _Z(indentsys) == count + _P(indentsys)
195198

199+
200+
def test_nyquist_indent_dont(indentsys):
196201
# first value of default omega vector was 0.1, replaced by 0. for contour
197202
# 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,
203+
count, contour = ct.nyquist_plot(indentsys,
204+
plot=False,
205+
indent_radius=.1007,
199206
return_contour=True)
200207
np.testing.assert_allclose(contour[0], .1007+0.j)
201208
# second value of omega_vector is larger than indent_radius: not indented
202209
assert np.all(contour.real[2:] == 0.)
203210

211+
212+
def test_nyquist_indent_do(indentsys):
204213
plt.figure();
205-
count, contour = ct.nyquist_plot(sys, indent_radius=0.01,
214+
count, contour = ct.nyquist_plot(indentsys,
215+
indent_radius=0.01,
206216
return_contour=True)
207217
plt.title("Pole at origin; indent_radius=0.01; encirclements = %d" % count)
208-
assert _Z(sys) == count + _P(sys)
218+
assert _Z(indentsys) == count + _P(indentsys)
209219
# indent radius is smaller than the start of the default omega vector
210220
# check that a quarter circle around the pole at origin has been added.
211221
np.testing.assert_allclose(contour[:50].real**2 + contour[:50].imag**2,
212222
0.01**2)
213223

224+
225+
def test_nyquist_indent_left(indentsys):
214226
plt.figure();
215-
count = ct.nyquist_plot(sys, indent_direction='left')
227+
count = ct.nyquist_plot(indentsys, indent_direction='left')
216228
plt.title(
217229
"Pole at origin; indent_direction='left'; encirclements = %d" % count)
218-
assert _Z(sys) == count + _P(sys, indent='left')
230+
assert _Z(indentsys) == count + _P(indentsys, indent='left')
231+
219232

220-
# System with poles on the imaginary axis
233+
def test_nyquist_indent_im():
234+
"""Test system with poles on the imaginary axis."""
221235
sys = ct.tf([1, 1], [1, 0, 1])
222236

223237
# Imaginary poles with standard indentation

‎control/tests/sisotool_test.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -102,8 +102,8 @@ def test_sisotool(self, tsys):
102102

103103
# Check if the bode_mag line has moved
104104
bode_mag_moved = np.array(
105-
[674.0242, 667.8354, 661.7033, 655.6275, 649.6074, 643.6426,
106-
637.7324, 631.8765, 626.0742, 620.3252])
105+
[69.0065, 68.6749, 68.3448, 68.0161, 67.6889, 67.3631, 67.0388,
106+
66.7159, 66.3944, 66.0743])
107107
assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],
108108
bode_mag_moved, 4)
109109

0 commit comments

Comments
 (0)
0