0

我想使用涉及贝塞尔函数的 JuMP 最小化目标函数。完整的函数相当复杂,我没有这个函数的导数。似乎自动推导期间的 JuMP 使用等于 NaN 的参数调用此函数,这会导致函数失败。

这是完整的示例:

using JuMP
using SpecialFunctions
using Ipopt

xdata = linspace(0.1,5,100);

data = 3 * besselk.(2,2 * xdata)

function fit4(p1,p2)
    @show p1,p2
    return sum((data - p1 * besselk.(2,p2 * xdata)).^2)
end

# is 0 as it should
@show fit4(3,2)

m = Model(solver=IpoptSolver())
@variable(m, 0.01 <= x[1:2] <= 5)
JuMP.register(m,:fit4,2,fit4; autodiff = true)

@NLobjective(m, Min, fit4(x[1],x[2]))
status = solve(m)

它产生一个错误:

 ERROR: AmosException with id 2: overflow.
 Stacktrace:
  [1] _besselk(::Float64, ::Complex{Float64}, ::Int32) at /home/abarth/.julia/v0.6/SpecialFunctions/src/bessel.jl:239
  [2] besselk(::Int64, ::Float64) at /home/abarth/.julia/v0.6/SpecialFunctions/src/bessel.jl:450
  [3] besselk at /home/abarth/.julia/v0.6/ForwardDiff/src/dual.jl:202 [inlined]
  [4] (::##8#10)(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}) at ./<missing>:0
  [5] macro expansion at ./broadcast.jl:155 [inlined]
  [6] macro expansion at ./simdloop.jl:73 [inlined]
  [7] macro expansion at ./broadcast.jl:149 [inlined]
  [8] _broadcast!(::##8#10, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2},1}, ::Tuple{Tuple{Bool}}, ::Tuple{Tuple{Int64}}, ::StepRangeLen{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}}}, ::Tuple{}, ::Type{Val{0}}, ::CartesianRange{CartesianIndex{1}}) at ./broadcast.jl:141
  [9] broadcast_t(::Function, ::Type{T} where T, ::Tuple{Base.OneTo{Int64}}, ::CartesianRange{CartesianIndex{1}}, ::StepRangeLen{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}}}) at ./broadcast.jl:270
[...]

besselk 似乎是 Fortran Amos 库中的一个函数。JuMP(和 ForewardDiff 等人)是否要求所有函数都在纯 Julia 中定义?

4

0 回答 0