Examples with PoissonDist umontreal.ssj.probdist.PoissonDist used on opensource projects

Example 1 with PoissonDist

use of umontreal.ssj.probdist.PoissonDist in project Stochastic-Inventory by RobinChen121.

the class LevelFitsS method main.

public static void main(String[] args) {
    double truncationQuantile = 0.9999;
    double stepSize = 1;
    double minInventory = -500;
    double maxInventory = 1000;
    double fixedOrderingCost = 250;
    double variOrderingCost = 0;
    double penaltyCost = 26;
    // double[] meanDemand = { 20, 40, 60, 40 };
    double holdingCost = 1;
    int maxOrderQuantity = 41;
    boolean isForDrawGy = true;
    // // get demand possibilities for each period
    // int T = meanDemand.length;
    // Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T)
    // //.mapToObj(i -> new PoissonDist(meanDemand[i]))
    // .mapToObj(i -> new NormalDist(meanDemand[i], coeValue * meanDemand[i]))
    // .toArray(Distribution[]::new);
    // double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    // int T = 8;
    // double[] values = {6, 7};
    // double[] probs = {0.95, 0.05};
    // Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T)
    // .mapToObj(i -> new DiscreteDistribution(values, probs, values.length)) // can be changed to other distributions
    // .toArray(DiscreteDistribution[]::new);
    // double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    int T = 4;
    double[][] values = { { 34, 159, 281, 286 }, { 14, 223, 225, 232 }, { 5, 64, 115, 171 }, { 35, 48, 145, 210 } };
    double[][] probs = { { 0.018, 0.888, 0.046, 0.048 }, { 0.028, 0.271, 0.17, 0.531 }, { 0.041, 0.027, 0.889, 0.043 }, { 0.069, 0.008, 0.019, 0.904 } };
    Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T).mapToObj(// can be changed to other distributions
    i -> new DiscreteDistribution(values[i], probs[i], values[i].length)).toArray(DiscreteDistribution[]::new);
    double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    // feasible actions
    Function<State, double[]> getFeasibleAction = s -> {
        double[] feasibleActions = new double[(int) (maxOrderQuantity / stepSize) + 1];
        int index = 0;
        for (double i = 0; i <= maxOrderQuantity; i = i + stepSize) {
            feasibleActions[index] = i;
            index++;
        }
        return feasibleActions;
    };
    // state transition function
    StateTransitionFunction<State, Double, Double, State> stateTransition = (state, action, randomDemand) -> {
        double nextInventory = state.getIniInventory() + action - randomDemand;
        nextInventory = nextInventory > maxInventory ? maxInventory : nextInventory;
        nextInventory = nextInventory < minInventory ? minInventory : nextInventory;
        return new State(state.getPeriod() + 1, nextInventory);
    };
    // immediate value
    ImmediateValueFunction<State, Double, Double, Double> immediateValue = (state, action, randomDemand) -> {
        double fixedCost = 0, variableCost = 0, inventoryLevel = 0, holdingCosts = 0, penaltyCosts = 0;
        fixedCost = action > 0 ? fixedOrderingCost : 0;
        variableCost = variOrderingCost * action;
        inventoryLevel = state.getIniInventory() + action - randomDemand;
        holdingCosts = holdingCost * Math.max(inventoryLevel, 0);
        penaltyCosts = penaltyCost * Math.max(-inventoryLevel, 0);
        double totalCosts = fixedCost + variableCost + holdingCosts + penaltyCosts;
        return totalCosts;
    };
    /**
     *****************************************************************
     * Solve
     */
    Recursion recursion = new Recursion(OptDirection.MIN, pmf, getFeasibleAction, stateTransition, immediateValue);
    int period = 1;
    double iniInventory = 500;
    State initialState = new State(period, iniInventory);
    long currTime = System.currentTimeMillis();
    double opt = recursion.getExpectedValue(initialState);
    System.out.println("final optimal expected value is: " + opt);
    System.out.println("optimal order quantity in the first priod is : " + recursion.getAction(initialState));
    double time = (System.currentTimeMillis() - currTime) / 1000;
    System.out.println("running time is " + time + "s");
    /**
     *****************************************************************
     * Simulating sdp results
     */
    int sampleNum = 10000;
    SimulateFitsS simuation = new SimulateFitsS(distributions, sampleNum, recursion);
    simuation.simulateSDPGivenSamplNum(initialState);
    double error = 0.0001;
    double confidence = 0.95;
    simuation.simulateSDPwithErrorConfidence(initialState, error, confidence);
    /**
     *****************************************************************
     * Fit (s, S) levels
     */
    System.out.println("");
    double[][] optTable = recursion.getOptTable();
    MIPFitsS findsS = new MIPFitsS(maxOrderQuantity, T);
    double[][] optsS = findsS.getSinglesS(optTable);
    System.out.println("single s, S level: " + Arrays.deepToString(optsS));
    optsS = findsS.getTwosS(optTable);
    System.out.println("two s, S level: " + Arrays.deepToString(optsS));
    optsS = findsS.getThreesS(optTable);
    System.out.println("three s, S level: " + Arrays.deepToString(optsS));
    double sim1 = simuation.simulateSinglesS(initialState, optsS, maxOrderQuantity);
    System.out.printf("one level gap is %.2f%%\n", (sim1 - opt) * 100 / opt);
    double sim2 = simuation.simulateTwosS(initialState, optsS, maxOrderQuantity);
    System.out.printf("two level gap is %.2f%%\n", (sim2 - opt) * 100 / opt);
    double sim3 = simuation.simulateThreesS(initialState, optsS, maxOrderQuantity);
    System.out.printf("three level gap is %.2f%%\n", (sim3 - opt) * 100 / opt);
}

