Einleitung
Am Institut für Systemdynamik wird ein numerischer Optimierungsalgorithmus zur Lösung von nichtlinearen/quadratischen Optimalsteuerungsproblemen entwickelt. Dabei wird besonderer Wert auf eine Ausnutzung der besonderen Struktur von Optimalsteuerungsproblemen gelegt, sodass eine schnelle Lösung möglich ist.
Der Algorithmus wird in ANSI-C geschrieben und ist somit plattformunabhängig. Als Bibliothek für Vektor-/Matrixoperationen kommt die ebenfalls in C implementierte MESCHACH-Bibliothek zum Einsatz. Da für die meisten Embedded-Systeme C-Compiler verfügbar sind, ist eine Anwendung auf diesen möglich.
Algorithmus
Die nichtlinearen Optimalsteuerungsprobleme, auf welche hier besonderes Augenmerk gelegt wird, haben eine Stufenstruktur der Form
$$\begin{align} & \rlap{\min_{\lbrace x_{k+1}, u_k \rbrace_{k=0}^{N-1}} \sum_{k=0}^{N-1} J_k( x_k, u_k ) + J_N( x_N ) } & & & & & & \sf{Gütefunktional} \\ {\sf{u.Nb.\quad}} x_{k+1} &= F_k(x_k, u_k) & & & & k = 0,\ldots,N-1 & & \sf{Systemzustandsgleichung} \\ x_0 &= x_I & & & & & & \sf{Anfangsbedingung} \\ 0 &\leq c_k( x_k, u_k ) & & & & k = 0,\ldots,N-1 & & \sf{Pfadbeschränkungen} \\ 0 &\leq c_N( x_N ) & & & & & & \sf{Endzustandsbeschränkung} \end{align} $$Die Lösung des o.g. Optimalsteuerungsproblems wird mithilfe einer Kombination eines Vollschritt-SQP Verfahrens [2] und eines unterlagerten, primal-dual Aktive-Menge-Verfahrens [6] berechnet. Die Schrittweite wird im SQP-Verfahren durch einen Vertrauensbereich begrenzt, welcher über einen Filtermechanismus angepasst wird [5]. Bei Auftreten eines unzulässigen QPs wird eine sog. Restaurationsphase [5] angewendet.
Die Modellfunktionen \(F_k\), \(J_k\) und \(c_k\) müssen zweifach stetig differenzierbar bezüglich der Zustände \(x_k\) und Stellgrößen \(u_k\) sein. Die Hessematrix des Gütefunktionals darf dabei indefinit sein, da bei der Lösung der unterlagerten QPs Richtungen mit negativer Krümmung bzw. mit unendlich-linearem Abstieg explizit berücksichtigt werden [3,4].
Die Struktur eines linear-quadratischen Optimalsteuerungsproblems ergibt sich zu
\begin{align} & \rlap{\min_{\lbrace x_{k+1}, u_k \rbrace_{k=0}^{N-1}} \sum_{k=0}^{N-1} \frac{1}{2} \begin{pmatrix} x_k \\ u_k \end{pmatrix}^T \begin{pmatrix} H_{xx,k} & H_{xu,k} \\ H_{xu,k}^T & H_{uu,k} \end{pmatrix} \begin{pmatrix} x_k \\ u_k \end{pmatrix} + \begin{pmatrix} f_{x,k} \\ f_{u,k} \end{pmatrix}^T \begin{pmatrix} x_k \\ u_k \end{pmatrix} + j_k} \\ & \rlap{\qquad\qquad\qquad +\frac{1}{2}x_N^T H_{xx,N} x_N + f_{x,N}^T \cdot x_N + j_N } & & & & & & \\ {\sf{u.Nb.\quad}} x_{k+1} &= A_k \cdot x_k + B_k \cdot u_k + b_k & & & & k = 0,\ldots,N-1 & & \sf{Systemzustandsgleichung} \\ x_0 &= x_I & & & & & & \sf{Anfangsbedingung} \\ 0 &\leq C_{x,k} x_k + C_{u,k} u_k + c_k & & & & k = 0,\ldots,N-1 & & \sf{Pfadbeschränkungen} \\ 0 &\leq C_{x,N} x_N + c_N & & & & & & \sf{Endzustandsbeschränkung} \end{align}wobei die Besetztheitsstruktur des KKT-Systems typischerweise ähnlich zu der in Abbildung 1 aussieht.
Abbildung 1: Besetztheitsstruktur des KKT-Systems
Bei der Lösung des QPs wird die Struktur der schwach besetzten Beschränkungs- und Hessematrizen explizit ausgenutzt und nicht eine Kondensierung durch Elimination der Systemzustände durchgeführt. Dadurch ergibt sich ein effizientes, rekursives Lösungsverfahren [1].
Es kommt ein Nullraumverfahren zum Einsatz, in dem die Beschränkungen der im Zuge des Aktive-Menge-Verfahren autretenden gleichungsbeschränkten QPs mithilfe einer QR-Zerlegung aufgelöst werden. Bei Aktivierung und Deaktivierung einer Beschränkung wird ausgenutzt, dass nachfolgende Stufen nicht neu berechnet werden müssen. Dadurch verringert sich der Rechenaufwand.
Neue SQP-Iterationen werden mit der vorigen QP-Lösung initialisiert. Falls die alte Steuertrajektorie nicht zulässig ist, wird ein Phase-1 LP gelöst, welches die Steuertrajektorie so modifiziert, dass ein möglichst guter Warmstart mit einer zulässigen Trajektorie erfolgen kann.
MATLAB®-Interface
Zur komfortablen prototypischen Anwendung des Solvers existiert ein MATLAB® MEX-Interface, welches eine Implementierung der Modellfunktionen direkt in MATLAB® ermöglicht. Dabei ist auch eine Eingabe der Zustandstransitionsgleichungen \(F_k\) als Differentialgleichungssystem zur automatischen numerischen Integration möglich, wobei hier eine Anzahl an Integrationsverfahren wählbar ist (Euler-Vorwärts, RK2-HEUN, RK4, DOPRI5).
Download
Es ist geplant in Zukunft zunächst das Matlab-Interface als Binärdatei und anschließend den Quelltext hier zu veröffentlichen.Literatur
| [1] |
Arnold, E., Tatjewski, P. and Wołochowicz, P. Two methods for large-scale nonlinear optimization and their comparison on a case study of hydropower optimization Journal of Optimization Theory and Applications, Springer, 1994, Vol. 81(2), pp. 221-248 |
| [2] |
Boggs, P. and Tolle, J. Sequential quadratic programming Acta numerica, Cambridge Univ Press, 1995, Vol. 4(1), pp. 51 |
| [3] |
Casas, E. and Pola, C. An algorithm for indefinite quadratic programming based on a partial Cholesky factorization RAIRO. Recherche opérationnelle, EDP Sciences, 1993, Vol. 27(4), pp. 401-426 |
| [4] |
Conn, A. R. and Gould, N. I. On the location of directions of infinite descent for nonlinear programming algorithms SIAM Journal on Numerical Analysis, SIAM, 1984, Vol. 21(6), pp. 1162-1179 |
| [5] |
Fletcher, R. and Leyffer, S. Nonlinear programming without a penalty function Mathematical Programming, Springer, 1997, Vol. 91(2), pp. 239-269 |
| [6] |
Goldfarb, D. and Idnani, A. A numerically stable dual method for solving strictly convex quadratic programs Mathematical programming, Springer, 1983, Vol. 27(1), pp. 1-33 |