Revenue Management Problem

import hexaly.optimizer
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 hexaly.optimizer.HexalyOptimizer() as optimizer:
        #
        # Declare the optimization model
        #
        model = optimizer.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 optimizer
        if time_limit is not None:
            optimizer.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"])

        optimizer.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)

#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <vector>

using namespace hexaly;
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 HxExternalFunction<hxdouble> {
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; }

    hxdouble call(const HxExternalArgumentValues& 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:
    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Hexaly Program variables
    vector<HxExpression> variables;
    HxExpression funcCall;

    void solve(int timeLimit, int evaluationLimit) {
        // Declare the optimization model
        HxModel model = optimizer.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
        HxExpression 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
        HxExternalContext context = func.getExternalContext();
        HxSurrogateParameters surrogateParams = context.enableSurrogateModeling();

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

        model.close();

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

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

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

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

using System;
using System.IO;
using System.Collections.Generic;
using Hexaly.Optimizer;

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

    // Hexaly Optimizer
    private HexalyOptimizer optimizer;

    // Hexaly Program variables
    private HxExpression[] variables;
    private HxExpression funcCall;

    public RevenueManagement()
    {
        optimizer = new HexalyOptimizer();
    }

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

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

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

        // Create the function
        HxDoubleExternalFunction func = new HxDoubleExternalFunction(revenueManagement.Call);
        HxExpression 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
        HxExternalContext context = funcExpr.GetExternalContext();
        HxSurrogateParameters surrogateParams = context.EnableSurrogateModeling();

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

        model.Close();

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

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

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

        optimizer.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);
        }
    }
}

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

import com.hexaly.optimizer.*;

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 HxDoubleExternalFunction {

        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(HxExternalArgumentValues 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);
        }
    }

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Hexaly Program variables
    private HxExpression[] variables;
    private HxExpression funcCall;

    private RevenueManagement(HexalyOptimizer optimizer) {
        this.optimizer = optimizer;
    }

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

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

        // Create the function
        HxExpression 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
        HxExternalContext context = func.getExternalContext();
        HxSurrogateParameters surrogateParams = context.enableSurrogateModeling();

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

        model.close();

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

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

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

        optimizer.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 (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            RevenueManagement model = new RevenueManagement(optimizer);
            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);
        }
    }
}