Revenue Management¶

Principles learned¶

• Add constraints to a black-box model
• Call an expensive black-box function
• Set bounds for the black-box function
• Add evaluation points for the black-box function

Problem¶ A businessman wants to maximize the income from the sale of a product over a certain time horizon, split into several periods. At the beginning of the time horizon, he has to decide the total amount of product to buy. Then in each period he has to choose the number of units to sell, provided that the total number of units sold during the time horizon does not exceed the initial amount purchased. The price of the product increases over the periods, this means the businessman has to determine the number of units to reserve for customers who arrive later, because they will pay more. The later demand must therefore be considered to make a wise decision at each period. Since the demand for the product is stochastic, the businessman simulates each repartition of units a large number of times to get a robust estimate of his income.

For more details, see: revenue_management.html

Program¶

Parameters are fixed to the recommended values: 3 periods are used and the initial cost of the product is set to \$80.

The demand at each period t is defined by: Dₜ=μₜXYₜ, where:

• Yₜ has an exponential distribution with a rate parameter λ=1.
• X has a gamma distribution with a shape parameter k=1 and a scale parameter θ=1, which is equivalent to the standard exponential distribution.
• μₜ is the mean demand for this period.

The distributions are implemented in each language. The prices and mean demands at each period are in the table below:

Period 1 2 3
Price 100 300 400
Mean demand 50 20 30

To have a robust estimate of the income, the simulation is realized a large number of times (1.000.000) using a Monte Carlo method. Each simulation takes several seconds to run, that’s why the next point to be evaluated has to be wisely chosen. A black-box function is thus used to compute the average income from the simulations.

Three integer decision variables are declared. The first one corresponds to the initial quantity purchased at the beginning of the time horizon. The second determines the amount of product that has to be reserved for periods 2 and 3, and the third the amount of product available for period 3. The domains of these variables are all [0, 100]. To ensure the feasibility, each variable is constrained to be lesser or equal than the previous one.

As the black-box function is provided by the user, LocalSolver cannot compute anything about it. It is then useful to parametrize it via the LSBlackBoxContext. In this example, the simulation will never return a negative value, because the prices at any periods are above the initial cost, and all decision variables are positive. The lower bound of the function is thus set to 0. The maximum number of evaluations is set to 30.

For this simulation, an evaluation point was previously computed and its value was saved: the variables [100, 50, 30] generates a mean revenue of 4740.99. This point is added with a LSBlackBoxEvaluationPoint to warm start the solver.

Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python revenue_management.py
Execution (Linux)
export PYTHONPATH=/opt/localsolver_10_5/bin/python
python revenue_management.py
########## revenue_management.py ##########

import localsolver
import sys
import math
import random

class RevenueManagementFunction:

def __init__(self, seed):
self.nb_periods = 3
self.prices = [100, 300, 400]
self.mean_demands = [50, 20, 30]
self.purchase_price = 80
self.evaluated_points = [{
"point": [100, 50, 30],
"value": 4740.99
}]
self.nb_simulations = int(1e6)
self.seed = seed

# Black-box function
def evaluate(self, argument_values):
variables = [argument_values.get(i) for i in range(argument_values.count())]
# Initial quantity purchased
nb_units_purchased = variables
# Number of units that should be left for future periods
nb_units_reserved = variables[1:] + 

# Sets seed for reproducibility
random.seed(self.seed)
# Creates distribution
X = [gamma_sample() for i in range(self.nb_simulations)]
Y = [[exponential_sample() for i in range(self.nb_periods)]
for j in range(self.nb_simulations)]

# Runs simulations
sum_profit = 0.0
for i in range(self.nb_simulations):
remaining_capacity = nb_units_purchased
for j in range(self.nb_periods):
# Generates demand for period j
demand_j = int(self.mean_demands[j] * X[i] * Y[i][j])
nb_units_sold = min(max(remaining_capacity - nb_units_reserved[j], 0),
demand_j)
remaining_capacity = remaining_capacity - nb_units_sold
sum_profit += self.prices[j] * nb_units_sold

