Skip to content

Commit 3b3c005

Browse files
RoboX2020RoboX2020
committed
Fix discrete-time frequency handling in describing_function_response
For a discrete-time system, describing_function_response found limit cycles with about the right amplitude but reported the wrong frequency: it took the oscillation frequency as the imaginary part of the Nyquist contour, which for a discrete system lies on the unit circle exp(1j*w*dt), so it reported about sin(w*dt) - the same number whatever dt is. The refinement cost function and the returned positions also evaluated H at z = 1j*w rather than on the unit circle, and refinement typically failed with "not able to refine result; returning estimate". For example, H(z) = 0.5/(z^2 - 0.5z) with dt = 1 crosses the negative real axis at w = 1.3181, but the reported frequency was 0.9703, and dt = 0.5 gave the same 0.9703 instead of 2.6362. When H.isdtime(strict=True) - matching the timebase semantics used by nyquist_response - take the frequencies as np.angle(contour) / dt and evaluate H at exp(1j*w*dt) in the refinement cost function and for the returned positions. The strict check keeps dt=None (unspecified timebase) on the continuous-time path; dt=True (discrete, unspecified sample time) reports rad/sample. Adds regression tests: a known unit-circle crossing for dt = 1 and dt = 0.5, plus dt=None continuous-equivalence and dt=True coverage. The discrete tests fail at baseline; continuous-time behavior is unchanged. Fixes #1269
1 parent bfe9380 commit 3b3c005

2 files changed

Lines changed: 63 additions & 3 deletions

File tree

‎control/descfcn.py‎

Lines changed: 15 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -346,7 +346,19 @@ def describing_function_response(
346346
response = nyquist_response(
347347
H, omega, warn_encirclements=warn_nyquist, warn_nyquist=warn_nyquist,
348348
_check_kwargs=_check_kwargs, **kwargs)
349-
H_omega, H_vals = response.contour.imag, H(response.contour)
349+
# For a discrete-time system the contour lies on the unit circle,
350+
# exp(1j * omega * dt): the frequency is its angle divided by dt,
351+
# not its imaginary part, and H must be evaluated on the unit
352+
# circle rather than at 1j * omega
353+
if H.isdtime(strict=True):
354+
H_omega = np.angle(response.contour) / H.dt
355+
def _H_eval(w):
356+
return H(np.exp(1j * w * H.dt))
357+
else:
358+
H_omega = response.contour.imag
359+
def _H_eval(w):
360+
return H(1j * w)
361+
H_vals = H(response.contour)
350362

351363
# Compute the describing function
352364
df = describing_function(F, A)
@@ -375,7 +387,7 @@ def _cost(x):
375387
# Note: imposing bounds messed up the optimization (?)
376388
if x[0] < 0 or x[1] < 0:
377389
return 1
378-
return abs(1 + H(1j * x[1]) *
390+
return abs(1 + _H_eval(x[1]) *
379391
describing_function(F, x[0]))**2
380392
res = scipy.optimize.minimize(
381393
_cost, [a_guess, omega_guess])
@@ -386,7 +398,7 @@ def _cost(x):
386398
else:
387399
a_final, omega_final = res.x[0], res.x[1]
388400

389-
pos = H(1j * omega_final)
401+
pos = _H_eval(omega_final)
390402

391403
# Save the final estimate
392404
positions.append(pos)

‎control/tests/descfcn_test.py‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,54 @@ def test_describing_function_response():
172172
H_multiple(1j*w), decimal=5)
173173

174174

175+
@pytest.mark.parametrize("dt,omega_ref", [(1, 1.3181), (0.5, 2.6362)])
176+
def test_describing_function_response_discrete(dt, omega_ref):
177+
"""Discrete-time systems: frequency on the unit circle, not sin(w*dt).
178+
179+
gh-1269: describing_function_response reported about sin(w*dt) for
180+
the oscillation frequency and evaluated H at z = 1j*w instead of
181+
exp(1j*w*dt). The crossing amplitude and frequency for this loop
182+
are known from a dense unit-circle scan: H(exp(j*w*dt)) crosses the
183+
negative real axis at w = omega_ref with H = -0.5, giving
184+
a = 4/pi * 0.5 = 0.6366 for the sign nonlinearity.
185+
"""
186+
H = ct.tf([0.5], [1, -0.5, 0], dt=dt)
187+
response = ct.describing_function_response(
188+
H, np.sign, np.linspace(0.01, 5, 400),
189+
omega=np.linspace(1e-3, np.pi / dt, 2000))
190+
assert len(response.intersections) == 1
191+
a, w = response.intersections[0]
192+
assert w == pytest.approx(omega_ref, rel=2e-2)
193+
assert a == pytest.approx(0.6366, rel=2e-2)
194+
# the returned position is the response on the unit circle
195+
np.testing.assert_almost_equal(
196+
response.positions[0], H(np.exp(1j * w * dt)), decimal=5)
197+
198+
199+
def test_describing_function_response_timebase():
200+
"""Timebase edge cases (gh-1269 review).
201+
202+
dt=None means an unspecified (continuous-compatible) timebase and
203+
must behave exactly like continuous time; dt=True is discrete time
204+
with unspecified sample time, with frequency in rad/sample.
205+
"""
206+
amp = np.linspace(1, 4, 10)
207+
F = ct.saturation_nonlinearity(1)
208+
xsects_ref = ct.describing_function_response(
209+
ct.tf([8], [1, 2, 2, 1]), F, amp)
210+
xsects_none = ct.describing_function_response(
211+
ct.tf([8], [1, 2, 2, 1], dt=None), F, amp)
212+
np.testing.assert_allclose(xsects_ref.intersections,
213+
xsects_none.intersections)
214+
215+
H_true = ct.tf([0.5], [1, -0.5, 0], dt=True)
216+
xsects_true = ct.describing_function_response(
217+
H_true, np.sign, np.linspace(0.01, 5, 400),
218+
omega=np.linspace(1e-3, np.pi, 2000))
219+
a, w = xsects_true.intersections[0]
220+
assert w == pytest.approx(1.3181, rel=2e-2)
221+
222+
175223
def test_describing_function_plot():
176224
# Simple linear system with at most 1 intersection
177225
H_larger = ct.tf([1], [1, 2, 2, 1]) * 8

0 commit comments

Comments
 (0)