ECON457 R lab 2 Optimization in R

http://cran.r-project.org/web/views/Optimization.html

https://www.rdocumentation.org/packages/ROI/versions/0.2-1

The R Optimization Infrastructure (ROI)

The R Optimization Infrastructure (ROI) package promotes the development and use of interoperable (open source) optimization problem solvers for R.

ROI handle LP up to MILP and MIQCP problems using the following supported solvers :

1. lpSolve

2. quadprog

3. Rcplex

4. Rglpk (default)

5. Rsymphony

Solvers

CPLEX is a commercial solver. Students can get academic license.

GLPK is free. R has the high level interface Rglpk.

http://winglpk.sourceforge.net/

https://cran.r-project.org/web/packages/Rglpk/index.html

Linear Progarmming

#install.packages("ROI")
#install.packages(c( "PortfolioAnalytics","ECOSolveR","Rsymphony","Rglpk","ROI.plugin.symphony","ROI.plugin.glpk","fPortfolio","ROI.plugin.quadprog","ROI.plugin.nloptr","ROI.plugin.ipop","ROI.plugin.ecos"))

library("ROI")

ROI: R Optimization Infrastructure

Registered solver plugins: nlminb, ecos, glpk, ipop, nloptr, quadprog, symphony.

Default solver: auto.

#library("ROI.plugin.glpk")
#library("ROI.plugin.symphony")
#library(ROI.plugin.ipop)

#library(ROI.plugin.nloptr)
#library(ROI.plugin.quadprog)
#library(ROI.plugin.ecos)

ROI_registered_solvers()

               nlminb                  ecos                  glpk 
  "ROI.plugin.nlminb"     "ROI.plugin.ecos"     "ROI.plugin.glpk" 
                 ipop                nloptr              quadprog 
    "ROI.plugin.ipop"   "ROI.plugin.nloptr" "ROI.plugin.quadprog" 
             symphony 
"ROI.plugin.symphony" 

ROI.plugin.nlminb

ROI.plugin.ecos

ROI.plugin.glpk

ROI.plugin.ipop

ROI.plugin.nloptr

ROI.plugin.quadprog

ROI.plugin.symphony

## Simple linear program.
## maximize:   2 x_1 + 4 x_2 + 3 x_3
## subject to: 3 x_1 + 4 x_2 + 2 x_3 <= 60
##             2 x_1 +   x_2 +   x_3 <= 40
##               x_1 + 3 x_2 + 2 x_3 <= 80
##               x_1, x_2, x_3 are non-negative real numbers

LP <- OP( c(2, 4, 3),
          L_constraint(L = matrix(c(3, 2, 1, 4, 1, 3, 2, 2, 2), nrow = 3),
                       dir = c("<=", "<=", "<="),
                       rhs = c(60, 40, 80)),
          max = TRUE )
LP

ROI Optimization Problem:

Maximize a linear objective function with
- 3 objective variables,

subject to
- 3 constraints of type linear.

Solve

sol <- ROI_solve(LP)# , solver = "glpk")
sol

Optimal solution found.
The objective value is: 7.666667e+01

Solutions

solution(sol, type = c("primal"))

[1]  0.000000  6.666667 16.666667

solution(sol, type = c("dual"))

[1] -1.833333  0.000000  0.000000

Extract solutions

sol$solution

[1]  0.000000  6.666667 16.666667

sol$objval

[1] 76.66667

sol$status

$code
[1] 0

$msg
  solver glpk
    code 5
  symbol GLP_OPT
 message Solution is optimal.
roi_code 0

sol$message

$optimum
[1] 76.66667

$solution
[1]  0.000000  6.666667 16.666667

$status
[1] 5

$solution_dual
[1] -1.833333  0.000000  0.000000

$auxiliary
$auxiliary$primal
[1] 60.00000 40.00000 53.33333

$auxiliary$dual
[1] 0.8333333 0.6666667 0.0000000

Quadratic Progarmming

It will require putting problem in matrix form.

## Simple quadratic program.
## minimize: - 5 x_2 + 1/2 (x_1^2 + x_2^2 + x_3^2)
## subject to: -4 x_1 - 3 x_2       >= -8
##              2 x_1 +   x_2       >=  2
##                    - 2 x_2 + x_3 >=  0

QP <- OP( Q_objective (Q = diag(1, 3), L = c(0, -5, 0)),
          L_constraint(L = matrix(c(-4,-3,0,2,1,0,0,-2,1),
                                  ncol = 3, byrow = TRUE),
                       dir = rep(">=", 3),
                       rhs = c(-8,2,0)) )
