{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solve a linear-quadratic ODE Optimal Control Problem\n", "\n", "## Notebook to demonstrate how to solve an linear quadratic optimal control problem via oppy\n", "\n", "We consider the linear quadratic optimal control problem \n", "\\begin{align*}\n", "&\\min_{y,u}J(y,u),\\\\\n", "\\mathrm{s.t.}& \\quad\\dot{y}(t)=Ay(t)+Bu(t)\\quad \\mathrm{in}\\,t\\in(0,T)\\quad \\mathrm{with}\\quad y(0)=y_0\n", "\\end{align*}\n", "\n", "\n", "where\n", "$J(y,u)=\\frac{1}{2}||y(T)-y_T||^2+\\frac{\\mu}{2}||y(t)-y_d(t)||^2+\\frac{\\nu}{2}||u(t)||^2$\n", "is the costfunctional, $u\\in L^2(0,T;\\mathbb{R^m})$ is the control function, $y$ is the solution of the state equation, $y_T=[0,1]^T$ is the desired final state, $y_d\\equiv[1,1]^T$ desired trajectory and $\\mu,\\nu$ regularization parameters.\n", "\n", "To approach this problem we define the reduced cost functional:\n", "\\begin{align*}\n", "\\hat{J}(u)=J(\\mathcal{S}(u),u),\n", "\\end{align*}\n", "where $\\mathcal{S}$ is the solution operator, which means $y=\\mathcal{S}(u)$ solves the differntial equation in the constraints with the corresponding control $u$. One can show that solving the reduced problem is equivalent to the original problem (see [1], Chapter 3.1). \n", "\n", "The idea is to use an unconstrained optimization method on the reduced problem. \n", "Therefore, following the Lagrange method (see [1], Chapter 4.1) we derive the gradient condition for computing the optimal solution $(\\hat{y},\\hat{u})$ of the linear quadratic optimal control problem:\n", "\\begin{align*}\n", "\\nabla \\hat{J}(\\hat{u})= \\nu \\hat{u}+B^T\\hat{p}=0,\n", "\\end{align*}\n", "where $\\hat{p}$ solves the adjoint equation:\n", "\\begin{align*}\n", "\\begin{array}{r l r}\n", "-\\dot{p}(t)&=A^Tp(t)+\\mu(\\hat{y}-y_d),&\\mathrm{for}\\,t\\in[0,T],\\\\\n", "p(T)&=\\hat{y}(T)-y_T.&\\end{array}\n", "\\end{align*}\n", "Moreover for evaluating the reduced cost functional $\\hat{J}(u)$ we solve the differential equation in the constraints by integrating forward and for computing the reduced gradient (if required) $\\nabla\\hat{J}(u)$ we integrate the adjoint equation backwards. These algorithms are written in a seperate file called $\\small{\\textbf{functions.py}}$. \n", "\n", "The parameter data for the linear quadratic problem are written in the $\\small{\\textbf{problem_data.py}}$ and thus have to be imported.\n", "\n", "[1] Gabriele Ciaramella, Lecture Notes: $\\scriptsize{\\textit{Optimal Control of Ordinary Differential Equations}}$, July 21, 2019.\n", "\n", "First let's import all dependencies needed." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "from oppy.options import Options\n", "from oppy.unconOpt import gradmethod\n", "from oppy.unconOpt import armijo\n", "\n", "import functions_linear_quadratic as fun\n", "import problem_data_linear_quadratic as data\n", "\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we define the linesearch methods Armijo and Wolfe for the Gradient method, Quasi Newton method and Fletcher Reeves method. Therefore the reduced cost functional and its first derivation is required which is given by $\\texttt{fun.fAndG}$.\n", "\n", "Note: The parameters $y_0, y_T$ are in $\\mathbb{R}^M$ and the matrix $y_d\\in\\mathbb{R}^{M\\times N}$ are the values of a time dependent function on an equidistant grid with mesh size $dt$, where $M$ is the dimension number in space and $N$ is the number of discretization points. \n", "\n", "Moreover, take care that oppy requires one dimensional arrays as input functions!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define cost function\n", "def fhandle(x):\n", " # Check the shape of x\n", " if x.ndim == 1:\n", " x = np.reshape(x,(data.M, data.N))\n", " f = fun.fAndG(x, data.data, 1)\n", " return f\n", "\n", "def dfhandle(x):\n", " if x.ndim == 1:\n", " x = np.reshape(x,(data.M, data.N))\n", " f,grad = fun.fAndG(x, data.data,2)\n", " grad = np.reshape(grad, (data.M* data.N))\n", " return grad\n", "\n", "# Armijo and Wolfe linesearch options\n", "initstep = 1\n", "beta = 0.5\n", "amax = 30\n", "c1 = 1e-4\n", "c2 = 0.4\n", "\n", "options_arm = Options(initstep=initstep, c1=c1, beta=beta, nmax=amax)\n", "options_wol = Options(initstep=initstep, c1=c1, c2=c2, beta=beta, nmax=amax)\n", "\n", "def lineNorm(x,d):\n", " nor = data.dt*np.sum(np.sum(dfhandle(x)*d))\n", " return nor \n", " \n", "# Define Armijo linesearch\n", "def linesearch_armijo(x, d):\n", " t, flag = armijo(fhandle, dfhandle, x, d, options_arm)\n", " return t, flag" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have to set the options for the optimization methods. The standard discretized $L_2$-norm has been chosen for this problem and for all optimization methods. Furthermore the options with armijo or wolfe line search has to be set separately and the Symmetric Rank 1 update method for the Quasi Newton method is used." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Options for Gradient scheme with Armijo linesearch\n", "opt_grad_armijo = Options(disp=True, linesearch=linesearch_armijo)\n", "\n", "# Initial guess \n", "u0 = np.random.rand(data.M*data.N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we can run the optimization methods with an random initial control $u\\in\\mathbb{R}^{M\\times N}$. For plotting purpose we reshape the final solution." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Gradient scheme with Armijo\n", "result_grad_arm = gradmethod(fhandle, dfhandle, u0, opt_grad_armijo)\n", "res_grad_arm = result_grad_arm.res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The gradient method with armijo linesearch method stops with a total iteration number of 40. For less iterations we try another linesearch e.g. Wolfe linesearch. Therefore, we import the dependencies, set the options for the linesearch and finally execute the gradient algorithm." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from oppy.unconOpt import wolfe\n", "\n", "# Define Wolfe linesearch\n", "def linesearch_wolfe(x, d):\n", " t, flag = wolfe(fhandle, dfhandle, x, d, options_wol)\n", " return t, flag\n", "opt_grad_wolfe = Options(disp=True, linesearch=linesearch_wolfe)\n", "\n", "# Options for Gradient scheme with Wolfe linesearch\n", "result_grad_wol = gradmethod(fhandle, dfhandle, u0, opt_grad_wolfe)\n", "residual_grad_wol = result_grad_wol.res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Wolfe linesearch improves the gradient algorithm, since the iteration number is 27 and thus less than the gradient method with Armijo linesearch. We try another optimization algorithm which might be faster than the previous examples. Therefore, we import the nonlinear CG algorithm with default settings (except for $\\textit{disp}$ option) and use the same initial conditions as before." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "from oppy.unconOpt import nonlinear_cg\n", "\n", "# Nonlinear_CG scheme\n", "result_nonlin_cg = nonlinear_cg(fhandle, dfhandle, u0, Options(disp=True))\n", "residual_nonlin_cg = result_nonlin_cg.res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The nonlinear CG method is faster than the previous methods with a total iteration number of 22. One can try a non-default linesearch which will mostly lead to similar results. The next method will be the fastest of all methods before. We import and execute the Quasi Newton method with the default setting." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "from oppy.unconOpt import quasinewton\n", "\n", "# Quasi Newton scheme\n", "result_quasi_newton = quasinewton(fhandle, dfhandle, u0, Options(disp=True))\n", "residual_quasi_newton = result_quasi_newton.res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This method only needs 6 iterations steps to terminate. Similar to the nonlinear CG, a different linesearch does not change the result much. An error plot with respect to iteration numbers of the different optimization methods with a semi y-axis log scale is used to compare the different schemes." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "resisdual_grad_arm = result_grad_arm.res\n", "resisdual_grad_wol = result_grad_wol.res\n", "resisdual_quasi_newton = result_quasi_newton.res\n", "resisdual_nonlin_cg = result_nonlin_cg.res\n", "\n", "# Plot convergence rate\n", "plt.semilogy(resisdual_grad_arm, label='Gradient + Armijo')\n", "plt.semilogy(resisdual_grad_wol, label='Gradient + Wolfe')\n", "plt.semilogy(resisdual_nonlin_cg, label='Nonlinear CG')\n", "plt.semilogy(resisdual_quasi_newton, label='Quasi Newton')\n", "plt.title('Norm of the gradient of different schemes')\n", "plt.grid(True)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All methods converge with an absolute and relative residual norm tolerance of $10^{-4}$. In this linear quadratic problem the Quasi Newton method with the Wolfe linesearch provides the best solution with less than 10 total number of iterations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For better understanding of the problem, we plot the states $y$ without a control function $u$ with respect to the time $t$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Plot uncontrolled state\n", "u0 = np.zeros((data.M,data.N))\n", "y = fun.forward(u0,data)\n", "\n", "xAxe = np.linspace(0,data.T, num=data.N)\n", "plt.plot(xAxe,y[0,:], label='$y_1(t)$')\n", "plt.plot(xAxe,y[1,:], label='$y_2(t)$')\n", "plt.title('Solution of the state with no control')\n", "plt.grid(True)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now extract the optimal solution of the results of the quasi Newton method and solve the differential equation, which is the constraint of the optimal control problem. Therefore we use the forward solver." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot optimal solution\n", "u_new_sol = np.reshape(result_quasi_newton.x, (data.M, data.N))\n", "yfinal = fun.forward(u_new_sol,data.data)\n", "\n", "xAxe = np.linspace(0,data.T, num=data.N)\n", "plt.plot(xAxe,yfinal[0,:], label='$y_1(t)$')\n", "plt.plot(xAxe,yfinal[1,:], label='$y_2(t)$')\n", "plt.plot(xAxe,u_new_sol[0,:], label='$u_1(t)$')\n", "plt.plot(xAxe,u_new_sol[1,:], label='$u_2(t)$')\n", "plt.title('Optimal solution $y_1$ and $y_2$ with the optimal controls $u_1$ and $u_2$')\n", "plt.grid(True)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The state outcome with optimal control aligns with the desired state and desired trajectory in the cost function, since they were chosen as $y_T=[0,1]^T$ respectively $y_d\\equiv[1,1]^T$. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" } }, "nbformat": 4, "nbformat_minor": 4 }