Asymptotic Stabilization

Control Engineering with Python

Symbols

🐍	Code	🔍	Worked Example
📈	Graph	🧩	Exercise
🏷️	Definition	💻	Numerical Method
💎	Theorem	🧮	Analytical Method
📝	Remark	🧠	Theory
ℹ️	Information	🗝️	Hint
⚠️	Warning	🔓	Solution

🐍 Imports

from numpy import *
from numpy.linalg import *
from numpy.testing import *
from scipy.integrate import *
from scipy.linalg import *
from matplotlib.pyplot import *

🧭 Asymptotic Stabilization

When the system

$\dot{x} = A x, x \in R^{n}$

is not asymptotically stable at the origin, maybe there are some inputs $u \in R^{m}$ such that

$\dot{x} = A x + B u$

that we can use to stabilize asymptotically the system?

🏷️ Linear Feedback

We search for $u$ as a linear feedback:

$u (t) = - K x (t)$

for some $K \in R^{m \times n}$ .

📝 Note

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.

🏷️ Closed-Loop Diagram

💎 Closed-Loop Dynamics

When

$\begin{array}{rcl} \dot{x} & = & A x + B u \\ u & = & - K x \end{array}$

the state $x \in R^{n}$ evolves according to:

$\dot{x} = (A - B K) x$

💎 Reminder

The closed-loop system is asymptotically stable iff every eigenvalue of the matrix

$A - B K$

is in the open left-hand plane.

🏷️ Spectrum as a Multiset

Multisets remember the multiplicity of their elements. It’s convenient to describe the spectrum of matrices:

$A := [\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2 \end{array}] \Rightarrow σ (A) = {1, 1, 2}$

$0 \notin σ (A), 1 \in σ (A), 1 \in^{2} σ (A), 1 \notin^{3} σ (A)$

$2 \in σ (A), 2 \notin^{2} σ (A)$

💎 Pole Assignment

Assumptions

The system

$\dot{x} = A x + B u, x \in R^{n}, u \in R^{p}$

is controllable.
$Λ$ is a symmetric multiset of $n$ complex numbers:

$Λ = {λ_{1}, \dots, λ_{n}} \subset C and λ \in^{k} Λ \Rightarrow \overset{―}{λ} \in^{k} Λ .$

💎 Pole Assignment

Conclusion

$\Rightarrow$ There is a matrix $K \in R^{n \times m}$ such that

$σ (A - B K) = Λ .$

🔍 Pole Assignment

Consider the double integrator $\ddot{x} = u$

$\frac{d}{d t} [\begin{matrix} x \\ \dot{x} \end{matrix}] = [\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}] [\begin{matrix} x \\ \dot{x} \end{matrix}] + [\begin{matrix} 0 \\ 1 \end{matrix}] 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])

🐍 📊 Spectrum

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)

⚠️ Limitations

⛔ The place_poles function rejects eigenvalues whose multiplicity is higher than the rank of $B$ .

In the previous example, $rank B = 1$ , so

❌ poles = [-1, -1] won’t work.
✔️ poles = [-1, -2] will.

🧩 Pole Assignment

Consider the system with dynamics

$\begin{array}{ccr} {\dot{x}}_{1} & = & x_{1} - x_{2} + u \\ {\dot{x}}_{2} & = & - x_{1} + x_{2} + u \end{array}$

We apply the control law

$u = - k_{1} x_{1} - k_{2} x_{2} .$

1. 🧮

Can we assign the poles of the closed-loop system freely by a suitable choice of $k_{1}$ and $k_{2}$ ?

2. 🧠

Explain this result.

🔓 Pole Assignment

1. 🔓

$\begin{aligned} A - B K & = [\begin{array}{rr} 1 & - 1 \\ - 1 & 1 \end{array}] - [\begin{array}{r} 1 \\ 1 \end{array}] [\begin{array}{rr} k_{1} & k_{2} \end{array}] \\ = [\begin{array}{rr} 1 - k_{1} & - 1 - k_{2} \\ - 1 - k_{1} & 1 - k_{2} \end{array}] \end{aligned}$

$\begin{aligned} det A - B K & = det ([\begin{array}{rr} s - 1 + k_{1} & 1 + k_{2} \\ 1 + k_{1} & s - 1 + k_{2} \end{array}]) \\ = (s - 1 + k_{1}) (s - 1 + k_{2}) - (1 + k_{1}) (1 + k_{2}) \\ = s^{2} + (k_{1} + k_{2}) s - 2 (k_{1} + k_{2}) \end{aligned}$

$\begin{aligned} σ (A - B K) & = {λ_{1}, λ_{2}} \\ = {λ \in C | s^{2} + (k_{1} + k_{2}) s - 2 (k_{1} + k_{2}) = 0} . \end{aligned}$

Since the characteristic polynomial is also $(s - λ_{1}) (s - λ_{2})$ we get $k_{1} + k_{2} = - λ_{1} - λ_{2}, - 2 (k_{1} + k_{2}) = λ_{1} λ_{2}$

