Example of a simple (very routine) problem in economics

[azev77]

Hi & thank you for this package!
I'd like to use it to solve a very simple (routine) problem in economics.

This is a consumption-saving problem.
The household lives during time: t ∈ [0, 10]
The household's wealth @ time t is: B(t)
The household's choice is consumption: c(t)
The household is endowed w/ $100: B(0) = 100
The household cannot die in debt: B(10) = 0

I have the closed form solution for comparison.
Is it possible to solve this problem w/ your package?

[pulsipher]

Hi there, thanks for reaching out. The principle limitation of InfiniteOpt (in its current form) is support for non-quadratic nonlinear functions of decision variables (e.g., log(c(t))). This is due to our dependence on JuMP which does allow general nonlinear functions/expressions when modeling directly in JuMP but this nonlinear interface is unextendible for packages like ours. This limitation will be lifted once jump-dev/MathOptInterface.jl#846 is addressed.
Until then, our currently accepted expression types are described here.

There are plans to make a temporary patch to address this problem (#16). I hope to address this sometime this summer, but we'll see what my schedule allows.

It might be possible to use a change of variables to adapt your formulation. Otherwise, you can transcribe the problem manually into a traditional NLP form that can be used by JuMP.

[odow]

something like:

using JuMP, Ipopt
Δt = 0.01
T = 10
time = 0.0:Δt:T
model = Model(Ipopt.Optimizer)
@variable(model, c[time] >= 0.0001)  # To avoid log(0).
@variable(model, B[time])
fix(B[0.0], 100.0)
fix(B[T], 0.0)
@NLobjective(model, Max, sum(exp(-0.01 * t) * log(c[t]) for t in time))
@constraint(
    model,
    [i = 2:length(T)],
    B[t] == B[t - 1] + Δt * (0.05 * B[t - 1] - c[t]),
)
optimize!(model)

[azev77]

@pulsipher
I repeated it w/ quadratic utility (-1)*(c - 100)^2 and it works beautifully except for at the very corners:

using InfiniteOpt, JuMP, Ipopt, Plots
m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t in [0, 10],num_supports = 10_001)
@infinite_variable(m, B(t))
@infinite_variable(m, 0 <= c(t) <= 99)  #m, c(t) >= 0
@objective(m, Max, integral( (-1)*(c - 100)^2, t) )
@BDconstraint(m, (t == 0), B == 100)      # B[t=0] = 100
@BDconstraint(m, (t == 10), B == 0)
@constraint(m, c1, @deriv(B, t) == 0.01*B - c)

m
optimize!(m)
termination_status(m)
opt_obj = objective_value(m) #
#c, B, value.(c), value.(B)
c_opt, B_opt = value.(c), value.(B)
u_ts = supports.(c)
u_ts =[u_ts[j][1] for j in eachindex(u_ts)]

plot(u_ts,   B_opt, lab = "B: wealth balance")
plot!(u_ts,  c_opt, lab = "c: consumption")

Here are the control/state vars from your package:

From Mathematica:

1. This can be a nice example of a (very) simple econ problem for your package!
2. Is it yet possible to allow exponentially discounted objective? exp(-0.01*t)*(-1)*(c(t) - 100)^2

[azev77]

@odow thanks!
I tried your method applied to the quadratic example (above) & it doesn't seem to work:

using JuMP, Ipopt, Plots
Δt = 0.01; T = 10.0; time = 0.0:Δt:T;
m = Model(Ipopt.Optimizer)
@variable(m, 0.0001 <= c[time] <= 99.0) 
@variable(m, B[time])
fix(B[0.0], 100.0)
fix(B[T], 0.0)
@NLobjective(m, Max, sum((-1.0)*(c[t] - 100.0)^2.0 for t in time) )
@constraint(m, 
    [i = 2:length(T)], 
    B[t] == B[t - 1] + Δt * (0.01 * B[t - 1] - c[t])
    )
optimize!(m)
objective_value(m)
c_opt, B_opt = value.(c), value.(B)
plot(time,   B_opt.data, lab = "B: wealth balance")
plot!(time,  c_opt.data, lab = "c: consumption")

Do you have any idea what it might be?

[odow]

The perils of writing code without checking that it runs:

Δt = 0.01; T = 10.0; time = 0.0:Δt:T;
m = Model(Ipopt.Optimizer)
@variable(m, 0.0001 <= c[time] <= 99.0) 
@variable(m, B[time])
fix(B[0.0], 100.0)
fix(B[T], 0.0)
@NLobjective(m, Max, sum((-1.0)*(c[t] - 100.0)^2.0 for t in time) )
for i in 2:length(time)
    t_1, t = time[i - 1], time[i]
    @constraint(m, B[t] == B[t_1] + Δt * (0.01 * B[t_1] - c[t]))
