Review of the past two weeks

During the last two weeks, I focused on the better integration of Optimization solvers for BVP problem solving in BoundaryValueDiffEq.jl. The basic funtionalities are done and the related PR, Optimization based solvers, is now mergable. What left in issue Onto Dynamic Optimization: Allow for solving via OptimizationProblem and allow inequality constraints is the interface for optimal control problems.

A classical optimal control problem is the rocket launching problem(aka Goddard problem). Say we have a rocket with limited fuel and is launched vertically. And we want to control the final altitude of this rocket so that we can make the best of the fuel in rocket to get to the highest altitude. The state variables are:

Velocity of the rocket:
Altitude of the rocket:
Mass of the rocket and the fuel:

The control variable is

Thrust of the rocket:

The dynamics of the launching can be formulated with three differential equations:

where the drag is a function of altitude and velocity:

gravity is a function of altitude:

is a constant. Suppose the final time is , we here want to maximize the . Similar solving for such optimal control problem can be found on JuMP.jl and InfiniteOpt.jl. The detailed parameters are taken from COPS.

So to model optimal control problems and solve them with BoundaryValueDiffEq.jl is like:

using BoundaryValueDiffEqMIRK, OptimizationMOI, Ipopt
h_0 = 1                      # Initial height
v_0 = 0                      # Initial velocity
m_0 = 1.0                    # Initial mass
m_T = 0.6                    # Final mass
g_0 = 1                      # Gravity at the surface
h_c = 500                    # Used for drag
c = 0.5 * sqrt(g_0 * h_0)    # Thrust-to-fuel mass
D_c = 0.5 * 620 * m_0 / g_0  # Drag scaling
u_t_max = 3.5 * g_0 * m_0    # Maximum thrust
T_max = 0.2                  # Number of seconds
T = 1_000                    # Number of time steps
Δt = 0.2 / T;                # Time per discretized step

tspan = (0.0, 0.2)
drag(x_h, x_v) = D_c * x_v^2 * exp(-h_c * (x_h - h_0) / h_0)
g(x_h) = g_0 * (h_0 / x_h)^2
function rocket_launch!(du, u, p, t)
    # u_t is the control variable (thrust)
    x_v, x_h, x_m, u_t = u[1], u[2], u[3], u[4]
    du[1] = (u_t-drag(x_h, x_v))/x_m - g(x_h)
    du[2] = x_v
    du[3] = -u_t/c
end
function rocket_launch_bc!(res, u, p, t)
    res[1] = u(0.0)[1] - v_0
    res[2] = u(0.0)[2] - h_0
    res[3] = u(0.0)[3] - m_0
    res[4] = u(0.2)[4] - 0.0
end
function constraints!(res, u, p)
    res[1] = u[1]
    res[2] = u[2]
    res[3] = u[3]
    res[4] = u[4]
end
cost_fun(u, p) = -u[end-2] # To minimize, only temporary, need to use temporary solution here similar to what we do in boundary condition evaluations.
#cost_fun(sol, p) = -sol(0.2)[4]
u0 = [v_0, h_0, m_0, 0.0]
rocket_launch_fun = BVPFunction(rocket_launch!, rocket_launch_bc!; cost = cost_fun, inequality = constraints!, bcresid_prototype = zeros(4))
rocket_launch_prob = BVProblem(rocket_launch_fun, u0, tspan; lcons = [0.0, 0.0, m_T, 0.0], ucons = [Inf, Inf, m_0, u_t_max])
sol = solve(rocket_launch_prob, MIRK4(; optimize = Ipopt.Optimizer()); dt = 0.002)

The basic idea is to formulate the optimal control problem as an unknown parameter estimation problem like we did in PR Add unknown parameters estimation. The control variables do not have dynamics attached to them, so their values can only be inferred by estimating with several constraints, e.g. boundary conditions, inequality constraints, equality constraints etc.

TODOs in next two weeks

In the following two weeks, there are a few TODOs:

Merge PR Optimization based solvers ASAP.
Finalize the interface for optimal control problems.
Finalize a suitable convention for the cost function in formulating optimal control problems, mainly about the interpolating and integral of control variables.
The generalization to heterogeneous GPU backends for Shooting methods is relatively easy, the remain issue is the evaluation of boundary conditions routine, for example, if we have a boundary condition of:

function bc!(residual, u, p, t)
    residual[1] = u(pi / 4)[1] - 1.0
    residual[2] = u(pi / 2)[1] - 1.0
end

Which would always involves scalar indexing which is always disallowed by default. We need to figure out a better way to handling this.