phase2: Function to solve programs such as (24)
In cmottet/GLP: GLP

Description Usage Arguments Value Examples

View source: R/OptimizationFunctions.R

Solves programs such as (24) using the generalized linear programming approach described in Algorithm 1. We slightly generalize this algorithm so that inequalities in (24) can be replaced by lower inequalities or equalities if desired.

1
2
3

phase2(initBFS, objFun, constFun, constRHS, constDir, constLambda = rep(0,
  length(constRHS)), objLambda = 0, C = 10000, IterMax = 100,
  err = 1e-06)

`initBFS`	List containing a boolean feasible indicating wether an initial feasible was found for (24), and if so, initBFS should also contain two vectors p and x representing a feasible distribution function.
`objFun`	Function H
`constFun`	List containing the functions G_j
`constRHS`	Vector containing the bounds gamma_j
`constDir`	Vector containing the direction of the constraints (e.g. '=', '<=', ">=" )
`constLambda`	Vector containing the the values of the parameters λ_{j,M} (use default vector of 0's to solve classic generalized linear programs)
`objLambda`	Scalar containing the value of λ_M (use default 0 to solve classic generalized linear programs)
`C`	Upper bound of the support of feasible distribution functions (default 1e4)
`IterMax`	Maximum number of iterations for the procedure (default = 100)
`err`	Tolerance of the algorithm (default 1e-6)

A list containing

`p`	Vector containing the point masses of the optimal distribution function when the program is feasible and such distribution exists
`x`	Vector containing the point supports of the optimal distribution function when the program is feasible and such distribution exists
`s`	Optimal value of the variable s when the program is feasible and an optimal solution exists
`lpdual`	Vector containing the optimal dual multipliers when the program is feasible and such an optimal solution exists
`bound`	Optimal objective value
`status`	Integer describing the final status of the procedure (0 => Reached a solution, 1 => the algorithm terminated by reaching the maximum number of iterations, 2 => the algorithm entered a cycle)
`nIter`	Number of iterations reached when the procedure terminated
`eps`	Opposite value of the inner optimization program when the procedure terminated
`lastx`	Value minimizing the inner optimization program when the procedure terminated

####
#### Finding a the optimal probability measure (p,x) to the problem
####
#### max P(X > d)
#### s.t. sum(p)  = 1
####      sum(px) = 1
####      sum(px^2) = 2
####
#### where c is the 90th percentile of a standard exponential distribution
#### Note that the solution to this problem is known (see Theorem 3.3 of Bertsimas and Popescu)

# Function and parameters for the integral of the objective
rate <- 1
d <- qexp(0.9,rate)
objFun <- function(x) return(as.numeric(d <= x))

# Function for the integrals of the constraints inequality
constFun = rep(list(
function(x) 1,
function(x) x,
function(x) x^2
),2)

# Direction of the inequality constraints
constDir <- rep(c("<=", ">="), each = 3)

# Bounds on the inequalities
constRHS <- rep(c(1,1,2), 2)

# Values on the RHS of each inequality
# here we choose the moment of order 0, 1, and 2 of an exponential distribution
rate <- 1
mu0 <- 1
mu1 <- 1/rate
mu2 <- 2/rate^2

# Lambdas for the objective function and the constraints functions
constLambda <- rep(c(0,0,0),2)
objLambda <- 0

# Get a basic feasible solution
initBFS  <-  phase1(constFun, constRHS,constDir, constLambda)

# Check feasibility
with(initBFS, c(sum(p), sum(p*x), sum(p*x^2)))

# Solve the optimization program of interest
output <- phase2(initBFS,
                 objFun,
                 constFun,
                 constRHS,
                 constDir,
                 constLambda,
                 objLambda,
                 C = 1000,
                 err = 5e-10)

# Check that the output matches the analytical solution
CMsquare <- (mu2 - mu1^2)/mu1^2
delta <-  (d/mu1-1)

data.frame(Algorithm = output$bound, Analytical = CMsquare/(CMsquare + delta^2))