Example 2 with PoissonDist

use of umontreal.ssj.probdist.PoissonDist in project Stochastic-Inventory by RobinChen121.

the class TwoLevelFitsSTest method main.

public static void main(String[] args) {
    String headString = "K, v, h, I0, pai, Qmax, DemandPatt, OpValue, Time(sec), simValue, error";
    WriteToCsv.writeToFile("./" + "test_results.csv", headString);
    double[][] demands = { { 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30 }, { 6.6, 9.3, 11.1, 12.9, 16.8, 21.6, 24, 26.4, 31.5, 33.9, 36.3, 40.8, 44.4, 47.1, 48.3, 50.1, 50.1, 48.9, 47.1, 44.4 }, { 45.9, 48.6, 49.8, 49.8, 48.6, 45.9, 42.3, 37.8, 35.4, 33, 30.3, 27.9, 25.5, 23.1, 20.7, 18.3, 14.4, 10.8, 8.1, 6.3 }, { 36.3, 30, 23.7, 21, 23.7, 30, 36.3, 39, 36.3, 30, 23.7, 21, 23.7, 30, 36.3, 30.9, 24.3, 21.3, 26.4, 33 }, { 47.1, 30, 12.9, 6, 12.9, 30, 47.1, 54, 47.1, 30, 12.9, 6, 12.9, 30, 47.1, 29.7, 15, 7.5, 11.4, 29.7 }, { 62.7, 27.3, 9.9, 23.7, 1, 22.8, 32.7, 34.5, 67.2, 6.6, 14.4, 40.8, 3.9, 63.3, 25.5, 45, 53.1, 24.6, 10.2, 49.8 }, { 1, 15.3, 45.6, 140.1, 80.4, 146.7, 133.8, 74.4, 84.3, 108.9, 46.5, 87.9, 66, 27.9, 32.1, 88.5, 161.7, 36, 32.4, 49.8 }, { 14.1, 24.3, 70.8, 118.2, 49.2, 86.1, 152.4, 117.3, 226.2, 208.2, 78.3, 58.5, 96, 33.3, 57.3, 115.5, 17.7, 134.7, 127.8, 179.7 }, { 13.2, 34.8, 79.2, 43.2, 43.8, 59.4, 22.2, 54.9, 61.2, 34.2, 49.5, 95.4, 35.7, 144.6, 160.2, 103.8, 150.6, 86.4, 122.7, 63.6 }, { 14.7, 56.4, 19.2, 83.7, 135.9, 67.2, 66.9, 155.1, 87.3, 164.1, 193.8, 67.2, 64.5, 132, 34.8, 131.4, 132.9, 35.7, 172.5, 152.1 } };
    // 1500, 1000, 500
    double[] K = { 500, 800, 200 };
    double[] v = { 0, 5, 10 };
    // 15, 10, 5
    double[] pai = { 15, 10, 5 };
    // 3, 5, 7
    double[] capacity = { 2, 3, 4 };
    double truncationQuantile = 0.9999;
    double stepSize = 1;
    double minInventory = -300;
    double maxInventory = 800;
    double holdingCost = 1;
    for (int iK = 0; iK < K.length; iK++) {
        for (int iv = 0; iv < v.length; iv++) {
            for (int ipai = 0; ipai < pai.length; ipai++) {
                for (int idemand = 0; idemand < demands.length; idemand++) for (int icapacity = 0; icapacity < capacity.length; icapacity++) {
                    // double[] meanDemand = demands[idemand];
                    double[] meanDemand = { 9, 23, 53, 29 };
                    double fixedOrderingCost = K[iK];
                    double proportionalOrderingCost = v[iv];
                    double penaltyCost = pai[ipai];
                    int maxOrderQuantity = (int) (Math.round(Arrays.stream(meanDemand).sum() / meanDemand.length) * capacity[icapacity]);
                    // get demand possibilities for each period
                    int T = meanDemand.length;
                    Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T).mapToObj(// can be changed to other distributions
                    i -> new PoissonDist(meanDemand[i])).toArray(PoissonDist[]::new);
                    double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
                    // feasible actions
                    Function<State, double[]> getFeasibleAction = s -> {
                        double[] feasibleActions = new double[(int) (maxOrderQuantity / stepSize) + 1];
                        int index = 0;
                        for (double i = 0; i <= maxOrderQuantity; i = i + stepSize) {
                            feasibleActions[index] = i;
                            index++;
                        }
                        return feasibleActions;
                    };
                    // state transition function
                    StateTransitionFunction<State, Double, Double, State> stateTransition = (state, action, randomDemand) -> {
                        double nextInventory = state.getIniInventory() + action - randomDemand;
                        nextInventory = nextInventory > maxInventory ? maxInventory : nextInventory;
                        nextInventory = nextInventory < minInventory ? minInventory : nextInventory;
                        return new State(state.getPeriod() + 1, nextInventory);
                    };
                    // immediate value
                    ImmediateValueFunction<State, Double, Double, Double> immediateValue = (state, action, randomDemand) -> {
                        double fixedCost = 0, variableCost = 0, inventoryLevel = 0, holdingCosts = 0, penaltyCosts = 0;
                        fixedCost = action > 0 ? fixedOrderingCost : 0;
                        variableCost = proportionalOrderingCost * action;
                        inventoryLevel = state.getIniInventory() + action - randomDemand;
                        holdingCosts = holdingCost * Math.max(inventoryLevel, 0);
                        penaltyCosts = penaltyCost * Math.max(-inventoryLevel, 0);
                        double totalCosts = fixedCost + variableCost + holdingCosts + penaltyCosts;
                        return totalCosts;
                    };
                    /**
                     *****************************************************************
                     * Solve
                     */
                    Recursion recursion = new Recursion(OptDirection.MIN, pmf, getFeasibleAction, stateTransition, immediateValue);
                    // using tree map when finding levels for s and S
                    recursion.setTreeMapCacheAction();
                    int period = 1;
                    double iniInventory = 0;
                    State initialState = new State(period, iniInventory);
                    long currTime = System.currentTimeMillis();
                    double finalValue = recursion.getExpectedValue(initialState);
                    System.out.println("final optimal expected value is: " + finalValue);
                    System.out.println("optimal order quantity in the first priod is : " + recursion.getAction(initialState));
                    double time = (System.currentTimeMillis() - currTime) / 1000;
                    System.out.println("running time is " + time + "s");
                    /**
                     *****************************************************************
                     * Simulating two level s S results
                     */
                    int sampleNum = 10000;
                    SimulateFitsS simuation = new SimulateFitsS(distributions, sampleNum, recursion);
                    double[][] optTable = recursion.getOptTable();
                    MIPFitsS findsS = new MIPFitsS(maxOrderQuantity, T);
                    double[][] optsS = findsS.getTwosS(optTable);
                    System.out.println(Arrays.deepToString(optsS));
                    simuation.simulateSDPGivenSamplNum(initialState);
                    double simFinalValue = simuation.simulateTwosS(initialState, optsS, maxOrderQuantity);
                    System.out.printf("Optimality gap is: %.2f%%\n", (simFinalValue - finalValue) / finalValue * 100);
                    String out = fixedOrderingCost + ",\t" + proportionalOrderingCost + ",\t" + holdingCost + ",\t" + iniInventory + ",\t" + penaltyCost + ",\t" + maxOrderQuantity + ",\t" + (idemand + 1) + ",\t" + finalValue + ",\t" + time + ",\t" + simFinalValue + ",\t" + (simFinalValue - finalValue) / finalValue;
                    WriteToCsv.writeToFile("./" + "test_results.csv", out);
                }
            }
        }
    }
}

