Skip to content

Commit

Permalink
test
Browse files Browse the repository at this point in the history
  • Loading branch information
RickGelhausen committed Jun 25, 2024
1 parent 59fff1b commit 52d6c0e
Showing 1 changed file with 0 additions and 111 deletions.
111 changes: 0 additions & 111 deletions exercise-sheet-5.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -162,114 +162,3 @@ mit $ r_i(R) = R $ für $ i \in \{1, 2, 3\} $ und $ r_i(R) = 2R $ für $ i \in \
```{r, echo=FALSE, out.width='50%', fig.align='center', fig.show='hold', fig.cap='**Abbildung 4** - Links: Skizze des Pendels. Rechts: Trajektorie des Pendels mit Initialzustand $x_0 = \\left[\\frac{3\\pi}{4} \\quad 0\\right]^T$ und ohne Steuerung, $u_k = 0 \\forall k$.'}
knitr::include_graphics("figures/sheet-5/p5.png")
```

Wir betrachten ein Pendel.
Dessen Position ist durch den Winkel $\theta$ eindeutig bestimmt. Dabei entspricht $\theta = \pi$ der Position, in der es gerade nach unten hängt.
Eine Illustration finden Sie in Abbildung 4.
In der Aufhängung des Pendels sitzt ein Motor, sodass es anhand eines Drehmoments $u$ gesteuert werden kann.
Die Winkelgeschwindigkeit ist $\omega = \dot{\theta}$.
Fassen wir Position und Geschwindigkeit im Zustandsvektor $x = \begin{bmatrix} \theta \\ \omega \end{bmatrix}^T$ zusammen, können wir die Dynamik des Pendels durch die gewöhnliche Differentialgleichung (ordinary differential equation, ODE)

\[
\dot{x} = \begin{bmatrix} \dot{\theta} \\
\dot{\omega} \end{bmatrix} = f(x, u) := \begin{bmatrix} \omega \\
\sin \theta + u \end{bmatrix}
\]

beschreiben. (Zur Vereinfachung haben wir hier alle Einheiten ignoriert. Eigentlich müsste vor den Sinus ein Faktor mit Einheit,
da dieser dann mit dem Drehmoment summiert wird.)


Wir betrachten das Pendel über die Zeitdauer $ T $ und diskretisieren diese in $ N $ Zeitschritte. Zur Simulation verwenden wir das Runge-Kutta-Verfahren 4. Ordnung (RK4) und erhalten dadurch die diskretisierte Dynamik

\[
x_{k+1} = F(x_k, u_k),
\]

wobei $ x_k $ der Zustand zum diskreten Zeitpunkt $ k $, $ k = 0, \ldots, N $, und $ h = \frac{T}{N} $ der Integrationschritt ist. Wir fassen die Zustände und die Steuerungsinputs zu allein Zeitpunkten in den Matrizen

\[
X := \begin{bmatrix} x_0 & \ldots & x_N \end{bmatrix} \in \mathbb{R}^{2 \times (N+1)} \quad \text{und} \quad U := \begin{bmatrix} u_0 & \ldots & u_{N-1} \end{bmatrix} \in \mathbb{R}^{1 \times N}
\]

zusammen.

Unser Ziel ist es nun, das Pendel aus seiner herabhängenden Ruhelage $x_0 = \begin{bmatrix} \pi & 0 \end{bmatrix}^T$ zum Zeitpunkt $k = 0$ in die aufrechte stehende Position $x_N = \begin{bmatrix} 0 & 0 \end{bmatrix}^T$ zum Zeitpunkt $k = N$ zu schwingen.
Dabei wollen wir den Steuerungsaufwand $ L(U) := \sum_{k=0}^{N-1} u_k^2 $ minimieren.

**Aufgaben:**

