Exakte Eingangs- Ausgangslinearisierung

Wir wollen hier eine Exakte Eingangs-Ausgangslinearisierung im SISO-Fall mithilfe differentialgeometrischen Methoden durchführen. Jürgen Adamy gibt mit seinem Buch Nichtlineare Systeme und Regelungen eine sehr gute Übersicht über die nichtlineare Regelungstechnik. Wer eine tiefe Einführung in die differentialgeometrischen Methoden der Physik will, kann diese in dem herausragenden Buch The Geometry of Physics finden.

Auf eine theoretische Einführung wollen wir hier erstmal verzichten. Diese soll in einem späteren Blogeintrag nachgereicht werden.

Es soll nur gezeigt werden, dass Python / Sympy auch für nichtlineare Methoden gut geeignet ist.

Eingang/Aungangs - Linearisierung

Ein nichtlineares SISO-System ist in der Form

\[ \begin{split} \dot{x} &= f (x) + g(x)u \\ y &= h(x) \end{split} \]

gegeben.

Durch das Anwenden der Lie Ableitung können wir den relativen Grad ermitteln:

\[ \begin{split} y &= L^0_f h \\ \dot{y} &= L^1_f h \\ \ddot{y} &= L^2_f h \\ \vdots \\ y^{(r)} &= L^r_f h + L_g (L^{r-1}_f h)u. \end{split} \]

In der letzten Gleichung tritt zum ersten Mal der Eingang \(u\) in Erscheinung. Die Anzahl der notwendigen Ableitungen wird als relativer Grad \(r\) bezeichnet. Wir wollen hier einfachheitshalber annehmen, dass der relative Grad \(r\), gleich der Anzahl der Zustände \(x_i\) ist (\(r=n\)).

Oft werden die Abkürzungen \(\alpha(x) = L^r_f h\) und \(\beta(x) = L_g (L^{r-1}_f g)\) für \(y^{(r)} = \alpha(x)+\beta(x) u\) eingeführt.

Wir suchen nun einen Eingang, der die Nichtlinearitäten im System kompensiert, d. h. \(y^{(r)} = \alpha(x) + \beta(x) u \equiv v\). Daraus folgt:

\[ \begin{split} v = \alpha(x) + \beta(x) u \\ u = \frac{1}{\beta(x)}[-\alpha(x) + v] \end{split} \]

Wir können neue Zustände \(z\) einführen:

\[ \begin{align} z_1 &= y = L^0_f h \\ z_2 &= \dot{y} = L^1_f h \\ \vdots \\ z_r &= y^{(r-1)} = L^{r-1}_f h \\ z_{r+1} &= y^{(r)} = \alpha + \beta u = v \end{align} \]

In diesen neuen Zuständen ist das System nun linear.

Das System kann also in Form eines Zustandsraums wie folgt dargestellt werden

\[ \begin{split} \begin{bmatrix} \dot{z}_1 \\ \dot{z}_2 \\ \vdots \\ \dot{z}_r \\ \end{bmatrix} &= \begin{bmatrix} 0 & 1 & \cdots & 0 \\ 0 & 0 & \ddots & 0 \\ 0 & 0 & \cdots & 1 \\ 0 & 0 & \cdots & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \\ \vdots \\ z_r \\ \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 1 \\ \end{bmatrix} \\ y &= \begin{bmatrix} 1 & 0 & \cdots 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \\ \vdots \\ z_r \\ \end{bmatrix} \end{split} \]

und entspricht eine Kaskade von r Integratoren. Jetzt können wir jede Art von Regler für das lineare System mit \(r\) Integratoren entwerfen. Am besten geeignet wäre ein Zustandsregler. Die Transformation zwischen den beiden Zustandsräumen, \(x\) und \(z\) wollen wir mit \(\Phi(x)\) bezeichnen:

\[ z = \Phi(x) = \begin{bmatrix} \varphi_1(x) \\ \varphi_2(x) \\ \vdots \\ \varphi_r(x) \end{bmatrix} = \begin{bmatrix} L^0_f h \\ L^1_f h \\ \vdots \\ L^{r-1}_f h \end{bmatrix} \]

Hinweis

Algorithmus für die Eingangs-/Ausgangslinearisierung:

1. Ermitteln des relativen Grads \(r\), des nichtlinearen Systems

\[ \begin{split} L_g [L^k_f h(x)] &= 0 \\ L_g [L^{r-1}_f h(x)] &\neq 0 \end{split} \]

1. Ermitteln der Zustandstransformation

