Control systems exercises#
Three problems built on the control systems gallery. Try each one before opening its solution. The solutions run when the documentation is built, and the test suite checks every number they print.
1. The largest stable gain#
A plant \(G(s) = 1/((s+1)(s+2)(s+3))\) is under proportional control with gain \(K\). Use Routh’s criterion to find the range of \(K\) for which the closed loop is stable, and the frequency at which it oscillates at the limit. Check with the closed-loop poles.
Builds on: Routh’s stability criterion.
Solution
import numpy as np
from engineeringkit.control_systems import routh_first_column, routh_stable
den = np.polymul(np.polymul([1, 1], [1, 2]), [1, 3]) # s^3 + 6 s^2 + 11 s + 6
print("characteristic polynomial without K:", den)
# characteristic polynomial without K: [ 1 6 11 6]
# The s^1 row of the Routh array is (6*11 - (6 + K)) / 6, which vanishes at K = 60.
for K in (59.0, 61.0):
print(f"K = {K:.0f}: first column {np.round(routh_first_column(den + [0, 0, 0, K]), 3)}, stable {routh_stable(den + [0, 0, 0, K])}")
# K = 59: first column [ 1. 6. 0.167 65. ], stable True
# K = 61: first column [ 1. 6. -0.167 67. ], stable False
poles = np.roots(den + [0, 0, 0, 60.0])
print("poles at K = 60:", np.round(np.sort_complex(poles), 4))
# poles at K = 60: [-6.+0.j 0.-3.3166j 0.+3.3166j]
assert np.isclose(np.abs(poles.imag).max(), np.sqrt(11))
The loop is stable for \(0 < K < 60\). At \(K = 60\) the auxiliary polynomial \(6s^2 + 66\) gives a pair of poles at \(\pm j\sqrt{11}\), an oscillation at 3.317 rad/s.
2. A servomechanism to specification#
A position servo behaves as a standard second-order system. It must overshoot by no more than 10% and settle (2%) within 2 s. Find the smallest damping ratio and natural frequency that do it, then simulate the step response and measure both.
Builds on: Hazen’s theory of servomechanisms.
Solution
import numpy as np
from engineeringkit.control_systems import second_order_metrics, step_response, transfer_function
ln_os = np.log(0.10)
zeta = -ln_os / np.sqrt(np.pi**2 + ln_os**2)
wn = 4 / (zeta * 2.0)
metrics = second_order_metrics(zeta, wn)
print(f"zeta = {zeta:.4f}, wn = {wn:.3f} rad/s")
print(f"predicted: overshoot {metrics['percent_overshoot']:.2f}%, 2% settling {metrics['settling_time_2pct']:.2f} s")
# zeta = 0.5912, wn = 3.383 rad/s
# predicted: overshoot 10.00%, 2% settling 2.00 s
t, y = step_response(transfer_function([wn**2], [1, 2 * zeta * wn, wn**2]), T=np.linspace(0, 5, 5001))
outside = np.flatnonzero(abs(y - 1) > 0.02)
print(f"simulated: overshoot {100 * (y.max() - 1):.2f}%, settles at {t[outside[-1] + 1]:.2f} s")
# simulated: overshoot 10.00%, settles at 1.75 s
assert np.isclose(100 * (y.max() - 1), 10.0, atol=0.01)
The overshoot formula is exact. The settling time \(4/(\zeta\omega_n)\) is an envelope bound, so the simulated response settles a little sooner.
3. An LQR for a double integrator#
A mass on a frictionless track, \(\ddot x = u\), has state \((x, \dot x)\). Find the LQR gain for \(Q = I\) and \(R = 1\), and check it against the hand solution \(K = [1, \sqrt 3]\). Use Lyapunov’s test to confirm that the closed loop is stable.
Builds on: Lyapunov’s direct method.
Solution
import numpy as np
from engineeringkit.control_systems import lqr, lyapunov_stable
A = np.array([[0.0, 1.0], [0.0, 0.0]])
B = np.array([[0.0], [1.0]])
K, P, closed_poles = lqr(A, B, np.eye(2), np.eye(1))
print("K =", np.round(K, 6), " closed-loop poles:", np.round(closed_poles, 4))
print("Lyapunov stable:", lyapunov_stable(A - B @ K), "(open loop:", lyapunov_stable(A), ")")
# K = [[1. 1.732051]] closed-loop poles: [-0.866+0.5j -0.866-0.5j]
# Lyapunov stable: True (open loop: False )
assert np.allclose(K, [[1.0, np.sqrt(3)]])
The closed-loop poles are \((-\sqrt3 \pm j)/2\), a damping ratio of 0.866: the LQR with equal weights gives a well-damped loop without any pole placement by hand.