Models

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 matplotlib.pyplot import *

🏷️ Ordinary Differential Equation (ODE)

The “simple” version:

$\dot{x} = f (x)$

where:

State: $x \in R^{n}$
State space: $R^{n}$
Vector field: $f : R^{n} \to R^{n}$ .

More general versions:

Time-dependent vector-field: $\dot{x} = f (t, x), t \in I \subset R,$
$x \in X$ , open subset of $R^{n}$ ,
$x \in X$ , $n$ -dimensional manifold.

🏷️ Vector Field

Visualize $f (x)$ as an arrow with origin the point $x$ .
Visualize $f$ as a field of such arrows.
In the plane ( $n = 2$ ), use quiver from Matplotlib.

🐍 Helper

We define a Q function helper whose arguments are

f: the vector field (a function)
xs, ys: the coordinates (two 1d arrays)

and which returns:

the tuple of arguments expected by quiver.

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)

🔍 Rotation Vector Field

Consider $f (x, y) = (- y, x) .$

def f(xy):
    x, y = xy
    return array([-y, x])

📈 Vector Field

figure()
x = y = linspace(-1.0, 1.0, 20)
ticks = [-1.0, 0.0, 1.0]
xticks(ticks); yticks(ticks)
gca().set_aspect(1.0)
quiver(*Q(f, x, y))

🏷️ ODE Solution

A solution of $\dot{x} = f (x)$ is

a (continuously) differentiable function $x : I \to R^{n},$
defined on a (possibly unbounded) interval $I$ of $R$ ,
such that for every $t \in I,$

$\dot{x} (t) = d x (t) / d t = f (x (t)) .$

📈 Stream Plot

When $n = 2$ , represent a diverse set of solutions in the state space with streamplot

figure()
x = y = linspace(-1.0, 1.0, 20)
gca().set_aspect(1.0)
streamplot(*Q(f, x, y), color="k")

🏷️ Initial Value Problem (IVP)

Solutions $x (t)$ , for $t \geq t_{0}$ , of

$\dot{x} = f (x)$

such that

$x (t_{0}) = x_{0} \in R^{n} .$

🏷️

The initial condition $(t_{0}, x_{0})$ is made of

the initial time $t_{0} \in R$ and
the initial value or initial state $x_{0} \in R^{n}$ .

The point $x (t)$ is the state at time $t$ .

🏷️ Higher-Order ODEs

(Scalar) differential equations whose structure is

$y^{(n)} (t) = g (y, \dot{y}, \ddot{y}, \dots, y^{(n - 1)})$

where $n > 1$ .

💎 Higher-Order ODEs

The previous $n$ -th order ODE is equivalent to the first-order ODE

$\dot{x} = f (x), x \in R^{n}$

with

$f (y_{0}, \dots, y_{n - 2}, y_{n - 1}) := (y_{1}, \dots, y_{n - 1}, g (y_{0}, \dots, y_{n - 1})) .$

🗝️

The result is more obvious if we expand the first-order equation:

$\begin{array}{ccl} {\dot{y}}_{0} & = & y_{1} \\ {\dot{y}}_{1} & = & y_{2} \\ ⋮ & ⋮ & ⋮ \\ {\dot{y}}_{n} & = & g (y_{0}, y_{1}, \dots, y_{n - 1}) \end{array}$

🧩 Pendulum

1. 🧠 🧮

Establish the equations governing the pendulum dynamics.

2. 🧠 🧮

Generalize the dynamics when there is a friction torque $c = - b \dot{θ}$ for some $b \geq 0$ .

We denote $ω$ the pendulum angular velocity:

$ω := \dot{θ} .$

3. 🧠 🧮

Transform the dynamics into a first-order ODE with state $x = (θ, ω)$ .

4. 📈

Draw the system stream plot when $m = 1$ , $ℓ = 1$ , $g = 9.81$ and $b = 0$ .

5. 🧠 🧮

Determine least possible angular velocity $ω_{0} > 0$ such that when $θ (0) = 0$ and $\dot{θ} (0) = ω_{0}$ , the pendulum reaches (or overshoots) $θ (t) = π$ for some $t > 0$ .

🔓 Pendulum

1. 🔓

The pendulum total mechanical energy $E$ is the sum of its kinetic energy $K$ and its potential energy $V$ :

$E = K + V .$

The kinetic energy depends on the mass velocity $v$ :

$K = \frac{1}{2} m v^{2} = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2}$

The potential energy mass depends on the pendulum elevation $y$ . If we set the reference $y = 0$ when the pendulum is horizontal, we have

$V = m g y = - m g ℓ \cos θ$

$\Rightarrow E = K + V = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} - m g ℓ \cos θ .$

If the system evolves without any energy dissipation,

$\begin{aligned} \dot{E} & = \frac{d}{d t} (\frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} - m g ℓ \cos θ) \\ = m ℓ^{2} \dot{θ} \ddot{θ} + m g ℓ (\sin θ) \dot{θ} \\ = 0 \end{aligned}$

$\Rightarrow m ℓ^{2} \ddot{θ} + m g ℓ \sin θ = 0.$

2. 🔓

When there is an additional dissipative torque $c = - b θ$ , we have instead

$\dot{E} = c \dot{θ} = - b {\dot{θ}}^{2}$

and thus

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

3. 🔓

With $ω := \dot{θ}$ , the dynamics becomes

$\begin{array}{lll} \dot{θ} & = & ω \\ \dot{ω} & = & - (b / m ℓ^{2}) ω - (g / ℓ) \sin θ \end{array}$

4. 🔓

m=1.0; b=0.0; l=1.0; g=9.81
def f(theta_d_theta):
    theta, d_theta = theta_d_theta
    J = m * l * l
    d2_theta  = - b / J * d_theta
    d2_theta += - g / l * sin(theta)
    return array([d_theta, d2_theta])

📈

figure()
theta = linspace(-1.5 * pi, 1.5 * pi, 100)
d_theta = linspace(-5.0, 5.0, 100)
labels =  [r"$-\pi$", "$0$", r"$\pi$"]
xticks([-pi, 0, pi], labels)
yticks([-5, 0, 5])
streamplot(*Q(f, theta, d_theta), color="k")

5. 🔓

In the top vertical configuration, the total mechanical energy of the pendulum is

$E_{⊤} = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} - m g ℓ \cos π = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} + m g ℓ .$

Hence we have at least $E_{⊤} \geq m g ℓ$ .

On the other hand, in the bottom configuration,

$E_{⊥} = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} - m g ℓ \cos 0 = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} - m g ℓ .$

Hence, without any loss of energy, the initial velocity must satisfy $E_{⊥} \geq E_{⊤}$ for the mass to reach the top position.

That is

$E_{⊥} = \frac{1}{2} m ℓ^{2} {\dot{θ}}^{2} - m g ℓ \geq m g ℓ = E_{⊤}$

which leads to:

$| \dot{θ} | \geq 2 \sqrt{\frac{g}{ℓ}} .$