QP

ROI Optimization Problem:

Minimize a quadratic objective function with
- 3 objective variables,

subject to
- 3 constraints of type linear.

Solve

sol <- ROI_solve(QP, solver = "quadprog")
sol

Optimal solution found.
The objective value is: -2.380952e+00

Extract solutions

sol$solution

[1] 0.4761905 1.0476190 2.0952381

sol$objval

[1] -2.380952

sol$status

$code
[1] 0

$msg
  solver quadprog
    code 0
  symbol OPTIMAL
 message Solution is optimal
roi_code 0

sol$message

NULL

Solutions

solution(sol, type = c("primal"))

[1] 0.4761905 1.0476190 2.0952381

solution(sol, type = c("dual"))

[1] NA

Portfolio optimization - minimum variance

## Portfolio optimization - minimum variance
## -----------------------------------------
## get monthly returns of 30 US stocks
data( US30 )
r <- na.omit( US30 )
## objective function to minimize
obj <- Q_objective( 2*cov(r) )
## full investment constraint
full_invest <- L_constraint( rep(1, ncol(US30)), "==", 1 )
## create optimization problem / long-only

(op <- OP( objective = obj, constraints = full_invest ))

ROI Optimization Problem:

Minimize a quadratic objective function with
- 30 objective variables,

subject to
- 1 constraints of type linear.

Solution

## solve the problem - only works if a QP solver is registered
##  
res <- ROI_solve( op )
res

Optimal solution found.
The objective value is: 9.326635e-04

sol <- solution( res )
names( sol ) <- colnames( US30 )
round( sol[ which(sol > 1/10^6) ], 3 )

  IBM   MCD   MRK    PG     T    VZ   WMT   XOM 
0.165 0.301 0.005 0.138 0.031 0.091 0.169 0.101 

LP with boundary

suppressMessages(library(ROI))
lp_obj <- L_objective(c(1, 2))
 lp_con <- L_constraint(c(1, 1), dir="==", rhs=2)
 lp_bound <- V_bound(ui=1:2, ub=c(3, 3))
 #ui an integer vector specifying the indices of non-standard (i.e., values != Inf) upper bounds.
 # ub up bound
 lp <- OP(objective=lp_obj, constraints=lp_con, bounds=lp_bound, maximum=FALSE)
 bounds(lp)

ROI Variable Bounds:

0 lower and 2 upper non-standard variable bounds.

lp

ROI Optimization Problem:

Minimize a linear objective function with
- 2 objective variables,

subject to
- 1 constraints of type linear.

Solutions

 x <- ROI_solve(lp)
 x$objval

[1] 2

 x$solution

[1] 2 0

Change boundary

# change boundary 
 bounds(lp) <- V_bound(ui=1:2, ub=c(1, 1))

Solutions

 y <- ROI_solve(lp)
 y$objval

[1] 3

 y$solution

[1] 1 1

---

Other optimization package in R

We also can solve optimization problem using other packages.

#install.packages(c("lpSolve", "lpSolveAPI"))
#install.packages(c("quadprog","Rglpk", "nlme"))
suppressMessages(library(lpSolve))

library(lpSolveAPI)
suppressMessages(library(quadprog))
#suppressMessages(library(nlme))

suppressMessages(library(Rglpk))

args(lp)

NULL

args(solve.QP)

function (Dmat, dvec, Amat, bvec, meq = 0, factorized = FALSE) 
NULL

args(Rglpk_solve_LP)

function (obj, mat, dir, rhs, bounds = NULL, types = NULL, max = FALSE, 
    control = list(), ...) 
NULL

#args(nlminb)
#args(optim)
#args(Rcplex)

Linear Programming by lpsolve

Here is a list of some key features of lp_solve:

• Mixed Integer Linear Programming (MILP) solver

• Basically no limit on model size

• It is free and with sources

• Supports Integer variables, Semi-continuous variables and Special Ordered Sets

Example by lpSolveAPI

straightforward

Set up model

library(lpSolveAPI)
#create an empty model
lprec <- make.lp(nrow =0, ncol= 4)
# nrow  : # a nonnegative integer value specifying the number of constaints in the linear program.
# ncol  : # a nonnegative integer value specifying the number of decision variables in the linear program.

Setup objective function

# object
set.objfn(lprec, c(1, 3, 6.24, 0.1))

Add constraints