# Calculates mean revenue
mean_profit = sum_profit / self.nb_simulations
mean_revenue = mean_profit - self.purchase_price * nb_units_purchased

return mean_revenue

def exponential_sample(rate_param=1.0):
u = random.random()
return math.log(1 - u) / (-rate_param)

def gamma_sample(scale_param=1.0):
return exponential_sample(scale_param)

def solve(evaluation_limit, time_limit, output_file):
with localsolver.LocalSolver() as ls:
model = ls.model

# Generates data
revenue_management = RevenueManagementFunction(1)
nb_periods = revenue_management.nb_periods
# Declares decision variables
variables = [model.int(0, 100) for i in range(nb_periods)]

# Creates blackbox function
func_expr = model.create_double_blackbox_function(revenue_management.evaluate)
# Calls function
func_call = model.call(func_expr)

# Declares constraints
for i in range(1, nb_periods):
model.constraint(variables[i] <= variables[i-1])

# Maximizes function call
model.maximize(func_call)

# Sets lower bound
context = func_expr.blackbox_context
context.lower_bound = 0.0

model.close()

# Parametrizes the solver
if time_limit is not None:
ls.param.time_limit = time_limit

# Sets the maximum number of evaluations
context.evaluation_limit = evaluation_limit

for evaluated_point in revenue_management.evaluated_points:
evaluation_point = context.create_evaluation_point()
for i in range(nb_periods):
evaluation_point.set_return_value(evaluated_point["value"])

ls.solve()

# Writes the solution in a file
if output_file is not None:
with open(output_file, 'w') as f:
f.write("obj=%f\n" % func_call.value)
f.write("b=%f\n" % variables.value)
for i in range(1, nb_periods):
f.write("r%f=%f\n" % (i+1, variables[i].value))

if __name__ == '__main__':
output_file = sys.argv if len(sys.argv) > 1 else None
time_limit = int(sys.argv) if len(sys.argv) > 2 else None
evaluation_limit = int(sys.argv) if len(sys.argv) > 3 else 30

solve(evaluation_limit, time_limit, output_file)
Compilation / Execution (Windows)
cl /EHsc revenue_management.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver105.lib
revenue_management
Compilation / Execution (Linux)
g++ revenue_management.cpp -I/opt/localsolver_10_5/include -llocalsolver105 -lpthread -o revenue_management
revenue_management
//********* revenue_management.cpp *********

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <cmath>
#include <vector>
#include "localsolver.h"

using namespace localsolver;
using namespace std;

struct EvaluatedPoint {
public:
EvaluatedPoint(vector<int> point, double value) : point(point), value(value) {}
const int getPoint(int index) { return point[index]; }
const double getValue() { return value; }
private:
vector<int> point;
double value;
};

