airsp2.gms : Aircraft Allocation - stochastic optimization with DECIS

Description

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 DECIS:
  - fixed demand
  - expected value
  - deterministic equivalent

DECIS formulation with sampling
  - exact (deterministic equivalent)
  - expected value

Also see models AIRCRAFT and AIRSP.


Large Model of Types : DECIS lp


Category : GAMS Model library


Main file : airsp2.gms

$title Aircraft Allocation - Stochastic Optimization with DECIS (AIRSP2,SEQ=196)

$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 DECIS:
  - fixed demand
  - expected value
  - deterministic equivalent

DECIS formulation with sampling
  - exact (deterministic equivalent)
  - expected value

Also see models AIRCRAFT and AIRSP.


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

Keywords: linear programming, stochastic programming, aircraft allocation, routing
$offText

$if not set decisalg $set decisalg decism

Set
   i 'aircraft types and unassigned passengers' / a, b, c, d      /
   j 'assigned and unassigned routes'           / route-1*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 /
   fixeddemand(j) 'fixed demand (passengers in hundreds)'
                  / route-1 250
                    route-2 120
                    route-3 180
                    route-4  90
                    route-5 600 /
   costbumped(j)  'costs associated with bumping passengers'
                  / route-1 13
                    route-2 13
                    route-3  7
                    route-4  7
                    route-5  1 /;

*---------------------------------------------------------------------
* First we solve the "fixed demand" (deterministic demand) case.
* Notice that Dantzig reports a slightly different optimal solution
* for this problem in table 28-2-II on page 583. The solution
* reported there is slightly infeasible.
*---------------------------------------------------------------------
Positive Variable
   x(i,j)    'number of aircraft type i assigned to route j'
   bumped(j) 'passengers bumped';

Variable z   'objective variable';

Equation
   avail(i)  'aircraft availability constraints'
   demand(j) 'demand constraints'
   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= fixeddemand(j);

Model fixed  / cost, avail, demand /;

solve fixed using lp minimizing z;

Parameter
   results(*,i,j) 'first stage results'
   obj(*);

results('Fixed demand',i,j) = x.l(i,j);
obj('Fixed demand')         = z.l;

*---------------------------------------------------------------------
* Now introduce stochastic demand
*---------------------------------------------------------------------
option limCol = 0, limRow = 0;

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         ;

*---------------------------------------------------------------------
* Solve expected value problem directly. Expected value of demand is
* slightly different from fixed demand.
*---------------------------------------------------------------------
Parameter expdemand(j) 'expected value demand';
expdemand(j) = sum(k, probability(j,k)*stochasticdemand(j,k));

display expdemand, fixeddemand;

Equation expvaldemand(j);

expvaldemand(j).. sum(i, pcap(i,j)*x(i,j)) + bumped(j) =g= expdemand(j);

Model expval / cost, avail, expvaldemand /;

solve expval using lp minimizing z;

results('Expected value problem',i,j) = x.l(i,j);
obj('Expected value problem')         = z.l;

display results;

*---------------------------------------------------------------------
* Solve deterministic equivalent problem
* There are 750 scenario's.
*---------------------------------------------------------------------
Set s 'scenarios' / scen1*scen750 /;

Alias (k,k1,k2,k3,k4,k5);

Set ss(s) 'auxiliary set';

Scalar
   jprob 'joint probabilities'
   cnt   'counter';

Parameter
   scenprob(s)     'scenario probabilities'
   scendemand(s,j) 'demand for route j under scenario s';

cnt = 0;
loop((k1,k2,k3,k4,k5),
   jprob = probability('route-1',k1)*probability('route-2',k2)
         * probability('route-3',k3)*probability('route-4',k4)
         * probability('route-5',k5);
   if(jprob,
      cnt = cnt + 1;

* set ss to current element
      ss(s) = no;
      ss(s)$(ord(s) = cnt) = yes;

      scenprob(ss) = jprob;
      scendemand(ss,'route-1') =  stochasticdemand('route-1',k1);
      scendemand(ss,'route-2') =  stochasticdemand('route-2',k2);
      scendemand(ss,'route-3') =  stochasticdemand('route-3',k3);
      scendemand(ss,'route-4') =  stochasticdemand('route-4',k4);
      scendemand(ss,'route-5') =  stochasticdemand('route-5',k5);
   );
);