Example 3 with PoissonDist

use of umontreal.ssj.probdist.PoissonDist in project Stochastic-Inventory by RobinChen121.

the class MultiItemCashLookPolicy method main.

public static void main(String[] args) {
    double[] price = { 10, 5 };
    // higher margin vs lower margin
    double[] variCost = { 4, 2 };
    // initial cash
    double iniCash = 25;
    // initial inventory
    int iniInventory1 = 0;
    int iniInventory2 = 0;
    // higher average demand vs lower average demand
    double[][] demand = { { 6, 6 }, { 8, 8 } };
    // higher variance vs lower variance
    double[] coe = { 0.5, 0.25 };
    double[] salPrice = { 2, 1 };
    // horizon length
    int T = demand[0].length;
    double truncationQuantile = 0.999;
    int stepSize = 1;
    double minCashState = 0;
    double maxCashState = 10000;
    int minInventoryState = 0;
    int maxInventoryState = 200;
    int Qbound = 100;
    double discountFactor = 1;
    double Rmin = 25;
    double Rmax = 80;
    int incre = 2;
    int rowNum = (int) ((Rmax - Rmin) / incre) + 2;
    int row = 0;
    double[][] optResults = new double[rowNum][5];
    for (iniCash = Rmin; iniCash <= Rmax; iniCash = iniCash + incre) {
        // get demand possibilities for each period
        Distribution[][] distributions = new Distribution[demand.length][T];
        for (int t = 0; t < T; t++) {
            for (int i = 0; i < demand.length; i++) {
                distributions[t][i] = new PoissonDist(demand[i][t]);
            }
        }
        double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmfMulti();
        // build action list for two items
        Function<CashStateMulti, ArrayList<Actions>> buildActionList = s -> {
            ArrayList<Actions> actions = new ArrayList<>();
            for (int i = 0; i < Qbound; i++) for (int j = 0; j < Qbound; j++) {
                if (variCost[0] * i + variCost[1] * j < s.getIniCash() + 0.1) {
                    Actions thisAction = new Actions(i, j);
                    actions.add(thisAction);
                }
            }
            return actions;
        };
        // Immediate Value Function
        ImmediateValueFunction<CashStateMulti, Actions, Demands, Double> immediateValue = (IniState, Actions, RandomDemands) -> {
            double action1 = Actions.getFirstAction();
            double action2 = Actions.getSecondAction();
            double demand1 = RandomDemands.getFirstDemand();
            double demand2 = RandomDemands.getSecondDemand();
            double endInventory1 = Math.max(0, IniState.getIniInventory1() + action1 - demand1);
            double endInventory2 = Math.max(0, IniState.getIniInventory2() + action2 - demand2);
            double revenue = price[0] * (IniState.getIniInventory1() + action1 - endInventory1) + price[1] * (IniState.getIniInventory2() + action2 - endInventory2);
            double orderingCosts = variCost[0] * action1 + variCost[1] * action2;
            double salValue = 0;
            if (IniState.getPeriod() == T - 1) {
                salValue = salPrice[0] * endInventory1 + salPrice[1] * endInventory2;
            }
            return revenue - orderingCosts + salValue;
        };
        // State Transition Function
        StateTransitionFunction<CashStateMulti, Actions, Demands, CashStateMulti> stateTransition = (IniState, Actions, RandomDemands) -> {
            double endInventory1 = IniState.getIniInventory1() + Actions.getFirstAction() - RandomDemands.getFirstDemand();
            endInventory1 = Math.max(0, endInventory1);
            double endInventory2 = IniState.getIniInventory2() + Actions.getSecondAction() - RandomDemands.getSecondDemand();
            endInventory2 = Math.max(0, endInventory2);
            double nextCash = IniState.getIniCash() + immediateValue.apply(IniState, Actions, RandomDemands);
            nextCash = nextCash > maxCashState ? maxCashState : nextCash;
            nextCash = nextCash < minCashState ? minCashState : nextCash;
            endInventory1 = endInventory1 > maxInventoryState ? maxInventoryState : endInventory1;
            endInventory2 = endInventory2 < minInventoryState ? minInventoryState : endInventory2;
            // rounding states to save computing time
            nextCash = (int) nextCash;
            endInventory1 = (int) endInventory1;
            endInventory2 = (int) endInventory2;
            return new CashStateMulti(IniState.getPeriod() + 1, endInventory1, endInventory2, nextCash);
        };
        /**
         *****************************************************************
         * Solve
         */
        CashRecursionMulti recursion = new CashRecursionMulti(discountFactor, pmf, buildActionList, stateTransition, immediateValue, T);
        int period = 1;
        CashStateMulti iniState = new CashStateMulti(period, iniInventory1, iniInventory2, iniCash);
        long currTime = System.currentTimeMillis();
        double finalValue = iniCash + recursion.getExpectedValueMulti(iniState);
        System.out.println("final optimal cash  is " + finalValue);
        System.out.println("optimal order quantity in the first priod is :  Q1 = " + recursion.getAction(iniState).getFirstAction() + ", Q2 = " + recursion.getAction(iniState).getSecondAction());
        double time = (System.currentTimeMillis() - currTime) / 1000;
        System.out.println("running time is " + time + "s");
        /**
         *****************************************************************
         * Simulating sdp results
         *
         * simulating results a little lower than SDP
         */
        int sampleNum = 10000;
        CashSimulationMulti simuation = new CashSimulationMulti(sampleNum, distributions, discountFactor, recursion, stateTransition, immediateValue);
        double simFinalValue = simuation.simulateSDPGivenSamplNumMulti(iniState);
        System.out.println(simFinalValue);
        /**
         *****************************************************************
         * try to find some ordering patters from optTable
         *
         * output results to excel
         */
        double Q1 = recursion.getAction(iniState).getFirstAction();
        double Q2 = recursion.getAction(iniState).getSecondAction();
        System.out.println("");
        optResults[row][0] = iniInventory1;
        optResults[row][1] = iniInventory2;
        optResults[row][2] = iniCash;
        optResults[row][3] = Q1;
        optResults[row][4] = Q2;
        row++;
    }
    System.out.println("**************************************************");
    WriteToExcel wr = new WriteToExcel();
    String fileName = "optTable2.xls";
    String headString = "x1" + "\t" + "x2" + "\t" + "R" + "\t" + "Q1" + "\t" + "Q2";
    wr.writeArrayToExcel(optResults, fileName, headString);
}