#constraint
add.constraint(lprec, c(0, 78.26, 0, 2.9), ">=", 92.3)
add.constraint(lprec, c(0.24, 0, 11.31, 0), "<=", 14.8)
add.constraint(lprec, c(12.68, 0, 0.08, 0.9), ">=", 4)

Set bounds

set.bounds(lprec, lower = c(28.6, 18), columns = c(1, 4))
set.bounds(lprec, upper = 48.98, columns = 4)

Setup names

RowNames <- c("THISROW", "THATROW", "LASTROW")
ColNames <- c("COLONE", "COLTWO", "COLTHREE", "COLFOUR")
dimnames(lprec) <- list(RowNames, ColNames)
lprec

Model name: 
            COLONE    COLTWO  COLTHREE   COLFOUR          
Minimize         1         3      6.24       0.1          
THISROW          0     78.26         0       2.9  >=  92.3
THATROW       0.24         0     11.31         0  <=  14.8
LASTROW      12.68         0      0.08       0.9  >=     4
Kind           Std       Std       Std       Std          
Type          Real      Real      Real      Real          
Upper          Inf       Inf       Inf     48.98          
Lower         28.6         0         0        18          

Solve the model

solve(lprec)

[1] 0

  # 
  # 

Extract solutions

get.objective(lprec)

[1] 31.78276

  # [1] 31.78276
  # 
get.variables(lprec)

[1] 28.60000  0.00000  0.00000 31.82759

  # [1] 28.60000  0.00000  0.00000 31.82759
  # 
get.constraints(lprec)

[1]  92.3000   6.8640 391.2928

  # [1]  92.3000   6.8640 391.2928

Example by lpSolve

provide a lot of information including sensitivities analysis.

Setup objective function

library(lpSolve)
#
# Set up problem: maximize
#   x1 + 9 x2 +   x3 subject to
#   x1 + 2 x2 + 3 x3  <= 9
# 3 x1 + 2 x2 + 2 x3 <= 15
#
f.obj <- c(1, 9, 1)

Setup constraints

f.con <- matrix (c(1, 2, 3, 3, 2, 2), nrow=2, byrow=TRUE)
f.dir <- c("<=", "<=")
f.rhs <- c(9, 15)

Solve the model

#
lp ("max", f.obj, f.con, f.dir, f.rhs)

Success: the objective function is 40.5 

##  Success: the objective function is 40.5

Extract solutions

lp ("max", f.obj, f.con, f.dir, f.rhs)$solution

[1] 0.0 4.5 0.0

Get sensitivities

# Get sensitivities
#
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$sens.coef.from

[1] -1e+30  2e+00 -1e+30

##  [1] -1e+30  2e+00 -1e+30
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$sens.coef.to  

[1] 4.50e+00 1.00e+30 1.35e+01

##  [1] 4.50e+00 1.00e+30 1.35e+01
#
# Right now the dual values for the constraints and the variables are
# combined, constraints coming first. So in this example...
#
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals     

[1]   4.5   0.0  -3.5   0.0 -12.5

##  [1]   4.5   0.0  -3.5   0.0 -10.5
#
# ...the duals of the constraints are 4.5 and 0, and of the variables,
# -3.5, 0.0, -10.5. Here are the lower and upper limits on these:
#
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals.from

[1] -1.776357e-15 -1.000000e+30 -1.000000e+30 -1.000000e+30 -6.000000e+00

##  [1]  0e+00 -1e+30 -1e+30 -1e+30 -6e+00
lp ("max", f.obj, f.con, f.dir, f.rhs, compute.sens=TRUE)$duals.to  

[1] 1.5e+01 1.0e+30 3.0e+00 1.0e+30 3.0e+00

##  [1] 1.5e+01 1.0e+30 3.0e+00 1.0e+30 3.0e+00
#

Intergr Programming

# Run again, this time requiring that all three variables be integer
#
lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec=1:3)

Success: the objective function is 37 

##  Success: the objective function is 37
lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec=1:3)$solution

[1] 1 4 0

##  [1] 1 4 0
#
# You can get sensitivities in the integer case, but they're harder to
# interpret.
#
lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec=1:3, compute.sens=TRUE)$duals

[1]  1  0  0  7 -2

##  [1] 1 0 0 7 0
#

Quadratic Programming by quadprog

solving quadratic programming problems of the form min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0.

D is diagnal matix

library(quadprog)

