Optimal Control

Optimal Controls

The solution to optimal control problems for ordinary differential equations can be obtained by applying Pontryagin’s minimum principle. This usually yields in general a non-linear boundary value problem which has to be solved numerically. Such methods are called indirect methods as they are solving first order necessary conditions. The boundary value problems to be solved are of the form with the set of boundary conditions expressed in the form where is the vector function of the state and co-state variables, is the scalar or vector function of control variables, for some values of where each vector functions are independent. The boundary value problem also requires the satisfaction of two-point or multi-point boundary conditions. Of special interest are optimal control problems with constraints either for the control or the state variables. For such problems, the right hand side of the differential equation may be piecewise smooth, that is, there are points at which the right hand jumps as the control variable may show discontinuities. The non-smooth behaviour of the right hand side inhibits a reliable convergence of the numerical approximations towards the exact solution. A way around this problem is a transformation of the points with non-smooth behaviour to known, fixed locations . Then it is obvious that the numerical approximations converge with reliable speed of convergence. The aim of this project is to derive the required transformation in a systematic way and solve a number of typical problems.

1. Importing libraries:

       import numpy as np
       import scikits.bvp1lg.colnew as colnew
       from scipy.integrate import simps, trapz
   

2. Try Example 1: Opt Example 5 of optimal control tutorial. minimize J = int_0^T (u^2 + 3u - 2x) dt, with T = 2 subject to: x’ = x + u, p’ = 2-p, u \in [0,2], and q’ = 1 together with conditions x(0) = 5, p(q) + 7 = 0, p(q) + 3 = 0 and p(T) = 0
   • Setting up the initial parameters and guessed solution:

       T = 2.0
       S = [2 - np.log(4.5), 2 - np.log(2.5)]
     

   • Definition of

       def X(t):
     return 7*np.exp(t)-2
       def P(t):
     return 2*(1-np.exp(2-t))
       def fsub(t, z):
     x, p, q1, q2, j = z # ODE's for states and costates
     if t <= S[0]:
       u = 2
       w = q1/S[0]
     elif t>= S[1]:
       u = 0
         w = (T-q2)/(T-S[1])
     else:
       u = -.5*(p+3)
       w = (q2-q1)/(S[1]-S[0])
     return [w*(x+u), w*(2-p), 0, 0, w*(u**2 + 3*u - 2*x)]