Example 4 with PoissonDist

use of umontreal.ssj.probdist.PoissonDist in project Stochastic-Inventory by RobinChen121.

the class MipCashConstraint method findsCSPieceWise.

/**
 * this mip model is for lost sale, lost sale quantity is not a decision variable;
 * order-up-to level is a decision variable in this model.
 *
 * I^+ is represented by piecewise model
 *
 * revise the above model, but it seems the results of the two linear models are very similar
 *
 * @return s, C, S
 */
public double[][] findsCSPieceWise() {
    int T = distributions.length;
    double[] varz = null;
    double[] varS;
    double[] varx = null;
    double[] varxplus = null;
    double[] varR = null;
    int M = 1000000;
    try {
        IloCplex cplex = new IloCplex();
        // no cplex logging information
        cplex.setOut(null);
        // parameter values in array
        double[] K = new double[T];
        double[] h = new double[T];
        double[] v = new double[T];
        double[] p = new double[T];
        Arrays.fill(K, fixOrderCost);
        Arrays.fill(h, holdingCost);
        Arrays.fill(v, variCost);
        Arrays.fill(p, price);
        // decision variables
        // whether ordering in period t
        IloIntVar[] z = cplex.boolVarArray(T);
        // whether the ordering cycle is from period i to period j
        IloNumVar[][] P = new IloNumVar[T][T];
        for (int i = 0; i < P.length; i++) for (int j = 0; j < P.length; j++) P[i][j] = cplex.boolVar();
        // whether it is lost sale in period t
        IloIntVar[] delta = cplex.boolVarArray(T);
        // order-up-to level in period t
        IloNumVar[] S = cplex.numVarArray(T, 0.0, Double.MAX_VALUE);
        // inventory
        IloNumVar[] x = cplex.numVarArray(T, -Double.MAX_VALUE, Double.MAX_VALUE);
        IloNumVar[] xplus = cplex.numVarArray(T, 0.0, Double.MAX_VALUE);
        // may be useless
        IloNumVar[] xminus = cplex.numVarArray(T, 0.0, Double.MAX_VALUE);
        // end-of-period cash in each period
        IloNumVar[] R = cplex.numVarArray(T, 0.0, Double.MAX_VALUE);
        // objective function
        IloNumExpr finalCash = cplex.sum(cplex.prod(salvageValue, xplus[T - 1]), R[T - 1]);
        cplex.addMaximize(finalCash);
        IloNumExpr realDemand = cplex.numExpr();
        IloNumExpr tRevenue = cplex.numExpr();
        IloNumExpr tHoldCost = cplex.numExpr();
        IloNumExpr tPurchaseCost = cplex.numExpr();
        IloNumExpr tFixCost = cplex.numExpr();
        IloNumExpr tTotalCost = cplex.numExpr();
        IloNumExpr tTotalOrderCost = cplex.numExpr();
        int segmentNum = 5;
        double[] probs = new double[segmentNum];
        double segmentNumFloat = (float) segmentNum;
        double prob = (float) 1.0 / segmentNumFloat;
        Arrays.fill(probs, prob);
        int nbSamples = 10000;
        // sum lambdas from period i to period j
        double[][] sumLambdas = new double[T][T];
        for (int i = 0; i < T; i++) for (int j = 0; j < T; j++) {
            if (i > j)
                sumLambdas[i][j] = 0;
            else {
                for (int k = i; k <= j; k++) sumLambdas[i][j] += distributions[k].getMean();
            }
        }
        for (int t = 0; t < T; t++) {
            realDemand = cplex.diff(S[t], xplus[t]);
            tRevenue = cplex.prod(p[t], realDemand);
            tHoldCost = cplex.prod(h[t], xplus[t]);
            tFixCost = cplex.prod(K[t], z[t]);
            // sum Pjt == 1
            IloLinearNumExpr sumPjt;
            IloLinearNumExpr sumzjt;
            IloLinearNumExpr sumzjt2;
            sumPjt = cplex.linearNumExpr();
            for (int j = 0; j <= t; j++) sumPjt.addTerm(1, P[j][t]);
            // upper triangle
            cplex.addEq(sumPjt, 1);
            for (int j = t + 1; j < T; j++) // other Pjt = 0, or else cannot output values
            cplex.addEq(P[j][t], 0);
            // Pjt >= z_j - sum_{k = j+1}^{t}z_k
            // sum_{j=1}{t}z_j == 0 => P[0][t] == 1, this constraints are important for the extra piecewise constraints
            sumzjt2 = cplex.linearNumExpr();
            for (int j = 0; j <= t; j++) {
                sumzjt = cplex.linearNumExpr();
                for (int k = j + 1; k <= t; k++) sumzjt.addTerm(z[k], 1);
                cplex.addGe(P[j][t], cplex.diff(z[j], sumzjt));
                sumzjt2.addTerm(z[j], 1);
            }
            cplex.addGe(sumzjt2, cplex.diff(1, P[0][t]));
            // piecewise constraints
            // in each period t, get its conditional expectations for S_{jt}, j = 1, 2, ..., t
            long[] seed = { 1, 2, 3, 4, 5, 6 };
            double[][] conditionExpects = new double[t + 1][];
            IloLinearNumExpr sumdP = cplex.linearNumExpr();
            for (int j = 0; j <= t; j++) {
                Distribution[] distributions2 = new Distribution[1];
                distributions2[0] = new PoissonDist(sumLambdas[j][t]);
                PiecewiseComplementaryFirstOrderLossFunction pwcfolf = new PiecewiseComplementaryFirstOrderLossFunction(distributions2, seed);
                conditionExpects[j] = pwcfolf.getConditionalExpectations(probs, nbSamples);
                // System.out.println("max error is: " + pwcfolf.getMaxApproximationError(probs, nbSamples));
                sumdP.addTerm(P[j][t], sumLambdas[j][t]);
            }
            IloLinearNumExpr sumPrdP = cplex.linearNumExpr();
            IloNumExpr xsumdP = cplex.linearNumExpr();
            for (int i = 0; i < segmentNum; i++) {
                double slope = 0;
                double interval = 0;
                for (int k = 0; k <= i; k++) {
                    slope += probs[k];
                }
                xsumdP = cplex.prod(cplex.sum(x[t], sumdP), slope);
                for (int j = 0; j <= t; j++) {
                    for (int k = 0; k <= i; k++) {
                        interval -= probs[k] * conditionExpects[j][k];
                    }
                    sumPrdP.addTerm(interval, P[j][t]);
                }
                cplex.addGe(xplus[t], cplex.sum(xsumdP, sumPrdP));
            }
            // inventory flow: d_t = S_t - I_t
            cplex.addEq(distributions[t].getMean(), cplex.diff(S[t], x[t]));
            // relation of inventory: x_t = x_t^+ - x_t^-
            cplex.addEq(x[t], cplex.diff(xplus[t], xminus[t]));
            if (t == 0) {
                tPurchaseCost = cplex.prod(v[t], cplex.diff(S[t], iniInventory));
                tTotalOrderCost = cplex.sum(tPurchaseCost, tFixCost);
                tTotalCost = cplex.sum(tHoldCost, cplex.sum(tTotalOrderCost, overheadCost));
                // cash flow: R_{t} = R_{0} + p_t(S_t - x^+_t) - hx^+_t - v(S_t - x^+_{t-1}) - Kz_t
                cplex.addEq(cplex.diff(R[t], iniCash), cplex.diff(tRevenue, tTotalCost));
                // cash constraint: R_{0} >= Kz_t + v(S_t - x^+_{t-1})
                cplex.addLe(cplex.sum(overheadCost, tTotalOrderCost), iniCash);
                // z_t = 1 if S_1 - x^+_0 > 0
                // S_1 - x^+_0 >= 0
                // S_1 - x^+_0 <= z_t*M
                cplex.addGe(cplex.diff(S[t], iniInventory), 0);
                cplex.addLe(cplex.diff(S[t], iniInventory), cplex.prod(z[t], M));
            } else {
                tPurchaseCost = cplex.prod(v[t], cplex.diff(S[t], xplus[t - 1]));
                tTotalOrderCost = cplex.sum(tPurchaseCost, tFixCost);
                tTotalCost = cplex.sum(tHoldCost, cplex.sum(tTotalOrderCost, overheadCost));
                // cash flow: R_{t} = R_{0} + p_t(S_t - x^+_t) - hx^+_t - v(S_t - x^+_{t-1}) - Kz_t
                cplex.addEq(cplex.diff(R[t], R[t - 1]), cplex.diff(tRevenue, tTotalCost));
                // cash constraint: R_{0} >= Kz_t + v(S_t - x^+_{t-1})
                cplex.addLe(cplex.sum(overheadCost, tTotalOrderCost), R[t - 1]);
                // z_t = 1 if S_t - x^+_{t-1} > 0
                // S_t - x^+_{t-1} >= 0
                // S_t - x^+_{t-1} <= z_t*M
                cplex.addGe(cplex.diff(S[t], xplus[t - 1]), 0);
                cplex.addLe(cplex.diff(S[t], xplus[t - 1]), cplex.prod(z[t], 10000));
            }
            // x^+_t = max{S_t - d_t, 0}
            // x^+_t <= delta_t * M
            // S_t - d_t <= x^+_t + (1-delta_t)*M
            // S_t - d_t >= x^+_t - (1-delta_t)*M
            // s_t-I_{t-1} <= x_t*10000
            cplex.addLe(xplus[t], cplex.prod(delta[t], 10000));
            cplex.addLe(cplex.diff(S[t], distributions[t].getMean()), cplex.sum(xplus[t], cplex.prod(10000, cplex.diff(1, delta[t]))));
            cplex.addGe(cplex.diff(S[t], distributions[t].getMean()), cplex.diff(xplus[t], cplex.prod(10000, cplex.diff(1, delta[t]))));
        }
        if (cplex.solve()) {
            cplex.output().println("Solution status = " + cplex.getStatus());
            cplex.output().println("Solution value = " + cplex.getObjValue());
            System.out.println(cplex.getStatus());
            System.out.println("Solution value = " + cplex.getObjValue());
            varz = cplex.getValues(z);
            varS = cplex.getValues(S);
            varxplus = cplex.getValues(xplus);
            varx = cplex.getValues(x);
            varR = cplex.getValues(R);
            double[][] varP = new double[T][T];
            for (int i = 0; i < T; i++) for (int j = 0; j < T; j++) varP[i][j] = cplex.getValue(P[i][j]);
            System.out.println("z = ");
            System.out.println(Arrays.toString(varz));
            System.out.println("order-up-to level = ");
            System.out.println(Arrays.toString(varS));
            System.out.println("xplus = ");
            System.out.println(Arrays.toString(varxplus));
            System.out.println("x = ");
            System.out.println(Arrays.toString(varx));
            System.out.println("R = ");
            System.out.println(Arrays.toString(varR));
            System.out.println("P = ");
            System.out.println(Arrays.deepToString(varP));
        }
        cplex.end();
    } catch (IloException e) {
        System.err.println("Concert exception '" + e + "' caught");
    }
    return heuristicFindsCS(varz, varxplus, varR);
}

