Lindo Systems

! A very simple stochastic programming problem  (NewsVendScen);
! The newsvendor problem. Decide how much to stock when we
  have uncertain demand;
! Keywords: Charthisto, Newsboy problem, Normal distribution,
            Stochastic optimization, QRAND function, Quasi-random ;
sets:
 scene: u, v, x, cost, d, unif; ! Set of scenarios;
endsets
data:
 p = .6;  ! Penalty cost/unit for stocking too little;
 h = .3;  ! Cost/unit of holding too much;
 mu = 90; ! Mean demand;
 sd = 20; ! Standard deviation. To avoid negative demand,
             should have mu >> 3*sd;
 scene = 1..250;  ! Number of scenarios;
! Generate a vector of quasi-random uniforms with arbitrary seed;
 unif = @QRAND(1923436); 
 nbins = 15;  ! Number of bins to use in the histogram,
               0 means solver decides;
enddata

submodel findopt:
! Random parameters:
    unif( s) = uniform r.v, used to generate..
    d(s) = demand in scenario s,
  Decision variables,
    stage 0:
      x0 = amount to stock,
    stage 1:
      u(s) = amount under ~ short,
      v(s) = amount over;

 ! minimize average cost of under + over;
 min = obj;
   obj = @sum(scene(s): cost(s))/ns;
! For each demand scenario, compute Over and Under;
 @for(scene(s):
     cost(s) = (p*u(s) + h*v(s));
   ! Under - Over = demand - amount stocked;
     u(s) - v(s) = d(s) - x0;
      );
endsubmodel

calc:
  @SET( 'TERSEO',2);    ! Output level (0:verb, 1:terse, 2:only errors, 3:none);

  ns = @size( scene); ! Number of scenarios;

 !  Generate a vector of random variables with expected mean and sd;
 ! Click on: Edit | Paste Function | Distributions to
   see all available distributions;
  @FOR( scene(s):
     d( s) =  @PNORMINV( mu, sd, unif( s));
       );

  @solve( findopt);


  @write( 'Recommended amount to stock= ', @format( x0, '8.2f'),
             ' (Based on mean demand of: ', @format( mu, '8.2f'),')',  @newline(1)); 
  @write( 'Expected shortage + overage cost= ', @format( obj, '8.2f'), @newline(1)); 

!  Create a histogram with a specified number of bins. To see all available chart types
  click on:  Edit | Paste Function | Charting;
  @CHARTHISTO( 'Histogram of objective', 'Objective value', 'Frequency', 'legend', nbins, cost);
endcalc