🐍 | Code | 🔍 | Worked Example |
📈 | Graph | 🧩 | Exercise |
🏷️ | Definition | 💻 | Numerical Method |
💎 | Theorem | 🧮 | Analytical Method |
📝 | Remark | 🧠 | Theory |
ℹ️ | Information | 🗝️ | Hint |
⚠️ | Warning | 🔓 | Solution |
from numpy import *
import matplotlib; matplotlib.use("nbAgg")
%matplotlib notebook
from matplotlib.pyplot import *
When the system
is not asymptotically stable at the origin, maybe there are some inputs $u \in \mathbb{R}^m$ such that
that we can use to stabilize asymptotically the system?
We search for $u$ as a linear feedback:
for some $K \in \mathbb{R}^{m \times n}$.
In this scheme
⚠️ The full system state $x(t)$ must be measured.
🏷️ This information is then fed back into the system.
🏷️ The feedback closes the loop.
When
the state $x \in \mathbb{R}^n$ evolves according to:
The closed-loop system is asymptotically stable iff every eigenvalue of the matrix
is in the open left-hand plane.
Multisets remember the multiplicity of their elements. It’s convenient to describe the spectrum of matrices:
The system
$$ \dot{x} = A x + B u, \, x \in \mathbb{R}^n, \,u\in\mathbb{R}^p $$is controllable.
$\Lambda$ is a symmetric multiset of $n$ complex numbers:
$$ \Lambda = \{\lambda_1, \dots, \lambda_n\} \subset \mathbb{C} \; \mbox{ and } \; \lambda \in^k \Lambda \, \Rightarrow \, \overline{\lambda} \in^k \Lambda. $$$\Rightarrow$ There is a matrix $K \in \mathbb{R}^{n \times m}$ such that
Consider the double integrator $\ddot{x} = u$
(in standard form)
from scipy.signal import place_poles
A = array([[0, 1], [0, 0]])
B = array([[0], [1]])
poles = [-1, -2]
K = place_poles(A, B, poles).gain_matrix
assert_almost_equal(K, [[2.0, 3.0]])
eigenvalues, _ = eig(A - B @ K)
assert_almost_equal(eigenvalues, [-1, -2])
figure()
x = [real(s) for s in eigenvalues]
y = [imag(s) for s in eigenvalues]
plot(x, y, "kx")
xticks([-3, -2,-1, 0,1, 2,3])
yticks([-3, -2,-1, 0,1, 2,3])
plot([0, 0], [-3, 3], "k")
plot([-3, 3], [0, 0], "k")
title("Spectrum of $A-BK$"); grid(True)
⛔ The place_poles
function rejects eigenvalues whose multiplicity is
higher than the rank of $B$.
In the previous example, $\mbox{rank}\, B = 1$, so
❌ poles = [-1, -1]
won’t work.
✔️ poles = [-1, -2]
will.
Consider the system with dynamics
We apply the control law
Can we assign the poles of the closed-loop system freely by a suitable choice of $k_1$ and $k_2$?
Explain this result.
Since the characteristic polynomial is also $$ (s-\lambda_1)(s-\lambda_2) $$ we get $$ k_1 + k_2 = - \lambda_1 - \lambda_2, \; -2 (k_1 + k_2) = \lambda_1 \lambda_2 $$
Thus we have
and both poles cannot be assigned freely; for example if we select $\lambda_1 =1$, we end up with $\lambda_2 = -2$.
We have not checked the assumptions of 💎 Pole Assignment yet.
The commandability matrix is $$ \left[B, AB \right] = \left[ \begin{array}{rr} 1 & 0 \\ 1 & 0 \end{array} \right] $$ whose rank is $1 < 2$.
Since the system is not controllable, pole assignment may fail and it does here.
Consider the pendulum with dynamics:
Numerical Values:
Compute the linearized dynamics of the system around the equilibrium $\theta=\pi$ and $\dot{\theta} = 0$ ($u=0$).
Design a control law $$ u = -k_{1} (\theta - \pi) - k_{2} \dot{\theta} $$ such that the closed-loop linear system is asymptotically stable, with a time constant equal to $10$ sec.
Simulate this control law on the nonlinear systems when $\theta(0) = 0.9 \pi$ and $\dot{\theta}(0) = 0$.
Let $\Delta \theta = \theta - \pi$, $\omega = \dot{\theta}$, $\Delta \omega = \omega$, $\Delta u = u$.
We notice that $$ \begin{split} \sin \theta &= \sin (\pi + \Delta \theta) \\ &= -\sin \Delta \theta \\ &\approx -\Delta \theta \end{split} $$
The system dynamics can be approximated around $(\theta,\omega) = (\pi , 0)$ by
and $$ m \ell^2 (d/dt) \Delta \omega + b \Delta \omega - mg \ell \Delta \theta = \Delta u. $$
or in standard form
m = 1.0
l = 1.0
b = 0.1
g = 9.81
A = array([[ 0, 1],
[g/l , - b/(m*l*l)]])
B = array([[ 0],
[1/(m*l*l)]])
t1, t2 = 10.0, 5.0
poles = [-1/t1, -1/t2]
K = place_poles(A, B, poles).gain_matrix
def fun(t, theta_omega):
theta, omega = theta_omega
Δtheta, Δomega = theta - pi, omega
Δu = - K @ [Δtheta, Δomega]
u = Δu[0] # Δu has a (1,) shape
dtheta = omega
domega = - (g/l)*sin(theta) - b/(m*l*l)*omega \
+ 1.0/(m*l*l)*u
return array([dtheta, domega])
t_span = [0.0, 30.0]
y0 = [0.9*pi, 0.0]
r = solve_ivp(fun, t_span, y0, dense_output=True)
t = linspace(t_span[0], t_span[-1], 1000)
thetat, omega_t = r["sol"](t)
figure()
plot(t, thetat, label=r"$\theta(t)$")
xlabel("$t$")
yticks([0.9*pi, pi, 1.1*pi],
[r"$0.9\pi$", r"$\pi$", r"$1.1\pi$"])
grid(True); legend()
Consider the dynamics:
Numerical values: $$ m_1 = m_2 = 1, \; k_1 = 1, k_2 = 100, \; b_1 = 2, \; b_2 = 20 $$
Compute the poles of the system.
Is the origin asymptotically stable?
Use a linear feedback to:
kill the oscillatory behavior of the solutions,
“speed up” the dynamics.
Let $v_1 = \dot{x}_1$, $v_2 = \dot{x}_2$. With the state $(x_1, v_1, x_2, v_2)$:
m1 = m2 = 1
k1 = 1; k2 = 100
b1 = 2; b2 = 20
A = array([
[ 0, 1, 0, 0],
[-(k1+k2)/m1, -b1/m1, k2/m1, 0],
[ 0, 0, 0, 1],
[ k2/m2, 0, -k2/m2, -b2/m2]
])
B = array([[0.0], [0.0], [0.0], [1/m2]])
figure()
x = [real(s) for s in eigenvalues]
y = [imag(s) for s in eigenvalues]
plot(x, y, "kx")
plot(0.0, 0.0, "k.")
gca().set_aspect(1.0)
title("Spectrum of $A$"); grid(True)
y0 = [0.0, 0.0, 1.0, 0.0]
t = linspace(0.0, 20.0, 1000)
yt = array([expm(A * t_) for t_ in t]) @ y0
x1t, x2t = yt[:, 0], yt[:, 2]
figure()
plot(t, x1t, label="$x_1$")
plot(t, x2t, label="$x_2$")
xlabel("$t$")
grid(True); legend()
eigenvalues[3] = - 1 / 0.1
K = place_poles(A, B, eigenvalues).gain_matrix
print(repr(eig(A - B @ K)[0]))
figure()
x = [real(s) for s in eigenvalues]
y = [imag(s) for s in eigenvalues]
plot(x, y, "kx")
plot(0.0, 0.0, "k.")
gca().set_aspect(1.0)
title("Spectrum of $A - B K$"); grid(True)
y0 = [0.0, 0.0, 1.0, 0.0]
t = linspace(0.0, 20.0, 1000)
yt = array([expm((A-B@K) * t_) for t_ in t]) @ y0
x1t, x2t = yt[:, 0], yt[:, 2]
figure()
plot(t, x1t, label="$x_1$")
plot(t, x2t, label="$x_2$")
xlabel("$t$")
grid(True); legend()
eigenvalues[0] = - 1 / 0.09
eigenvalues[1] = - 1 / 0.08
K = place_poles(A, B, eigenvalues).gain_matrix
print(repr(eig(A - B @ K)[0]))
figure()
x = [real(s) for s in eigenvalues]
y = [imag(s) for s in eigenvalues]
plot(x, y, "kx")
plot(0.0, 0.0, "k.")
ylim(-12, 12)
gca().set_aspect(1.0)
title("Spectrum of $A - B K$"); grid(True)
y0 = [0.0, 0.0, 1.0, 0.0]
t = linspace(0.0, 20.0, 1000)
yt = array([expm((A-B@K) * t_) for t_ in t]) @ y0
x1t, x2t = yt[:, 0], yt[:, 2]
figure()
plot(t, x1t, label="$x_1$")
plot(t, x2t, label="$x_2$")
xlabel("$t$")
grid(True); legend()