Thus we have

$λ_{1} λ_{2} = 2 (λ_{1} + λ_{2}) \Rightarrow λ_{2} = \frac{2 λ_{1}}{λ_{1} - 2}$

and both poles cannot be assigned freely; for example if we select $λ_{1} = 1$ , we end up with $λ_{2} = - 2$ .

2. 🔓

We have not checked the assumptions of 💎 Pole Assignment yet.

The commandability matrix is $[B, A B] = [\begin{array}{rr} 1 & 0 \\ 1 & 0 \end{array}]$ whose rank is $1 < 2$ .

Since the system is not controllable, pole assignment may fail and it does here.

🧩 Pendulum

Consider the pendulum with dynamics:

$m ℓ^{2} \ddot{θ} + b \dot{θ} + m g ℓ \sin θ = u$

Numerical Values:

$m = 1.0, ℓ = 1.0, b = 0.1, g = 9.81$

1. 🧮

Compute the linearized dynamics of the system around the equilibrium $θ = π$ and $\dot{θ} = 0$ ( $u = 0$ ).

2. 💻

Design a control law $u = - k_{1} (θ - π) - k_{2} \dot{θ}$ such that the closed-loop linear system is asymptotically stable, with a time constant equal to $10$ sec.

3. 💻 🧮

Simulate this control law on the nonlinear systems when $θ (0) = 0.9 π$ and $\dot{θ} (0) = 0$ .

🔓 Pendulum

1. 🔓

Let $Δ θ = θ - π$ , $ω = \dot{θ}$ , $Δ ω = ω$ , $Δ u = u$ .

We notice that $\begin{aligned} \sin θ & = \sin (π + Δ θ) \\ = - \sin Δ θ \\ \approx - Δ θ \end{aligned}$

The system dynamics can be approximated around $(θ, ω) = (π, 0)$ by

$(d / d t) Δ θ = Δ ω$ and $m ℓ^{2} (d / d t) Δ ω + b Δ ω - m g ℓ Δ θ = Δ u .$

or in standard form

$\frac{d}{d t} [\begin{matrix} Δ θ \\ Δ ω \end{matrix}] = [\begin{array}{cc} 0 & 1 \\ g / ℓ & - b / (m ℓ^{2}) \end{array}] [\begin{matrix} Δ θ \\ Δ ω \end{matrix}] + [\begin{matrix} 0 \\ 1 / (m ℓ^{2}) \end{matrix}] Δ u$

2. 🔓

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

3. 🔓

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()

🧩 Double Spring

Consider the dynamics:

$\begin{array}{rcl} m_{1} {\ddot{x}}_{1} & = & - k_{1} x_{1} - k_{2} (x_{1} - x_{2}) - b_{1} {\dot{x}}_{1} \\ m_{2} {\ddot{x}}_{2} & = & - k_{2} (x_{2} - x_{1}) - b_{2} {\dot{x}}_{2} + u \end{array}$

Numerical values: $m_{1} = m_{2} = 1, k_{1} = 1, k_{2} = 100, b_{1} = 2, b_{2} = 20$

1. 🧮 💻

Compute the poles of the system.

Is the origin asymptotically stable?

2. 💻

Use a linear feedback to:

kill the oscillatory behavior of the solutions,
“speed up” the dynamics.

🔓 Double Spring System

1. 🔓

Let $v_{1} = {\dot{x}}_{1}$ , $v_{2} = {\dot{x}}_{2}$ . With the state $(x_{1}, v_{1}, x_{2}, v_{2})$ :

$A = [\begin{array}{cccc} 0 & 1 & 0 & 0 \\ - (k_{1} + k_{2}) / m_{1} & - b_{1} / m_{1} & k_{2} / m_{1} & 0 \\ 0 & 0 & 0 & 1 \\ k_{2} / m_{2} & 0 & - k_{2} / m_{2} & - b_{2} / m_{2} \end{array}]$

$B = [\begin{matrix} 0 \\ 0 \\ 0 \\ 1 / m_{2} \end{matrix}]$

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]])

eigenvalues, _ = eig(A)

>>> eigenvalues
array([-15.64029062 +0.j        ,  
        -3.15722141+11.45767938j,
        -3.15722141-11.45767938j,  
        -0.04526657 +0.j        ])

Since all eigenvalues have a negative real part, the double-spring system is asymptotically stable.

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()

2. 🔓

eigenvalues[3] = - 1 / 0.1
K = place_poles(A, B, eigenvalues).gain_matrix
print(repr(eig(A - B @ K)[0]))

eigenvalues, _ = eig(A - B @ K)

>>> eigenvalues
array([-15.64029062 +0.j        ,  
        -3.15722141+11.45767938j,
        -3.15722141-11.45767938j,  
        -1.         +0.j        ])

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]))

eigenvalues, _ = eig(A - B @ K)

>>> eigenvalues
array([-15.64029062+0.j, 
       -12.5       +0.j, 
       -11.11111111+0.j,
       -10.        +0.j])

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()