@lrennels , Hey!
Thanks a lot for these suggestions. I solved my previous problem. The issue was NLopt was giving me “FTOL_REACHED” as the output. As I am new to this optimization routine, I was expecting “NLOPT_SUCCESS”. I, therefore, got confused. I appreciate all the help.
I now have a new issue that I would appreciate a little help on. I am trying to optimize both the abatement rate and the savings rate. The way I am doing this is to pass one long parameter vector to the objective function. Is there a way I can pass two separate parameter vectors of length equivalent to the time dimension of the DICE2013 model? For instance, I would like to pass these to `max_objective’ function in NLopt. I am posting the code for the long parameter vector approach and also what I would ideally like to do below.
Long parameter vector approach:
# set directory
cd("C:/Users/singh/Dropbox/ResearchWork/Potential/mimiDICE2016R2")
# load packages
using NLopt, Mimi, MimiDICE2013, CSV, DataFrames
# write the function which takes as input an instance of DICE model
function optimise_model(m::Model=MimiDICE2013.get_model())
# define periods for which the optimization takes place
n_objectives = length(m.md.dim_dict[:time])
# define time by which the optimization should stop
stop_time = 640
# define the tolerance parameter for the functional value
tolerance = 1e-6
# determine the optimization algorithm
optimization_algorithm = :LN_SBPLX
# Create lower bound
lower_bound = [0.039;zeros(n_objectives-1);zeros(n_objectives)]
# Create upper bound
upper_bound = [0.039;ones(28);1.2*ones(31);ones(n_objectives)]
# Create initial condition for algorithm (this is a fraction of upper bound)
starting_point = [ones(n_objectives) ./ 1000;ones(n_objectives) ./ 1000]
# generate an optimization problem of NLopt type
opt = Opt(optimization_algorithm, n_objectives*2)
# Set the bounds.
lower_bounds!(opt, lower_bound)
upper_bounds!(opt, upper_bound)
# Assign the objective function to maximize.
max_objective!(opt, (x, grad) -> construct_objective(m, x))
# Set termination time.
maxtime!(opt, stop_time)
# Set optimizatoin tolerance (will stop if |Δf| / |f| < tolerance from one iteration to the next).
ftol_rel!(opt, tolerance)
# Optimize model.
maximum_objective_value, optimised_policy_vector, convergence_result = optimize(opt, starting_point)
# return the object: model and dictionary object
return (m = m, opt_result = Dict([("Maximum objective value", maximum_objective_value),
("Optimised policy vector", optimised_policy_vector),
("Convergence result", convergence_result)]))
end
# generate function that creates an objective function with the model and optimized abatement as the argument
function construct_objective(m::Model, optimised_abatement::Array{Float64,1})
# update abatement rate and run the model
update_param!(m, :MIU, optimised_abatement[1:n_objectives])
update_param!(m, :S, optimised_abatement[n_objectives+1:2*n_objectives])
run(m)
return m[:welfare, :UTILITY]
end
# run the model
model_opt = optimise_model()
# extract the policy vector
policy_opt = model_opt.opt_result[:"Optimised policy vector"]
# update the parameter and run model
m = MimiDICE2013.get_model()
update_param!(m, :MIU, policy_opt[1:60])
update_param!(m, :S, policy_opt[61:120])
run(m)
# export to CSV
CSV.write("optimized_model.csv", DataFrame((years = 2010:5:2305,
abatement = m[:emissions, :MIU],
cumulative_welfare = m[:welfare, :CUMCEMUTOTPER],
savings = m[:neteconomy, :S])))
What I would ideally like to do:
# set directory
cd("C:/Users/singh/Dropbox/ResearchWork/Potential/mimiDICE2016R2")
# load packages
using NLopt, Mimi, MimiDICE2013, CSV, DataFrames
# write the function which takes as input an instance of DICE model
function optimise_model(m::Model=MimiDICE2013.get_model())
# define periods for which the optimization takes place
n_objectives = length(m.md.dim_dict[:time])
# define time by which the optimization should stop
stop_time = 640
# define the tolerance parameter for the functional value
tolerance = 1e-6
# determine the optimization algorithm
optimization_algorithm = :LN_SBPLX
# Create lower bound
lower_bound = [0.039;zeros(n_objectives-1);zeros(n_objectives)]
# Create upper bound
upper_bound = [0.039;ones(28);1.2*ones(31);ones(n_objectives)]
# Create initial condition for algorithm (this is a fraction of upper bound)
starting_point = [ones(n_objectives) ./ 1000;ones(n_objectives) ./ 1000]
# generate an optimization problem of NLopt type
opt = Opt(optimization_algorithm, n_objectives*2)
# Set the bounds.
lower_bounds!(opt, lower_bound)
upper_bounds!(opt, upper_bound)
# Assign the objective function to maximize.
max_objective!(opt, (x, grad) -> construct_objective(m, x))
# Set termination time.
maxtime!(opt, stop_time)
# Set optimizatoin tolerance (will stop if |Δf| / |f| < tolerance from one iteration to the next).
ftol_rel!(opt, tolerance)
# Optimize model.
maximum_objective_value, optimised_policy_vector, convergence_result = optimize(opt, starting_point)
# return the object: model and dictionary object
return (m = m, opt_result = Dict([("Maximum objective value", maximum_objective_value),
("Optimised policy vector", optimised_policy_vector),
("Convergence result", convergence_result)]))
end
# generate function that creates an objective function with the model and optimized abatement as the argument
function construct_objective(m::Model, optimised_abatement::Array{Float64,1}, optimised_savings::Array{Float64,1})
# update abatement rate and run the model
update_param!(m, :MIU, optimised_abatement)
update_param!(m, :S, optimised_savings)
run(m)
return m[:welfare, :UTILITY]
end
# run the model
model_opt = optimise_model()
