Revenue Management

Principles learned

  • Add constraints to a model with surrogate modeling

  • Call an expensive external function

  • Set bounds for an external function

  • Add evaluation points to a function

Problem

../_images/revenue_management_bb.svg

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.

Download the example


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. Thus, an external function is used to compute the average income from the simulations. To use the surrogate modeling feature, the method enableSurrogateModeling available on the LSExternalContext of the function is called.

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 external function is provided by the user, LocalSolver cannot compute anything about it. It is then useful to parametrize it via the LSExternalContext. 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 thanks to the LSSurrogateParameters.

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 LSEvaluationPoint to warm start the solver.

Execution:
localsolver revenue_management.lsp [evaluationLimit=] [solFileName=]
use io;
use random;

/* Define input data */
function input() {
    nbPeriods = 3;
    prices = { 100, 300, 400 };
    meanDemands = { 50, 20, 30 };
    purchasePrice = 80;
    evaluatedPoints = {{
        "point": { 100, 50, 30 },
        "value": 4740.99
    }};
    nbSimulations = round(1e6);
    seed = 1;

    // Create random module
    rdm = random.create();
}

// External function
function revenueManagement(v1, v2, v3) {
    // Initial quantity purchased
    nbUnitsPurchased = v1;
    // Number of units that should be left for future periods
    nbUnitsReserved = { v2, v3, 0 };

    // Set seed for reproducibility
    rdm.init(seed);
    // Create distribution
    rateParam = 1.0;
    scaleParam = 1.0;
    X[0..nbSimulations-1] = gammaSample(scaleParam);
    Y[0..nbSimulations-1][0..nbPeriods-1] = exponentialSample(rateParam);

    // Run simulations
    sumProfit = 0.0;
    for [i in 0..nbSimulations-1] {
        remainingCapacity = nbUnitsPurchased;
        for [j in 0..nbPeriods-1] {
            // Generate demand for period j
            demandJ = round(meanDemands[j] * X[i] * Y[i][j]);
            nbUnitsSold = min(max(remainingCapacity - nbUnitsReserved[j], 0),
                    demandJ);
            remainingCapacity -= nbUnitsSold;
            sumProfit += prices[j] * nbUnitsSold;
        }
    }

    // Calculate mean revenue
    meanProfit = sumProfit / nbSimulations;
    meanRevenue = meanProfit - purchasePrice * nbUnitsPurchased;

    return meanRevenue;
}

function exponentialSample(rateParam) {
    u = rdm.nextUniform(0, 1);
    return log(1 - u) / (-rateParam);
}

function gammaSample(scaleParam) {
    return exponentialSample(scaleParam);
}

/* Declare the optimization model */
function model() {
    // Declare decision variables
    variables[0..nbPeriods-1] <- int(0, 100);

    // Create the function
    funcExpr <- doubleExternalFunction(revenueManagement);
    // Call function
    funcArgs[0] <- funcExpr;
    funcArgs[i in 1..nbPeriods] <- variables[i-1];
    funcCall <- call[arg in funcArgs](arg);

    // Declare constraints
    for [i in 1..nbPeriods-1]
        constraint variables[i] <= variables[i - 1];

    // Maximize function call
    maximize funcCall;

    // Enable surrogate modeling
    context = funcExpr.context;
    surrogateParams = context.enableSurrogateModeling();

    // Set lower bound
    context.lowerBound = 0.0;
}

function param() {
    // Set the maximum number of evaluations
    if (evaluationLimit == nil) surrogateParams.evaluationLimit = 30;
    else surrogateParams.evaluationLimit = evaluationLimit;

    // Add evaluation points
    for [k in 0..count(evaluatedPoints)-1] {
        evaluationPoint = surrogateParams.createEvaluationPoint();
        for [i in 0..nbPeriods-1]
            evaluationPoint.addArgument(evaluatedPoints[k]["point"][i]);
        evaluationPoint.returnValue = evaluatedPoints[k]["value"];
    }
}

