# Solving an optimal control problem with JuMP

**URL:** https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570
**Category:** Optimization (Mathematical)
**Tags:** question, jump, optimization, forwarddiff
**Created:** [March 26, 2020, 10:18pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570 "2020-03-26T22:18:43Z")
**Posts on this page:** 1
**Showing post:** 17

<div class="post-metadata">

### Author: ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [March 28, 2020, 4:50pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/17 "2020-03-28T16:50:55Z")

</div>

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:

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

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](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-optim-jl/36635) 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](https://github.com/JuliaDiff/ForwardDiff.jl)) returns `Nan`

My optimization function needs to calculate [elliptic integrals](https://en.wikipedia.org/wiki/Elliptic_integral). After converting the code in [Elliptic.jl](https://github.com/nolta/Elliptic.jl) to be able to work with automatic differentiation, I get that the derivative of the [complete elliptic integral of the 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_first_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](https://github.com/nolta/Elliptic.jl)):

```julia-auto
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`).

---

_[View the full topic](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570)._