OptiMimi: Optimization of Mimi models
OptiMimi provides a simplified interface for finding optimal parameter
values for Mimi models (https://github.com/anthofflab/Mimi.jl). The
core interface consists of
problem to define the optimization
problem, and calls to
solution to solve it.
The package provides two major approaches for performing optimization:
General and Linear programming. The General approach allows takes
full models and applies non-linear optimization techniques to optimize
parameters within them. The Linear programming approach allows models
to define linear operation matrices which represent the computations
they perform in Mimi's
run_timestep function. In addition, OptiMimi
offers a way to automatically generate these matrices.
OptiMimi supports autodifferentiation using ForwardDiff. To use it,
Model object must be created with the optional argument
Number, to specify a general-enough type to handle dual-numbers:
m = Model(Number) If this is used, when the optimization
problem is created, OptiMimi use gradient algorithms.
The General approach in OptiMimi can use algorithms in both NLopt and BlackBoxOptim. The Linear programming approach uses any solver supported by MathProgBase.
Constructing an optimization problem (General approach)
Setup an optimization problem with the
problem(model, names, lowers, uppers, objective; constraints, algorithm)
modelis a Mimi model, with some parameters intended for optimization.
namesare lists of the parameter names to be optimized, and must be the same length.
uppersare a list of lower and upper bounds and must be the same length as
names; the same bounds are used for all values within a parameter array.
objectiveis a function that takes a Mimi Model, with parameters set and fully executed, and returns a value to be maximized.
constraints(optional) is a vector of inequality constraint functions; each takes a Mimi Model, with parameters set (but not necessarily executed), and should return < 0 when the constraint is satisfied.
algorithm(optional) is a symbol, currently chosen from the NLopt algorithms.
The return value is an object of the OptimizationProblem type, to be passed to
Start by creating a Mimi model and ensuring that it runs with all
parameters set. In the example below,
my_model is a model with an
agriculture component, in which N regions are evaluated in a single
timestep to consume energy and produce corn.
The optimization maximizes economic value, trading off the value of the corn against the cost of the energy for fertilizer. We also add a constraint that the total fertilizer cannot be more than 1 million kg, to reduce environmental impacts.
using OptiMimi # Prices of goods p_F = 0.25 # the global price of food (per kg of corn) p_E = 0.4 # the global price of fuel (per kWh) # Objective to maximize economic output function objective(model::Model) sum(my_model[:agriculture, :cornproduction] * p_F - my_model[:agriculture, :cornenergyuse] * p_E) end constraints = [model -> sum(model.components[:agriculture].Parameters.fertilizer) - 1e6] # Setup the optimization optprob = problem(my_model, [:agriculture], [:fertilizer], [0.], [1e6], objective, constraints=constraints)
Note that (1) the objective function is provided with the prepared
model, not with the raw initialization values, and (2) even though
there are N values to be set and optimized over in the
parameter, the lower and upper bounds are only specified once.
Solving the optimization problem
The optimization problem, returned by
problem is solved by
solution(optprob, generator; maxiter, verbose)
optprobis the result of the
generatoris a function of no arguments, which returns a full set of parameter values, with values concatenated across parameters in the order of
namesabove. This should generally be stochastic, and if the specified model fails the constraints then
generatorwill be called again until it succeeds.
maxiter(optional) is the maximum number of iterations for the optimization; currently it only is used for the maximum number of times that
generatorwill be called.
verbose(optional) is a boolean to specify if status messages should be printed.
The return value is a tuple of the maximum found objective value, and the concatenated collection of model parameters that produced it.
Continuing the example above, we solve the optimization problem:
(maxf, maxx) = solution(optprob, () -> [0. for i in 1:5]) println(maxf) println(maxx)
Our generator function can only generate a single initial condition: all 0's.
Once you have a solution, you can initialize a model with it usnig the
setparameters(my_model, [:agriculture], [:fertilizer], maxx)
This is most useful when there are multiple parameters being optimized, or the parameters have multiple dimensions.
Constructing an optimization problem (Linear programming approach)
Linear programming allows for vastly faster optimizations, so long as the constraints and objective can be translated into linear algebra relationships. See https://en.wikipedia.org/wiki/Linear_programming for details.
Within OptiMimi, the large matrix and vectors which define the linear programming constraints and objective are developed in segments. In line with the model organization structure of Mimi, the computations which relate variables to parameters are kept separate, and organized with each component. These computations can be used directly, if variables of a component are constrained and the parameters of that same component are optimized over; or, they can be connected across multiple components. Rather than defining the entire matrix at once, OptiMimi allows segments of the matrix specific to each component to be defined separately.
Segments of the constraint matrix are encapsulated in
LinearProgrammingRoom objects, which combine a sparse matrix with
information about a single model parameter and variable. A column
vector, combined with information about a single model parameter or
variable, is encapsulated in a
LinearProgrammingShaft (its transpose) object. A
LinearProgrammingHouse contains the set of all matrices for the
There are a variety of functions available which create these objects or manipulate them. Some of the most commonly used operations are below:
roomdiagonal: creates a room for a variable which is a direct scaling of a parameter.
*: Allows two rooms to be multiplied together, which corresponds to "connecting" the variable of the first as the parameter of the second; or allows a room and a hall or shaft to be multiplied so which results in a weighted sum of the variables in the room (e.g., as an objective).
room_relabel: Relabel the variable of a room so that it corresponds to the parameter name of another room which it is to connected to (multiplied with).
The usual process for setting up a linear programming problem is as follows:
Mimi components are written as normal, with
@defcompcalls, and the Mimi model is constructed and external parameters are initialized.
Individual functions are specified for each component describing the gradient of a variable with respect to a parameter, using the
roomfunctions or the automated option in
makeroom.jl. The naming typically is as follows:
LinearProgrammingHosueis constructed, specifying the optimization parameters and constraint variables, like so:
parameters, constcomponents, variables)```
The objective is specified with a
setobjective!call, often as a sum over variable specified by a gradient function, e.g.,
Constraints are specified with
setconstraintoffset!calls. In all cases, the relationship must be specified so that
variable < offset. This looks like,
grad_COMPONENT_VARIABLE_PARAMETER(model))``` ```setconstraintoffset!(house, constraintoffset_COMPONENT_VARIABLE(model))```
The optimization is performed, using any solver supported by
houseoptimizefunction. For example:
using MathProgBase using Gurobi solver = GurobiSolver() sol = houseoptimize(house, solver)
- The result is studied using the
summarizeparameters, or, if the optimzation failed,