Decision Optimization 101

Never take a wrong decision again - simply compute your optimal decisions by modeling the reality correctly

Version 1-dev - last update: 2018-07-22

Click here to browse back to the list of tutorials.

Optimization 101

Never take a wrong decision again - simply compute your optimal decisions by modeling the reality correctly.

Summary

In this tutorial we are going to learn about all four important parts of an optimization model and how to formulate a quantitative decision problem such that it can be solved by an optimization solver. Furthermore four different ways to implement and solve optimization models in R are shown - two using a manual matrix-based approach and two using contemporary optimization modeling environments.

Libraries

We need the following R libraries to follow all the code fragments in this tutorial, while ROI and modopt.matlab are used for matrix-based modeling the packages CVXR and ompr are optimization modeling frameworks.

### Libraries
library(magrittr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)
library(CVXR)
library(modopt.matlab)

Optimization 101

In general, books on optimization put the main focus mainly on algorithms and solution techniques. This whole tutorial Decision Optimization with R has been written such that an edited excerpt of these tutorials will be used as an appendix to the book Become a Quant with R which applies Simulation, Optimization and Data Science to the field of Quantitative, Computational and Algorithmic Finance. Especially from the perspective of using optimization basically to solve our Quantitative Finance tasks (but this applies to all other areas of application too) we should take a slightly different view on Optimization. We should look at this beautiful set of methods solely from a decision perspective. This is what the scientific field of Operations Research (OR) is mainly doing and this is why they rebranded their scientific field as The Science of Better, i.e. the scientific approach of taking better decisions. This is exactly what we need to be able to and ever since the bible of optimization in OR by Hillier and Lieberman (1982) there were lots of books related to applications, e.g. we will use many examples from Bertsimas and Freund (2000).

Most books from the field of OR start looking at the problems by starting with linear programming and the simplex algorithm. In Finance on the other hand people start by learning optimization by solving a Markowitz portfolio optimization problem using quadratic programming techniques - sometimes even solving the problem manually using Lagrange multipliers and seeking the extremum of the Lagrangian. Luckily this restrictive starting point (i.e. linear programming and/or quadratic programming) has been changed to a rather general look on convex optimization especially since the publication of a book convincingly entitled Convex Optimization S. P. Boyd and Vandenberghe (2004) which is freely available from the authors homepage and highly recommended to anyone interested in the field of Optimization.

In this regard for those that are even more interested in the basics of numerical algebra and related topics may also take a look at the new book Introduction to Applied Linear Algebra – Vectors, Matrices, and Least Squares S. Boyd and Vandenberghe (2018) which is also freely available from the authors homepage

The four parts of an optimization model

So whenever we want to solve a decision optimization problem with a computer, we need to build a model of the problem we want to be solved. This is a no brainer, but while the solution methods are pretty complex, the craft of modeling a perfect decision optimization problem is actually straightforward. There are four imporatant parts of the model to be identified and defined:

  1. Decision Variables
  2. Objective Function
  3. Constraints
  4. Parameters

Part One: Decision Variables

Decision variables are all values, i.e. numbers that we want to be determined by our optimizer. These values are often represented as a vector of real-values or integer numbers. Sometimes it is clear which decision variables have to be specified for a certain decision problem but sometimes figuring out which variables are required can be tricky. The best way to learn this specification is by taking a look at various use-cases and application examples.

Part Two: Objective function

The objective function captures the specification of our goal. While in engineering and other disciplines of science it is often clear what the goal is because there are laws of nature and similar given rules to consider. However in decision optimization things may be different. There are cases where the goal is not clear in advance or where different ideas of the goal exist.

The objective function is mainly either a minimization or a maximization. In some cases more than one goal is combined - then we talk about multi-criteria optimization. The latter sub-field using more than one objective function simultaneously is extremely important for the field of Finance as actually all portfolio optimization problems in the spirit of Markowitz risk/return analysis are bi-criteria optimization problems where the investor is trying to maximize the (expected) return/profit as well as aiming at minimizing the risk of the investment, i.e. the meta-problem can be written down as follows:

maximize (expected) return minimize risk

More about different ways to formulate bi-criteria optimization models especially in regard to portfolio selection will be covered in another tutorial.

Part Three: Constraints

Constraints define all restrictions determined by nature, our company, the customers, the regulator and everything that hinders us e.g. from maximizing our profit to infinity or lowering our costs to zero. Specifying constraints allows us to focus on the integration of objective real-world constraints and dynamics, i.e. the general event handling.

Part Four: Parameters

Parameters are values, which are determined by nature, our company, the customers, the regulator and are exogeneously given. These parameters are expecially important if uncertainty is introduced because then we need a certain way to specify the uncertainty e.g. by some probability distribution function. In the latter case it is possible to integrate subjective views into our decision optimization model. We will get back to this later on.

The first Use-Case: Business Optimization

We consider an optimization model from the field of Business Optimization. Business optimization can teach a lot about decision optimization in general. We consider an example from the book by Bertsimas and Freund (2000) and want to find the optimal production plan for the Gemstone Tool Company (GTC).

The Winnipeg plant of the Gemstone Tool Company (GTC) produces wrenches and pliers. Both tools are made from steel, and the production requires (1) molding the tools on a molding machine and (2) assembling the tools on an assembly machine.

GTC wants to determine the daily production of wrenches and pliers so as to maximize their contribution to earnings, subject to the daily demand limits.

Wrenches Pliers Availability
Required steel (lbs) 1.5 1.0 27000 lbs/day
Required time on molding machine (h) 1.0 1.0 21000 h/day
Required time on assembly machine (h) 0.3 0.5 9000 h/day
Demand limit (units/day) 15000 16000
Contribution to earnings ($/unit) $0.13 $0.10

So e.g. one wrench requires 1.5 lbs of steel, one hour each on the molding machine as well as the assembly machine and we know that we cannot sell more than 15000 wrenches per day.

This is clearly a deterministic decision optimization problem as all values are directly specified. From a decision optimization perspective the extension to optimization under undercertainty is often crucial.

In any way the questions which the Gemstone Tool Company (GTC) raises are:

  1. What is the optimal daily production plan?
  2. What are the daily contributions to earnings from this plan?
  3. Which resources would be most critical in this plan?

Whe need to take the following approach towards building our optimization model, which includes the following steps:

  • Step 1: Identify the decision variables.
  • Step 2: Identify the objective function.
  • Step 3: Identify the constraints.

Step 1: Identify the decision variables

What decision(s) must be made in order to solve the problem? Actually we need to determine the optimal number of wrenches and the optimal number of pliers which have to be produced per day, so we end up in a decision vector of length two, i.e. \((W, P)\) where

W = number of *wrenches* produced per day
P = number of *pliers* produced per day

Step 2: Identify the objective function

For a production company the objective is typically either

  1. to maximize revenue, or
  2. to minimize cost.

In this specific case we aim at maximize the contribution to earnings, so our objective function looks as follows:

maximize 0.13 W + 0.10 P

Step 3: Identify the constraints

What is restricting GTC from making infinite revenues or incurring zero costs?

First GTC is restricted by the availability of steel per day and both wrenches and pliers need a different amount of steel per unit. We do have 27000 lbs per day, a wrench required 1.5 lbs and one pair of pliers 1 lbs, such that we obtain the following constraint:

1.5 W + 1 P <= 27000

Next we add the constraint regarding the availability of molding machines where we have 21000h in total per day and both products require 1h of molding, i.e.

1 W + 1 P <= 21000
Availability of assembly machine 0.3W + 0.5P <= 9000

Now we have to care about the bounds on our decision variables. First, we have to ensure not to exceed our daily demand limits for each product.

W <= 15000, P <= 16000

