Setting up free final time in JuMP optimization

Specific Domains Optimization (Mathematical)
1

I am trying to set up a nonlinear optimization problem using JuMP. I have been referencing some examples included in its documentation. The problem I am trying to set up is similar to the Rocket Control example: Rocket Control · JuMP

Instead of fixing time and trying to maximize altitude, what would this example look like if the final time was free and the goal was to hit a specific altitude while minimizing thrust? My main inquiry is about how to incorporate a free final time.

Thanks in advance!

2

In general, you cannot have the number of JuMP variables be a decision variable in the optimization.

You probably need to choose a fixed time horizon (and therefore, the number of discretized decision variables) that is an upper bound on the true time required.

1 Like
3

Instead of defining Δt = 0.2 / T you can define Δt as a variable.

@variable(model, eps() <= Δt)

Note that the final time is tf = T*Δt
Since you want a specific final altitude this can be added as a constraint:

@constraint(model, x_h[T] == specific_altitude)

The objective to maximize thrust is a bit unclear. Since thrust is a trajectory what do you want to maximize? The 2 norm? Maybe what you want is to reach the specific_altitude in minimum time. In that case the objective function becomes:

@objective(model,  Min, Δt)

You can also have a look at the RK4 code over here. Note that it uses function tracing and solves the problem of Dubins vehicle travelling from one orientation to another in minimum time. With function tracing JuMP is nearly as powerful as CasADi to solve optimal control problems.

2 Likes
4

you can define Δt as a variable.

This was how we had it in an earlier version of the tutorial, but doesn’t it mean the solver can choose large step sizes that might introduce inaccuracy? (I’m not an expert in these optimal control problems.)

Any suggestions for things that we could improve?

5

Yes, I do remember because I had translated the code directly for my use. So was a bit confused by seeing the example again today.

Yes, In final time free problems the final time is not known. So the step size depends on final_time/N where N is the number of steps. One way to keep step size small is to select a large value of N. But this increases the number of optimization variables, constraints etc. Another approach is to use mesh refinement where one can start with a small value of N and increase the number of grid points where the solution is having more error. It is an iterative process. I am trying to develop a package but it is in preliminary stage.

So actually I have not compared the formulation time with respect to CasADi (Since the major time is spent in the solver). I was comparing the flexibility in formulating problems. I could even translate the direct collocation example

Direct Collocation 
using JuMP
using Ipopt
using FastGaussQuadrature
using Polynomials
using GLMakie

# d: Number of points for interpolation. Degree of interpolating polynomial = d-1
d = 4
x, w = gausslegendre(d-1)

# Collocation points
tau_root = [0; (1.0 .+ x) ./ 2]

# Coefficients of the collocation equation
C = zeros(d, d)

# Coefficients of continuity equation
D = zeros(d)

# Coefficients of quadrature function
B = zeros(d)

# Construct polynomial basis
Lj::Polynomial = Polynomial([1])
for j in 1:d
    # Construct lagrange polynomials to get the polynomial basis at the collocation point
    global tau_root, C, D, B, Lj
    Lj = Polynomial([1])
    for r in 1:d
        if r != j
            Lj = Lj*Polynomial([-tau_root[r], 1]) / (tau_root[j] - tau_root[r])
        end
    end

    # Evaluate the polynomial at the end of interval
    D[j] = Lj(1.0)

    # Evaluate the time derivatives at the collocation point
    Ljd = derivative(Lj)
    for r = 1:d
        C[j, r] = Ljd(tau_root[r])
    end

    # Evaluate the integral of the polynomial to get the coefficients of the quadrature function
    LjI = integrate(Lj)
    B[j] = LjI(1.0)
end

# Time Horizon
T = 10
N = 1000
h = T/N
time = range(start = 0, step = 0.5, length = N)

ns = 2
nu = 1

# Declare Model variables
dc = Model(Ipopt.Optimizer)
@variable(dc, x[1:ns, 1:d, 1:N], start = 0)
@variable(dc, -1 <= u[1:nu, 1:N] <= 1, start = 0)

# State constraints
@constraint(dc, x[1, 1:d, 1:N] .>= -0.25)

# Initial state
@constraint(dc, x[1:ns, 1, 1] .== [0, 1])

# System model
dyn(x, u) = [((1 - x[2]^2)*x[1] - x[2] + u[1]), x[1]]

