# Setting up free final time in JuMP optimization

**URL:** <https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [December 7, 2023, 8:59pm UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302 "2023-12-07T20:59:54Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![pluie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pluie/32/35949_2.png) [@pluie](https://discourse.julialang.org/u/pluie)\
**Post date:** [December 7, 2023, 8:59pm UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/1 "2023-12-07T20:59:54Z")

</div>

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](https://jump.dev/JuMP.jl/stable/tutorials/nonlinear/rocket_control/#Rocket-Control)

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!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 7, 2023, 9:30pm UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/2 "2023-12-07T21:30:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![A\_C](https://avatars.discourse-cdn.com/v4/letter/a/da6949/32.png) [@A\_C](https://discourse.julialang.org/u/A_C)\
**Post date:** [December 8, 2023, 3:44am UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/3 "2023-12-08T03:44:03Z")

</div>

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

```julia
@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:

```julia
@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:

```julia
@objective(model, Min, Δt)

```

You can also have a look at the RK4 code over [here](https://discourse.julialang.org/t/jump-user-defined-functions-with-vector-inputs-and-outputs/106025). 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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 8, 2023, 4:27am UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/4 "2023-12-08T04:27:16Z")

</div>

> 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.)

> [@A\_C](#):
>
> With function tracing JuMP is nearly as powerful as CasADi to solve optimal control problems.

Any suggestions for things that we could improve?

---

<div class="post-metadata">

**Author:** ![A\_C](https://avatars.discourse-cdn.com/v4/letter/a/da6949/32.png) [@A\_C](https://discourse.julialang.org/u/A_C)\
**Post date:** [December 8, 2023, 5:13am UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/5 "2023-12-08T05:13:02Z")

</div>

> [@odow](#):
>
> This was how we had it in an earlier version of the tutorial

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.

> [@odow](#):
>
> doesn’t it mean the solver can choose large step sizes

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](https://github.com/A-C1/NumericalOptimalControl/blob/master/src/meshrecomputation.jl) but it is in preliminary stage.

> [@odow](#):
>
> Any suggestions for things that we could improve?

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](https://github.com/casadi/casadi/blob/main/docs/examples/matlab/direct_collocation.m)

> **Direct Collocation**
>
> ```julia
> 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?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 8, 2023, 5:21am UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/6 "2023-12-08T05:21:31Z")

</div>

> [@A\_C](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![pluie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pluie/32/35949_2.png) [@pluie](https://discourse.julialang.org/u/pluie)\
**Post date:** [December 8, 2023, 2:07pm UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/7 "2023-12-08T14:07:25Z")

</div>

> [@A\_C](#):
>
> The objective to maximize thrust is a bit unclear.

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

> [@A\_C](#):
>
> Instead of defining `Δt = 0.2 / T` you can define `Δt ` as a variable.
> 
> ```julia
> @variable(model, eps() <= Δt)
> 
> ```
> 
> Note that the final time is `tf = T*Δt`

Thank you! This worked!

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [December 8, 2023, 2:29pm UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/8 "2023-12-08T14:29:01Z")

</div>

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](https://doi.org/10.1137/1.9781611976199) , section 4.5.1 (variable phase length).

I hope it helps.

---

<div class="post-metadata">

**Author:** ![pulsipher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulsipher/32/10461_2.png) [@pulsipher](https://discourse.julialang.org/u/pulsipher)\
**Post date:** [December 12, 2023, 2:11pm UTC](https://discourse.julialang.org/t/setting-up-free-final-time-in-jump-optimization/107302/9 "2023-12-12T14:11:07Z")

</div>

A convenient way to solve optimal control problems like this is with [InfiniteOpt](https://github.com/infiniteopt/InfiniteOpt.jl) 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:

- [Free terminal time problems · infiniteopt/InfiniteOpt.jl · Discussion #289 · GitHub](https://github.com/infiniteopt/InfiniteOpt.jl/discussions/289)
- [Minimum time manouevering problem + boundaries · infiniteopt/InfiniteOpt.jl · Discussion #219 · GitHub](https://github.com/infiniteopt/InfiniteOpt.jl/discussions/219#discussioncomment-2160386)
- [Minimum time objective ensuring waypoints touched in hovercraft example · infiniteopt/InfiniteOpt.jl · Discussion #192 · GitHub](https://github.com/infiniteopt/InfiniteOpt.jl/discussions/192)
- [Minimum time problems · infiniteopt/InfiniteOpt.jl · Discussion #188 · GitHub](https://github.com/infiniteopt/InfiniteOpt.jl/discussions/188)