Scalar sumprob;
sumprob = sum(s,scenprob(s));

display sumprob;

abort$(abs(sumprob - 1) > 0.001) "Error in probabilities";

Positive Variable scenbumped(s, j) 'passengers bumped under scenario s';

Variable          capacityuse(j)   'intermediate variables to reduce number of nonzeroes';

Equation
   deteqobj
   deteqdemand(s,j)
   defcap(j);

deteqobj .. z =e= sum((i,j), c(i,j)*x(i,j))
               +  sum(s, scenprob(s)*sum(j, costbumped(j)*scenbumped(s,j)));

* deteqdemand(s,j).. sum(i, pcap(i,j)*x(i,j)) + scenbumped(s,j) =g= scendemand(s,j);
* the above formulation is intuitive but repeats sum(i, pcap(i,j)*x(i,j)) for
* every s, we introduce the intermediate variable capacityuse instead.

deteqdemand(s,j).. capacityuse(j) + scenbumped(s,j) =g= scendemand(s,j);

defcap(j)..        capacityuse(j) =e= sum(i, pcap(i,j)*x(i,j));

Model deteq / deteqobj, avail, deteqdemand, defcap /;

solve deteq using lp minimizing z;

results('Deterministic equivalent',i,j) = x.l(i,j);
obj('Deterministic equivalent')         = z.l;

*---------------------------------------------------------------------
* Output stochastic file
*---------------------------------------------------------------------
File stg / MODEL.STG /;
put stg "INDEP DISCRETE"/;
loop((j,k)$probability(j,k), put "RHS demand ",j.tl," ",stochasticdemand(j,k)," PERIOD2 ",probability(j,k)/;);
putClose;

$onText
The file MODEL.STG will look like:

INDEP DISCRETE
RHS demand route-1            200.00 PERIOD2         0.20
RHS demand route-1            220.00 PERIOD2         0.05
.
.
RHS demand route-5            600.00 PERIOD2         0.80
RHS demand route-5            620.00 PERIOD2         0.10
$offText

*---------------------------------------------------------------------
* Define stages
*---------------------------------------------------------------------
x.stage(i,j)    = 1;
bumped.stage(j) = 2;
avail.stage(i)  = 1;
demand.stage(j) = 2;

*---------------------------------------------------------------------
* Output a MINOS option file
*---------------------------------------------------------------------
File mopt / MINOS.SPC /;
put  mopt;
put "begin"/;
put "rows 250"/;
put "columns 250"/;
put "elements 10000"/;
put "end"/;
putClose;

*---------------------------------------------------------------------
* solve the stochastic model
*---------------------------------------------------------------------
option lp = %decisalg%;

solve fixed using lp minimizing z;
results('DECIS default sampling',i,j) = x.l(i,j);
obj('DECIS default sampling')         = z.l;

*---------------------------------------------------------------------
* Let DECIS solve the model exactly
* Stochastic Universe option: 4 "ISTRAT"
*---------------------------------------------------------------------
File decopt / %decisalg%.opt /;
put decopt;
put '4 "ISTRAT"'/;
putClose;

fixed.optFile = 1;
solve fixed using lp minimizing z;

results('DECIS universe problem',i,j) = x.l(i,j);
obj('DECIS universe problem')         = z.l;

*---------------------------------------------------------------------
* Solve expected value problem using DECIS
*---------------------------------------------------------------------
put decopt;
put '1 "ISTRAT"'/;
putClose;

option lp = %decisalg%;

solve fixed using lp minimizing z;

results('DECIS expected value problem',i,j) = x.l(i,j);
obj('DECIS expected value problem')         = z.l;

option  obj:2:0:1;
display obj;

option  results:2:1:2;
display results;