\[ z = \Phi(x) = \begin{bmatrix} \varphi_1(x) \\ \varphi_2(x) \\ \vdots \\ \varphi_r(x) \end{bmatrix} = \begin{bmatrix} L^0_f h \\ L^1_f h \\ \vdots \\ L^{r-1}_f h \end{bmatrix} \]

1. Definieren des neuen Eingangs, sodass die Zustände linearisiert werden

\[ \begin{split} \alpha(x) &= L^{r}_{f} h \\ \beta(x) &= L_g (L^{(r-1)}_{f} h) \\ u &= \frac{1}{\beta(x)}\left[-\alpha(x) + v \right] \end{split} \]

Beispiel

Gegeben sei das nichtlineare Eingrößensystem:

\[ \begin{split} \dot{x} &= \begin{bmatrix} 0 \\ x_1 + x_2^2 \\ x_1 - x_2 \end{bmatrix} + \begin{bmatrix} e^{x_2} \\ e^{x_2} \\ 0 \end{bmatrix} u \\ y &= h(x) = x_3 \end{split} \]

CODE: Beispiel

Umsetzung der Eingangs-Ausgangslinearisierung in Python / Sympy.

from sympy import init_printing
init_printing()

from sympy import *

Jacobi Matrix

def jac(f,x):
    return f.jacobian(x)

Lie Ableitung

def lie_deriv(f, h, x, n=1):
    if n==0:
        Lfh = h
    else:
        Lfh = h.jacobian(x)@f
        for i in range(1,n):
            Lfh = Lfh.jacobian(x)@f
    return Lfh

Zustände und Eingänge

x1, x2, x3 = symbols('x1 x2 x3', real=True)
x_vec =  Matrix([[x1], [x2], [x3]])

u = symbols('u', real=True)
u_vec =  Matrix([[u]])

z1, z2, z3 = symbols('z1 z2 z3', real=True)
z_vec =  Matrix([[z1], [z2], [z3]])


v = symbols('v', real=True)
v_vec =  Matrix([[v]])

Anlegen des nichtlinearen Systems

f = Matrix([
    [0],
    [x1+x2**2],
    [x1-x2],
])

f

\(\displaystyle \left[\begin{matrix}0\\x_{1} + x_{2}^{2}\\x_{1} - x_{2}\end{matrix}\right]\)

g = Matrix([
    [exp(x2)],
    [exp(x2)],
    [0],
])

g

\(\displaystyle \left[\begin{matrix}e^{x_{2}}\\e^{x_{2}}\\0\end{matrix}\right]\)

h = Matrix([
    [x3],
])

h

\(\displaystyle \left[\begin{matrix}x_{3}\end{matrix}\right]\)

Ermitteln des relativen Grades

Lf0h = lie_deriv(f, h, x_vec, 0)
display(Lf0h)
LgLf0h = lie_deriv(g, Lf0h, x_vec, 1)
display(LgLf0h)

\(\displaystyle \left[\begin{matrix}x_{3}\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}0\end{matrix}\right]\)

Lf1h = lie_deriv(f, h, x_vec, 1)
display(Lf1h)
LgLf1h = lie_deriv(g, Lf1h, x_vec, 1)
display(LgLf1h)

\(\displaystyle \left[\begin{matrix}x_{1} - x_{2}\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}0\end{matrix}\right]\)

Lf2h = lie_deriv(f, h, x_vec, 2)
display(Lf2h)
LgLf2h = lie_deriv(g, Lf2h, x_vec, 1)
display(LgLf2h)

\(\displaystyle \left[\begin{matrix}- x_{1} - x_{2}^{2}\end{matrix}\right]\)

\(\displaystyle \left[\begin{matrix}- 2 x_{2} e^{x_{2}} - e^{x_{2}}\end{matrix}\right]\)

Zustandstransformationen

Lf0h

\(\displaystyle \left[\begin{matrix}x_{3}\end{matrix}\right]\)

Lf1h

\(\displaystyle \left[\begin{matrix}x_{1} - x_{2}\end{matrix}\right]\)

Lf2h

\(\displaystyle \left[\begin{matrix}- x_{1} - x_{2}^{2}\end{matrix}\right]\)

phi_x = Matrix([[Lf0h],[Lf1h],[Lf2h]])
phi_x

\(\displaystyle \left[\begin{matrix}x_{3}\\x_{1} - x_{2}\\- x_{1} - x_{2}^{2}\end{matrix}\right]\)

z2x = dict(zip(z_vec, phi_x))
z2x

\(\displaystyle \left\{ z_{1} : x_{3}, \ z_{2} : x_{1} - x_{2}, \ z_{3} : - x_{1} - x_{2}^{2}\right\}\)

