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

Example 1 with LognormalDist

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

the class ChanceCash method main.

public static void main(String[] args) {
    double iniCash = 150;
    double iniI = 0;
    double trunQuantile = 0.999;
    // the higher value results in slower running speed. maximum negative possibility rate is 1 - serviceRate.
    double serviceRate = 0.6;
    // sample number in each period, the number of samples in the first period can have big influence
    int[] sampleNums = { 5, 5, 5, 5 };
    double[] meanDemand = { 10, 5, 10, 5 };
    int T = sampleNums.length;
    int[] sampleNumsRolling = { 5, 5, 5, 5 };
    // rolling horizon length
    int rollingLength = 2;
    double meanDemandSum = Arrays.stream(meanDemand).sum();
    double rollingDemandSum = Arrays.stream(meanDemand).limit(rollingLength).sum();
    double portion = rollingDemandSum / meanDemandSum;
    double rollingServiceRate = Math.pow(serviceRate, portion);
    // simulating sample number in testing SAA and extended SAA
    int sampleNum = 100;
    double holdCostUnit = 0;
    double salvageValueUnit = 5;
    double[] prices = new double[T];
    double[] variCostUnits = new double[T];
    double[] overheadCosts = new double[T];
    double[] mus = new double[T];
    double sigmasCoe = 0.25;
    Arrays.fill(prices, 22);
    Arrays.fill(variCostUnits, 10);
    // overhead costs
    Arrays.fill(overheadCosts, 100);
    // Arrays.fill(mus, 3.6);
    // Arrays.fill(sigmas, 0.6);
    // maximum ordering quantity when having enough cash
    double maxOrderQuantity = 300;
    Distribution[] distributions = IntStream.iterate(0, i -> i + 1).limit(T).mapToObj(i -> new PoissonDist(meanDemand[i])).toArray(// Poisson demand
    Distribution[]::new);
    // .mapToObj(i -> new LognormalDist(mus[i], sigmas[i])).toArray(Distribution[]::new);
    /**
     * solve the problem by SAA
     */
    // generate scenarios, samples in each period form a scenario tree
    Sampling sampling = new Sampling();
    double[][] scenarios = sampling.generateLHSamples(distributions, sampleNums);
    int sampleNumTotal = IntStream.of(sampleNums).reduce(1, (a, b) -> a * b);
    int sampleNumTotalSimulate = IntStream.of(sampleNumsRolling).limit(rollingLength).reduce(1, (a, b) -> a * b);
    int negativeScenarioNumRequire = (int) (sampleNumTotal * (1 - serviceRate));
    LostSaleChance model = new LostSaleChance(distributions, sampleNums, iniCash, iniI, prices, variCostUnits, salvageValueUnit, holdCostUnit, overheadCosts, serviceRate, scenarios);
    long currTime = System.currentTimeMillis();
    double[] result;
    double time1;
    double positiveScenario;
    double survivalProb;
    double lostRate;
    NumberFormat nf = NumberFormat.getPercentInstance();
    nf.setMinimumFractionDigits(5);
    DecimalFormat df = new DecimalFormat("###, ###");
    result = model.solveMaxSurvival();
    time1 = (System.currentTimeMillis() - currTime) / 1000.00;
    currTime = System.currentTimeMillis();
    System.out.println("**********************************************");
    System.out.println("result of SAA-scenario tree before sorting scenarios: ");
    System.out.println("running time is " + time1 + "s");
    System.out.printf("first stage decison Q is: %.2f\n", result[0]);
    positiveScenario = result[1];
    System.out.printf("Objective value is: %.0f in %d scenarios\n", result[1], sampleNumTotal);
    survivalProb = 100 * result[1] / sampleNumTotal;
    System.out.printf("Survival probability is: %.5f%%\n", survivalProb);
    System.out.println("lost sale scenario number in the solution is : " + result[2]);
    System.out.println("maximum lost sale scenario number allowed is: " + negativeScenarioNumRequire);
    lostRate = result[2] / (double) sampleNumTotal;
    System.out.println("lost sale rate of SAA is: " + nf.format(lostRate));
    System.out.println("lost sale max required rate is: " + nf.format(1 - serviceRate));
    System.out.println();
    /**
     * Simulate the restult of SAA
     */
    int stepSize = 1;
    double fixOrderCost = 0;
    double depositeRate = 0;
    double minInventoryState = 0;
    double maxInventoryState = 500;
    // can affect results, should be smaller than minus fixedOrderCost
    double minCashState = -1000;
    double maxCashState = 2000;
    double discountFactor = 1;
    double[][][] pmf = new GetPmf(distributions, trunQuantile, stepSize).getpmf();
    // immediate value
    ImmediateValueFunction<CashState, Double, Double, Double> immediateValue = (state, action, randomDemand) -> {
        int t = state.getPeriod() - 1;
        double revenue = prices[t] * Math.min(state.getIniInventory() + action, randomDemand);
        double fixedCost = action > 0 ? fixOrderCost : 0;
        double variableCost = variCostUnits[t] * action;
        double deposite = (state.getIniCash() - fixedCost - variableCost) * (1 + depositeRate);
        double inventoryLevel = state.getIniInventory() + action - randomDemand;
        double holdCosts = holdCostUnit * Math.max(inventoryLevel, 0);
        double cashIncrement = revenue + deposite - holdCosts - overheadCosts[t] - state.getIniCash();
        double salValue = state.getPeriod() == T ? salvageValueUnit * Math.max(inventoryLevel, 0) : 0;
        cashIncrement += salValue;
        return cashIncrement;
    };
    // state transition function
    StateTransitionFunction<CashState, Double, Double, CashState> stateTransition = (state, action, randomDemand) -> {
        double nextInventory = Math.max(0, 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;
        // cash is integer or not
        // the right should be a decimal
        nextCash = Math.round(nextCash * 1) / 1;
        return new CashState(state.getPeriod() + 1, nextInventory, nextCash);
    };
    int period = 1;
    CashState initialState = new CashState(period, iniI, iniCash);
    double[] result1;
    // no need to add overheadCost in this class
    CashSimulation simulation1 = new CashSimulation(distributions, sampleNum, immediateValue, stateTransition);
    double error;
    double thisServiceRate;
    result1 = simulation1.simulateSAA(initialState, result[0], serviceRate, sampleNums, prices, variCostUnits, overheadCosts, salvageValueUnit, holdCostUnit, scenarios, sampleNum);
    System.out.println("final simulated survival probability of SAA in " + df.format(sampleNum) + " samples is: " + nf.format(result1[0]));
    error = 1.96 * Math.sqrt(result1[1] * (1 - result1[1]) / sampleNum);
    thisServiceRate = 1 - result1[1];
    System.out.println("final simulated service sale rate of SAA " + " is: " + nf.format(thisServiceRate) + " with error " + nf.format(error));
    /**
     * solve the problem by extended formulation of SAA,
     * sort scenarios in the whole planning horizon
     */
    currTime = System.currentTimeMillis();
    // same result with soveSort or solveSort2, but less computational time
    result = model.solveSortWhole();
    // former name is solveSortFurther()
    time1 = (System.currentTimeMillis() - currTime) / 1000.00;
    currTime = System.currentTimeMillis();
    System.out.println("**********************************************");
    System.out.println("after sorting scenarios in the whole planning horizon, result of SAA-scenario tree: ");
    System.out.println("running time is " + time1 + "s");
    System.out.printf("first stage decison Q is: %.2f\n", result[0]);
    positiveScenario = result[1];
    System.out.printf("Objective value is: %.0f in %d scenarios\n", positiveScenario, sampleNumTotal);
    survivalProb = 100 * result[1] / sampleNumTotal;
    System.out.printf("Survival probability is: %.5f%%\n", survivalProb);
    System.out.println("lost sale scenario number in the solution is : " + result[2]);
    System.out.println("maximum lost sale scenario number allowed is: " + negativeScenarioNumRequire);
    lostRate = result[2] / (double) sampleNumTotal;
    System.out.println("lost sale rate of SAA is: " + nf.format(lostRate));
    System.out.println("lost sale max required rate is: " + nf.format(1 - serviceRate));
    System.out.println();
    /**
     * Simulate the result of extended SAA
     */
    result1 = simulation1.simulateExtendSAAWhole(initialState, result[0], serviceRate, sampleNums, prices, variCostUnits, overheadCosts, salvageValueUnit, holdCostUnit, scenarios, sampleNum);
    System.out.println("final simulated survival probability of extended SAA(sort whole planning horizon) in " + df.format(sampleNum) + " samples is: " + nf.format(result1[0]));
    error = 1.96 * Math.sqrt(result1[1] * (1 - result1[1]) / sampleNum);
    thisServiceRate = 1 - result1[1];
    System.out.println("final simulated service rate of extended SAA(sort whole planning horizon) " + " is: " + nf.format(thisServiceRate) + " with error " + nf.format(error));
    /**
     * solve the problem by extended formulation of SAA,
     * sort scenarios in each period.
     */
    // currTime = System.currentTimeMillis();
    // result = model.solveSortEach();	// same result with soveSort or solveSort2, but less computational time
    // // former name is solveSortFurther()
    // 
    // time1 = (System.currentTimeMillis() - currTime) / 1000.00;
    // currTime = System.currentTimeMillis();
    // System.out.println("**********************************************");
    // System.out.println("after sorting scenarios in each period, the result of SAA-scenario tree: ");
    // System.out.println("running time is " + time1 + "s");
    // System.out.printf("first stage decison Q is: %.2f\n", result[0]);
    // positiveScenario = result[1];
    // System.out.printf("Objective value is: %.0f in %d scenarios\n", positiveScenario, sampleNumTotal);
    // survivalProb = 100 * result[1] / sampleNumTotal;
    // System.out.printf("Survival probability is: %.5f%%\n", survivalProb);
    // System.out.println("lost sale scenario number in the solution is : " + result[2]);
    // System.out.println("maximum lost sale scenario number allowed is: " + negativeScenarioNumRequire);
    // lostRate = result[2] / (double) sampleNumTotal;
    // System.out.println("lost sale rate of SAA is: " + nf.format(lostRate));
    // System.out.println("lost sale max required rate is: " + nf.format(1 - serviceRate));
    // System.out.println();
    // 
    // /**
    // * Simulate the result of extended SAA sorting each period
    // */
    // result1 = simulation1.simulateExtendSAAEach(initialState, result[0], serviceRate, sampleNums, prices, variCostUnits, overheadCosts, salvageValueUnit, holdCostUnit, scenarios, sampleNum);
    // System.out.println("final simulated survival probability of extended SAA sorting each period in " + df.format(sampleNum) + " samples is: " + nf.format(result1[0]));
    // System.out.println("final simulated lost sale rate of extended SAA sorting each period " + " is: " + nf.format(result1[1]));
    /**
     * solve the problem by SDP when there is no joint chance constraint
     */
    // feasible actions
    Function<CashState, double[]> getFeasibleAction1 = s -> {
        int t = s.getPeriod() - 1;
        double thisPeriodServRate = (1 - serviceRate) / T;
        // double minQ = Math.ceil(distributions[t].inverseF(1 - thisPeriodServRate)); // minimum ordering quantity in each period
        double minQ = 0;
        return DoubleStream.iterate(minQ, i -> i + stepSize).limit((int) maxOrderQuantity + 1).toArray();
    };
    /**
     *****************************************************************
     * Solve
     */
    CashRecursion recursion = new CashRecursion(OptDirection.MAX, pmf, getFeasibleAction1, stateTransition, immediateValue, discountFactor);
    currTime = System.currentTimeMillis();
    recursion.setTreeMapCacheAction();
    double finalValue = recursion.getSurvProb(initialState);
    System.out.println("**********************************************");
    System.out.println("result of SDP with no service rate constraint is: ");
    System.out.println("survival probability for this initial state is: " + nf.format(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 sdp results
     */
    sampleNum = 10000;
    // no need to add overheadCost in this class
    CashSimulation simulation = new CashSimulation(distributions, sampleNum, recursion, discountFactor);
    double[] result2 = simulation.simulateSDPGivenSamplNum(initialState, immediateValue);
    System.out.println("final simulated survival probability in " + df.format(sampleNum) + " samples is: " + nf.format(result2[0]));
    System.out.println("final simulated lost sale rate " + " is: " + nf.format(result2[1]));
    /**
     *****************************************************************
     * solve the problem by SDP when there is individual chance constraint approximation
     */
    // feasible actions 2
    // in fact, no cash constraint in this paper
    Function<CashState, double[]> getFeasibleAction2 = s -> {
        int t = s.getPeriod() - 1;
        double thisPeriodServRate = (1 - serviceRate) / T;
        // minimum ordering quantity in each period
        double minQ = Math.ceil(distributions[t].inverseF(1 - thisPeriodServRate));
        return DoubleStream.iterate(minQ, i -> i + stepSize).limit((int) maxOrderQuantity + 1).toArray();
    };
    /**
     *****************************************************************
     * Solve
     */
    recursion = new CashRecursion(OptDirection.MAX, pmf, getFeasibleAction2, stateTransition, immediateValue, discountFactor);
    period = 1;
    initialState = new CashState(period, iniI, iniCash);
    currTime = System.currentTimeMillis();
    recursion.setTreeMapCacheAction();
    finalValue = recursion.getSurvProb(initialState);
    System.out.println("**********************************************");
    System.out.println("result of SDP with service rate constraint is: ");
    finalValue = finalValue * 100;
    System.out.printf("survival probability for this initial state is: %.2f%%\n", finalValue);
    System.out.println("optimal order quantity in the first priod is : " + recursion.getAction(initialState));
    time = (System.currentTimeMillis() - currTime) / 1000;
    System.out.println("running time is " + time + "s");
    /**
     *****************************************************************
     * Simulating sdp results
     */
    sampleNum = 10000;
    // no need to add overheadCost in this class
    simulation = new CashSimulation(distributions, sampleNum, recursion, discountFactor);
    result2 = simulation.simulateSDPGivenSamplNum(initialState, immediateValue);
    System.out.println("final simulated survival probability in " + df.format(sampleNum) + " samples is: " + nf.format(result2[0]));
    System.out.println("final simulated lost sale rate " + " is: " + nf.format(result2[1]));
    /**
     * solve the problem by rolling horizon of further SAA
     */
    // number of scenarios for rolling SAA
    sampleNum = 100;
    currTime = System.currentTimeMillis();
    System.out.println("**********************************************");
    result1 = simulation1.rollingHoirzonFurtherExtendSAA(rollingLength, initialState, rollingServiceRate, sampleNumsRolling, prices, variCostUnits, overheadCosts, salvageValueUnit, holdCostUnit, scenarios, sampleNum);
    time1 = (System.currentTimeMillis() - currTime) / 1000.00;
    System.out.println("after rolling horizon for length " + rollingLength + ", result of SAA-scenario tree: ");
    System.out.println("running time is " + time1 + "s");
    System.out.println("final simulated survival probability of rolling further SAA in " + df.format(sampleNum) + " samples is: " + nf.format(result1[0]));
    System.out.println("final simulated lost sale rate of rolling further SAA " + " is: " + nf.format(result1[1]));
}