/* Write the solution in a file */
function output() {
    if (solFileName != nil) {
        local solFile = io.openWrite(solFileName);
        solFile.println("obj=", funcCall.value);
        solFile.println("b=", variables[0].value);
        for [i in 1..nbPeriods-1]
            solFile.println("r", i + 1, "=", variables[i].value);
    }
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python
python revenue_management.py
Execution (Linux)
export PYTHONPATH=/opt/localsolver_11_5/bin/python
python 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

    # External 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[0]
        # Number of units that should be left for future periods
        nb_units_reserved = variables[1:] + [0]

        # Set seed for reproducibility
        random.seed(self.seed)
        # Create 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)]

        # Run simulations
        sum_profit = 0.0
        for i in range(self.nb_simulations):
            remaining_capacity = nb_units_purchased
            for j in range(self.nb_periods):
                # Generate 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

        # Calculate 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:
        #
        # Declare the optimization model
        #
        model = ls.model

        # Generate data
        revenue_management = RevenueManagementFunction(1)
        nb_periods = revenue_management.nb_periods
        # Declare decision variables
        variables = [model.int(0, 100) for _ in range(nb_periods)]

        # Create the function
        func_expr = model.create_double_external_function(revenue_management.evaluate)
        # Call function
        func_call = model.call(func_expr)
        func_call.add_operands(variables)

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

        # Maximize function call
        model.maximize(func_call)

        # Enable surrogate modeling
        context = func_expr.external_context
        surrogate_params = context.enable_surrogate_modeling()

        # Set lower bound
        context.lower_bound = 0.0

        model.close()

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

        # Set the maximum number of evaluations
        surrogate_params.evaluation_limit = evaluation_limit

        # Add evaluation points
        for evaluated_point in revenue_management.evaluated_points:
            evaluation_point = surrogate_params.create_evaluation_point()
            for i in range(nb_periods):
                evaluation_point.add_argument(evaluated_point["point"][i])
            evaluation_point.set_return_value(evaluated_point["value"])

        ls.solve()

        # Write 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[0].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[1] if len(sys.argv) > 1 else None
    time_limit = int(sys.argv[2]) if len(sys.argv) > 2 else None
    evaluation_limit = int(sys.argv[3]) 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\localsolver115.lib
revenue_management
Compilation / Execution (Linux)
g++ revenue_management.cpp -I/opt/localsolver_11_5/include -llocalsolver115 -lpthread -o revenue_management
revenue_management
#include "localsolver.h"
#include <cmath>
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <vector>

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;
};

/* External function */
class RevenueManagementFunction : public LSExternalFunction<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:
    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 LSExternalArgumentValues& 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);
        }
        // Set seed for reproducibility
        srand(seed);
        // Create 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);
        }

        // Run simulations
        double sumProfit = 0;
        for (unsigned int i = 0; i < nbSimulations; ++i) {
            int remainingCapacity = nbUnitsPurchased;
            for (unsigned int j = 0; j < nbPeriods; ++j) {
                // Generate 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;
            }
        }
        // Calculate mean revenue
        double meanProfit = sumProfit / nbSimulations;
        double meanRevenue = meanProfit - purchasePrice * nbUnitsPurchased;

        return meanRevenue;
    }
};

class RevenueManagement {
public:
    // LocalSolver
    LocalSolver localsolver;

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

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

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

        // Create the function
        LSExpression func = model.createExternalFunction(&revenueManagement);
        // Call function
        funcCall = model.call(func);
        for (unsigned int i = 0; i < nbPeriods; ++i) {
            funcCall.addOperand(variables[i]);
        }

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

        // Maximize function call
        model.maximize(funcCall);

        // Enable surrogate modeling
        LSExternalContext context = func.getExternalContext();
        LSSurrogateParameters surrogateParams = context.enableSurrogateModeling();

        // Set lower bound
        context.setLowerBound(0.0);

        model.close();

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

        // Set the maximum number of evaluations
        surrogateParams.setEvaluationLimit(evaluationLimit);

        // Add evaluation points
        for (EvaluatedPoint evaluatedPoint : revenueManagement.getEvaluatedPoints()) {
            LSEvaluationPoint evaluationPoint = surrogateParams.createEvaluationPoint();
            for (int i = 0; i < nbPeriods; ++i) {
                evaluationPoint.addArgument((lsint)evaluatedPoint.getPoint(i));
            }
            evaluationPoint.setReturnValue(evaluatedPoint.getValue());
        }