_x2z = dict(zip(phi_x, z_vec))
_x2z

\(\displaystyle \left\{ x_{3} : z_{1}, \ - x_{1} - x_{2}^{2} : z_{3}, \ x_{1} - x_{2} : z_{2}\right\}\)

Neuer Eingang, der die Zustände linearisiert

alpha = Lf3h = lie_deriv(f, h, x_vec, 3)
alpha

\(\displaystyle \left[\begin{matrix}- 2 x_{2} \left(x_{1} + x_{2}^{2}\right)\end{matrix}\right]\)

beta = LgLf2h
beta

\(\displaystyle \left[\begin{matrix}- 2 x_{2} e^{x_{2}} - e^{x_{2}}\end{matrix}\right]\)

u = 1/beta[0]*(-alpha[0]+v)
u

\(\displaystyle \frac{v + 2 x_{2} \left(x_{1} + x_{2}^{2}\right)}{- 2 x_{2} e^{x_{2}} - e^{x_{2}}}\)

fc = simplify(f+g*u)
fc

\(\displaystyle \left[\begin{matrix}- \frac{v + 2 x_{2} \left(x_{1} + x_{2}^{2}\right)}{2 x_{2} + 1}\\\frac{- v + x_{1} + x_{2}^{2}}{2 x_{2} + 1}\\x_{1} - x_{2}\end{matrix}\right]\)

Global linearisiertes System

fg_z = simplify(jac(phi_vec,x_vec)@fc).subs(_x2z)
fg_z

\(\displaystyle \left[\begin{matrix}z_{2}\\z_{3}\\v\end{matrix}\right]\)

h_z = h.subs(_x2z)
h_z

\(\displaystyle \left[\begin{matrix}z_{1}\end{matrix}\right]\)

Az =  jac(fg_z, z_vec)
Az

\(\displaystyle \left[\begin{matrix}0 & 1 & 0\\0 & 0 & 1\\0 & 0 & 0\end{matrix}\right]\)

Bz =  jac(fg_z, v_vec)
Bz

\(\displaystyle \left[\begin{matrix}0\\0\\1\end{matrix}\right]\)

Cz = jac(h_z, z_vec)
Cz

\(\displaystyle \left[\begin{matrix}1 & 0 & 0\end{matrix}\right]\)

Wir haben also das System

\[ \begin{split} \dot{x} &= \begin{bmatrix} 0 \\ x_1 + x_2^2 \\ x_1 - x_2 \end{bmatrix} + \begin{bmatrix} e^{x_2} \\ e^{x_2} \\ 0 \end{bmatrix} u \\ y &= h(x) = x_3 \end{split} \]

mittels eines Reglers gobal linearisiert:

\[ \begin{split} \begin{bmatrix} \dot{z}_1 \\ \dot{z}_2 \\ \dot{z}_3 \\ \end{bmatrix} &= \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \\ z_3 \\ \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} v \\ y &= \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \\ z_3 \end{bmatrix} \end{split} \]

Fazit

Python / Sympy bietet viele Möglichkeiten für die nichtlineare Regelungtechnik. Dieser Blogeintrag bezieht sich auf differentialgeometrischen Methoden, genauer auf die Exakte Eingangs- Ausgangslinearisierung im Eingrößenfall. Gerade das sympy.diffgeom Paket scheint sich sehr gut für diese Methoden zu eignen. Jedoch sollte der theoretische Aufwand, um sich in solche Module einzuarbeiten, nicht unterschätzt werden. Deswegen haben wir hier nur Basis-Befehle benutzt, was für erste Schritte ausreichend ist.

Gut an diesen Methoden ist, dass diese mühelos in Steuerungen implementiert werden können. Selbst Arduino und Raspberry Pico Controller sind dafür geeignet.

Aber auch andere nichtlinearen Methoden können mit Python angegangen werden. Als Beispiel soll der nichtlineare MPC Regler genannt werden. Will man aber die Optimierung auf den Steuerungen in Echtzeit durchführen, braucht man auch dort geeignete Optimierer. Besonders interessant, könnte hier der Einsatz von Echtzeit-Linux-Systemen sein, da die Entwicklungsplattform dann große Ähnlichkeit mit der Produktivsystem hat.

Referenzen

https://docs.sympy.org/latest/modules/diffgeom.html
Nichtlineare Systeme und Regelungen, 3. Auflage (Jürgen Adamy, 2018)
The Geometry of Physics, Third Edition (Theodore Frankel, 2012)