/* Black-box function */
class RevenueManagementFunction : public LSBlackBoxFunction<lsdouble> {
private:
int seed;
const int nbPeriods = 3;
const int purchasePrice = 80;
const int nbSimulations = (int) 1e6;
vector<EvaluatedPoint> evaluatedPoints;

const int prices(int index) {
const int p[] = {100, 300, 400};
return p[index];
}

const int meanDemands(int index) {
const int d[] = {50, 20, 30};
return d[index];
}

double exponentialSample(double rateParam = 1.0) {
double u = (double) rand() / RAND_MAX;
return log(1 - u) / (-rateParam);
}

double gammaSample(double scaleParam = 1.0) {
return exponentialSample(scaleParam);
}

public:
// Constructor
RevenueManagementFunction(int seed) : seed(seed) {
evaluatedPoints.push_back(EvaluatedPoint({100, 50, 30}, 4740.99));
}

const unsigned int getNbPeriods() { return nbPeriods; }
const vector<EvaluatedPoint> getEvaluatedPoints() { return evaluatedPoints; }

lsdouble call(const LSBlackBoxArgumentValues& argumentValues) override {
// Initial quantity purchased
int nbUnitsPurchased = argumentValues.getIntValue(0);
// Number of units that should be left for future periods
vector<int> nbUnitsReserved(nbPeriods, 0);
for (unsigned int j = 0; j < nbPeriods - 1; j++) {
nbUnitsReserved[j] = argumentValues.getIntValue(j+1);
}
// Sets seed for reproducibility
srand(seed);
// Creates distribution
vector<double> X;
for (unsigned int i = 0; i < nbSimulations; i++) {
X.push_back(gammaSample());
}
vector<vector<double>> Y;
for (unsigned int i = 0; i < nbSimulations; i++) {
vector<double> yt;
for (unsigned int j = 0; j < nbPeriods; j++) {
yt.push_back(exponentialSample());
}
Y.push_back(yt);
}

// Runs simulations
double sumProfit = 0;
for (unsigned int i = 0; i < nbSimulations; i++) {
int remainingCapacity = nbUnitsPurchased;
for (unsigned int j = 0; j < nbPeriods; j++) {
// Generates demand for period j
int demand = (int) (meanDemands(j) * X[i] * Y[i][j]);
int nbUnitsSold = min(max(remainingCapacity - nbUnitsReserved[j], 0),
demand);
remainingCapacity = remainingCapacity - nbUnitsSold;
sumProfit += prices(j) * nbUnitsSold;
}
}
// Calculates mean revenue
double meanProfit = sumProfit / nbSimulations;
double meanRevenue = meanProfit - purchasePrice * nbUnitsPurchased;

return meanRevenue;
}
};

class RevenueManagement {
public:
// Solver
LocalSolver localsolver;

// LS Program variables
vector<LSExpression> variables;
LSExpression bbCall;

void solve(int timeLimit, int evaluationLimit) {
// Declares the optimization model
LSModel model = localsolver.getModel();

// Generates data
RevenueManagementFunction revenueManagement(1);
unsigned int nbPeriods = revenueManagement.getNbPeriods();
// Declares decision variables
for (unsigned int i = 0; i < nbPeriods; i++) {
variables.push_back(model.intVar(0, 100));
}

// Creates blackbox function
LSExpression bbFunc = model.createBlackBoxFunction(&revenueManagement);
// Calls function
bbCall = model.call(bbFunc);
for (unsigned int i = 0; i < nbPeriods; i++) {
}

// Declares constraints
for (unsigned int i = 1; i < nbPeriods; i++) {
model.constraint(variables[i] <= variables[i-1]);
}

// Maximizes function call
model.maximize(bbCall);

// Sets lower bound
LSBlackBoxContext context = bbFunc.getBlackBoxContext();
context.setLowerBound(0.0);

model.close();

// Parametrizes the solver
if (timeLimit != 0) {
localsolver.getParam().setTimeLimit(timeLimit);
}

// Sets the maximum number of evaluations
context.setEvaluationLimit(evaluationLimit);

for (EvaluatedPoint evaluatedPoint : revenueManagement.getEvaluatedPoints()) {
LSBlackBoxEvaluationPoint evaluationPoint = context.createEvaluationPoint();
for (int i = 0; i < nbPeriods; i++) {
}
evaluationPoint.setReturnValue(evaluatedPoint.getValue());
}

localsolver.solve();
}

// Writes the solution in a file
void writeSolution(const string& fileName) {
ofstream outfile;
outfile.open(fileName.c_str());
outfile << "obj=" << bbCall.getDoubleValue() << endl;
outfile << "b=" << variables.getIntValue() << endl;
for (unsigned int i = 1; i < variables.size(); i++) {
outfile << "r" << (i+1) << "=" << variables[i].getIntValue() << endl;
}
}
};

int main(int argc, char** argv) {
const char* solFile = argc > 1 ? argv : NULL;
const char* strTimeLimit = argc > 2 ? argv : "0";
const char* strEvaluationLimit = argc > 3 ? argv : "30";

try {
RevenueManagement model;
model.solve(atoi(strTimeLimit), atoi(strEvaluationLimit));
if (solFile != NULL) model.writeSolution(solFile);
} catch (const exception& e) {
cerr << "An error occurred: " << e.what() << endl;
return 1;
}
return 0;
}
Compilation / Execution (Windows)
copy %LS_HOME%\bin\localsolvernet.dll .
csc RevenueManagement.cs /reference:localsolvernet.dll
RevenueManagement
/********** RevenueManagement.cs **********/

