spbenders3.gms : Stochastic Benders - Sequential GamsModelInstance

Description

This example demonstrates a stochastic Benders implementation for the
simple transport example.

This is the third example of a sequence of stochastic Benders
implementations using various methods to solve the master and
subproblem.

This third example implements the stochastic Benders algorithm using
the sequential solves of master and subproblem. Here these models
are implemented as Python OO-API GamsModelInstance objects. The
advantage is that the models need to be generated only once and solved
with varying data. Since the model rim of a GamsModelInstance cannot
be changed, the master model includes all possible cuts with non-binding
constraint at the beginning: sum(j, eps*received(j)) =l= bigM. During
the cause of the algorithm the right hand side and the coefficients for
received(j) are updated with the real cut data.

Keywords: linear programming, stochastic Benders algorithm, transportation
          problem, GAMS embedded code facility, Python

  On the major platforms, GMSPYTHONHOME gets set automatically, otherwise the
  user has to set it. This condition can also be removed, if one has set up its
  Python environment appropriately


Reference

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

Small Model of Type : LP


Category : GAMS Model library


Main file : spbenders3.gms

$title Stochastic Benders - Sequential GamsModelInstance (SPBENDERS3,SEQ=420)

$onText
This example demonstrates a stochastic Benders implementation for the
simple transport example.

This is the third example of a sequence of stochastic Benders
implementations using various methods to solve the master and
subproblem.

This third example implements the stochastic Benders algorithm using
the sequential solves of master and subproblem. Here these models
are implemented as Python OO-API GamsModelInstance objects. The
advantage is that the models need to be generated only once and solved
with varying data. Since the model rim of a GamsModelInstance cannot
be changed, the master model includes all possible cuts with non-binding
constraint at the beginning: sum(j, eps*received(j)) =l= bigM. During
the cause of the algorithm the right hand side and the coefficients for
received(j) are updated with the real cut data.

Keywords: linear programming, stochastic Benders algorithm, transportation
          problem, GAMS embedded code facility, Python
$offText

* On the major platforms, GMSPYTHONHOME gets set automatically, otherwise the
* user has to set it. This condition can also be removed, if one has set up its
* Python environment appropriately
$if not setEnv GMSPYTHONHOME
$if %gams.pySetup%==1 $abort.noError Embedded code Python not ready to be used
$if not setEnv GMSPYTHONHOME $log --- Using external Python

Set
   i 'factories'            / f1*f3 /
   j 'distribution centers' / d1*d5 /;

Parameter
   capacity(i) 'unit capacity at factories'
               / f1 500, f2 450, f3 650 /
   demand(j)   'unit demand at distribution centers'
               / d1 160, d2 120, d3 270, d4 325, d5 700 /
   prodcost    'unit production cost'                    / 14 /
   price       'sales price'                             / 24 /
   wastecost   'cost of removal of overstocked products' /  4 /;

Table transcost(i,j) 'unit transportation cost'
          d1    d2    d3    d4    d5
   f1   2.49  5.21  3.76  4.85  2.07
   f2   1.46  2.54  1.83  1.86  4.76
   f3   3.26  3.08  2.60  3.76  4.45;

$ifThen not set useBig
   Set s 'scenarios' / lo, mid, hi /;

   Table ScenarioData(s,*) 'possible outcomes for demand plus probabilities'
             d1   d2   d3   d4   d5  prob
      lo    150  100  250  300  600  0.25
      mid   160  120  270  325  700  0.50
      hi    170  135  300  350  800  0.25;
$else
$  if not set nrScen $set nrScen 10
   Set s 'scenarios' / s1*s%nrScen% /;
   Parameter ScenarioData(s,*) 'possible outcomes for demand plus probabilities';
   option seed = 1234;
   ScenarioData(s,'prob') = 1/card(s);
   ScenarioData(s,j)      = demand(j)*uniform(0.6,1.4);
$endIf

* Benders master problem
$if not set maxiter $set maxiter 25
Set iter            'max Benders iterations' / 1*%maxiter% /;

Alias (iter,it);

Parameter
   cutconst(iter)   'constants in optimality cuts'    / #iter    0 /
   cutcoeff(iter,j) 'coefficients in optimality cuts' / #iter.#j 0 /;

