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