# Solving an optimal control problem with JuMP

I am trying to solve an optimal control problem using `JuMP`.

My problem involves a set of ODEs and an objective function that needs the final conditions of the solution of the ODEs. The input on which I want to optimize is actually a function (called the control function). This is the kind of problem that one would usually solve with NLOptControl.jl.

I am using the `DifferentialEquations` package to solve the ODEs. I give the control function as an argument to the ODE function.

Let me illustrate this using an example from the `DifferentialEquations` docs. In this example you can see that `M` is a function of time that is passed to the `ODEProblem`, and it is named `p(t)` in `function pendulum!(du,u,p,t)`

``````using DifferentialEquations
using Plots

l = 1.0                             # length [m]
m = 1.0                             # mass[m]
g = 9.81                            # gravitational acceleration [m/s²]

function pendulum!(du,u,p,t)
du = u                    # θ'(t) = ω(t)
du = -3g/(2l)*sin(u) + 3/(m*l^2)*p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end

θ₀ = 0.01                           # initial angular deflection [rad]
ω₀ = 0.0                            # initial angular velocity [rad/s]
u₀ = [θ₀, ω₀]                       # initial state vector
tspan = (0.0,10.0)                  # time interval

M = t->0.1sin(t)                    # external torque [Nm]

prob = ODEProblem(pendulum!,u₀,tspan,M)
sol = solve(prob)

``````

Since the optimizer cannot get a function as an input, I use a discrete array of N points, which are constructed into a function using linear interpolation.

Here is an example of doing it with just 2 points `g1, g2`:

``````function opt_func_g(g1, g2)
tspan = (0.0, 10.0)

g_func(t) = g1 + (g2-g1)*t/10

u0 = [0.0, 1.0]
prob = ODEProblem(pendulum, u0, tspan, g_func)
sol = DifferentialEquations.solve(prob)

final_cond = sol[1, end]
return final_cond
end
``````

Then to solve with `JuMP`

``````model = Model()
@variable(model, g1)
@variable(model, g2)
@constraint(model, -10 <= g1 <= 10)
@constraint(model, -10 <= g2 <= 10)
register(model, :opt_func_g, 2, opt_func_g, autodiff=true)
@NLobjective(model, Min, opt_func_g(g1, g2))
set_optimizer(model, Ipopt.Optimizer)
optimize!(model)
``````

I get an error related to the forward differentiation:

``````This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        4
Number of nonzeros in Lagrangian Hessian.............:        0

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#107#109"{typeof(opt_func_g)},Float64},Float64,2})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
Float64(::T) where T<:Number at boot.jl:718
Float64(!Matched::Int8) at float.jl:60
...

Stacktrace:
 convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#107#109"{typeof(opt_func_g)},Float64},Float64,2}) at ./number.jl:7
``````

I see that the problem is giving the `ODEProblem` a function `g(t)` based on the parameters `g1, g2`. Is there a way to solve this issue?

One way I have tried is to around this is to give the `ODEProblem` the two parameters `g1, g2` and then define the function inside the ODE function (i.e. `pendulum` in my example above). Meaning:

``````function pendulum!(du,u,p,t)
g1, g2 = p
M(t) = g1 + (g2-g1)*t/10
du = u                    # θ'(t) = ω(t)
du = -3g/(2l)*sin(u) + 3/(m*l^2)*p(t) # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)
end
``````

But this does not work as well (I get the same error).

Is there a way to solve such an optimal control problem using this package?

(Please forgive any ignorance on my side - this is my first experience solving an optimization problem of this sort).

2 Likes

Is the reason you cant use NLOptControl.jl, that they haven’t updated to Julia 1.0+?

what if you do `typeof(g1)[0.0,1.0]` to make the state propagate the dual numbers?

I found this part in the documentation of `DifferentialEquations.jl` that discusses the issue of forward differentiation that helped me understand what you are referring to.

This indeed helps!

I have now another issue that hopefully can be solved in a similar way.

My target function needs to compute an elliptic integral. For that I use the Elliptic.jl package. Say I have a function