Example 5 with PoissonDist

use of umontreal.ssj.probdist.PoissonDist in project Stochastic-Inventory by RobinChen121.

the class CheckFG method main.

public static void main(String[] args) {
    double[] meanDemand = { 8, 8, 8 };
    double iniCash = 13;
    double iniInventory = 0;
    double fixOrderCost = 10;
    double variCost = 1;
    double price = 8;
    double salvageValue = 0.5;
    double holdingCost = 0;
    // minimum cash balance the retailer can withstand
    double overheadCost = 0;
    // maximum ordering quantity when having enough cash
    double maxOrderQuantity = 200;
    double truncationQuantile = 0.9999;
    int stepSize = 1;
    double minInventoryState = 0;
    double maxInventoryState = 500;
    // can affect results, should be smaller than minus fixedOrderCost
    double minCashState = -100;
    double maxCashState = 2000;
    double discountFactor = 1;
    boolean isForDrawGy = true;
    // get demand possibilities for each period
    int T = meanDemand.length;
    Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T).mapToObj(i -> new PoissonDist(meanDemand[i])).toArray(Distribution[]::new);
    double[][][] pmf = new GetPmf(distributions, truncationQuantile, stepSize).getpmf();
    // feasible actions
    Function<CashState, double[]> getFeasibleAction = s -> {
        double maxQ = (int) Math.min(maxOrderQuantity, Math.max(0, (s.getIniCash() - overheadCost - fixOrderCost) / variCost));
        return DoubleStream.iterate(0, i -> i + stepSize).limit((int) maxQ + 1).toArray();
    };
    // immediate value
    ImmediateValueFunction<CashState, Double, Double, Double> immediateValue = (state, action, randomDemand) -> {
        double revenue = 0;
        double fixedCost = 0;
        double variableCost = 0;
        double inventoryLevel = 0;
        revenue = price * Math.min(state.getIniInventory() + action, randomDemand);
        fixedCost = action > 0 ? fixOrderCost : 0;
        variableCost = variCost * action;
        inventoryLevel = state.getIniInventory() + action - randomDemand;
        double holdCosts = holdingCost * Math.max(inventoryLevel, 0);
        double cashIncrement = revenue - fixedCost - variableCost - holdCosts;
        double salValue = state.getPeriod() == T ? salvageValue * Math.max(inventoryLevel, 0) : 0;
        cashIncrement += salValue;
        return cashIncrement;
    };
    // state transition function
    StateTransitionFunction<CashState, Double, Double, CashState> stateTransition = (state, action, randomDemand) -> {
        double nextInventory = state.getIniInventory() + action - randomDemand;
        double nextCash = state.getIniCash() + immediateValue.apply(state, action, randomDemand);
        nextCash = nextCash > maxCashState ? maxCashState : nextCash;
        nextCash = nextCash < minCashState ? minCashState : nextCash;
        nextInventory = nextInventory > maxInventoryState ? maxInventoryState : nextInventory;
        nextInventory = nextInventory < minInventoryState ? minInventoryState : nextInventory;
        // nextCash = Math.round(nextCash * 100) / 100.00;
        return new CashState(state.getPeriod() + 1, nextInventory, nextCash);
    };
    /**
     *****************************************************************
     * Solve F(x, R)
     */
    int minInventorys = 0;
    // for drawing pictures
    int maxInventorys = 100;
    int minCash = 0;
    int maxCash = (int) fixOrderCost + 80;
    int RLength = maxCash - minCash + 1;
    int xLength = maxInventorys - minInventorys + 1;
    int period = 1;
    int index = 0;
    int rowIndex = 0;
    int columnIndex = 0;
    long currTime = System.currentTimeMillis();
    CashRecursion recursion = new CashRecursion(OptDirection.MAX, pmf, getFeasibleAction, stateTransition, immediateValue, discountFactor);
    double[][] yG = new double[xLength * RLength][3];
    double[][] resultTableF = new double[RLength][xLength];
    double[][] resultTableQ = new double[RLength][xLength];
    for (double initialCash = minCash; initialCash <= maxCash; initialCash++) {
        columnIndex = 0;
        for (double initialInventory = minInventorys; initialInventory <= maxInventorys; initialInventory++) {
            // initialInventory
            yG[index][0] = initialCash;
            yG[index][1] = initialInventory;
            // iniCash
            yG[index][2] = recursion.getExpectedValue(new CashState(period, initialInventory, initialCash));
            resultTableF[rowIndex][columnIndex] = yG[index][2];
            resultTableQ[rowIndex][columnIndex] = recursion.getAction(new CashState(period, initialInventory, initialCash));
            index++;
            columnIndex++;
        }
        rowIndex++;
    }
    // immediate value for GB and GA
    ImmediateValueFunction<CashState, Double, Double, Double> immediateValue2 = (state, action, randomDemand) -> {
        double revenue = 0;
        double fixedCost = 0;
        double variableCost = 0;
        double inventoryLevel = 0;
        if (isForDrawGy == true && state.getPeriod() == 1) {
            revenue = price * Math.min(state.getIniInventory(), randomDemand);
            fixedCost = 0;
            // variCost * state.getIniInventory(); // be careful
            variableCost = 0;
            inventoryLevel = state.getIniInventory() - randomDemand;
        } else {
            revenue = price * Math.min(state.getIniInventory() + action, randomDemand);
            fixedCost = action > 0 ? fixOrderCost : 0;
            variableCost = variCost * action;
            inventoryLevel = state.getIniInventory() + action - randomDemand;
        }
        double holdCosts = holdingCost * Math.max(inventoryLevel, 0);
        double cashIncrement = revenue - fixedCost - variableCost - holdCosts - overheadCost;
        double salValue = state.getPeriod() == T ? salvageValue * Math.max(inventoryLevel, 0) : 0;
        cashIncrement += salValue;
        return cashIncrement;
    };
    // state transition function for GA
    StateTransitionFunction<CashState, Double, Double, CashState> stateTransition3 = (state, action, randomDemand) -> {
        double nextInventory = isForDrawGy && state.getPeriod() == 1 ? state.getIniInventory() - randomDemand : state.getIniInventory() + action - randomDemand;
        double nextCash = state.getIniCash() + immediateValue2.apply(state, action, randomDemand);
        nextCash = nextCash > maxCashState ? maxCashState : nextCash;
        nextCash = nextCash < minCashState ? minCashState : nextCash;
        nextInventory = nextInventory > maxInventoryState ? maxInventoryState : nextInventory;
        nextInventory = nextInventory < minInventoryState ? minInventoryState : nextInventory;
        return new CashState(state.getPeriod() + 1, nextInventory, nextCash);
    };
    CashRecursion recursion3 = new CashRecursion(OptDirection.MAX, pmf, getFeasibleAction, stateTransition3, immediateValue2, discountFactor);
    double[][] yG3 = new double[xLength * RLength][3];
    index = 0;
    double[][] resultTableGA = new double[xLength][RLength];
    rowIndex = 0;
    for (double initialInventory = minInventoryState; initialInventory <= maxInventorys; initialInventory++) {
        columnIndex = 0;
        for (double initialCash = minCash; initialCash <= maxCash; initialCash++) {
            // initialInventory
            yG3[index][0] = initialInventory;
            // initialCash
            yG3[index][1] = initialCash;
            yG3[index][2] = recursion3.getExpectedValue(new CashState(period, initialInventory, initialCash)) - variCost * initialInventory;
            // careful, minus cy
            resultTableGA[rowIndex][columnIndex] = yG3[index][2];
            index++;
            columnIndex++;
        }
        rowIndex++;
    }
    double[][] resultH = recursion3.getH(resultTableGA, minCash, fixOrderCost, variCost);
    System.out.println("Single-crossing property: " + recursion3.checkSingleCrossing(resultH));
    WriteToCsv wr = new WriteToCsv();
    wr.writeArrayCSVLabel(resultTableQ, minCash, minInventorys, "Q.csv");
    wr.writeArrayCSVLabel(resultTableF, minCash, minInventorys, "F.csv");
    wr.writeArrayCSVLabel(resultTableGA, minCash, minInventorys, "G.csv");
    wr.writeArrayCSVLabel(resultH, minCash, minInventorys, "H.csv");
    // wr.writeArrayCSV(recursion3.getH3Column(resultTableGA, minCash, fixOrderCost, variCost), "G3column.csv");
    double time = (System.currentTimeMillis() - currTime) / 1000;
    System.out.println("running time is " + time + "s");
}