Skip to content

Commit 8fe71bc

Browse files
dongxuelian2dongxuelian2
committed
Fix Nyquist indentation geometry and automatic direction selection
Assisted-by: OpenAI Codex
1 parent bfe9380 commit 8fe71bc

3 files changed

Lines changed: 245 additions & 82 deletions

File tree

‎control/freqplot.py‎

Lines changed: 114 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -1104,7 +1104,7 @@ def gen_zero_centered_series(val_min, val_max, period):
11041104
'nyquist.arrow_size': 8, # pixel size for arrows
11051105
'nyquist.encirclement_threshold': 0.05, # warning threshold
11061106
'nyquist.indent_radius': 1e-4, # indentation radius
1107-
'nyquist.indent_direction': 'right', # indentation direction
1107+
'nyquist.indent_direction': None, # automatic indentation
11081108
'nyquist.indent_points': 200, # number of points to insert
11091109
'nyquist.max_curve_magnitude': 15, # rescale large values
11101110
'nyquist.blend_fraction': 0.15, # when to start scaling
@@ -1240,10 +1240,13 @@ def nyquist_response(
12401240
Define the threshold for generating a warning if the number of net
12411241
encirclements is a non-integer value. Default value is 0.05 and can
12421242
be set using `config.defaults['nyquist.encirclement_threshold']`.
1243-
indent_direction : str, optional
1244-
For poles on the imaginary axis, set the direction of indentation to
1245-
be 'right' (default), 'left', or 'none'. The default value can
1246-
be set using `config.defaults['nyquist.indent_direction']`.
1243+
indent_direction : str or None, optional
1244+
Direction to indent around poles within `indent_radius` of the
1245+
imaginary axis: 'right', 'left', or 'none' (disable indentation).
1246+
The default, None, indents to the right of stable poles and poles
1247+
exactly on the imaginary axis, and to the left of unstable poles.
1248+
The default value can be set using
1249+
`config.defaults['nyquist.indent_direction']`.
12471250
indent_points : int, optional
12481251
Number of points to insert in the Nyquist contour around poles that
12491252
are at or near the imaginary axis.
@@ -1277,10 +1280,18 @@ def nyquist_response(
12771280
If a continuous-time system contains poles on or near the imaginary
12781281
axis, a small indentation will be used to avoid the pole. The radius
12791282
of the indentation is given by `indent_radius` and it is taken to the
1280-
right of stable poles and the left of unstable poles. If a pole is
1281-
exactly on the imaginary axis, the `indent_direction` parameter can be
1282-
used to set the direction of indentation. Setting `indent_direction`
1283-
to 'none' will turn off indentation.
1283+
right of stable poles and the left of unstable poles by default. Poles
1284+
exactly on the imaginary axis are indented to the right. Setting
1285+
`indent_direction` to 'left' or 'right' overrides this choice for all
1286+
poles within `indent_radius` of the axis; 'none' turns off indentation.
1287+
Explicit directions can change which open-loop poles are enclosed by
1288+
the contour, and hence the number of encirclements.
1289+
1290+
When `omega` is supplied, no extra points are added for indentations.
1291+
The frequencies must resolve each detour sufficiently. Automatic
1292+
frequency selection adds `indent_points` around nearby poles. A forced
1293+
detour around the far side of an off-axis pole can change both the real
1294+
and imaginary parts of the corresponding contour points.
12841295
12851296
For those portions of the Nyquist plot in which the contour is indented
12861297
to avoid poles, resulting in a scaling of the Nyquist plot, the line
@@ -1312,8 +1323,6 @@ def nyquist_response(
13121323
# Use existing dictionary, to keep track of processed keywords
13131324
_kwargs |= kwargs
13141325

1315-
indent_direction_arg = 'indent_direction' in _kwargs
1316-
13171326
# Get values for params
13181327
omega_num_given = omega_num is not None
13191328
omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)
@@ -1322,11 +1331,16 @@ def nyquist_response(
13221331
encirclement_threshold = config._get_param(
13231332
'nyquist', 'encirclement_threshold', _kwargs,
13241333
_nyquist_defaults, pop=True)
1325-
indent_direction = config._get_param(
1326-
'nyquist', 'indent_direction', _kwargs, _nyquist_defaults, pop=True)
1334+
# Unlike most configuration parameters, None explicitly selects automatic
1335+
# indentation, even if the configured default specifies a direction.
1336+
indent_direction = _kwargs.pop('indent_direction', config._get_param(
1337+
'nyquist', 'indent_direction', defval=_nyquist_defaults))
13271338
indent_points = config._get_param(
13281339
'nyquist', 'indent_points', _kwargs, _nyquist_defaults, pop=True)
13291340

1341+
if indent_direction not in (None, 'left', 'right', 'none'):
1342+
raise ValueError("unknown value for indent_direction")
1343+
13301344
if _check_kwargs and _kwargs:
13311345
raise TypeError("unrecognized keywords: ", str(_kwargs))
13321346

@@ -1401,43 +1415,53 @@ def nyquist_response(
14011415
#
14021416
# Check to make sure indent radius is small enough
14031417
#
1404-
# If there is a closed loop pole that is near the imaginary axis
1405-
# at a point that is near an open loop pole, it is possible that
1418+
# If there is a closed loop pole inside an indentation circle,
14061419
# indentation might skip or create an extraneous encirclement.
1420+
# Forced detours can reach past the imaginary-axis strip of
1421+
# width indent_radius, so check proximity to the indented poles.
14071422
# We check for that situation here and generate a warning if that
14081423
# could happen.
14091424
#
1425+
indent_poles = splane_poles[abs(splane_poles.real) < indent_radius]
14101426
for p_cl in splane_cl_poles:
1411-
# See if any closed loop poles are near the imaginary axis
1412-
if abs(p_cl.real) <= indent_radius:
1427+
if len(indent_poles) > 0:
14131428
# See if any open loop poles are close to closed loop poles
1414-
if len(splane_poles) > 0:
1415-
p_ol = splane_poles[
1416-
(np.abs(splane_poles - p_cl)).argmin()]
1429+
p_ol = indent_poles[
1430+
(np.abs(indent_poles - p_cl)).argmin()]
14171431

1418-
if abs(p_ol - p_cl) <= indent_radius and \
1419-
warn_encirclements:
1420-
warnings.warn(
1421-
"indented contour may miss closed loop pole; "
1422-
"consider reducing indent_radius to below "
1423-
f"{abs(p_ol - p_cl):5.2g}", stacklevel=2)
1432+
if abs(p_ol - p_cl) <= indent_radius and \
1433+
warn_encirclements:
1434+
warnings.warn(
1435+
"indented contour may miss closed loop pole; "
1436+
"consider reducing indent_radius to below "
1437+
f"{abs(p_ol - p_cl):5.2g}", stacklevel=2)
14241438

14251439
#
14261440
# See if we should add some frequency points near imaginary poles
14271441
#
14281442
for p in splane_poles:
14291443
# See if we need to process this pole (skip if on the negative
14301444
# imaginary axis or not near imaginary axis + user override)
1431-
if p.imag < 0 or abs(p.real) > indent_radius or \
1445+
if p.imag < 0 or abs(p.real) >= indent_radius or \
14321446
omega_range_given:
14331447
continue
14341448

1449+
# Find the circle's intersections with the imaginary axis.
1450+
half_height = np.sqrt(indent_radius**2 - p.real**2)
1451+
start_freq = max(0, p.imag - half_height)
1452+
end_freq = p.imag + half_height
1453+
side = 1 if indent_direction == 'right' or (
1454+
indent_direction is None and p.real <= 0) else -1
1455+
if side * p.real > 0 and p.imag < indent_radius:
1456+
# A major arc crossing the real axis starts there, so that
1457+
# the primary contour still joins its mirror image.
1458+
start_freq = 0
1459+
14351460
# Find the frequencies before the pole frequency
14361461
below_points = np.argwhere(
1437-
splane_contour.imag - abs(p.imag) < -indent_radius)
1462+
splane_contour.imag < start_freq)
14381463
if below_points.size > 0:
14391464
first_point = below_points[-1].item()
1440-
start_freq = p.imag - indent_radius
14411465
else:
14421466
# Add the points starting at the beginning of the contour
14431467
assert splane_contour[0] == 0
@@ -1446,57 +1470,61 @@ def nyquist_response(
14461470

14471471
# Find the frequencies after the pole frequency
14481472
above_points = np.argwhere(
1449-
splane_contour.imag - abs(p.imag) > indent_radius)
1473+
splane_contour.imag > end_freq)
14501474
last_point = above_points[0].item()
14511475

14521476
# Add points for half/quarter circle around pole frequency
14531477
# (these will get indented left or right below)
14541478
splane_contour = np.concatenate((
14551479
splane_contour[0:first_point+1],
14561480
(1j * np.linspace(
1457-
start_freq, p.imag + indent_radius, indent_points)),
1481+
start_freq, end_freq, indent_points)),
14581482
splane_contour[last_point:]))
14591483

14601484
# Indent points that are too close to a pole
14611485
if len(splane_poles) > 0: # accommodate no splane poles if dtime sys
14621486
for i, s in enumerate(splane_contour):
14631487
# Find the nearest pole
14641488
p = splane_poles[(np.abs(splane_poles - s)).argmin()]
1489+
if s.imag == 0 and p.imag < 0:
1490+
p = p.conjugate()
14651491

1466-
# See if we need to indent around it
1467-
if abs(s - p) < indent_radius:
1468-
# Figure out how much to offset (simple trigonometry)
1469-
offset = np.sqrt(
1470-
indent_radius ** 2 - (s - p).imag ** 2) \
1471-
- (s - p).real
1472-
1473-
# Figure out which way to offset the contour point
1474-
if indent_direction_arg:
1475-
if indent_direction == 'right':
1476-
# Indent to the right
1477-
splane_contour[i] += offset
1478-
1479-
elif indent_direction == 'left':
1480-
# Indent to the left
1481-
splane_contour[i] -= offset
1482-
1483-
else:
1484-
raise ValueError(
1485-
"unknown value for indent_direction")
1486-
1487-
elif p.real < 0 or (p.real == 0 and
1488-
indent_direction == 'right'):
1489-
# Indent to the right
1490-
splane_contour[i] += offset
1491-
1492-
elif p.real > 0 or (p.real == 0 and
1493-
indent_direction == 'left'):
1494-
# Indent to the left
1495-
splane_contour[i] -= offset
1492+
side = 1 if indent_direction == 'right' or (
1493+
indent_direction is None and p.real <= 0) else -1
1494+
major_arc = side * p.real > 0
1495+
crosses_origin = major_arc and 0 <= p.imag < indent_radius
14961496

1497+
# See if we need to indent around it
1498+
if abs(p.real) >= indent_radius:
1499+
continue
1500+
half_height = np.sqrt(indent_radius**2 - p.real**2)
1501+
if abs(s - p) < indent_radius or (
1502+
crosses_origin and
1503+
0 <= s.imag <= p.imag + half_height):
1504+
if major_arc:
1505+
# The far side of an off-axis pole requires a major
1506+
# arc to join the imaginary axis continuously.
1507+
angle = np.arccos(-side * p.real / indent_radius)
1508+
start_freq = max(0, p.imag - half_height)
1509+
start_angle = -angle
1510+
if crosses_origin:
1511+
start_freq = 0
1512+
start_angle = -np.arcsin(
1513+
p.imag / indent_radius)
1514+
fraction = (s.imag - start_freq) / (
1515+
p.imag + half_height - start_freq)
1516+
theta = start_angle + fraction * (
1517+
angle - start_angle)
1518+
splane_contour[i] = p + indent_radius * (
1519+
side * np.cos(theta) + 1j * np.sin(theta))
1520+
if s.imag == 0:
1521+
splane_contour[i] = splane_contour[i].real
14971522
else:
1498-
raise ValueError(
1499-
"unknown value for indent_direction")
1523+
# Preserve the frequencies on the near-side arc.
1524+
offset = np.sqrt(
1525+
indent_radius**2 - (s - p).imag**2)
1526+
splane_contour[i] = (
1527+
p.real + side * offset + 1j * s.imag)
15001528

15011529
# change contour to z-plane if necessary
15021530
if sys.isctime():
@@ -1536,18 +1564,22 @@ def nyquist_response(
15361564
if isinstance(sys, (StateSpace, TransferFunction)):
15371565
# Count the number of open/closed loop RHP poles
15381566
if sys.isctime():
1539-
if indent_direction == 'right':
1540-
P = (sys.poles().real > 0).sum()
1541-
else:
1542-
P = (sys.poles().real >= 0).sum()
1567+
P = (sys.poles().real >= 0).sum()
15431568
Z = (sys.feedback().poles().real >= 0).sum()
15441569
else:
1545-
if indent_direction == 'right':
1546-
P = (np.abs(sys.poles()) > 1).sum()
1547-
else:
1548-
P = (np.abs(sys.poles()) >= 1).sum()
1570+
P = (np.abs(sys.poles()) >= 1).sum()
15491571
Z = (np.abs(sys.feedback().poles()) >= 1).sum()
15501572

1573+
if indent_direction != 'none':
1574+
# Explicit directions can include stable poles or exclude
1575+
# unstable poles. Use the mapped poles for discrete systems
1576+
# so the pole count and the contour use the same boundary.
1577+
enclosed = splane_poles.real > 0
1578+
if indent_direction is not None:
1579+
enclosed[abs(splane_poles.real) < indent_radius] = (
1580+
indent_direction == 'left')
1581+
P = enclosed.sum()
1582+
15511583
# Check to make sure the results make sense; warn if not
15521584
if Z != count + P and warn_encirclements:
15531585
warnings.warn(
@@ -1662,9 +1694,10 @@ def nyquist_plot(
16621694
Define the threshold for generating a warning if the number of net
16631695
encirclements is a non-integer value. Default value is 0.05 and can
16641696
be set using `config.defaults['nyquist.encirclement_threshold']`.
1665-
indent_direction : str, optional
1666-
For poles on the imaginary axis, set the direction of indentation to
1667-
be 'right' (default), 'left', or 'none'.
1697+
indent_direction : str or None, optional
1698+
Direction to indent around poles on or near the imaginary axis:
1699+
'right', 'left', or 'none' (disable indentation). The default, None,
1700+
selects the side automatically; see `nyquist_response`.
16681701
indent_points : int, optional
16691702
Number of points to insert in the Nyquist contour around poles that
16701703
are at or near the imaginary axis.
@@ -1761,11 +1794,13 @@ def nyquist_plot(
17611794
If a continuous-time system contains poles on or near the imaginary
17621795
axis, a small indentation will be used to avoid the pole. The radius
17631796
of the indentation is given by `indent_radius` and it is taken to the
1764-
right of stable poles and the left of unstable poles. If a pole is
1765-
exactly on the imaginary axis, the `indent_direction` parameter can be
1766-
used to set the direction of indentation. Setting `indent_direction`
1767-
to 'none' will turn off indentation. If `return_contour` is True,
1768-
the exact contour used for evaluation is returned.
1797+
right of stable poles and the left of unstable poles by default. Poles
1798+
exactly on the imaginary axis are indented to the right. Setting
1799+
`indent_direction` to 'left' or 'right' overrides this choice for all
1800+
poles within `indent_radius` of the axis; 'none' turns off indentation.
1801+
Explicit directions can change which open-loop poles are enclosed by
1802+
the contour, and hence the number of encirclements. If `return_contour`
1803+
is True, the exact contour used for evaluation is returned.
17691804
17701805
For those portions of the Nyquist plot in which the contour is indented
17711806
to avoid poles, resulting in a scaling of the Nyquist plot, the line

0 commit comments

Comments
 (0)