``````function target(final_cond, params)
# do some stuff
# among them calculate the elliptic integral
# I put here some operations in the argument to indicate
# that I do something to the parameters that go into Elliptic.
res = Elliptic.F(final_cond/params, params + params)
return res
end
``````

And this function is placed in the optimization function

``````function opt_func_g(g1, g2)
tspan = (0.0, 10.0)

g_func(t) = g1 + (g2-g1)*t/10

u0 = typeof(g1)[0.0, 1.0]  # Thanks to @ChrisRackauckas
prob = ODEProblem(pendulum, u0, tspan, g_func)
sol = DifferentialEquations.solve(prob)

final_cond = sol[1, end]

# Here is the target function
result = target(final_cond, params)
return result
end
``````

I guess I need to perform the same `typeof` assignment here. But I can’t seem to do it properly. Do you know how?

I have installed Julia 0.6 in order to use NLOptControl.jl and tested it. But unfortunately there is no way there to access the value of time when specifying the equations of motion (EOM). For example, writing an equation of the form `x2[j] * sin(t)` is impossible.

I saw some examples of solving optimal control problems with `JuMP` but they were all specifying the EOM as constraints to the solver (e.g. the rocket example). I think that this kind of implementation will not give the proper accuracy in the solution of the EOM. Therefore, I was trying to use the `DifferentialEquations` solver and optimize over my function this way. But unfortunately it seems the `ForwardDiff` cannot pass the solver.

Is the issue with this one that Elliptic.jl is not AD-differentiable?

I am not sure exactly what it means to be AD-differentiable.

The error I get is

``````This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        4
Number of nonzeros in Lagrangian Hessian.............:        0

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#107#109"{typeof(opt_func_g_arr_elliptic)},Float64},Float64,4})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:200
Float64(::T) where T<:Number at boot.jl:718
Float64(!Matched::Int8) at float.jl:60
...

