$title Aircraft allocation - stochastic optimization (AIRSP3,SEQ=69)

$onText

The objective of this model is to allocate aircrafts to routes to maximize
the expected profit when traffic demand is uncertain. Several problems
are solved with standard LP solvers and an EMP SP solver:

Also see models AIRCRAFT, AIRSP and AIRSP2.

Dantzig, G B, Chapter 28. In Linear Programming and Extensions.
Princeton University Press, Princeton, New Jersey, 1963.

$offText

set i     aircraft types and unassigned passengers / a, b, c, d /
    j     assigned and unassigned routes           / route-1, route-2, route-3, route-4, route-5 /
;

table  c(i,j)       costs per aircraft (1000s)

          route-1    route-2    route-3     route-4     route-5
 a          18         21         18          16           10
 b                     15         16          14            9
 c                     10                      9            6
 d          17         16         17          15           10
;
table  pcap(i,j)    passenger capacity of aircraft i on route j

         route-1     route-2     route-3     route-4     route-5
 a          16         15          28          23          81
 b                     10          14          15          57
 c                      5                       7          29
 d           9         11          22          17          55
;
parameter aircraft(i)  aircraft availability  / a   10
                                                b   19
                                                c   25
                                                d   15 /;

parameter costbumped(j) costs associated with bumping passengers
       / route-1 13
         route-2 13
         route-3  7
         route-4  7
         route-5  1 /;

set k possible outcomes /k1*k5/;

table  stochasticdemand(j,k)     demand distribution on route j
             k1      k2     k3     k4     k5
 route-1    200     220    250    270    300
 route-2     50     150
 route-3    140     160    180    200    220
 route-4     10      50     80    100    340
 route-5    580     600    620
;

table  probability(j,k)  "probabilities corresponding to sd(j,k)"

              k1     k2     k3     k4     k5
 route-1      .2     .05    .35    .2     .2
 route-2      .3     .7
 route-3      .1     .2     .4     .2     .1
 route-4      .2     .2     .3     .2     .1
 route-5      .1     .8     .1
;

Parameter rdemand(j) random demand initialize with mean;
rdemand(j) = sum(k$probability(j,k), probability(j,k)*stochasticdemand(j,k));

positive variable x(i,j)    number of aircraft type i assigned to route j;
positive variable bumped(j) passengers bumped;
variable          z         objective variable;

equation          avail(i)  aircraft availability constraints;
equation          demand(j) demand constraints;
equation          cost      objective;

cost..       z =e= sum((i,j), c(i,j)*x(i,j)) + sum(j, costbumped(j)*bumped(j));
avail(i)..   sum(j, x(i,j)) =l= aircraft(i);
demand(j)..  sum(i, pcap(i,j)*x(i,j)) + bumped(j) =g= rdemand(j);

model fixed /cost,avail,demand/;

file emp / '%emp.info%' /; put emp '* problem %gams.i%';
loop(j,
  put / 'randvar ' rdemand.tn(j) ' discrete ';
  loop(k$probability(j,k),
    put / probability(j,k) stochasticdemand(j,k));
);
putclose / 'stage 2 demand bumped rdemand';

$eval scenmax card(k)**card(j)

Set s            scenarios / s1*s%scenmax% /;
Parameter
    srep(s,*)       scenario probability / #s.prob 0 /
    s_demand(s,j)   demand realization by scenario
    rep             solution metric report;

Set dictSP / s      .scenario.''
             ''     .opt.     srep
             rdemand.randvar. s_demand /;

* Solve stochastic recourse problem or here-and-now
solve fixed using emp min z scenario dictSP;
rep('here-and-now') = z.l;

* Solve expected value (EV) problem
solve fixed using lp min z;

* Compute expected result of using the EV solution by evaluating first stage decision of EV solution under all scenarios:
x.fx(i,j) = x.l(i,j);
Parameter
    s_cost(s)       cost by scenario;

Set dictMS / s      .scenario.''
             rdemand.param.   s_demand
             z      .level.   s_cost /;

solve fixed using lp min z scenario dictMS;
rep('Expected Value of EV solution') = sum(s, s_cost(s)*srep(s,'prob'));

* Now assume perfect forecast and solve the wait-and-see model
x.lo(i,j) = 0; x.up(i,j) = inf;
solve fixed using lp min z scenario dictMS;
rep('wait-and-see') = sum(s, s_cost(s)*srep(s,'prob'));

* Now calculate the different measure
rep('Expected Value of Perfect Information (EVPI)') = rep('here-and-now') - rep('wait-and-see');
rep('Value of Stochastic Solution (VSS)')           = rep('Expected Value of EV solution') - rep('here-and-now');

option rep:2:0:1 display rep;