Nonlinear objective

[azev77]

Hi, I'm trying to test out the nonlinear capabilities.
I repeated the example in the docs (except w/ log(c) as the objective)

using InfiniteOpt, Ipopt
ρ = 0.025  # discount rate
T = 10.0   # life horizon
r = 0.05   # interest rate
B0 = 100.0 # endowment
#u(c) = -(c - 100)^2       # utility function
u(c) = log(c)                    # utility function

discount(t; ρ=ρ) = exp(-ρ*t) # discount function
BC(B, c; r=r) = r*B - c      # budget constraint

opt = Ipopt.Optimizer   # desired solver
ns = 1_000;             # number of points in the time grid

m = InfiniteModel(opt)
@infinite_parameter(m, t in [0, T], num_supports = ns)
@variable(m, B, Infinite(t)) ## state variables
@variable(m, c, Infinite(t)) ## control variables

@objective(m, Max, integral(u(c), t, weight_func = discount))

@constraint(m, B(0) == B0)
@constraint(m, B(T) == 0)
@constraint(m, c1, deriv(B, t) == BC(B, c; r=r))
optimize!(m)

The message:

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.............:     1000


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

EXIT: Invalid number in NLP function or derivative detected.

[notesofdabbler]

I think the default initial guess is 0 for ipopt which likely causes invalid number for log(c). If you use log(c + 1.0e-5), then the code works.

[pulsipher]

@notesofdabbler is correct, the default guess for c invokes an invalid evaluation for Ipopt (i.e., log(0)). We can correct this by simply specifying the initial guess value of c:

@variable(m, c, Infinite(t), start = 1)

This modification resolves the problem.

[azev77]

@notesofdabbler @pulsipher Thank you!
The current example works using @variable(m, c, Infinite(t), start = 1)

1. the prev error message didn't give me any guidance. Can this be made more explicit either in InfiniteOpt.jl/JuMP.jl/ or elsewhere? @odow
2. I get a strange instability for the objective u(c) = ( c^(1- γ) ) / (1- γ) when the parameter γ=2

using InfiniteOpt, Ipopt
ρ = 0.025  # discount rate
T = 10.0   # life horizon
r = 0.05   # interest rate
B0 = 100.0 # endowment

function f(γ)
    u(c) = ( c^(1- γ) ) / (1- γ)
    discount(t; ρ=ρ) = exp(-ρ*t) # discount function
    BC(B, c; r=r) = r*B - c      # budget constraint

    opt = Ipopt.Optimizer   # desired solver
    ns = 1_000;             # number of points in the time grid

    m = InfiniteModel(opt)
    @infinite_parameter(m, t in [0, T], num_supports = ns)
    @variable(m, B, Infinite(t)) ## state variables
    @variable(m, c, Infinite(t), start = 1) ## control variables. Initial guess c(t)=1. Ipopt default initial guess c(t)=1.

    @objective(m, Max, integral(u(c), t, weight_func = discount))

    @constraint(m, B(0) == B0)
    @constraint(m, B(T) == 0)
    @constraint(m, c1, deriv(B, t) == BC(B, c; r=r))
    optimize!(m)
    termination_status(m)
    c_opt = value(c)
    B_opt = value(B)
    ts = supports(t)
    opt_obj = objective_value(m) # V(B0, 0)

    ix = 2:(length(ts)-1) # index for plotting
    plot!(ts[ix],   B_opt[ix], lab = "B: wealth balance")
    plot!(ts[ix],  c_opt[ix], lab = "c: consumption")
end

plot()
f(2.0 - .1)  # correct
f(2.0 + .1) # correct
f(2.0 + 0) # wrong

[pulsipher]

1. the prev error message didn't give me any guidance. Can this be made more explicit either in InfiniteOpt.jl/JuMP.jl/ or elsewhere?

I agree that a more useful error message would be good. The tricky thing is that InfiniteOpt.jl and JuMP.jl don't check for feasibility and/or numerical instability of a model before passing it on to the solver, since this is nontrivial to implement in general and would likely add some overhead. Hence, the error isn't triggered until Ipopt is called. Perhaps we could do something there...

2. I get a strange instability for the objective u(c) = ( c^(1- γ) ) / (1- γ) when the parameter γ=2

It is not immediately obvious to me why this is, I'll need to think about it so more.