Simple Firm Investment Problem giving errors

[azev77]

using InfiniteOpt, Ipopt, Plots
δ=0.10; θ=2.00; r=0.15;
z= r + δ -.05;

T = 120.0   # horizon
k0 = 1.0 # capital endowment
u(i,k; z=z,θ=θ,δ=δ) = z*k -i -(θ/2)*k*(i/k - δ)^2 # dividend = return function
discount(t; r=r) = exp(-r*t) # discount function
BC(k, i; δ=δ) = i - δ*k        # LOM 
opt = Ipopt.Optimizer   # desired solver
ns = 1_000;             # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables

@objective(m, Max, integral(u(k,i), t, weight_func = discount))
ERROR: LoadError: /(::GeneralVariableRef,::GeneralVariableRef) is not defined. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective.
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:33
 [2] /(#unused#::GeneralVariableRef, #unused#::GeneralVariableRef)
   @ JuMP C:\Users\azevelev\.julia\packages\JuMP\b3hGi\src\operators.jl:418
 [3] u(k::GeneralVariableRef, i::GeneralVariableRef; z::Float64, θ::Float64, δ::Float64)
   @ Main c:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\4prob\3a_Firm_Invest\Investment.jl:417
 [4] u(k::GeneralVariableRef, i::GeneralVariableRef)
   @ Main c:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\4prob\3a_Firm_Invest\Investment.jl:417
 [5] macro expansion
   @ C:\Users\azevelev\.julia\packages\MutableArithmetics\8xkW3\src\rewrite.jl:279 [inlined]
 [6] macro expansion
   @ C:\Users\azevelev\.julia\packages\JuMP\b3hGi\src\macros.jl:1260 [inlined]
 [7] top-level scope
   @ c:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\4prob\3a_Firm_Invest\Investment.jl:432
in expression starting at c:\Users\azevelev\Dropbox\Computation\Julia\DE_OC_VFI\4prob\3a_Firm_Invest\Investment.jl:432

Also, can this package solve infinite horizon problems yet?

[pulsipher]

What version of InfiniteOpt are you running? You'll need to update to the latest version for nonlinear expressions to work.

Infinite horizon problems are still on the todo list, I put most of my development time recently toward adding general nonlinear support. Admittedly, I don't foresee myself having much time in the immediate future for adding full support for infinite-horizon problems, so we'll probably need someone else to lead that addition. We can currently represent infinite-horizon problems, we just need to add a transformation/solution approach to solve them in a general manner. I opened #190 to track this feature request.

[azev77]

Just updated. Still get errors:

using InfiniteOpt, Ipopt, Plots
p=1.00; δ=0.10; θ=2.00; r=0.15;
z= r + δ -.02   
T = 120.0   # life horizon
k0 = 1.0 # endowment
u(k,i;z=z,θ=θ,δ=δ) = z*k -i -(θ/2)*k*(i/k - δ)^2 # utility function
discount(t; r=r) = exp(-r*t) # discount function
BC(k, i; δ=δ) = i - δ*k        # LOM 

opt = Ipopt.Optimizer   # desired solver
ns = 1_000;             # number of gridpoints

m = InfiniteModel(opt)
@infinite_parameter(m, t ∈ [0, T], num_supports = ns)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@objective(m, Max, integral(u(k,i), t, weight_func = discount))
@constraint(m, c1, deriv(k, t) == BC(k, i; δ=δ))
@constraint(m, k == k0, DomainRestrictions(t => 0))
@constraint(m, k == 0, DomainRestrictions(t => T))

optimize!(m)

This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.

Number of nonzeros in equality constraint Jacobian...:     5999
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     3000


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          = 0
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 0.011

EXIT: Invalid number in NLP function or derivative detected.

Btw, I'm not sure if this is part of the issue, but there should be a salvage value at T, equal to e^(-r*T) * (1-\delta)K_T

[pulsipher]

It looks like the update resolved the nonlinear error and now InfiniteOpt is able to formulate the problem as you specify it.

Now Ipopt is throwing an error since the objective function is not well posed numerically since it divides by zero if k(t) has any points that approach 0. This is what occurs and induces the EXIT: Invalid number in NLP function or derivative detected. error from Ipopt. It appears that constraining k(t) to be strictly positive makes the problem infeasible. Adding the salvage value is a constant (it has not state/control variables in it), so adding it as an extra term to the objective wouldn't affect the solution.

Hence, the objective function will need to be reformulated into a better posed function. We can accomplish this by defining an auxiliary variable a == i / k - δ which we then pose as a quadratic constraint k * (a + δ) == i. Then we can use a when defining the objective function:

using InfiniteOpt, Ipopt
p = 1.00; δ = 0.10; θ = 2.00; r = 0.15
z = r + δ -.02   
T = 120.0   # life horizon
k0 = 1.0 # endowment

m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t ∈ [0, T], num_supports = 1000)
@variable(m, k, Infinite(t)) ## state variables
@variable(m, i, Infinite(t)) ## control variables
@variable(m, a, Infinite(t))
@constraint(m, k * (a + δ) == i) # defines a == i / k - δ
@objective(m, Max, ∫(z * k - i - (θ / 2) * k * a^2, t, weight_func = t -> exp(-r * t)))
@constraint(m, c1, ∂(k, t) == i - δ * k)
@constraint(m, k(0) == k0)
@constraint(m, k(T) == 0)

optimize!(m)

This is able to then solve on my machine.

Since this doesn't constitute a bug in InfiniteOpt.jl, I please ask that further discussion on how to formulate this problem be done via the discussion forum (https://github.com/pulsipher/InfiniteOpt.jl/discussions).

[pulsipher]

This discussion is continued in #194.