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[1] = u[2] # θ'(t) = ω(t)
du[2] = -3g/(2l)*sin(u[1]) + 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)
plot(sol,linewidth=2,xaxis="t",label=["θ [rad]" "ω [rad/s]"],layout=(2,1))
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:
[1] 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[1] = u[2] # θ'(t) = ω(t)
du[2] = -3g/(2l)*sin(u[1]) + 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).