% Observers % 👤 Sébastien Boisgérault
🐍 | Code | 🔍 | Worked Example |
📈 | Graph | 🧩 | Exercise |
🏷️ | Definition | 💻 | Numerical Method |
💎 | Theorem | 🧮 | Analytical Method |
📝 | Remark | 🧠 | Theory |
ℹ️ | Information | 🗝️ | Hint |
Warning | 🔓 | Solution |
from numpy import *
from numpy.linalg import *
from scipy.integrate import solve_ivp
from scipy.linalg import solve_continuous_are
from matplotlib.pyplot import *
from numpy.testing 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)
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
def Q(f, xs, ys):
X, Y = meshgrid(xs, ys)
fx = vectorize(lambda x, y: f([x, y])[0])
fy = vectorize(lambda x, y: f([x, y])[1])
return X, Y, fx(X, Y), fy(X, Y)
Controling a system generally requires the knowledge of the state
Can we reduce the amount of physical sensors and still be able to compute the state with "virtual" or "software" sensors ?
Control engineers call these software devices observers.
First we adress their mathematical feasibility.
The system
is observable if the knowledge of
-
The knowledge of
$x(0)$ determines uniquely$x(t)$ via the system dynamics. -
Later, observers will provide merely asymptotically exact estimates
$\hat{x}(t)$ of$x(t)$ , that satisfy$\hat{x}(t) - x(t) \to 0$ when$t \to +\infty.$
The definition of observability may be extended to systems with (known) inputs
In general, the input
But for linear systems, the choice of
Indeed, if
and we can deduce
then in the general case, when we measure
we can compute
and deduce
The position
where
-
we don't know where the car is at
$t=0$ , -
we don't know what its initial speed is,
-
we do know that the car doesn't accelerate (
$u=0$ ).
If we measure the position
-
$x(0) = y(0)$ is known, -
$\dot{x}(0) = \dot{y}(0)$ is also computable.
Thus the system is observable.
What if we measure the speed instead of the location ?
The system dynamics
We can't deduce the position of the car from the measure of its speed; the system is not observable.
The system
-
"$,$" row concatenation of matrices.
-
"$;$" column concatenation of matrices.
We have
The system
The system
Consider
with
Is the system observable ?
Yes! The rank of its observability matrix
$$
\left[
\begin{array}{c}
C \
CA \
\vdots \
C A^{n-1}
\end{array}
\right]
$$
is at most
$$\dot{x}n = 0, , \dot{x}{n-1} = x_n, , \cdots ,, \dot{x}_1 = x_2, , y=x_1$$
Show that the system is observable.
The standard form of the dynamics associated to the state
Thus,
$$
\begin{array}{rl}
C &=& \left[1, 0, 0, \dots, 0 \right] \
CA &=& \left[0, 1,0, \dots, 0 \right] \
\vdots &=& \vdots \
CA^{n-1} &=& \left[0, 0, 0, \dots, 1 \right]
\end{array}
$$
The observability matrix has rank
-
$d T_1/dt = 0 + (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)$ -
$y = T_4$
Show that the system is observable.
Is it still true if the four cells are organized as a square and the temperature sensor is in any of the corners ?
Can you make the system observable with two (adequatly located) sensors?
The standard form of the dynamics associated to the state
Therefore the observability matrix is $$ \left[ \begin{array}{rrrr} 0 & 0 & 0 & 1 \ 0 & 0 & 1 & -1 \ 0 & 1 & -3 & 2 \ 1 & -5 & 9 & 5 \end{array} \right] $$ whose rank is 4: the system is observable.
In the square configuration, there is no way to get observability with a single thermometer.
Indeed the new
\right]
$$
and the new
The corresponding observability matrices are
All the possible observability matrices in this case have rank 2 or 3 < 4.
With 2 sensors, "it depends" (on the location of the sensors). For example:
The first case corresponds to
the second one to
The observability matrices are
The first has rank 3, the second rank 4.
Simulate the system behavior
and since we don't know better,
Is
The dynamics of the state estimate error
The state estimator error
doesn't satisfy in general
We have
for every value of
Change the observer dynamics to account for differences between
for some observer gain matrix
The new dynamics of
The system
The system
In this case, we can perform arbitrary pole assignment:
-
for any conjugate set
$\Lambda$ of eigenvalues, -
there is a matrix
$K \in \mathbb{R}^{p \times n}$ such that$$ \sigma(A^t - C^t K) = \Lambda $$
Since
Thus, if we set
we have solved the pole assignment problem for observers:
Consider the double integrator
\left[\begin{array}{cx} 0 & 1 \ 0 & 0\end{array}\right] \left[\begin{array}{c} x_1 \ x_2 \end{array}\right] + \left[\begin{array}{c} 0 \ 1 \end{array}\right] u $$
(in standard form)
from scipy.signal import place_poles
A = array([[0, 1], [0, 0]])
C = array([[1, 0]])
poles = [-1, -2]
K = place_poles(A.T, C.T, poles).gain_matrix
L = K.T
assert_almost_equal(K, [[3.0, 2.0]])
\left[\begin{array}{cx} 0 & 1 \ 0 & 0\end{array}\right] \left[\begin{array}{c} \hat{x}_1 \ \hat{x}_2 \end{array}\right] + \left[\begin{array}{c} 0 \ 1 \end{array}\right] u
- \left[\begin{array}{c} 3 \2 \end{array}\right] (\hat{y} - y) $$
def fun(t, X_Xhat):
x, x_hat = X_Xhat[0:2], X_Xhat[2:4]
y, y_hat = C.dot(x), C.dot(x_hat)
dx = A.dot(x)
dx_hat = A.dot(x_hat) - L.dot(y_hat - y)
return r_[dx, dx_hat]
y0 = [-2.0, 1.0, 0.0, 0.0]
result = solve_ivp(
fun=fun,
t_span=[0.0, 5.0],
y0=y0,
max_step=0.1
)
figure()
t = result["t"]
y = result["y"]
plot(t, y[0], "C0", label="$x_1$")
plot(t, y[2], "C0--", label=r"$\hat{x}_1$")
plot(t, y[1], "C1", label="$x_2$")
plot(t, y[3], "C1--", label=r"$\hat{x}_2$")
xlabel("$t$"); grid(); legend()
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
pp.gcf().subplots_adjust(bottom=0.2)
save("images/observer-trajectories")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
Consider
-
the state
$x(t)$ is unknown ($x(0)$ is unknown), -
only (a noisy version of)
$y(t)$ is available.
We want a sensible estimation
We now assume the existence of state and output disturbances
Thes disturbances (or "noises") are unknown; we are searching for the
estimate
For a known
where:
-
$Q \in \mathbb{R}^{n \times n}$ and$R \in \mathbb{R}^{p\times p}$ , -
(to be continued ...)
-
$Q$ and$R$ are symmetric ($R^t = R$ and$Q^t = Q$ ), -
$Q$ and$R$ are positive definite (denoted "$>0$")
If it is known that there is a large state disturbance but small output
disturbance, it makes sense to reduce the impact of the state disturbance
in the composition of
Assume that
There is a state estimation
The dynamics of the corresponding estimation error
The gain matrix
where
Solve the Ricatti equation for optimal control with
then define
Consider the system $$ \begin{split} \dot{x} &= v\ y &= x + w \end{split} $$
If we believe that the state and output perturbation are of the same scale, we may try
With
whose only positive solution is
With
Thus, the optimal filter is
Consider the double integrator
\left[\begin{array}{cc} 0 & 1 \ 0 & 0\end{array}\right] \left[\begin{array}{c} x \ \dot{x} \end{array}\right] + \left[\begin{array}{c} v_1 \ v_2 \end{array}\right], y = \left[\begin{array}{c} 1 & 0 \end{array} \right] \left[\begin{array}{c} x \ \dot{x} \end{array}\right]
- w $$
(in standard form)
A = array([[0, 1], [0, 0]])
B = array([[0], [1]])
Q = array([[1, 0], [0, 1]])
R = array([[1]])
sca = solve_continuous_are
Sigma = sca(A.T, C.T, inv(Q), inv(R))
L = Sigma @ C.T @ R
eigenvalues, _ = eig(A - L @ C)
assert all([real(s) < 0 for s in eigenvalues])
figure()
x = [real(s) for s in eigenvalues]
y = [imag(s) for s in eigenvalues]
plot(x, y, "kx"); xlim(-2, 2); ylim(-2, 2)
plot([0, 0], [-2, 2], "k");
plot([-2, 2], [0, 0], "k")
grid(True); axis("square")
title("Eigenvalues")
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
axis("square")
axis([-2, 2, -2, 2])
save("images/poles-Kalman")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
🐍
def fun(t, X_Xhat):
x, x_hat = X_Xhat[0:2], X_Xhat[2:4]
y, y_hat = C.dot(x), C.dot(x_hat)
dx = A.dot(x)
dx_hat = A.dot(x_hat) - L.dot(y_hat - y)
return r_[dx, dx_hat]
y0 = [-2.0, 1.0, 0.0, 0.0]
result = solve_ivp(
fun=fun,
t_span=[0.0, 5.0],
y0=y0,
max_step=0.1
)
figure()
t = result["t"]; y = result["y"]
plot(t, y[0], "C0", label="$x_1$")
plot(t, y[2], "C0--", label=r"$\hat{x}_1$")
plot(t, y[1], "C1", label="$x_2$")
plot(t, y[3], "C1--", label=r"$\hat{x}_2$")
xlabel("$t$")
grid(); legend()
::: hidden :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
pp.gcf().subplots_adjust(bottom=0.2)
save("images/observer-Kalman-trajectories")
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::: slides :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
<style> .reveal p { text-align: left; } .reveal section img { text-align: center; border:0; height:50vh; width:auto; display:block; margin: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; } section p.author { text-align: center; margin: auto; } </style>