Variable
   ship(i,j)        'shipments'
   product(i)       'production'
   received(j)      'quantity sent to market'
   zmaster          'objective variable of master problem'
   theta            'future profit';

Positive Variable ship;

Equation
   masterobj        'master objective function'
   production(i)    'calculate production in each factory'
   receive(j)       'calculate quantity to be send to markets'
   optcut(iter)     'Benders optimality cuts';

masterobj..
   zmaster =e= theta - sum((i,j), transcost(i,j)*ship(i,j))
                     - sum(i, prodcost*product(i));

receive(j)..    received(j) =e= sum(i, ship(i,j));

production(i).. product(i)  =e= sum(j, ship(i,j));

optcut(iter)..
   theta =l= cutconst(iter) + sum(j, cutcoeff(iter,j)*received(j));

product.up(i) = capacity(i);

Model masterproblem / all /;

* Benders' subproblem
Parameter r(j) 'received from master';

Variable
   sales(j)    'sales (actually sold)'
   waste(j)    'overstocked products'
   zsub        'objective variable of sub problem';

Positive Variable sales, waste;

Equation
   subobj      'subproblem objective function'
   selling(j)  'part of received is sold'
   market(j)   'upperbound on sales';

subobj..     zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j));

selling(j).. sales(j) + waste(j) =e= r(j);

market(j)..  sales(j) =l= demand(j);

Model subproblem / subobj, selling, market /;

* Benders loop
Scalar
   rgap       'relative gap'       /    0 /
   lowerBound 'global lower bound' / -inf /
   upperBound 'global upper bound' / +inf /
   objMaster                       /    0 /
   objSub                          /    0 /;

$set solverlog
$if set useSolverLog $set solverlog output=sys.stdout

embeddedCode Python:
def solveMI(mi, symIn=[], symOut=[]):
  for sym in symIn:
    gams.db[sym].copy_symbol(mi.sync_db[sym])
  mi.solve(%solverlog%)
  for sym in symOut:
    mi.sync_db[sym].copy_symbol(gams.db[sym])
pauseEmbeddedCode
abort$execerror 'Python error. Check the log';

* Initialize cut to be non-binding
cutconst(iter)   = 1e15;
cutcoeff(iter,j) = eps;
r(j)      = 0;
objMaster = 0;

$libInclude pyEmbMI miSub    'subproblem max zsub using lp'       -all_model_types=cplexd  r.Zero demand.Zero
$libInclude pyEmbMI miMaster 'masterproblem max zmaster using lp' -all_model_types=cplexd  cutconst.Accumulate cutcoeff.Accumulate

loop(it,
*  The clear of the cut data below goes together with the Accumulate updateType
*  of the update symbols. It also works without the clear and updateType BaseCase
*  but requires much more data exchange because we communicate in every iteration
*  the data of all cuts generated so far
   option clear = cutconst, clear = cutcoeff;
   objSub = 0;
   cutconst(it) = eps;
   continueEmbeddedCode:
   gams.db['r'].copy_symbol(miSub.sync_db['r'])
   pauseEmbeddedCode
   loop(s,
      demand(j) = scenarioData(s,j);
      continueEmbeddedCode:
      solveMI(miSub,['demand'],['market','selling','zsub'])
      pauseEmbeddedCode market, selling, zsub
      objSub = objSub + ScenarioData(s,'prob')*zsub.l;
      cutconst(it)   = cutconst(it)   + ScenarioData(s,'prob')*sum(j, market.m(j)*demand(j));
      cutcoeff(it,j) = cutcoeff(it,j) + ScenarioData(s,'prob')*selling.m(j);
   );
   if(lowerBound < objMaster + objSub, lowerBound = objMaster + objSub);
   rgap = (upperBound - lowerBound)/(1 + abs(upperBound));
   break$(rgap < 0.001);
   continueEmbeddedCode:
   solveMI(miMaster,['cutconst','cutcoeff'],['received','zmaster','theta'])
   pauseEmbeddedCode received, zmaster, theta
   upperBound = zmaster.l;
   objMaster  = zmaster.l - theta.l;
   r(j) = received.l(j);
);
abort$(rgap >= 0.001) 'need more iterations', lowerbound, upperbound;
display 'optimal solution', lowerbound, upperbound;