Ipopt and JuMP - Invalid number in NLP function or derivative

Specific Domains Optimization (Mathematical)
,
1

We are generalising the following julia code

to two or more dimensional grids.

The following is the output (for each grid point):

"Number of Iterations…: 0

Number of objective function evaluations = 0
Number of objective gradient evaluations = 1
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 1
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 0
Total CPU secs in IPOPT (w/o function evaluations) = 0.026
Total CPU secs in NLP function evaluations = 0.000

EXIT: Invalid number in NLP function or derivative detected."

The following is the full working example:

using MathOptInterface, LinearAlgebra, Statistics, Plots, QuantEcon, Interpolations, NLsolve, Optim, Random, IterTools, JuMP, Ipopt, MathOptInterface;

u(x) = sum(log(x[i]) for i in 1:numSectors)
w(x) = sum(5*log(x[i]) for i in 1:numSectors)

β = 0.96
α = 0.4
f(x) = x^α

numSectors = 2

numPoints1D = 2

grid = Vector{Vector{Float64}}(undef,(numPoints1D)^numSectors)

gridMax = 2
gridMin = 1

iter=1
for p in product(LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D))
    grid[iter] = collect(p)
    global iter += 1
end

wVal = w.(grid)
   
function T(wVal, grid, β, f ; compute_policy = false)

    wVals = zeros(numPoints1D,numPoints1D);
    for j in numPoints1D
        for i in numPoints1D
            wVals[i,j] = wVal[i+(i-1)*j]
        end
    end
    function wFunc(x,y)
        return interpolate((LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D)), wVals, Gridded(Linear()))(x,y)
    end
    Tw = zeros(length(grid))
    σ = similar(grid)
    for n in 1:length(grid)
        y = grid[n]
        modTrial = Model(Ipopt.Optimizer);
        @variable(modTrial,  c[1:numSectors])
        @variable(modTrial, k[1:numSectors], start=1.5)
        for i in 1:numSectors
            @constraint(modTrial, 0.99*gridMin <= c[i] <= 0.999*y[i])
	    @constraint(modTrial, k[i] == y[i] - c[i])
        end
	register(modTrial, :wFunc, 2, wFunc; autodiff = true)
	register(modTrial, :f, 1, f; autodiff = true)
	@NLobjective(modTrial, Max, sum(log(c[i]) for i in 1:numSectors) + β*wFunc(f(k[1]),f(k[2])))
        optimize!(modTrial)
        Tw[n] = JuMP.objective_value(modTrial)
        if compute_policy
            σ[n] = value.(c)
            return Tw, σ
        end
    end
    return Tw 
end
wVal = T(wVal, grid, β, f)

Would anyone be able to help with what the problem is here, and a way to possibly fix it?
Thanks

2

log(c[i])

Your c variable has free bounds, but you’re calling log on it, which isn’t defined for c <= 0. Try adding

@variable(modTrial,  c[1:numSectors] >= 0.0001)
1 Like
3

Thanks @odow , that makes good sense. We are now, however getting the following error in relation to our call of the interpolate function:

“ERROR: LoadError: BoundsError: attempt to access 2×2 interpolate((range(1.0, stop=2.0, length=2),range(1.0,
stop=2.0, length=2)), ::Matrix{Float64}, Gridded(Linear())) with element type Float64 at index [0.18688702907494092, 0.18688702907494092]”

Not sure why it is trying to access the index when we want that pair of numbers as inputs to the function.

Our function wFunc(x,y) behaves well in the following minimal working example, but not within the optimisation problem.

julia> testGrid = [1.0,2.0]
julia> testGridA = [1.0,2.0]
julia> testRand = [1.0 2.0; 3.0 4.0]
julia> interpolate((testGrid,testGridA), testRand, Gridded(Linear()))
2×2 interpolate((::Vector{Float64},::Vector{Float64}), ::Matrix{Float64}, Gridded(Linear())) with element type Float64:
 1.0  2.0
 3.0  4.0

julia> f(x) = x^0.4
f (generic function with 1 method)

julia> interpolate((testGrid,testGridA), testRand, Gridded(Linear()))(f(1.0),f(1.5))
1.1760790225246736
4

You have more domain issues:

  • k is a free variable, so f(k[1]) isn’t defined for k < 0
  • For small k, f(k) < 1, and the point [0.18688702907494092, 0.18688702907494092] is outside your interpolation domain.
1 Like
5

@odow Thankyou so much, that’s running now.

1 Like