Furthermore it is clear that we cannot produce negative amounts of tools. While this should be obvious it is always a good idea to add all possible bounds to all decision variables.

W >= 0, P >= 0

Once we indentified all the important things we can set up one meta-model of the optimization problem, e.g.

variables W, P
maximize 0.13W + 0.10P
subject to
  1.5 W + 1 P <= 27000
  1 W + 1 P <= 21000
  0.3 W + 0.5 P <= 9000
  W <= 15000, P <= 16000
  W >= 0, P >= 0

Solving our optimization model using R

Once we specified our optimization model (see the meta model above) we need to tell our computer about it such that we may obtain a quantitative solution.

Manual solution - Formulation

The classical approach to model an optimization problem is the manual way, i.e. specifying the problem as vectors and matrices. In this case we have to keep track of the ordering of our variables, so let’s fix the ordering of first the wrenches and then the pliers \(W, P\) - in matrix format we write down complete the model as follows:

\[ \begin{array}{lrcrcclcrl} & & & & W & P & & & & \vert \text{ decision variables}\\ \text{maximize} & & & ( & 0.13, & 0.1 & ) & & & \vert \text{ objective function} \\ \text{subject to} & & & ( & 1.5, & 1 & ) & \leq & 27000 & \vert \text{ availability of steel } \\ \text{(constraints)} & & & ( & 1, & 1 & ) & \leq & 21000 & \vert \text{ molding machine} \\ & & & ( & 0.3, & 0.5 & ) & \leq & 9000 & \vert \text{ assembly machine} \\ & & & ( & 0, & 0 & ) & & & \vert \text{ lower bound } \\ & & & ( & 15000, & 16000 & ) & & & \vert \text{ upper bound } \\ \end{array} \] In this case the objective function is a vector of the length of decision variables, i.e. two, so we obtain \(f = ( 0.13, 0.1 )\) which we need to maximize. The constraints consist of a matrix \(A\) and a vector \(b\). In optimization theory the matrix \(A\) is sometimes called the left-hand side (of the constraints) while \(b\) is likewise called the right-hand side.

The size of \(A\) is the number of variables (columns) times the number of constraints (rows) while the length of \(b\) equals the number of constraints. These two parts can be extracted as follows:

\[ A = \begin{pmatrix} 1.5 & 1 \\ 1 & 1 \\ 0.3 & 0.5 \\ \end{pmatrix}, b = \begin{pmatrix} 27000 \\ 21000 \\ 9000 \end{pmatrix} \]

It should be added that sometimes the matrix \(A\) and vector \(b\) are split into two parts, one inequality \(A \leq b\) and one equality part \(A = b\). In our example we only definied inequality constraints.

Furthermore we explicitely specify lower and upper bounds on the decision variables as two vectors each of the size of the number of decision variables. In our case we obtain \(l = (0, 0)\) and \(b = (15000, 16000)\).

Manual solution - using modopt.matlab

Using the R package modopt.matlab (disclaimer: this package has been created by the author of this tutorial) we just have to enter all the details of the above matrix structure directly to obtain the solution:

### Manual solution - using modopt.matlab
sol <- linprog(f=-c(0.13, 0.1),
               A=matrix(c(1.5, 1, 1, 1, 0.3, 0.5), nrow = 3, byrow=TRUE),
               b=c(27000, 21000, 9000),
               lb=c(0, 0),
               ub=c(15000, 16000))

solution <- round(sol$x)
print(solution)
## [1] 12000  9000

We obtain an optional solution to produce 1.210^{4} wrenches and 9000 pliers.

Manual solution - using ROI

The package modopt.matlab provides functions to rewrite optimization models into the rather complete optimization package ROI, i.e. the underlying optimization functions are based on this package. If one aims to employ ROI directly the model looks as follows:

