Hello, I am trying to numerically solve the dynamic problem of optimal natural resource depletion (aka “eating the cake”) with InfiniteOpt.jl.
The script below works fine, but for the time horizon (i.e. the upper bound of the domain of the infinite dimension “time”) I need to specify it from the analytical resolution that I already have.
Is there a way to let the time horizon “free”, i.e. let it to be endogenously retrieved in the optimisation ?
[Side question] Is is possible to use an objective version that is itself an integral computed with
welfare_auto ?). I am having errors with it, so I had to compute it manually…
# "Eating the cake" model of depletion of natural resources cd(@__DIR__) using Pkg Pkg.activate(".") # Pkg.add(["QuadGK", "InfiniteOpt", "Ipopt", "Plots"]) # only once to install packages using QuadGK, InfiniteOpt, Ipopt, Plots # Parameters k=100; a = 10 # Parameters of the demand function S₀ = 1000 # Initial stock ρ = 0.01 # Discount rate Sτ = 0 # Final stocks test = sqrt(2*S₀*a/ρ) # this comes from the analytical resolution that I already have !! T = test # time horizon # Demand function at each point in time demand_price(r;k=k,a=a) = k*exp(-a*r) # OBJECTIVE FUNCTION # Welfare (utility) at each point in time welfare_auto(r;k=k,a=a) = quadgk(q -> demand_price(q,k=k,a=a), 0, r, rtol=1e-4) # quadgk first output is the value of the integral, second value is an estimated error welfare_manual(r;k=k,a=a) = - k*(exp(-a*r) - 1) / a discount(t; ρ=ρ) = exp(-ρ*t) # discount function # Equation of motion of the state variable var_stocks(r) = -r 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, S, Infinite(t)) ## state variables @variable(m, r, Infinite(t)) ## control variables @objective(m, Max, integral(welfare_manual(r), t, weight_func = discount)) @constraint(m, S(0) == S₀) @constraint(m, S(test) == Sτ) @constraint(m, r >= 0) @constraint(m, c1, deriv(S, t) == var_stocks(r)) optimize!(m) termination_status(m) r_opt = value(r) S_opt = value(S) ts = supports(t) opt_obj = objective_value(m) ix = 2:(length(ts)-1) # index for plotting plot(ts[ix], S_opt[ix], lab = "S: resource stock") plot(ts[ix], r_opt[ix], lab = "r: resource extraction") plot(ts[ix], demand_price.(r_opt[ix]), lab = "p: price") # For checking... price_manual(t;k=k,ρ=ρ,T=T) = k*exp(ρ*(t-T)) plot!(ts[ix], price_manual.(ts[ix]), lab = "p: price (manual")