Stacktrace:
 F(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#107#109"{typeof(opt_func_g_arr_elliptic)},Float64},Float64,4}, ::Float64) at /Users/roi/.julia/packages/Elliptic/WKWjT/src/Elliptic.jl:79
 θ_elliptic_opt(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#107#109"{typeof(opt_func_g_arr_elliptic)},Float64},Float64,4}, ::Float64, ::Float64) at ./In:3
``````

Yeah, it looks like it’s just that Elliptic.jl fails on Dual numbers. It’s explicitly trying to make things Float64 instead:

so that would need to get opened up / fixed.

Does this mean that if I take the code and instead of `Float64` I put in `Real` it will work?

``````function K(m::Float64)
if m < 1.
drf,ierr = SLATEC.DRF(0., 1. - m, 1.)
@assert ierr == 0
return drf
elseif m == 1.
return Inf
elseif isnan(m)
return NaN
else
throw(DomainError(m, "argument m not <= 1"))
end
end
K(x::Float32) = Float32(K(Float64(x)))
K(x::Real) = K(Float64(x))
``````

put in

``````function K(m::Real)
if m < 1.
drf,ierr = SLATEC.DRF(0., 1. - m, 1.)
@assert ierr == 0
return drf
elseif m == 1.
return Inf
elseif isnan(m)
return NaN
else
throw(DomainError(m, "argument m not <= 1"))
end
end
``````

You then probably need to change `0.` to `zero(m)` etc.

I don’t know if I am missing something or not.

When I do the following

``````julia> m = 1.0
1.0

julia> typeof(m)
Float64

julia> n = zero(m)
0.0

julia> typeof(n)
Float64

julia> k = convert(Real, m)
1.0

julia> typeof(k)
Float64
``````

Whatever I do, the type remains `Float64`. How can I figure out what I need to change in the function? Is there a systematic way to understand this?

I see that any number `number` I can replace with `typeof(m)(number)`.

What can I do with this?
`const D1MACH1 = floatmin(Float64)`

What Chris was trying to say is that the type constraint should be loosened from explicitly asking for `Float64` to the abstract type `Real`. Note that actual values can’t have type `Real` as it is an abstract type:

``````julia> isabstracttype(Real)
true

julia> Float64 <: Real
true

julia> supertype(Float64)
AbstractFloat

julia> supertype(AbstractFloat)
Real
``````

so if you write a function with `Real` inputs it will allow more than just `Float64`, namely all types `T` for which `T <: Real == true`:

``````julia> f(x::Float64) = "Looks like we've got a Float here!"
f (generic function with 1 method)

julia> f(x::Real) = "This could be any Real now"
f (generic function with 2 methods)

julia> f(2.0)
"Looks like we've got a Float here!"

julia> f(Float32(2.0))
"This could be any Real now"
``````

As AD relies on dual numbers, functions will have to accept dual numbers to work with AD. With a `Float64` costraint this doesn’t work, while with a `Real` constraint it will work provided that `DualNumber <: Real` (not sure which dual type is relevant here exactly).

1 Like

I see. Thank you!!

So from what you wrote I understand that I can replace any number `x` with `typeof(m)(x)` .
Say, take `sin(4)` and put in `sin(typeof(m)(4))`. Correct?

And, do you know what can I dot for these constant decelerations?

``````const D1MACH1 = floatmin(Float64)
const D1MACH2 = floatmax(Float64)
const D1MACH3 = eps(Float64)/2
const D1MACH4 = eps(Float64)
``````

Well ideally there wouldn’t be any hardcoded numbers in the code (other than maybe zeros or ones as Chris mentioned, which you could replace with `zero(x)` etc)

As for those constants that should just work if yo replace `Float64` with `typeof(x)` where `x` is your function input.

1 Like

Is it not possible to add an extra state `t` with state-derivative `tdot = 1`? This should serve as time.

1 Like

This is a great idea. I tried to define another state like this

``````dx = :(1.0)
``````

But I get the following error:

``````MethodError: Cannot `convert` an object of type Float64 to an object of type Expr
This may have arisen from a call to the constructor Expr(...),
since type constructors fall back to convert methods.

Stacktrace:
 setindex!(::Array{Expr,1}, ::Float64, ::Int64) at ./array.jl:578
``````

I tried e.g. converting to an Abstract type

``````dx = :(convert(Real, 1.0))
``````

and this if fine, but then when I do `configure` I get

``````julia> configure!(n;(:Nck=>),(:finalTimeDV=>false), (:tf=>tplan))

Unrecognized function "convert" used in nonlinear expression.

Stacktrace:
 macro expansion at /Users/roi/.julia/v0.6/JuMP/src/parsenlp.jl:96 [inlined]
 macro expansion at /Users/roi/.julia/v0.6/JuMP/src/macros.jl:1412 [inlined]
 macro expansion at /Users/roi/.julia/v0.6/NLOptControl/src/diffeq.jl:87 [inlined]
 anonymous at ./<missing>:?
 NLExpr(::NLOptControl.NLOpt, ::Expr, ::Array{Any,2}, ::Vararg{Any,N} where N) at /Users/roi/.julia/v0.6/NLOptControl/src/diffeq.jl:89
 DiffEq(::NLOptControl.NLOpt, ::Array{Any,2}, ::Array{Any,2}, ::Int64, ::Int64) at /Users/roi/.julia/v0.6/NLOptControl/src/diffeq.jl:36
 OCPdef!(::NLOptControl.NLOpt) at /Users/roi/.julia/v0.6/NLOptControl/src/setup.jl:305
 #configure!#52(::Array{Any,1}, ::Function, ::NLOptControl.NLOpt) at /Users/roi/.julia/v0.6/NLOptControl/src/setup.jl:474
 (::NLOptControl.#kw##configure!)(::Array{Any,1}, ::NLOptControl.#configure!, ::NLOptControl.NLOpt) at ./<missing>:0
 anonymous at ./<missing>:?
``````

So I am not sure how to implement this.

This post is meant to provide some summary about the things that have helped so far.

Indeed, the error of the forward differentiation I got was fixed by transforming things to abstract type. As @ChrisRackauckas wrote:

this helped me in changing the type of the initial conditions of the ODE solver from `Float64` to abstract type.

Then I had to do more changes in the code to make sure everything gets abstract types.

Now, when I run `optimize` I don’t get the error, but the differentiation takes a long time (the code is running for more than 15 hours, without any output so far).

In the meantime I am trying to solve the problem using `Optim.jl` and `BlackBoxOptim.jl`. See this post for more information.

There is still a chance that `NLOptControl.jl` will be helpful, but I am not sure yet about it.

# Update - automatic differentiation (ForwardDiff.jl) returns `Nan`

My optimization function needs to calculate elliptic integrals. After converting the code in Elliptic.jl to be able to work with automatic differentiation, I get that the derivative of the complete elliptic integral of the 1st kind is `NaN`. The derivative is actually given by other elliptic integrals (as can be seen in the last link at the end of the section).

Here’s the minimal code (note that the credit to the code of calculating the elliptic integral belongs to nolta’s package Elliptic.jl):

``````julia> function DRF(X::Real, Y::Real, Z::Real)
D1MACH1 = typeof(X)(1e-308)
D1MACH2 = typeof(X)(1e308)
D1MACH3 = typeof(X)(5e-16)/typeof(X)(2.0)
D1MACH4 = typeof(X)(5e-16)
D1MACH5 = log10(typeof(X)(2.))

ERRTOL = (typeof(X)(4.0)*D1MACH3)^typeof(X)(1.0/6.0)
LOLIM  = typeof(X)(5.0) * D1MACH1
UPLIM  = D1MACH2/typeof(X)(5.0)
C1 = typeof(X)(1.0/24.0)
C2 = typeof(X)(3.0/44.0)
C3 = typeof(X)(1.0/14.0)

ans = zero(X)
if min(X,Y,Z) < zero(X)
return ans, one(X)
end

if max(X,Y,Z) > UPLIM
return ans, typeof(X)(3)
end

if min(X+Y,X+Z,Y+Z) < LOLIM
return ans, typeof(X)(2.0)
end

XN = X
YN = Y
ZN = Z
MU = zero(X)
XNDEV = zero(X)
YNDEV = zero(X)
ZNDEV = zero(X)

while true
MU = (XN+YN+ZN)/typeof(X)(3.0)
XNDEV = typeof(X)(2.0) - (MU+XN)/MU
YNDEV = typeof(X)(2.0) - (MU+YN)/MU
ZNDEV = typeof(X)(2.0) - (MU+ZN)/MU
EPSLON = max(abs(XNDEV),abs(YNDEV),abs(ZNDEV))
if (EPSLON < ERRTOL) break end
XNROOT = sqrt(XN)
YNROOT = sqrt(YN)
ZNROOT = sqrt(ZN)
LAMDA = XNROOT*(YNROOT+ZNROOT) + YNROOT*ZNROOT
XN = (XN+LAMDA)*typeof(X)(0.250)
YN = (YN+LAMDA)*typeof(X)(0.250)
ZN = (ZN+LAMDA)*typeof(X)(0.250)
end

E2 = XNDEV*YNDEV - ZNDEV*ZNDEV
E3 = XNDEV*YNDEV*ZNDEV
S  = typeof(X)(1.0) + (C1*E2-typeof(X)(0.10)-C2*E3)*E2 + C3*E3
ans = S/sqrt(MU)

return ans, typeof(X)(0)
end
DRF (generic function with 1 method)

julia> function Kdiff(m::Real)
if m < one(m)
drf,ierr = DRF(zero(m), one(m) - m, one(m))
@assert ierr == zero(m)
return drf
elseif m == one(m)
return Inf
elseif isnan(m)
return NaN
else
throw(DomainError(m, "argument m not <= 1"))
end
end
Kdiff (generic function with 1 method)

julia> Kdiff(0.5)
1.854074677301372

julia> using ForwardDiff

julia> ForwardDiff.derivative(Kdiff, 0.5)
NaN

``````

Is there a way to fix this `Nan` issue? or to provide somehow the derivative only to this part of the code? (I mean that my fitness function needs to calculate the elliptic integral at some point, but I do not know the derivatives of other parts of the fitness function).

In the meantime, I gave this fitness function to `JuMP` and it runs. I get some results which I try to compare to other algorithm, and they seem for now to make sense.
But, shouldn’t it be a problem that the derivative is `Nan`? In what way does it affect the solver (I am using `Ipopt`).