# Objective function
L(x, u) = x[1]^2 + x[2]^2 + u[1]^2

pdot = zeros(ns)
J = 0
Xend = zeros(ns)

for i in 1:N-1 
    global pdot, J, Xend
    Xend = x[1:ns, 1, i]*D[1]
    for r = 2:d
        pdot = zeros(ns)
        for j in 1:d 
            pdot = pdot + C[j, r] * x[1:ns, j, i]
        end
        @constraint(dc, pdot .== h*dyn(x[1:ns, r, i], u[1:nu, i]))
        J = J + h*B[r]*L(x[1:ns, r, i], u[1:nu, i])
        Xend = Xend + x[1:ns, r, i]*D[r]
    end
    @constraint(dc, Xend .== x[1:ns, 1, i+1])
end

@objective(dc, Min, J)

optimize!(dc)

un = u[1:nu, :]
xn = x[1:ns,1,:]

f1, ax1, s1 = stairs(time, value.(un)[1, :], color = :blue, label = "u")
lines!(ax1, time, value.(xn)[1, :], color = :black, label = "x1")
lines!(ax1, time, value.(xn)[2, :], color = :red, label = "x2")
axislegend()

I said nearly because of the limitation of registered function limitation. But not sure if it would have that much impact for optimal control problems. One question I had was for user defined operator: “if we define the hessian of the user defined operator the matrix is considered as dense”. So will this cause loss of sparsity only for the variables involved in the user defined operator or the overall hessian matrix would be dense?

1 Like
6

Only the Hessian of the user-defined function is dense. If the function is nested in other terms, then whether this propagates to become an issue is another story. Where possible, we try to preserve sparsity.

1 Like
7

Minimize thrust (not maximize). Minimizing thrust minimizes fuel, so I am trying to formulate a fuel optimal problem.

Thank you! This worked!

1 Like
8

I am picking from my notes:

The general optimal control problem I am going to consider is this

\begin{aligned} \text{minimize} &\; \left[\phi(x(t_\mathrm{f}),t_\mathrm{f}) + \int_{t_\mathrm{i}}^{t_\mathrm{f}} L(x(t),u(t),t)\mathrm{d}t\right ]\\ \text{subject to} &\; \dot x(t) = f(x(t),u(t),t),\\ &\; x(t_\mathrm{i}) = r_\mathrm{i},\\ &\; h(x(t),u(t)) \leq 0,\\ &\; \nu (x(t_f)) = 0, \;\text{for example}\; x(t_\mathrm{f}) = r_\mathrm{f} \;\text{or undetermined}. \end{aligned}

If the final time t_\mathrm{f} , or the whole interval \sigma = t_\mathrm{f}-t_\mathrm{i} is free, it can be added to the optimization variables after the following transformation
t = t_\mathrm{i} + \underbrace{(t_\mathrm{f}-t_\mathrm{i})}_{\sigma}\tau, where the new (independent) time is \tau and satisfies \tau\in[0,1] .

Finding differential on both sides: \mathrm{d}t = \sigma\mathrm{d}\tau.

But then the model of the dynamics changes to \frac{\mathrm{d}x(t)}{\mathrm{d}t} = \frac{\mathrm{d}x(t(\tau))}{\sigma\mathrm{d}\tau} = \frac{\mathrm{d}\tilde x(\tau)}{\sigma\mathrm{d}\tau}, and therefore

\frac{\mathrm{d}\tilde x(\tau)}{\mathrm{d}\tau} = \sigma f(\tilde x(\tau),\tilde u(\tau),\tau).

The cost function transforms to
\text{minimize} \; \left[\phi(x(1),1) + \int_{0}^{1} L(x(\tau),u(\tau),\tau)\mathrm{d}\tau\right ]

Source: Betts, John T. Practical Methods for Optimal Control Using Nonlinear Programming. 3rd ed. Advances in Design and Control. Society for Industrial and Applied Mathematics, 2020. https://doi.org/10.1137/1.9781611976199 , section 4.5.1 (variable phase length).

I hope it helps.

9

A convenient way to solve optimal control problems like this is with InfiniteOpt which is JuMP extension for infinite-dimensional optimization problems (e.g., dynamic optimization). It automates more complex discretization schemes (e.g., orthogonal collocation over finite elements).

It is pretty easy to formulate final time problems like this in InfiniteOpt:

1 Like