# 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:** 17\
**Page:** 1

<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 26, 2020, 10:18pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/1 "2020-03-26T22:18:43Z")

</div>

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](https://github.com/JuliaMPC/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](https://docs.sciml.ai/stable/tutorials/ode_example/#Example-3:-Solving-Nonhomogeneous-Equations-using-Parameterized-Functions-1). 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)`

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

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

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

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

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

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [March 27, 2020, 12:45am UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/2 "2020-03-27T00:45:22Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 27, 2020, 1:00am UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/3 "2020-03-27T01:00:43Z")

</div>

> [@roi.holtzman](#):
>
> u0 = [0.0, 1.0]

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

---

<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 27, 2020, 7:53am UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/4 "2020-03-27T07:53:58Z")

</div>

> [@ChrisRackauckas](#):
>
> 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](https://docs.sciml.ai/stable/basics/faq/index.html#Autodifferentiation-and-Dual-Numbers-1) 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](https://github.com/nolta/Elliptic.jl) package. Say I have a function

```julia
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[1], params[1] + params[2])
  return res
end

```

And this function is placed in the optimization function

```julia
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?

---

<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 27, 2020, 8:11am UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/5 "2020-03-27T08:11:52Z")

</div>

> [@Albert\_Zevelev](#):
>
> Is the reason you cant use NLOptControl.jl, that they haven’t updated to Julia 1.0+?

I have installed Julia 0.6 in order to use [NLOptControl.jl](https://github.com/JuliaMPC/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](https://nbviewer.jupyter.org/github/JuliaOpt/JuMPTutorials.jl/blob/master/notebook/modelling/rocket_control.ipynb)). 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 27, 2020, 2:28pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/6 "2020-03-27T14:28:05Z")

</div>

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

---

<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 27, 2020, 2:35pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/7 "2020-03-27T14:35:09Z")

</div>

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

The error I get is

```julia
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:
 [1] 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
 [2] θ_elliptic_opt(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#107#109"{typeof(opt_func_g_arr_elliptic)},Float64},Float64,4}, ::Float64, ::Float64) at ./In[7]:3

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 27, 2020, 2:38pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/8 "2020-03-27T14:38:39Z")

</div>

> [@roi.holtzman](#):
>
> I am not sure exactly what it means to be AD-differentiable.
> 
> The error I get is

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

> <https://github.com/nolta/Elliptic.jl/blob/master/src/jacobi.jl#L52-L53>

so that would need to get opened up / fixed.

---

<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 27, 2020, 3:27pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/9 "2020-03-27T15:27:17Z")

</div>

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

E.g. instead of this

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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 27, 2020, 3:37pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/10 "2020-03-27T15:37:24Z")

</div>

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

---

<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 27, 2020, 3:59pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/11 "2020-03-27T15:59:29Z")

</div>

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

When I do the following

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

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [March 27, 2020, 4:18pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/12 "2020-03-27T16:18:55Z")

</div>

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

---

<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 27, 2020, 4:24pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/13 "2020-03-27T16:24:34Z")

</div>

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?

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

```

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [March 27, 2020, 4:27pm UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/14 "2020-03-27T16:27:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [March 28, 2020, 4:48am UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/15 "2020-03-28T04:48:09Z")

</div>

> [@roi.holtzman](#):
>
> writing an equation of the form `x2[j] * sin(t)` is impossible.

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

---

<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, 10:11am UTC](https://discourse.julialang.org/t/solving-an-optimal-control-problem-with-jump/36570/16 "2020-03-28T10:11:18Z")

</div>

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

```julia
dx[5] = :(1.0)

```

But I get the following error:

```julia
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:
 [1] setindex!(::Array{Expr,1}, ::Float64, ::Int64) at ./array.jl:578
 [2] include_string(::String, ::String) at ./loading.jl:522

```

I tried e.g. converting to an Abstract type

```julia
dx[5] = :(convert(Real, 1.0))

```

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

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

Unrecognized function "convert" used in nonlinear expression.

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

```

So I am not sure how to implement this.

---

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