using System;
using System.IO;
using System.Collections.Generic;
using localsolver;

public class RevenueManagement : IDisposable
{
public class EvaluatedPoint
{
private int[] point;
private double value;

public EvaluatedPoint(int[] point, double value)
{
this.point = point;
this.value = value;
}

public int GetPoint(int index) {
return point[index];
}

public double GetValue() {
return value;
}
}

/* Black-box function */
public class RevenueManagementFunction
{
private int seed;
private const int nbPeriods = 3;
private const int purchasePrice = 80;
private const int nbSimulations = (int) 1e6;
private readonly int[] prices = {100, 300, 400};
private readonly int[] meanDemands = {50, 20, 30};
private List<EvaluatedPoint> evaluatedPoints = new List<EvaluatedPoint>();

public RevenueManagementFunction(int seed)
{
this.seed = seed;
int[] point = {100, 50, 30};
}

public double Call(LSBlackBoxArgumentValues argumentValues)
{
// Initial quantity purchased
int nbUnitsPurchased = (int) argumentValues.GetIntValue(0);
// Number of units that should be left for future periods
int[] nbUnitsReserved = new int[nbPeriods];
for (int j = 0; j < nbPeriods - 1; j++)
{
nbUnitsReserved[j] = (int) argumentValues.GetIntValue(j+1);
}
nbUnitsReserved[nbPeriods - 1] = 0;
// Sets seed for reproducibility
Random rng = new Random(seed);
// Creates distribution
double[] X = new double[nbSimulations];
for (int i = 0; i < nbSimulations; i++)
{
X[i] = GammaSample(rng);
}
double[,] Y = new double[nbSimulations, nbPeriods];
for (int i = 0; i < nbSimulations; i++)
{
for (int j = 0; j < nbPeriods; j++)
{
Y[i,j] = ExponentialSample(rng);
}
}

// Runs simulations
double sumProfit = 0;
for (int i = 0; i < nbSimulations; i++)
{
int remainingCapacity = nbUnitsPurchased;
for (int j = 0; j < nbPeriods; j++)
{
// Generates demand for period j
int demand = (int) (meanDemands[j] * X[i] * Y[i,j]);
int nbUnitsSold = Math.Min(Math.Max(remainingCapacity - nbUnitsReserved[j],
0), demand);
remainingCapacity = remainingCapacity - nbUnitsSold;
sumProfit += prices[j] * nbUnitsSold;
}
}
// Calculates mean revenue
double meanProfit = sumProfit / nbSimulations;
double meanRevenue = meanProfit - purchasePrice * nbUnitsPurchased;

return meanRevenue;
}

private static double ExponentialSample(Random rng, double rateParam = 1.0)
{
double u = rng.NextDouble();
return Math.Log(1 - u) / (-rateParam);
}

private static double GammaSample(Random rng, double scaleParam = 1.0)
{
return ExponentialSample(rng, scaleParam);
}

public int GetNbPeriods()
{
return nbPeriods;
}

public List<EvaluatedPoint> GetEvaluatedPoints() {
return evaluatedPoints;
}
}

// Solver
private LocalSolver localsolver;

// LS Program variables
private LSExpression[] variables;
private LSExpression bbCall;

public RevenueManagement()
{
localsolver = new LocalSolver();
}

public void Dispose()
{
if (localsolver != null)
localsolver.Dispose();
}

public void Solve(int timeLimit, int evaluationLimit)
{
// Declares the optimization model
LSModel model = localsolver.GetModel();

// Generates data
RevenueManagementFunction revenueManagement = new RevenueManagementFunction(1);
int nbPeriods = revenueManagement.GetNbPeriods();
// Declares decision variables
variables = new LSExpression[nbPeriods];
for (int i = 0; i < nbPeriods; i++)
{
variables[i] = model.Int(0, 100);
}

// Creates blackbox function
LSDoubleBlackBoxFunction bbFunc = new LSDoubleBlackBoxFunction(revenueManagement.Call);
LSExpression bbFuncExpr = model.DoubleBlackBoxFunction(bbFunc);
// Calls function
bbCall = model.Call(bbFuncExpr);
for (int i = 0; i < nbPeriods; i++)
{
}

// Declares constraints
for (int i = 1; i < nbPeriods; i++)
{
model.Constraint(variables[i] <= variables[i-1]);
}

// Maximizes function call
model.Maximize(bbCall);

// Sets lower bound
LSBlackBoxContext context = bbFuncExpr.GetBlackBoxContext();
context.SetLowerBound(0.0);

model.Close();

// Parametrizes the solver
if (timeLimit != 0)
{
localsolver.GetParam().SetTimeLimit(timeLimit);
}

// Sets the maximum number of evaluations
context.SetEvaluationLimit(evaluationLimit);

foreach (EvaluatedPoint evaluatedPoint in revenueManagement.GetEvaluatedPoints())
{
LSBlackBoxEvaluationPoint evaluationPoint = context.CreateEvaluationPoint();
for (int i = 0; i < nbPeriods; i++)
{
}
evaluationPoint.SetReturnValue(evaluatedPoint.GetValue());
}

localsolver.Solve();
}

// Writes the solution in a file
public void WriteSolution(string fileName)
{
using (StreamWriter output = new StreamWriter(fileName))
{
output.WriteLine("obj=" + bbCall.GetDoubleValue());
output.WriteLine("b=" + variables.GetIntValue());
for (int i = 1; i < variables.Length; i++) {
output.WriteLine("r" + i + "=" + variables[i].GetIntValue());
}
}
}

public static void Main(string[] args)
{
string outputFile = args.Length > 0 ? args : null;
string strTimeLimit = args.Length > 1 ? args : "0";
string strEvaluationLimit = args.Length > 2 ? args : "30";

using (RevenueManagement model = new RevenueManagement())
{
model.Solve(int.Parse(strTimeLimit), int.Parse(strEvaluationLimit));
if (outputFile != null)
model.WriteSolution(outputFile);
}
}
}
Compilation / Execution (Windows)
javac RevenueManagement.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. RevenueManagement
Compilation / Execution (Linux)
javac RevenueManagement.java -cp /opt/localsolver_10_5/bin/localsolver.jar
java -cp /opt/localsolver_10_5/bin/localsolver.jar:. RevenueManagement
/********** RevenueManagement.java **********/

