% Controllability % ๐ค Sรฉbastien Boisgรฉrault
-
๐ Documents (GitHub)
-
ยฉ๏ธ License CC BY 4.0
๐ | Code | ๐ | Worked Example |
๐ | Graph | ๐งฉ | Exercise |
๐ท๏ธ | Definition | ๐ป | Numerical Method |
๐ | Theorem | ๐งฎ | Analytical Method |
๐ | Remark | ๐ง | Theory |
โน๏ธ | Information | ๐๏ธ | Hint |
Warning | ๐ | Solution |
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
from numpy import *
from numpy.linalg import *
from numpy.testing import *
from scipy.integrate import *
from scipy.linalg import *
from matplotlib.pyplot import *
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: notebook :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
from numpy import *
import matplotlib; matplotlib.use("nbAgg")
%matplotlib notebook
from matplotlib.pyplot import *
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Python 3.x Standard Library
import gc
import os
# Third-Party Packages
import numpy as np; np.seterr(all="ignore")
import numpy.linalg as la
import scipy.misc
import matplotlib as mpl; mpl.use("Agg")
import matplotlib.pyplot as pp
import matplotlib.axes as ax
import matplotlib.patches as pa
#
# Matplotlib Configuration & Helper Functions
# --------------------------------------------------------------------------
# TODO: also reconsider line width and markersize stuff "for the web
# settings".
fontsize = 10
width = 345 / 72.27
height = width / (16/9)
rc = {
"text.usetex": True,
"pgf.preamble": r"\usepackage{amsmath,amsfonts,amssymb}",
#"font.family": "serif",
"font.serif": [],
#"font.sans-serif": [],
"legend.fontsize": fontsize,
"axes.titlesize": fontsize,
"axes.labelsize": fontsize,
"xtick.labelsize": fontsize,
"ytick.labelsize": fontsize,
"figure.max_open_warning": 100,
#"savefig.dpi": 300,
#"figure.dpi": 300,
"figure.figsize": [width, height],
"lines.linewidth": 1.0,
}
mpl.rcParams.update(rc)
# Web target: 160 / 9 inches (that's ~45 cm, this is huge) at 90 dpi
# (the "standard" dpi for Web computations) gives 1600 px.
width_in = 160 / 9
def save(name, **options):
cwd = os.getcwd()
root = os.path.dirname(os.path.realpath(__file__))
os.chdir(root)
pp.savefig(name + ".svg", **options)
os.chdir(cwd)
def set_ratio(ratio=1.0, bottom=0.1, top=0.1, left=0.1, right=0.1):
height_in = (1.0 - left - right)/(1.0 - bottom - top) * width_in / ratio
pp.gcf().set_size_inches((width_in, height_in))
pp.gcf().subplots_adjust(bottom=bottom, top=1.0-top, left=left, right=1.0-right)
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
The system
-
for any
$t_0 \in \mathbb{R}$ ,$x_0 \in \mathbb{R}^n$ and$x_f \in \mathbb{R}^n$ , -
there are
$t_f > t_0$ and$u: [t_0, t_f] \to \mathbb{R}^m$ such that -
the solution
$x(t)$ such that$x(t_0)=x_0$ satisfies$$ x(t_f) = x_f. $$
Let
It is admissible if there is a function
The position
where
The car is initially at the origin of a road and motionless.
We would like to cross the end of the road (location
Numerical values:
-
$m=1500 , \mbox{kg}$ , -
$t_f=10 , \mbox{s}$ ,$d_f=100 , \mbox{m}$ and$v_f=100 , \mbox{km/h}$ .
We search for a reference trajectory for the state
such that:
-
$d_r(0)=0$ ,$\dot{d}_r(0) = 0$ , -
$d_r(t_f) = x_f$ ,$\dot{d}_r(t_f) = v_f$ .
We check that this reference trajectory is admissible,
i.e. that we can find a control
Here, if
Thus,
We can find
with
(equivalently, with
m = 1500.0
xf = 100.0
vf = 100.0 * 1000 / 3600 # m/s
tf = 10.0
alpha = vf/tf**2 - 2*xf/tf**3
beta = 3*xf/tf**2 - vf/tf
def x(t):
return alpha * t**3 + beta * t**2
def d2_x(t):
return 6 * alpha * t + 2 * beta
def u(t):
return m * d2_x(t)
y0 = [0.0, 0.0]
def fun(t, y):
x, d_x = y
d2_x = u(t) / m
return [d_x, d2_x]
result = solve_ivp(
fun, [0.0, tf], y0, dense_output=True
)
figure()
t = linspace(0, tf, 1000)
xt = result["sol"](t)[0]
plot(t, xt)
grid(True); xlabel("$t$"); title("$d(t)$")
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/car")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
figure()
vt = result["sol"](t)[1]
plot(t, 3.6 * vt)
grid(True); xlabel("$t$")
title(r"$\dot{d}(t)$ km/h")
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/car-speed")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Let
Find a smooth reference trajectory
The first line of the vector equation
Any trajectory
Consider the pendulum with dynamics:
Find a smooth reference trajectory that leads the pendulum from the bottom configuration
to the top configuration
Show that the reference trajectory is admissible and compute the
corresponding input
Simulate the result and visualize the solution.
What should theoretically happen at
Numerical Values:
Since the initial and final conditions form a set of 4 equations, we search for an admissible trajectory defined by a 3rd-order polynomial
with unknown coefficients
The initial conditions
The final condition
and
The latter equation can be transformed into
and the former becomes
and thus
Finally, the following trajectory
By construction our reference trajectory meets the initial condition. There is a unique input that makes the solution follows this reference; its is defined by
where
m = 1.0
l = 1.0
b = 0.1
g = 9.81
t0, tf = 0.0, 10.0
def theta_r(t):
s = t/tf
return -2*pi*s**3 + 3*pi*s**2
def dtheta_r(t):
s = t/tf
return -6*pi/tf*s**2 + 6*pi/tf*s
def d2theta_r(t):
s = t/tf
return -12*pi/(tf*tf)*s + 6*pi/(tf*tf)
def u_r(t):
return (m*l*l*d2theta_r(t) +
b*dtheta_r(t) +
m*g*l*sin(theta_r(t)))
def fun(t, theta_dtheta):
theta, dtheta = theta_dtheta
d2theta = ((-b*dtheta - m*g*l*sin(theta) + u_r(t))
/ (m * l * l))
return dtheta, d2theta
t_span = [t0, tf]
y0 = [0.0, 0.0]
r = solve_ivp(fun, t_span, y0, dense_output=True)
t = linspace(t0, tf, 1000)
solt = r["sol"](t)
thetat, dthetat = solt
figure()
plot(t, theta_r(t), "k--", label=r"$\theta_r(t)$")
plot(t, thetat, label=r"$\theta(t)$")
yticks([0, 0.5*pi, pi], ["$0$", r"$\pi/2$", r"$\pi$"])
title(r"Simulation of $\theta(t)$")
grid(True); legend()
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/pendulum-miss")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Theoretically, we should have
Since the top equilibrium is unstable, small (numerical) errors in the state may cause large deviations of the trajectory.
We may reduce the simulation error thresholds to alleviate this problem.
rhp = solve_ivp(fun, t_span, y0,
dense_output=True,
rtol=1e-6, atol=1e-9)
t = linspace(t0, tf, 1000)
solt_hp = rhp["sol"](t)
thetat_hp, dthetat_hp = solt_hp
figure()
plot(t, theta_r(t), "k--", label=r"$\theta_r(t)$")
plot(t, thetat, label=r"$\theta(t)$ (standard)")
plot(t, thetat_hp, label=r"$\theta(t)$ (low error)")
yticks([0, 0.5*pi, pi], ["$0$", r"$\pi/2$", r"$\pi$"])
title(r"Simulation of $\theta(t)$")
grid(True); legend()
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
tight_layout()
save("images/pendulum-miss-less")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
A system
-
from the origin
$x_0 = 0$ at$t_0=0$ , -
we can reach any state
$x_f \in \mathbb{R}^n$ .
The system
def KCM(A, B):
n = shape(A)[0]
mp = matrix_power
cs = column_stack
return cs([mp(A, k) @ B for k in range(n)])
n = 3
A = zeros((n, n))
for i in range(0, n-1):
A[i,i+1] = 1.0
B = zeros((n, 1))
B[n-1, 0] = 1.0
C = KCM(A, B)
C_expected = [[0, 0, 1], [0, 1, 0], [1, 0, 0]]
assert_almost_equal(C, C_expected)
assert matrix_rank(C) == n
-
This implementation of
KCM
is not optimized: fewer computations are achievable using the identity:$$ A^n B = A \times (A^{n-1}B). $$
-
Rank computations are subject to (catastrophic) numerical errors; sensitivity analysis or symbolic computations may alleviate the problem.
Consider
-
$x \in \mathbb{R}^n$ , -
$u \in\mathbb{R}^m$ and, -
$\mathrm{rank} , B = n$ .
Show that
Is the system controllable ?
Given
The shape of
By assumption
The system is controllable since
Since
Then by construction
If we define
To join
Indeed, that leads to
$$
\dot{x}(t) = A x(t) - A x_r(t) + \dot{x}_r(t), ; x(t_0) = x_r(t_0)
$$
or if we denote
$$\dot{x}n = u, , \dot{x}{n-1} = x_n, , \cdots ,, \dot{x}_1 = x_2.$$
Show that the system is controllable
We have
Thus,
The controllability matrix has full rank and the system is controllable.
-
$d T_1/dt = u + (T_2 - T_1)$ -
$d T_2/dt = (T_1 - T_2) + (T_3 - T_2)$ -
$d T_3/dt = (T_2 - T_3) + (T_4 - T_3)$ -
$d T_4/dt = (T_3 - T_4)$
Show that the system is controllable.
Assume now that the four cells are organized as a square.
-
Is the system still controllable?
-
Why?
-
How could you solve this problem?
We have
The controllability matrix
is full rank, thus the system is controllable.
If the system is organized as a square
The controllability matrix is
It is not full rank since the second and the third rows are equal. Thus the system is not controllable.
The lack of controllability is due to symmetry.
If the initial temperature in cell 2 and 3 are equal, since they play symmetric role in the system, their temperature will be equal for any subsequent time.
Hence, it will be impossible to reach a state with different temperature in cell 2 and 3, no matter what the input is.
To break this symmetry and restore controllability, we may for example try to add a second independent heat source sink in cell 2.
This leads to
and the controllability matrix
which has a full-rank.
<style> .reveal p { text-align: left; } .reveal section img { border:0; height:50vh; width:auto; } .reveal section img.medium { border:0; max-width:50vh; } .reveal section img.icon { display:inline; border:0; width:1em; margin:0em; box-shadow:none; vertical-align:-10%; } .reveal code { font-family: Inconsolata, monospace; } .reveal pre code { background-color: white; font-size: 1.5em; line-height: 1.5em; /_ max-height: 80wh; won't work, overriden _/ } /_ .reveal .slides .left { text-align: left; } _/ input { font-family: "Source Sans Pro", Helvetica, sans-serif; font-size: 42px; line-height: 54.6px; } code span.kw { color: inherit; font-weight: normal; } code span.cf { /_ return _/ color: inherit; font-weight: normal; } code span.fl { /_ floats _/ color: inherit; } code span.dv { /_ ints _/ color: inherit; } code span.co { /_ comments _/ font-style: normal; color: #adb5bd; /_ gray 5 _/} code span.st { /_ strings _/ color: inherit; } code span.op { /_ +, = _/ color: inherit; } /*** Details ******************************************************************/ details h1, details h2, details h3{ display: inline; } details summary { cursor: pointer; list-style: '๐ '; } details[open] summary { cursor: pointer; list-style: '๐ '; } summary::-webkit-details-marker { display: none } details[open] summary ~ * { animation: sweep .5s ease-in-out; } @keyframes sweep { 0% {opacity: 0} 100% {opacity: 1} } section p.author { text-align: center; margin: auto; } </style>