### Manual solution - using ROI
lp <- OP(objective = c(0.13, 0.1), 
         L_constraint(L = matrix(c(1.5, 1, 1, 1, 0.3, 0.5), nrow = 3, byrow=TRUE),
                      dir = c("<=", "<=", "<="),
                      rhs = c(27000, 21000, 9000)), 
         maximum = TRUE,
         bounds = V_bound(li=c(1:2), lb=c(0, 0), ui=c(1:2), ub=c(15000, 16000)))

sol <- ROI_solve(lp, solver = "glpk")

When we take a closer look at the ROI code then we see why for some cases it might be useful to use the much simpler modopt.matlab version. However, ROI is much more powerful and thus naturally needs more complicate structure. To specify an objective model, we use the function OP() and add the objective vector, the constraints: here we explicitely specify that the constraints are linear \(L\). Futhermore the bounds are specified using a special function V_bound() which allows for a sparse structure, i.e. by setting the indices (li, ui) and values (lb, ub) for each bound.

Using modeling languages - CVXR

For convenience reasons we may use an algebraic optimization modeling language. Since 2017 there are fortunately at least two modeling languages available in R, i.e. the packages CVXR and ompr. Both use a rather different syntax and also apply different optimization solvers. Let’s start with CVXR:

### Using modeling languages - CVXR
wrench <- Variable(1)
pliers <- Variable(1)

objective <- Maximize(0.13*wrench + 0.10*pliers)
constraints <- list(1.5*wrench + 1.0*pliers <= 27000,
                    1.0*wrench + 1.0*pliers <= 21000,
                    0.3*wrench + 0.5*pliers <= 9000,
                    wrench >= 0, pliers >= 0,
                    wrench <= 15000, pliers <= 16000)
gtc <- Problem(objective, constraints)
solution <- solve(gtc)

solution$getValue(wrench)
## [1] 12000
solution$getValue(pliers)
## [1] 9000

Using CVXR we start by defining out decision variables using the function Variable(). Then we specify the objective function using Maximize() and create a list with all the constraints, i.e. the budget constraint as well as the long/short limits. The objective function as well as the constraints are assembled together using the function Problem(). The resulting optimization model can be solved using solve().

Using modeling languages - ompr

First we need to remove the two previously optimized decision variables to avoid any confusion between the two optimization packages.

rm(wrench, pliers)

Then we may proceed with formulating the same problem with ompr and also making use of the magrittr pipe via the package dplyr.

### Using modeling languages - ompr
result <- MIPModel() %>%
  add_variable(wrench, type = "continuous", lb = 0, ub=15000) %>%
  add_variable(pliers, type = "continuous", lb = 0, ub=16000) %>%
  set_objective(0.13*wrench + 0.10*pliers, "max") %>%
  add_constraint(1.5*wrench + 1.0*pliers <= 27000) %>%
  add_constraint(1.0*wrench + 1.0*pliers <= 21000) %>%
  add_constraint(0.3*wrench + 0.5*pliers <= 9000) %>%
  solve_model(with_ROI(solver = "glpk")) 
get_solution(result, wrench)
## wrench 
##  12000
get_solution(result, pliers)
## pliers 
##   9000

With ompr we start creating the model issuing MIPModel() and apply different functions to this model. This fact allows for using the magrittr pipe (implicitly loaded through the package dpylr) conveniently. One difference is that we explicitly set the lower and upper bounds during the definition of the decision variables using add_variable(). The rest of the model should be self-explanatory.

References

Bertsimas, Dimitris, and Robert M. Freund. 2000. Data, Models, and Decisions: The Fundamentals of Management Science. 1st ed. South-Western College Publishing.

Boyd, Stephen P., and Lieven Vandenberghe. 2004. Convex Optimization. Cambridge University Press.

Boyd, Stephen, and Lieven Vandenberghe. 2018. Introduction to Applied Linear Algebra: Vectors, Matrices, and Least Squares. Cambridge University Press.

Hillier, Frederick, and Gerald J. Lieberman. 1982. Introduction to Operations Research. McGraw-Hill, USA.

Discussion