-
Notifications
You must be signed in to change notification settings - Fork 17
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Create example consumption_savings.jl (#141)
* Create consumption_savings.jl * Update consumption_savings.jl
- Loading branch information
Showing
1 changed file
with
113 additions
and
0 deletions.
There are no files selected for viewing
113 changes: 113 additions & 0 deletions
113
docs/src/examples/Optimal Control/consumption_savings.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,113 @@ | ||
# # Lifecycle consumption savings problem | ||
# In this case study, a household endowed with ``B_0`` dollars of wealth | ||
# must decide how much to consume and save | ||
# to maximize its (quadratic) utility over its finite lifecycle. | ||
|
||
# ## Formulation | ||
|
||
# The corresponding dynamic optimization problem is expressed: | ||
# ```math | ||
# \begin{aligned} | ||
# &V(B_0,0) = &&\underset{c(t), B(t)}{\text{max}} &&& \int_{t = 0}^{t=T} e^{-\rho t} u(c(t)) dt \\ | ||
# &\text{s.t.} \\ | ||
# &&& \frac{dB}{dt} = r \times B(t) - c(t), && t \in [0,T] \\ | ||
# &&& B(0) = B_0 \\ | ||
# &&& B(T) = 0 | ||
# \end{aligned} | ||
# ``` | ||
# where | ||
# the household lives during time ``t \in [0,T]``, | ||
# the state variable ``B(t)`` is the household's stock of wealth at time ``t``, | ||
# the choice variable ``c(t)`` is the household's consumption at time ``t``, | ||
# ``r`` is the interest rate, | ||
# ``B_0`` is the household's wealth endowment (initial condition), | ||
# ``B(T) = 0`` is the terminal condition, | ||
# and ``T`` is the time | ||
# horizon. | ||
|
||
# ## Model Definition | ||
|
||
# Let's implement this in `InfiniteOpt` and first import the packages we need: | ||
using InfiniteOpt, Ipopt | ||
|
||
# We set the preference and constraint parameters: | ||
ρ = 0.025 # discount rate | ||
k = 100.0 # utility bliss point | ||
T = 10.0 # life horizon | ||
r = 0.05 # interest rate | ||
B0 = 100.0 # endowment | ||
u(c; k=k) = -(c - k)^2 # utility function | ||
discount(t; ρ=ρ) = exp(-ρ*t) # discount function | ||
BC(B, c; r=r) = r*B - c # budget constraint | ||
|
||
# We set the hyperparameters: | ||
opt = Ipopt.Optimizer # desired solver | ||
ns = 1_000; # number of gridpoints | ||
|
||
# We initialize the infinite model and choose the Ipopt solver: | ||
m = InfiniteModel(opt) | ||
|
||
# Let's specify our infinite parameter which is time ``t \in [0, T]``: | ||
@infinite_parameter(m, t in [0, T], num_supports = ns) | ||
|
||
# Now let's specify the variables: | ||
@variable(m, B, Infinite(t)) ## state variables | ||
@variable(m, c, Infinite(t)) ## control variables | ||
|
||
# Specify the objective: | ||
@objective(m, Max, integral(u(c), t, weight_func = discount)) | ||
|
||
# Set the initial/terminal conditions: | ||
@constraint(m, B == B0, DomainRestrictions(t => 0)) | ||
@constraint(m, B == 0, DomainRestrictions(t => T)) | ||
|
||
# Set the budget constraint: | ||
@constraint(m, c1, deriv(B, t) == BC(B, c; r=r)) | ||
|
||
# ## Problem Solution | ||
|
||
# Optimize the model: | ||
optimize!(m) | ||
termination_status(m) | ||
|
||
# Extract the results: | ||
opt_obj = objective_value(m) | ||
c_opt, B_opt = value(c), value(B) | ||
ts = supports(t) | ||
|
||
# Plot the results: | ||
using Plots | ||
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") | ||
|
||
# That's it, now we have our optimal trajectory! | ||
|
||
# This very simple problem has a closed form solution: | ||
λ1 = exp((r)T) | ||
λ2 = exp(-(r-ρ)T) | ||
den = (λ1-λ2)r | ||
Ω1 = (k + (r*B0-k)λ2)/den | ||
Ω2 = (k + (r*B0-k)λ1)/den | ||
c0 = r*B0 + (r)Ω1 + (r-ρ)Ω2 | ||
BB(t) = (k/r) - Ω1*exp((r)t) + Ω2*exp(-(r-ρ)t) | ||
cc(t) = k + (c0-k)*exp(-(r-ρ)t) | ||
|
||
# Compare the solution given by `InfiniteOpt` with the closed form: | ||
plot(legend=:topright); | ||
plot!(ts[ix], c_opt[ix], color = 1, lab = "c: consumption, InfiniteOpt"); | ||
plot!(ts[ix], cc, color = 1, linestyle=:dash, lab = "c: consumption, closed form"); | ||
plot!(ts[ix], B_opt[ix], color = 4, lab = "B: wealth balance, InfiniteOpt"); | ||
plot!(ts[ix], BB, color = 4, linestyle=:dash, lab = "B: wealth balance, closed form") | ||
|
||
# Not bad! | ||
|
||
|
||
|
||
# ### Maintenance Tests | ||
# These are here to ensure this example stays up to date. | ||
using Test | ||
@test termination_status(m) == MOI.LOCALLY_SOLVED | ||
@test has_values(m) | ||
@test B_opt isa Vector{<:Real} | ||
@test c_opt isa Vector{<:Real} |