## Assume we want to minimize: -(0 5 0) %*% b + 1/2 b^T D b
## under the constraints:      A^T b >= b0
##          1 0 0
##      D = 0 1 0
##          0 0 1
## with b0 = (-8,2,0)^T
## and      (-4  2  0) 
##      A = (-3  1 -2)
##          ( 0  0  1)
## we can use solve.QP as follows:
##
Dmat       <- matrix(0,3,3)
diag(Dmat) <- 1
dvec       <- c(0,5,0)
Amat       <- matrix(c(-4,-3,0,2,1,0,0,-2,1),3,3)
bvec       <- c(-8,2,0)
solve.QP(Dmat,dvec,Amat,bvec=bvec)

$solution
[1] 0.4761905 1.0476190 2.0952381

$value
[1] -2.380952

$unconstrained.solution
[1] 0 5 0

$iterations
[1] 3 0

$Lagrangian
[1] 0.0000000 0.2380952 2.0952381

$iact
[1] 3 2

Football MIP by glpk

trying to pick the best possible fantasy football team given different constraints. goal is to pick the players that maximize the sum of their projected points.

The constraints are:

1. The team must include:

   • 1 QB

   • 2 RBs

   • 2 WRs

   • 1 TE

2. A player’s risk must not exceed 6

3. The sum of the players’ costs must not exceed 300.

#http://stackoverflow.com/questions/15147398/optimize-value-with-linear-or-non-linear-constraints-in-r



# We are going to solve:
# maximize f'x subject to A*x <dir> b
# where:
#   x is the variable to solve for: a vector of 0 or 1:
#     1 when the player is selected, 0 otherwise,
#   f is your objective vector,
#   A is a matrix, b a vector, and <dir> a vector of "<=", "==", or ">=",
#   defining your linear constraints.


name <- c("Aaron Rodgers","Tom Brady","Arian Foster","Ray Rice","LeSean McCoy","Calvin Johnson","Larry Fitzgerald","Wes Welker","Rob Gronkowski","Jimmy Graham")


# position
pos <- c("QB","QB","RB","RB","RB","WR","WR","WR","TE","TE")

# points
pts <- c(167, 136, 195, 174, 144, 135, 89, 81, 114, 111) 

risk <- c(2.9, 3.4, 0.7, 1.1, 3.5, 5.0, 6.7, 4.7, 3.7, 8.8) 

cost <- c(60, 47, 63, 62, 40, 60, 50, 35, 40, 40) 

#mydata <- data.frame(name, pos, pts, risk, cost) 

# number of variables
num.players <- length(name)



# objective:
f <- pts



# the variable are booleans
var.types <- rep("B", num.players)



# the constraints
A <- rbind(as.numeric(pos == "QB"), # num QB
           as.numeric(pos == "RB"), # num RB
           as.numeric(pos == "WR"), # num WR
           as.numeric(pos == "TE"), # num TE
           diag(risk),              # player's risk, a set of multiple constraints
           cost)                    # total cost

# diretion
dir <- c("==",
         "==",
         "==",
         "==",
         rep("<=", num.players),
         "<=")

# rhs b vector
b <- c(1,
       2,
       2,
       1,
       rep(6, num.players),
       300)



# solve
library(Rglpk)
sol <- Rglpk_solve_LP(obj = f, mat = A, dir = dir, rhs = b,
                      types = var.types, max = TRUE)
sol

$optimum
[1] 836

$solution
 [1] 1 0 1 0 1 1 0 1 1 0

$status
[1] 0

$solution_dual
[1] NA

$auxiliary
$auxiliary$primal
 [1]   1.0   2.0   2.0   1.0   2.9   0.0   0.7   0.0   3.5   5.0   0.0
[12]   4.7   3.7   0.0 298.0

$auxiliary$dual
[1] NA

# $optimum
# [1] 836                      ### <- the optimal total points



# $solution
#  [1] 1 0 1 0 1 1 0 1 1 0     ### <- a `1` for the selected players

# $status
# [1] 0                        ### <- an optimal solution has been found

# your dream team
name[sol$solution == 1]

[1] "Aaron Rodgers"  "Arian Foster"   "LeSean McCoy"   "Calvin Johnson"
[5] "Wes Welker"     "Rob Gronkowski"

Aaron Rodgers

Arian Foster

LeSean McCoy

Calvin Johnson

Wes Welker

Rob Gronkowski

# [1] "Aaron Rodgers"  "Arian Foster"   "LeSean McCoy"
# [4] "Calvin Johnson" "Wes Welker"     "Rob Gronkowski


# total cost
sum(cost * sol$solution)

[1] 298