import java.io.*;
import java.lang.Math;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;

import localsolver.*;

public class RevenueManagement {

private static class EvaluatedPoint {
private int[] point;
private double value;

public EvaluatedPoint(int[] point, double value) {
this.point = point;
this.value = value;
}
}

/* Black-box function */
private static class RevenueManagementFunction implements LSDoubleBlackBoxFunction {

private int seed;
private int nbPeriods = 3;
private int purchasePrice = 80;
private int nbSimulations = (int) 1e6;
private int[] prices = {100, 300, 400};
private int[] meanDemands = {50, 20, 30};
private List<EvaluatedPoint> evaluatedPoints = new ArrayList<EvaluatedPoint>();

public RevenueManagementFunction(int seed) {
this.seed = seed;
int[] point = {100, 50, 30};
}

@Override
public double call(LSBlackBoxArgumentValues argumentValues) {
// Initial quantity purchased
int nbUnitsPurchased = (int) argumentValues.getIntValue(0);
// Number of units that should be left for future periods
int[] nbUnitsReserved = new int[nbPeriods];
for (int j = 0; j < nbPeriods - 1; j++) {
nbUnitsReserved[j] = (int) argumentValues.getIntValue(j+1);
}
nbUnitsReserved[nbPeriods - 1] = 0;

// Sets seed for reproducibility
Random rng = new Random(seed);
// Creates distribution
double rateParam = 1.0;
double scaleParam = 1.0;
double[] X = new double[nbSimulations];
for (int i = 0; i < nbSimulations; i++) {
X[i] = gammaSample(rng, rateParam);
}
double[][] Y = new double[nbSimulations][nbPeriods];
for (int i = 0; i < nbSimulations; i++) {
for (int j = 0; j < nbPeriods; j++) {
Y[i][j] = exponentialSample(rng, scaleParam);
}
}

// Runs simulations
double sumProfit = 0;
for (int i = 0; i < nbSimulations; i++) {
int remainingCapacity = nbUnitsPurchased;
for (int j = 0; j < nbPeriods; j++) {
// Generates demand for period j
int demand = (int) (meanDemands[j] * X[i] * Y[i][j]);
int nbUnitsSold = Math.min(Math.max(remainingCapacity - nbUnitsReserved[j],
0), demand);
remainingCapacity = remainingCapacity - nbUnitsSold;
sumProfit += prices[j] * nbUnitsSold;
}
}
// Calculates mean revenue
double meanProfit = sumProfit / nbSimulations;
double meanRevenue = meanProfit - purchasePrice * nbUnitsPurchased;

return meanRevenue;
}

private static double exponentialSample(Random rng, double rateParam) {
double u = rng.nextDouble();
return Math.log(1 - u) / (-rateParam);
}

private static double gammaSample(Random rng, double scaleParam) {
return exponentialSample(rng, scaleParam);
}
}

// Solver
private final LocalSolver localsolver;

// LS Program variables
private LSExpression[] variables;
private LSExpression bbCall;

private RevenueManagement(LocalSolver localsolver) {
this.localsolver = localsolver;
}

// Declares the optimization model
private void solve(int timeLimit, int evaluationLimit) {
LSModel model = localsolver.getModel();

// Generates data
RevenueManagementFunction revenueManagement = new RevenueManagementFunction(1);
int nbPeriods = revenueManagement.nbPeriods;
// Declares decision variables
variables = new LSExpression[nbPeriods];
for (int i = 0; i < nbPeriods; i++) {
variables[i] = model.intVar(0, 100);
}

// Creates blackbox function
LSExpression bbFunc = model.doubleBlackBoxFunction(revenueManagement);
// Calls function with operands
bbCall = model.call(bbFunc);
for (int i = 0; i < nbPeriods; i++) {
}

// Declares constraints
for (int i = 1; i < nbPeriods; i++) {
model.constraint(model.leq(variables[i], variables[i-1]));
}

// Maximizes function call
model.maximize(bbCall);

// Sets lower bound
LSBlackBoxContext context = bbFunc.getBlackBoxContext();
context.setLowerBound(0.0);

model.close();

// Parametrizes the solver
if (timeLimit != 0) {
localsolver.getParam().setTimeLimit(timeLimit);
}

// Sets the maximum number of evaluations
context.setEvaluationLimit(evaluationLimit);

for (EvaluatedPoint evaluatedPoint : revenueManagement.evaluatedPoints) {
LSBlackBoxEvaluationPoint evaluationPoint = context.createEvaluationPoint();
for (int i = 0; i < nbPeriods; i++) {
}
evaluationPoint.setReturnValue(evaluatedPoint.value);
}

localsolver.solve();
}

// Writes the solution in a file
private void writeSolution(String fileName) throws IOException {
try (PrintWriter output = new PrintWriter(fileName)) {
output.println("obj=" + bbCall.getDoubleValue());
output.println("b=" + variables.getIntValue());
for (int i = 1; i < variables.length; i++) {
output.println("r" + i + "=" + variables[i].getIntValue());
}
}
}

public static void main(String[] args) {
String outputFile = args.length > 0 ? args : null;
String strTimeLimit = args.length > 1 ? args : "0";
String strEvaluationLimit = args.length > 2 ? args : "30";

try (LocalSolver localsolver = new LocalSolver()) {
RevenueManagement model = new RevenueManagement(localsolver);
model.solve(Integer.parseInt(strTimeLimit), Integer.parseInt(strEvaluationLimit));
if (outputFile != null) {
model.writeSolution(outputFile);
}
} catch (Exception ex) {
System.err.println(ex);
ex.printStackTrace();
System.exit(1);
}
}
}