end
optimize!(m)

[azev77]

@odow works beautifully. Thank you!
Do you know why there is an initial jump in the control variable @ first index:

[pulsipher]

@azev77 That's a cool example, nice job setting it up.

Yes you can add in the discount factor since it is only a function of the infinite parameter t which is not a decision variable. There are 2 ways to accomplish this.

Option 1: use a parameter function and embed it in the expression:

using InfiniteOpt, JuMP, Ipopt, Plots

discount(t) = exp(-0.01 * t) # this will work as a function for our parameter function

m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t in [0, 10], num_supports = 100)
@infinite_variable(m, B(t))
@infinite_variable(m, 0 <= c(t) <= 99)
pfunc = parameter_function(discount, t)
@objective(m, Min, integral(pfunc * (c - 100)^2, t))
@BDconstraint(m, (t == 0), B == 100)
@BDconstraint(m, (t == 10), B == 0)
@constraint(m, c1, deriv(B, t) == 0.01*B - c)
optimize!(m)
termination_status(m)
opt_obj = objective_value(m) #
c_opt, B_opt = value(c), value(B)
ts = supports(t)
plot(ts,   B_opt, lab = "B: wealth balance")
plot!(ts,  c_opt, lab = "c: consumption")

EDIT: See correction in the comments below.

Option 2: use the weight_func keyword in integral:

using InfiniteOpt, JuMP, Ipopt, Plots

discount(t) = exp(-0.01 * t) # this will work as a weight function in the integral

m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t in [0, 10], num_supports = 100)
@infinite_variable(m, B(t))
@infinite_variable(m, 0 <= c(t) <= 99)
@objective(m, Min, integral((c - 100)^2, t, weight_func = discount))
@BDconstraint(m, (t == 0), B == 100)
@BDconstraint(m, (t == 10), B == 0)
@constraint(m, c1, deriv(B, t) == 0.01*B - c)
optimize!(m)
termination_status(m)
opt_obj = objective_value(m) #
c_opt, B_opt = value(c), value(B)
ts = supports(t)
plot(ts,   B_opt, lab = "B: wealth balance")
plot!(ts,  c_opt, lab = "c: consumption")

[pulsipher]

Oh, my bad that temporarily introduces an expression with 3 variable reference multiplied together. This can be alleviated by introducing a placeholder variable:

using InfiniteOpt, JuMP, Ipopt, Plots

discount(t) = exp(-0.01 * t) # this will work as a function for our parameter function

m = InfiniteModel(Ipopt.Optimizer)
@infinite_parameter(m, t in [0, 10], num_supports = 100)
@infinite_variable(m, B(t))
@infinite_variable(m, 0 <= c(t) <= 99)
@infinite_variable(m, placeholder(t))
pfunc = parameter_function(discount, t)
@constraint(m, placeholder == (c - 100)^2)
@objective(m, Min, integral(pfunc * placeholder, t))
@BDconstraint(m, (t == 0), B == 100)
@BDconstraint(m, (t == 10), B == 0)
@constraint(m, c1, deriv(B, t) == 0.01*B - c)
optimize!(m)
termination_status(m)
opt_obj = objective_value(m) #
c_opt, B_opt = value(c), value(B)
ts = supports(t)
plot(ts,   B_opt, lab = "B: wealth balance")
plot!(ts,  c_opt, lab = "c: consumption")

Alternatively, option 2 provides a more compact avenue.

[azev77]

@odow is there a way to do this in JuMP.jl with:

1. derivative in the constraint instead of @constraint(m, B[t] == B[t_1] + Δt * (0.01 * B[t_1] - c[t]))
2. integral in the objective instead of @NLobjective(m, Max, sum((-1.0)*(c[t] - 100.0)^2.0 for t in time) )

[odow]

Do you know why there is an initial jump in the control variable @ first index:

c[0.0] isn't used, so it just picks a random value. You could fix it to zero or just plot from the second element onwards.

is there a way to do this in JuMP.jl

No. @pulsipher has done an excellent job adding these extensions to InfiniteOpt instead...

[azev77]

@pulsipher can I help you add this econ example to the docs?

[pulsipher]

@azev77 Ya that would be great. I think this would make a nice addition to the Examples page in the docs. Please refer the contribution guide (https://github.com/pulsipher/InfiniteOpt.jl/blob/master/CONTRIBUTING.md) for info on how to implement this and make a pull request. You can follow-up here with any questions you have.