So I tried to make a minimum example to ask questions based on a more complicated piece of code I have written:
- A HUGE common error I'm getting is expecting float64 and instead got ForwardDiff.Dual - can someone give me a tip how in general I always make sure I avoid this bug. I feel like every time I do a new optimization problem I have to reinvent the wheel to try to make this go away
- Apparently you cannot autodiff the julia exp() function? Does anyone know how to make it work?
- A workaround is I did a finite sum to approximate it via the taylor series. In my one function if I had 20 terms the autodiff worked, but it wasn't accurate enough - so I went to 40 terms but then julia told me to do factorial(big(k)) and then when I try to do that with autodiff it doesn't work now - anyone have a fix for this?
Any advice would be greatly appreciated!
using Cubature
using Juniper
using Ipopt
using JuMP
using LinearAlgebra
using Base.Threads
using Cbc
using DifferentialEquations
using Trapz
function mat_exp(x::AbstractVector{T},dim,num_terms,A) where T
sum = zeros(Complex{T},(dim,dim))
A[1,1] = A[1,1]*x[1]
A[2,2] = A[2,2]*x[2]
return exp(A)-1
end
function exp_approx_no_big(x::AbstractVector{T},dim,num_terms,A) where T
sum = zeros(Complex{T},(dim,dim))
A[1,1] = A[1,1]*x[1]
A[2,2] = A[2,2]*x[2]
for k=0:num_terms-1
sum = sum + (1.0/factorial(k))*A^k
end
return norm(sum)-1
end
function exp_approx_big(x::AbstractVector{T},dim,num_terms,A) where T
sum = zeros(Complex{T},(dim,dim))
A[1,1] = A[1,1]*x[1]
A[2,2] = A[2,2]*x[2]
for k=0:num_terms-1
sum = sum + (1.0/factorial(big(k)))*A^k
end
return norm(sum)-1
end
optimizer = Juniper.Optimizer
nl_solver= optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
mip_solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0, "threads"=>nthreads())
m = Model(optimizer_with_attributes(optimizer, "nl_solver"=>nl_solver, "mip_solver"=>mip_solver))
@variable(m, 0.0<=x[1:2]<=1.0)
dim=5
A=zeros(Complex,(dim,dim))
for k=1:dim
A[k,k]=1.0
end
println(A)
f(x...) = exp_approx_no_big(collect(x),dim,20,A)
g(x...) = exp_approx_big(collect(x),dim,40,A)
h(x...) = mat_exp(collect(x),dim,20,A)
register(m, :f, 2, f; autodiff = true)
@NLobjective(m, Min, f(x...))
optimize!(m)
println(JuMP.value.(x))
println(JuMP.objective_value(m))
println(JuMP.termination_status(m))