1. Formulieren Sie unser Optimalsteuerungsproblem als nichtlineares Programm. Dabei sollen $ x_0, \ldots, x_N $ und $ u_0, \ldots, u_{N-1} $ die Entscheidungsvariablen sein. Die Nebenbedingungen sind die Dynamik, sowie die Start- und Zielposition.
2. Diskutieren Sie kurz, ob das Problem konvex ist.
3. Benutzen Sie das bereitgestellte Template, um das NLP mit CasADi und IPOPT zu lösen. Erstellen Sie Plots der optimalen Trajektorien von $ \theta $, $ \omega $ und $ u $, mit der diskreten Zeit $ k $ auf der $ x $-Achse. Sie können Ihre Lösung außerdem mit der bereitgestellten Animation (pendulum.gif) vergleichen.
4. Wir führen nun eine zusätzliche Beschränkung der Steuerung ein. Zu allen Zeitpunkten soll gelten: $ |u_k| \leq u_{\text{max}} $. Diskutieren Sie kurz, wie sich der optimale Wert der Zielfunktion dadurch verändert.
5. Erweitern Sie Ihre NLP Formulierung um die zusätzliche Beschränkung. Beachten Sie, dass Sie hierbei die Betragsfunktion $ | \cdot | $ nicht verwenden sollten, da diese an der Stelle 0 nicht differenzierbar ist. Dies kann zu Problemen führen. Finden Sie stattdessen eine Umformulierung dieser Nebenbedingung.
6. Erweitern Sie Ihre Implementierung um die zusätzliche Nebenbedingung. Verwenden Sie $ u_{\text{max}} = 0.13 $.

---

# Aufgabe 3: Optimalsteuerung mit acados (Bonus)

In dieser Aufgabe lernen wir das open-source Softwarepaket acados kennen.
acados bietet eine Sammlung effizienter Algorithmen, die auf das Lösen von Optimalsteuerungsproblemen spezialisiert sind.
Dafür implementiert acados ein SQP-Verfahren sowie numerische Integratoren für Differentialgleichungen.
Zum Lösen der im SQP-Verfahren anfallenden QP wird auf moderne open-source QP-Löser zurückgegriffen, z.B. HP ipm, qpOASES, OSQP, DQAP.
Optimalsteuerungsprobleme werden mit

```{r, echo=FALSE, out.width='50%', fig.align='center', fig.show='hold', fig.cap='**Abbildung 4** - Links: Skizze des Pendels. Rechts: Trajektorie des Pendels mit Initialzustand $ x_0 = \left[\frac{3\pi}{4} \quad 0\right]^T $ und ohne Steuerung, $ u_k = 0 \forall k $.'}
knitr::include_graphics("figures/sheet-5/p6.png")
```

CasADi's symbolische Variablen definieren, welches auch zur Berechnung von Ableitungen verwendet wird. Für die grundlegenden Operationen der linearen Algebra (z.B. Matrix-Matrix-Multiplikationen) wird BLASFEO verwendet. Auf Grundlage der erwähnten Komponenten generiert acados dann C-Code, welcher ohne externe Abhängigkeiten auskommt. Dieser kann insbesondere auch auf eingebetteten Systemen ausgeführt werden, was die Optimalsteuerung technischer Systeme in Echtzeit ermöglicht. Unter folgenden Links finden Sie weitere Informationen:

- Dokumentation: [https://docs.acados.org/](https://docs.acados.org/)
- Installation: [https://docs.acados.org/installation/](https://docs.acados.org/installation/)
- Python-Interface Installation: [https://docs.acados.org/python_interface/](https://docs.acados.org/python_interface/)
- acados OCP-Formulierung: [https://github.com/acados/acados/blob/master/docs/problem_formulation/ocp_mex.pdf](https://github.com/acados/acados/blob/master/docs/problem_formulation/ocp_mex.pdf)
- acados Forum: [https://discourse.acados.org](https://discourse.acados.org)

### Pendel auf einem Wagen

Als Beispielproblem betrachten wir ein Pendel, das auf einem Wagen befestigt ist, siehe Abb. 5. Der Mittelpunkt des Wagens hat die horizontale Position $ p_x $, welche wir durch Ausüben einer Kraft $ F $ beeinflussen können. Die Auslenkung des Pendels, welches die Länge $ l $ hat, ist durch den Winkel $ \beta $ beschrieben. Der Wagen hat die Masse $ M $, und an der Spitze des Pendels ist ein Ball mit Masse $ m $ befestigt. Auf die Pendelmasse wirkt außerdem die Erdbeschleunigung $ g $. Die horizontale Geschwindigkeit des Wagens ist $ v_x $ und die Winkelgeschwindigkeit des Pendels ist $ \omega $. Durch Zusammenfassen der (Winkel)positionen und -geschwindigkeiten im Zustandsvektor $ x $ können wir das System durch folgende Differentialgleichung beschreiben:

\[
\dot{x} =
\begin{bmatrix}
v_x \\
\omega
\end{bmatrix} =
\begin{bmatrix}
-m \omega^2 \sin \theta + m \omega \cos \theta \sin \theta + F \\
-\frac{m \omega^2 \cos^2 \theta + M (1 - \cos^2 \theta) + \frac{(m + M) g \sin \theta}{l (1 - F \cos \theta)}}{M + \frac{F \cos^2 \theta}{m}}
\end{bmatrix}.
\]

Unser Ziel ist es nun das Pendel aus einer herabhängenden Position in eine aufrechte Position ($ \beta = 0 $) zu schwingen, während der Wagen am Ende die Position $ p_x = 0 $ haben soll. Unser Steuerungseingang ist hierbei $ u = F $. Dies soll innerhalb des Zeitintervalls $ t \in [0, T] $ passieren. Wir drücken dies als das folgende Optimalsteuerungsproblem aus,

\[
\min_{x(\cdot), u(\cdot)} \int_0^T \frac{1}{2} x(t)^T Q x(t) + \frac{1}{2} u(t)^T R u(t) dt + \frac{1}{2} x(T)^T Q_e x(T)
\]

unter den Nebenbedingungen

\[
\begin{aligned}
x(0) &= x_0, \\
\dot{x}(t) &= f(x(t), u(t)), \quad t \in [0, T], \\
-u_{\text{max}} &\leq u(t) \leq u_{\text{max}}, \quad t \in [0, T].
\end{aligned}
\]

wo bei $ \hat{x}_0 $ der gegebene initiale Zustand des Systems ist. Anders als wir es bisher in der Vorlesung gesehen haben, sind in dem obigen Optimalsteuerungsproblem die Entscheidungsvariablen $ x(\cdot) $ und $ u(\cdot) $ Funktionen der Zeit. Es handelt es sich deshalb nicht um ein NLP, und wir können es auch nicht ohne weiteres auf einem Computer repräsentieren. Hierfür muss es erst durch numerische Integration in der Zeit diskretisiert werden, wie wir es bereits in der vorherigen Aufgabe mit dem RK4-Verfahren gemacht haben. Da allerdings acados dies für uns übernimmt und eine Vielzahl effizienter Integrationsverfahren hierfür bereitstellt, übergeben wir das Optimalsteuerungsproblem in kontinuierlicher Zeit.

**Aufgaben:**

1. Installieren Sie acados sowie das zugehörige Python-Interface. Die Links dafür sind weiter oben gegeben. Versichern Sie sich, dass Ihre Installation funktioniert, indem Sie das Minimalbeispiel `minimal_example_ocp.py` ausführen (vgl. Installationsanleitung Python-Interface).
2. Das Optimalsteuerungsproblem ist für Sie bereits in `cartpole.py` implementiert. Machen Sie sich kurz mit dem Code vertraut und führen Sie ihn dann aus.
3. Wir wollen eine zusätzliche Nebenbedingung auf die Geschwindigkeit $ v_x $ einführen. Diese ist

\[
-u_{\text{max}} \leq u(t) \leq u_{\text{max}}, \quad t \in [0, T],
\]

mit $ v_x = 5 \, \text{m/s} $. Erweitern Sie `minimal_example_ocp.py` um diese Nebenbedingung, und lösen Sie das Optimalsteuerungsproblem erneut.

0 comments on commit 52d6c0e

Please sign in to comment.