Skip to content

Fix discrete-time frequency handling in describing_function_response #1270

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
RoboX2020 wants to merge 1 commit into from
Closed

Fix discrete-time frequency handling in describing_function_response #1270

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 15 additions & 3 deletions control/descfcn.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,19 @@ def describing_function_response(
response = nyquist_response(
H, omega, warn_encirclements=warn_nyquist, warn_nyquist=warn_nyquist,
_check_kwargs=_check_kwargs, **kwargs)
H_omega, H_vals = response.contour.imag, H(response.contour)
# For a discrete-time system the contour lies on the unit circle,
# exp(1j * omega * dt): the frequency is its angle divided by dt,
# not its imaginary part, and H must be evaluated on the unit
# circle rather than at 1j * omega
if H.isdtime(strict=True):
H_omega = np.angle(response.contour) / H.dt
def _H_eval(w):
return H(np.exp(1j * w * H.dt))
else:
H_omega = response.contour.imag
def _H_eval(w):
return H(1j * w)
H_vals = H(response.contour)

# Compute the describing function
df = describing_function(F, A)
Expand Down Expand Up @@ -375,7 +387,7 @@ def _cost(x):
# Note: imposing bounds messed up the optimization (?)
if x[0] < 0 or x[1] < 0:
return 1
return abs(1 + H(1j * x[1]) *
return abs(1 + _H_eval(x[1]) *
describing_function(F, x[0]))**2
res = scipy.optimize.minimize(
_cost, [a_guess, omega_guess])
Expand All @@ -386,7 +398,7 @@ def _cost(x):
else:
a_final, omega_final = res.x[0], res.x[1]

pos = H(1j * omega_final)
pos = _H_eval(omega_final)

# Save the final estimate
positions.append(pos)
Expand Down
48 changes: 48 additions & 0 deletions control/tests/descfcn_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,54 @@ def test_describing_function_response():
H_multiple(1j*w), decimal=5)


@pytest.mark.parametrize("dt,omega_ref", [(1, 1.3181), (0.5, 2.6362)])
def test_describing_function_response_discrete(dt, omega_ref):
"""Discrete-time systems: frequency on the unit circle, not sin(w*dt).

gh-1269: describing_function_response reported about sin(w*dt) for
the oscillation frequency and evaluated H at z = 1j*w instead of
exp(1j*w*dt). The crossing amplitude and frequency for this loop
are known from a dense unit-circle scan: H(exp(j*w*dt)) crosses the
negative real axis at w = omega_ref with H = -0.5, giving
a = 4/pi * 0.5 = 0.6366 for the sign nonlinearity.
"""
H = ct.tf([0.5], [1, -0.5, 0], dt=dt)
response = ct.describing_function_response(
H, np.sign, np.linspace(0.01, 5, 400),
omega=np.linspace(1e-3, np.pi / dt, 2000))
assert len(response.intersections) == 1
a, w = response.intersections[0]
assert w == pytest.approx(omega_ref, rel=2e-2)
assert a == pytest.approx(0.6366, rel=2e-2)
# the returned position is the response on the unit circle
np.testing.assert_almost_equal(
response.positions[0], H(np.exp(1j * w * dt)), decimal=5)


def test_describing_function_response_timebase():
"""Timebase edge cases (gh-1269 review).

dt=None means an unspecified (continuous-compatible) timebase and
must behave exactly like continuous time; dt=True is discrete time
with unspecified sample time, with frequency in rad/sample.
"""
amp = np.linspace(1, 4, 10)
F = ct.saturation_nonlinearity(1)
xsects_ref = ct.describing_function_response(
ct.tf([8], [1, 2, 2, 1]), F, amp)
xsects_none = ct.describing_function_response(
ct.tf([8], [1, 2, 2, 1], dt=None), F, amp)
np.testing.assert_allclose(xsects_ref.intersections,
xsects_none.intersections)

H_true = ct.tf([0.5], [1, -0.5, 0], dt=True)
xsects_true = ct.describing_function_response(
H_true, np.sign, np.linspace(0.01, 5, 400),
omega=np.linspace(1e-3, np.pi, 2000))
a, w = xsects_true.intersections[0]
assert w == pytest.approx(1.3181, rel=2e-2)


def test_describing_function_plot():
# Simple linear system with at most 1 intersection
H_larger = ct.tf([1], [1, 2, 2, 1]) * 8
Expand Down