        localsolver.solve();
    }

    /* Write the solution in a file */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());
        outfile << "obj=" << funcCall.getDoubleValue() << endl;
        outfile << "b=" << variables[0].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[1] : NULL;
    const char* strTimeLimit = argc > 2 ? argv[2] : "0";
    const char* strEvaluationLimit = argc > 3 ? argv[3] : "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
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;
        }
    }

    /* External 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 };
            evaluatedPoints.Add(new EvaluatedPoint(point, 4740.99));
        }

        public double Call(LSExternalArgumentValues 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;
            // Set seed for reproducibility
            Random rng = new Random(seed);
            // Create 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);
            }

            // Run simulations
            double sumProfit = 0;
            for (int i = 0; i < nbSimulations; ++i)
            {
                int remainingCapacity = nbUnitsPurchased;
                for (int j = 0; j < nbPeriods; ++j)
                {
                    // Generate 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;
                }
            }
            // Calculate 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;
        }
    }

    // LocalSolver
    private LocalSolver localsolver;

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

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

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

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

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

        // Create the function
        LSDoubleExternalFunction func = new LSDoubleExternalFunction(revenueManagement.Call);
        LSExpression funcExpr = model.DoubleExternalFunction(func);
        // Call function
        funcCall = model.Call(funcExpr);
        for (int i = 0; i < nbPeriods; ++i)
            funcCall.AddOperand(variables[i]);

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

        // Maximize function call
        model.Maximize(funcCall);

        // Enable surrogate modeling
        LSExternalContext context = funcExpr.GetExternalContext();
        LSSurrogateParameters surrogateParams = context.EnableSurrogateModeling();

        // Set lower bound
        context.SetLowerBound(0.0);

        model.Close();

        // Parametrize the solver
        if (timeLimit != 0)
            localsolver.GetParam().SetTimeLimit(timeLimit);

        // Set the maximum number of evaluations
        surrogateParams.SetEvaluationLimit(evaluationLimit);

        // Add evaluation points
        foreach (EvaluatedPoint evaluatedPoint in revenueManagement.GetEvaluatedPoints())
        {
            LSEvaluationPoint evaluationPoint = surrogateParams.CreateEvaluationPoint();
            for (int i = 0; i < nbPeriods; ++i)
                evaluationPoint.AddArgument(evaluatedPoint.GetPoint(i));
            evaluationPoint.SetReturnValue(evaluatedPoint.GetValue());
        }

        localsolver.Solve();
    }

    /* Write the solution in a file */
    public void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            output.WriteLine("obj=" + funcCall.GetDoubleValue());
            output.WriteLine("b=" + variables[0].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[0] : null;
        string strTimeLimit = args.Length > 1 ? args[1] : "0";
        string strEvaluationLimit = args.Length > 2 ? args[2] : "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_11_5/bin/localsolver.jar
java -cp /opt/localsolver_11_5/bin/localsolver.jar:. RevenueManagement
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;
        }
    }

    /* External function */
    private static class RevenueManagementFunction implements LSDoubleExternalFunction {

        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 };
            evaluatedPoints.add(new EvaluatedPoint(point, 4740.99));
        }

        @Override
        public double call(LSExternalArgumentValues 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;

            // Set seed for reproducibility
            Random rng = new Random(seed);
            // Create 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);
                }
            }

            // Run simulations
            double sumProfit = 0;
            for (int i = 0; i < nbSimulations; ++i) {
                int remainingCapacity = nbUnitsPurchased;
                for (int j = 0; j < nbPeriods; ++j) {
                    // Generate 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;
                }
            }
            // Calculate 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);
        }
    }

    // LocalSolver
    private final LocalSolver localsolver;

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

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

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

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

        // Create the function
        LSExpression func = model.doubleExternalFunction(revenueManagement);
        // Call function with operands
        funcCall = model.call(func);
        for (int i = 0; i < nbPeriods; ++i) {
            funcCall.addOperand(variables[i]);
        }

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

        // Maximize function call
        model.maximize(funcCall);

        // Enable surrogate modeling
        LSExternalContext context = func.getExternalContext();
        LSSurrogateParameters surrogateParams = context.enableSurrogateModeling();

        // Set lower bound
        context.setLowerBound(0.0);

        model.close();

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

        // Set the maximum number of evaluations
        surrogateParams.setEvaluationLimit(evaluationLimit);

        // Add evaluation points
        for (EvaluatedPoint evaluatedPoint : revenueManagement.evaluatedPoints) {
            LSEvaluationPoint evaluationPoint = surrogateParams.createEvaluationPoint();
            for (int i = 0; i < nbPeriods; ++i) {
                evaluationPoint.addArgument(evaluatedPoint.point[i]);
            }
            evaluationPoint.setReturnValue(evaluatedPoint.value);
        }

        localsolver.solve();
    }

    /* Write the solution in a file */
    private void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            output.println("obj=" + funcCall.getDoubleValue());
            output.println("b=" + variables[0].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[0] : null;
        String strTimeLimit = args.length > 1 ? args[1] : "0";
        String strEvaluationLimit = args.length > 2 ? args[2] : "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);
        }
    }
}