Example Problems and Code

A collection of sample problems

We describe a collection of 5 small multistage stochastic programming test problems that are available with the IBM Stochastic Solutions. Both the formulation of, and test data for, each problem is described. The SMPS input files for the problems are available in the Stochastic Solutions product directory .\oslstoch\samples\sampc\

Problem KandW3

Statement   Formulation   Data   Solution

This problem has been modified from the production planning example in the text book by Kall and Wallace.

Problem statement

A refinery can blend N raw materials into M different products. At present, the management is trying to decide how much of each of the raw materials to purchase and stock, so that they can be blended to satisfy the demand for the products in future T time periods. The demand has to be completely satisfied, and in case of raw material shortage the products can be outsourced at a higher cost. There is an inventory constraint on how much raw material can be stocked in total.

Formulation

Indices:
 
i  for raw materials (i = 1,...,N).
j  for products (j = 1,...,M).
t  for time periods (t = 1,...,T).
 

Variables:
xit  Quantity of raw material i purchased for use in period t .
yjt  Quantity of product j outsourced in period t .
 

Parameters:
ci  Present cost of raw material i .
aij  Raw material i required per unit of product j .
fjt  Cost of outsourcing product j in time period t .
b  Inventory capacity.
Random demand of product j in time period t .
 

KandW3 :
 
 
 

Test problem data

 
 

The scenario tree for the demand realizations is shown in Scenario tree for problem KandW3.

 
 
 
Figure 3.1: Scenario tree for problem KandW3
    

 

Solution

                                  
 

Problem PROD_MIX

Statement     Formulation     Data     Solution

This problem is taken from Examples from the Literature.

Problem statement

The problem consists of determining the product mix for a furniture shop with two workstations: carpentry and finishing. The availability of labor in man-hours at the two stations is limited. There are four product classes, each consuming a certain number of man-hours at the two stations. Each product earns a certain profit and the shop has the option to purchase labor from outside. The objective is to maximize the profit.

Formulation

Indices:
 
i  for product class (i = 1,...4).
j  for workstation (j=1,2 ).
 

Variables:
xi  Amount of product class j produced.
vj  Man-hours purchased from outside for workstation j .
 

Parameters:
ci  Profit per unit of product class i .
qj  Cost per man-hour for labor for workstation j .
Random man-hours at work station j required per unit of product class i .
Random man-hours available at workstation j .
 

PROD_MIX:

 
 
Note that the problem is in the form of a minimization problem in the data files.

Test problem data

 
 
Above   denotes a normal distribution with mean a and standard deviation b , and denotes a uniform distribution in the range [a,b] .

Solution

The solution of the above problem for a sample of 300 scenarios is as follows:

 

Note that the objective value is negative since the problem was supplied as a minimization problem.

Problem LandS

Statement     Formulation     Data     Solution

The following problem is taken from Louveaux and Smeers.

Problem statement

The following two-stage problem consists of determining the optimal capacity investment in various types of power plants so as to meet next period demands for electricity. Four power plants are considered and they can operate in three different modes. The next period demand for each of the three modes are to be met. There is a budget constraint and also a constraint on the minimum total capacity.

Formulation

Indices:
i  for power plant type (i = 1,...,4).
j  for operating mode (j = 1,...,3).
 

Variables:
xi  Installed capacity for power plant of type i .
yij  Operating level of power plant i in mode j .
 

Parameters:
ci  Cost per unit of capacity installed for plant type i .
fij  Cost per unit operating level of power plant i in mode j .
m  The minimum total capacity to be installed.
b  Budget available for capacity installment.
Random demand for power in mode j .
 

LandS:
      
 

Test problem data

                                                    
 
 

where    is a random variable taking the values 3.0,5.0 or 7.0 with probability 0.3, 0.4 and 0.3 respectively.

Solution

 

Problem HYDRO

Statement     Formulation     Data     Solution

The following problem is taken from Pereira and Pinto.

Problem statement

The objective of the problem is to determine a generation schedule for a hydrothermal generating system such that expected operating costs are minimized. These costs are composed of fuel costs for the thermal units and penalties for failure to meet power demand. The hydrothermal system consists of four hydroelectric reservoirs and a thermal generation unit. Constraints to the problem are the mass balance equations for the water supply across the reservoirs and across the stages; the demand constraints and the capacity of the operating units. The schedule has to be developed for three time periods.

Formulation

Indices:
i  for the reservoirs in the system (i = 1,...4).
i'  an alias for i .
t  for the time periods (t = 1,...3).
 

Variable:
GTt  Thermal generation in time period t .
Dt  Unmet power demand in time period t .
Vit  Water volume at reservoir i at the end of period t .
Qit  Water volume passed through reservoir i turbines in period t .
Sit  Water volume spilled from reservoir i in period t .
 

Parameters:
Mi  The set of reservoirs immediately upstream of reservoir i .
ct  Cost of thermal generation in period t .
ft  Penalty for unmet demand in period t .
Lt  Power demand in period t .
 ri Power generated per unit volume turbined at reservoir i .
Vi0  Initial volume in reservoir i .
Capacity of reservoir i .
Turbine capacity at reservoir i .
Random water inflow to reservoir i in period t .
 

HYDRO:
 

Test problem data

The reservoir configuration is shown in Reservoir configuration for problem HYDRO.

The scenario tree for the water inflow vectors is given in Scenario tree for problem HYDRO.

 
Figure 3.2: Reservoir configuration for problem HYDRO
                         
 
 
Figure 3.3: Scenario tree for problem HYDRO
    

 

Solution

The first stage solutions are reported below. Note that these solutions do not conform to those reported in the original paper Pereira and Pinto.


Problem CAP

Statement     Formulation     Data     Solution

Problem statement

The following problem consists of designing the capacities of the plants in a chemical processing network in order to maximize expected future profits. The network contains three processes and two chemical products, and a planning horizon of two periods is considered. The profit is the difference between the sales revenue of the chemicals and the capacity installation and operating costs for the plants.

formulation

Indices:
 i  for processes (i = 1,...,3).
 j  for chemicals (j=1,2 ).
 t  for time periods (t=1,2 ).
 t an alias for t .
 

Variables:
xi0  Design capacity of process i .
xit  Capacity installed for process i in time period t .
wit  Operating level of process i in period t .
 

Parameters:
ci0  Cost per unit of design capacity for process i .
cit  Cost per unit of capacity installed in process i in period t .
qit  Operating cost of process i in period t .
rjt  Sales price for chemical j in period t .
aij  Production of chemical j per unit operating level of process i .
Demand upper and lower bounds for chemical j in period t .
 

CAP:
 
 
 

 
Figure 3.4: Scenario tree for problem CAP
    

 

Test problem data

The demand lower bounds are assumed to be zero. The scenario tree for the demand upper bounds is given in Scenario tree for problem CAP.

Solution


Note that the objective value is negative since the problem was recast into the minimization form.
Shabbir Ahmed:
Operations Research Laboratory, Department of Mechanical and Industrial Engineering, University of Illinois, Urbana-Champaign.

Guide to the Driver Source Codes

The driver programs in the following sections illustrate the use of some of the options for input, solution and data access available in the IBM Stochastic Solution user callable modules. The C-language coding for the examples may be found in the ./examples directory. This code is reproduced in the Sample Driver Coding section.

All the input examples create an SPL-file <file>.spl, where the string <file> is input from the command line (except sgbd, which creates a file called sgbd.spl). The solver examples restore the stochastic program from <file>.spl, then proceed to solve the problem. For example, to solve the problem app0110 found in the ./data directory in SMPS format, execute the commands: 

    > exsmps data/app0110
    > exsolv data/app0110

Driver illustrating Tree Construction Subroutines

Driver exsgbdt is an example based on the aircraft allocation example of exsgbds. In this example, the distribution is sampled and the tree built up using ekks_AddToTree and ekks_GetScenarioFromTree subroutines. This is only a one-period example, nevertheless, multiperiod examples can easily be constructed on the basis of this code. The exsgbdt call produces the SP file sgbd_smpl.spl.

Driver Illustrating SMPS Input

The driver exsmps reads the SMPS input format. To use the driver, invoke it as follows: 
   > exsmps data/app0110
The application name app0110 is one of the applications found in the ./data directory. The smps call produces the SP file ../data/app0110.spl.

Direct Solution by exsolv

Driver exsolv restores the stochastic program from the SPL-file <file>.spl and solves it using decomposition. For instance, 
   > exsolv data/app0110
solves the problem generated by the invocation of exsmps in the example given above.

Post Solution Processing

Driver exsdst restores the stochastic program from the SP-file <file>.spl and solves the problem using decomposition. The optimal value along each scenario is then computed by using the user-access routines ekks_GetCoreData, ekks_GetScenarioTree and ekks_GetNodeData, and the results are displayed as a distribution. For instance, 
   > exsdst data/app0110
solves the problem generated by the invocation of smps in the example given above, and produces the distribution of objective values.

Sample Driver Coding

exsmps.c 

Read SMPS formatted files and save SP file. 
/*  
 *
 * exsmps.c:  Main driver for reading SMPS files.
 *
 *          Usage:   exsmps <file>
 *
 *          Reads CORE, TIME, and STOCH files from
 *          <file>.core, <file>.time, and <file>.stoch.
 *          Writes the SP-file as <file>.spl.
 *
 *          Implements INDEPENDENT DISCRETE, as well as SCENARIOS
 *          versions of the standard input format.
 *
 *          Demonstrates coding for taking SAMPLES from the independent
 *          distribution.
 *
 *         (compile w/ /DOSLMSDLL on INTEL platforms.
 */

#include <stdlib.h>
#include <string.h>
#include "ekkstoch.h"

#define SAMPLE 0                        /* (0,1) means (don't, do) sample */

#define MXSCN 400

void main( long argc, char *argv[] ) {

    Stoch_Data *s_data=NULL;                    /* Stoch_Data pointer */
    Tree *tree=NULL;
    long type;                                  /* distribution type */
    long replace;                               /* set by ekks_ReadMPS() */
    long rc=0;                                  /* return code */


/* filenames from argument list */
    char *kspc,*corefile,*timefile,*stocfile,*splfile;
    long nsize,nchar;
    long len;

/* get application name from command line */
    if ( argc != 2 ) {
      printf(" Usage:  smps <file>  , where "
             "{file.core,file.time,file.stoch} is a valid "
             "SMPS triple\n");
      return;
    }
    nchar = (long) strlen( argv[1] );

/*
 *  make file names:
 *      "app".core
 *      "app".time
 *      "app".stoch
 */
    nsize = 4*(nchar+1)+(5+5+6+3+5);
    kspc = (char *) malloc(sizeof(char)*nsize);
    if ( kspc==NULL ) {
      printf(" ERROR: smps is unable to store input file-names\n");
      return;
    }
    corefile = kspc; kspc += (nchar+1)+5;
     strcpy( corefile,argv[1] ); strcat( corefile,".core" );
    timefile = kspc; kspc += (nchar+1)+5;
     strcpy( timefile,argv[1] ); strcat( timefile,".time" );
    stocfile = kspc; kspc += (nchar+1)+6;
     strcpy( stocfile,argv[1] ); strcat( stocfile,".stoch" );
    splfile  = kspc; kspc += (nchar+1)+3;
     strcpy( splfile,argv[1] ); strcat( splfile,".spl" );

/* Initialize OSL environments */
    rc = ekks_InitializeOSL();           if(rc>0) return;

/* initialize Stoch_Data with maximum MXSCN scenarios */
    ekks_InitializeApplication( &rc,&s_data,MXSCN );               
        if ( rc>1 ) return;

/* Process MPS file + Process STOCH file + Process TIME file  */
    ekks_ReadMPS( &rc,s_data,&type,&replace, corefile, timefile, stocfile);
    if ( rc>1 ) return;

/* Process discrete distribution (if that's what it is) */
    if ( type==INDEP_DISCRETE || type==BLOCK_DISCRETE ){

        Independent *s_indp = NULL;
        long *indx;

/* ...get pointer to independent distribution structure */
        ekks_GetPointerToDistribution( &rc,s_data,&s_indp );          
                if ( rc>1 ) return;

/* ...malloc indexing arrays */
        nsize = (s_indp->nrvar) * sizeof(long);
        indx = (long *) malloc( nsize );        if (!indx) return;

/* ...loop for sampling independent distributions: */
        if (1==SAMPLE) {
                long jj,irand,ev_indx;
                long ns;
                double *dprobs;
                double dp,ds;
                Discrete *pdist;
                long nindp = s_indp->nrvar;

                ds = 1.0e0;
                ds /= MXSCN;
                for (ns=0;ns < MXSCN;ns++) {

/*   ...for each indep rv, process random number to get index */
                   for (jj=0;jj<nindp;jj++) {
                           ev_indx = 0;
                           pdist = s_indp->rvar[jj]->dist;
                           dprobs = pdist->dprb;
                           dp = dprobs[ev_indx];
                           irand = rand();
                           while ( irand > dp*RAND_MAX )
                                   dp += dprobs[++ev_indx];
                           indx[jj] = ev_indx;
                   }

/*   ...add event */
                   ekks_AddEvent( &rc,s_data,&tree,indx,ds,replace );
                        if ( rc>1 ) return;

                }

        } /* end if(SAMPLE) */

/* ...if NOT SAMPLE then we add all events: */
        if (0==SAMPLE) {
           long jj,ii,iii;
           long ns;
           double dp;
           long *incr;
           Discrete *pdist;

/* ... malloc increments array */
           incr = (long *) malloc( nsize );     if (!incr) return;

/* ...initialize the event increment to 1 and event index to 0 */
           for (jj=0;jj<s_indp->nrvar;jj++) incr[jj] = 1;
           memset(indx,0,nsize);

/* ...compute probability of first event */
           dp = 1.0;
           for (jj=0;jj<s_indp->nrvar;jj++) {
                pdist = s_indp->rvar[jj]->dist;
                dp *= pdist->dprb[ indx[jj] ];
           }

/* ...add first event */
                   ekks_AddEvent( &rc,s_data,&tree,indx,dp,replace );
                if ( rc>1 ) return;

/* ...main loop
        The logic of this loop is as follows:

        - If nevent[jj] divides the scenario index ii,
            reverse the increment direction incr[jj],
            and increase the random variable index jj by 1.
        - Increment the jj'th random variable by incr[jj]
            and call ekks_AddScenario with the new event.

        This works because the total number of events is
        the product of the number for each independent RV.
*/

           ns = 1;
           for (ii=0;ii<s_indp->nrvar;ii++) {
                pdist = s_indp->rvar[ii]->dist;
                ns *= pdist->nevt;
           }

           for (ii=1;ii<ns;ii++) {
                iii=ii; jj=0;
                pdist = s_indp->rvar[jj]->dist;
                while ( !(iii%pdist->nevt) ) {
                        iii /= pdist->nevt;
                        incr[jj] = -incr[jj];
                        jj++;
                        pdist = s_indp->rvar[jj]->dist;
                }

/* ...  adjust probability */
                dp /= pdist->dprb[ indx[jj] ];
                indx[jj] += incr[jj];
                dp *= pdist->dprb[ indx[jj] ];

/* ...  add event */
                   ekks_AddEvent( &rc,s_data,&tree,indx,dp,replace );
                                                if ( rc>1 ) return;

           } /* end for */
        } /* end if(!SAMPLE) */

    } /* end if ( type==INDEP_DISCRETE || type==BLOCK_DISCRETE ) */

/* Save SP-file */
    ekks_PutStochasticData( &rc,s_data,splfile );

/* Free Stoch_Data structures */
    ekks_FreeStructure(s_data);
}

exsgbds.c 

Generate aircraft allocation problem using ekks_AddScenario. 
/*  
 *
 * exsgbds.c:  Generic driver using ekks_AddScenario.
 *
 *         Implements G.B. Dantzig's Aircraft Allocation model,
 *         in: G.B. Dantzig, Linear Programming and Extensions
 *         (Princeton U. Press, 1963) pp. 572-597.
 *
 *         See also: A.J. King, "Stochastic programming problems:
 *         examples from the literature",
 *          in: Yu. Ermoliev and R.J-B Wets, eds, Numerical Techniques for
 *         Stochastic Optimization (Springer-Verlag, 1988) pp. 543-567.
 *         (compile w/  /DOSLMSDLL on INTEL platforms.
 */

#include  <stdlib.h>
#include "ekkstoch.h"

#define INF 1.0e31

void main( long argc, char *argv[] ) {

    Stoch_Data *s_data;                     /* Stoch_Data pointer */

/* Model dimensions */
    long ncol=27, nrow=9, nels=44;

/* Sparse matrix data...organized by row */
    long mrow[]={ 1, 1, 1, 1, 1,
            2, 2, 2, 2,
            3, 3, 3,
            4, 4, 4, 4, 4,
            5, 5, 5, 5,
            6, 6, 6, 6, 6, 6,
            7, 7, 7, 7, 7,
            8, 8, 8, 8, 8, 8,
            9, 9, 9, 9, 9, 9 };
    long mcol[]={ 1, 2, 3, 4, 5,
            6, 7, 8, 9,
            10, 11, 12,
            13, 14, 15, 16, 17,
            1,        13, 18, 19,
            2, 6, 10, 14, 20, 21,
            3, 7,     15, 22, 23,
            4, 8, 11, 16, 24, 25,
            5, 9, 12, 17, 26, 27 };
    double dels[] = { 1.0, 1.0, 1.0, 1.0, 1.0,
                1.0, 1.0, 1.0, 1.0,
                1.0, 1.0, 1.0,
                1.0, 1.0, 1.0, 1.0, 1.0,
                16.0,              9.0, -1.0, 1.0,
                15.0, 10.0,  5.0, 11.0, -1.0, 1.0,
                28.0, 14.0,       22.0, -1.0, 1.0,
                23.0, 15.0,  7.0, 17.0, -1.0, 1.0,
                81.0, 57.0, 29.0, 55.0, -1.0, 1.0 };

/* Objective */
    double dobj[]={ 18.0, 21.0, 18.0, 16.0, 10.0, 15.0, 16.0, 14.0, 9.0,
              10.0,  9.0,  6.0, 17.0, 16.0, 17.0, 15.0, 10.0, 0.0,
              13.0,  0.0, 13.0,  0.0,  7.0,  0.0,  7.0,  0.0, 1.0 };

/* Column bounds */
    double dclo[]={ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
              0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
              0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0 };
    double dcup[]={ INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,
              INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,
              INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF };

/* Row bounds */
    double drlo[]={ -INF, -INF, -INF, -INF,  0.0, 0.0, 0.0, 0.0, 0.0 };
    double drup[]={ 10.0, 19.0, 25.0, 15.0,  0.0, 0.0, 0.0, 0.0, 0.0 };

/* Stages */
    long nstg=2;
    long n_first_stg_rows=4;
    long n_first_stg_cols=17;
    long rstg[]={ 1,1,1,1,2,2,2,2,2 };
    long cstg[]={ 1,1,1,1,1,1,1,1,1,
            1,1,1,1,1,1,1,1,2,
            2,2,2,2,2,2,2,2,2 };

/* Stochastic data */
    long nindp=5;
    long nsamp[]={ 5, 2, 5, 5, 3 };
    double demand[]={ 200, 220, 250, 270, 300,
                50, 150,
                140, 160, 180, 200, 220,
                10, 50, 80, 100, 340,
                580, 600, 620 };
    double dprobs[]={ 0.2, 0.05, 0.35, 0.2, 0.2,
                0.3, 0.7,
                0.1, 0.2, 0.4, 0.2, 0.1,
                0.2, 0.2, 0.3, 0.2, 0.1,
                0.1, 0.8, 0.1 };

/* local variables */
    long rc,ns=1,ii,iii,jj,*indx,*incr;
    double dd,dp=1.0;

/* Initialize OSL environments */
   rc = ekks_InitializeOSL();
     if (rc > 0) return;

/* compute number of scenarios */
    for (ii=0;ii<nindp;ii++) ns *= nsamp[ii];
/* initialize Stoch_Data with ns scenarios */
    ekks_InitializeApplication( &rc,&s_data,ns );
    if ( rc>1 ) return;

/* Create Core */
    ekks_CreateCore( &rc,s_data,nstg,rstg,cstg,nrow,ncol,nels,
          dobj,drlo,drup,dclo,dcup,mrow,mcol,dels,0 );

/* Generate Scenarios */

/* ...zero everything */
    memset( dobj,0,ncol*sizeof(double) );
    memset( drlo,0,nrow*sizeof(double) );
    memset( drup,0,nrow*sizeof(double) );
    memset( dclo,0,ncol*sizeof(double) );
    memset( dcup,0,ncol*sizeof(double) );
    nels = 0;

/* ...initialization -- indx points to first sample of each rv */
    indx = (long *) malloc( (1+nindp)*sizeof(long) );
    memset( indx,0,(1+nindp)*sizeof(long));
    for (jj=0;jj<nindp;jj++) {
     indx[jj+1] += indx[jj] + nsamp[jj];
     dp *= dprobs[ indx[jj] ];
     drlo[ n_first_stg_rows + jj ] = demand[ indx[jj] ];
     drup[ n_first_stg_rows + jj ] = demand[ indx[jj] ];
    }

/* ...add first sample */
    ekks_AddScenario( &rc,s_data,1,0,2,dp,1,
          nrow,ncol,nels,dobj,drlo,drup,dclo,dcup,
          mrow,mcol,dels,ADD_TO_CORE_VALUES );

/* ...initialize the sample space increment for each rv to 1 */
    incr = (long *) malloc( nindp*sizeof(long) );
    for (jj=0;jj<nindp;jj++) incr[jj] = 1;

/* ...main loop
     The logic of this loop is as follows:
       While the sample size nsamp[jj] divides the scenario index ii,
         reverse the increment direction incr[jj],
          and increase the random variable index jj by 1.
       Increment the jj'th random variable by incr[jj]
         and call ekks_AddScenario with the new sample.
*/

    for (ii=1;ii<ns;ii++) {
     iii=ii; jj=0;
     while ( !(iii%nsamp[jj]) ) {
         iii /= nsamp[jj];
         incr[jj] = -incr[jj];
         jj++;
     }
     dd = dprobs[ indx[jj] ];
     indx[jj] += incr[jj];
     dp = dp* dprobs[ indx[jj] ]/dd;
     drlo[ n_first_stg_rows + jj ] = demand[ indx[jj] ];
     drup[ n_first_stg_rows + jj ] = demand[ indx[jj] ];
     ekks_AddScenario( &rc,s_data,ii+1,0,2,dp,1,
               nrow,ncol,nels,dobj,drlo,drup,dclo,dcup,
               mrow,mcol,dels,ADD_TO_CORE_VALUES );
    }

/* Save SP-file as "sgbd.spl" */
    ekks_PutStochasticData( &rc,s_data,"sgbd.spl" );

/* Free Stoch_Data structures */
    ekks_FreeStructure(s_data);
}

exsolv.c 

Read an SP file and solve using decomposition. 
/*
 * exsolv.c:  Main driver for reading SP-files and solving.
 *
 *         Reads problem from file .spl and solves using
 *         exsolv().
 *         (compile with  /DOSLMSDLL on INTEL platforms.)
 */

#include 
#include 
#include 
#include "ekkstoch.h"

#define MXSCN 750              /* maximum scenarios */

#define TARG_OPT_GAP  1.0e-8   /* default target optimality gap */

#define abs(a) (a>0?a:-(a))

void main( long argc, char *argv[] ) {

    Stoch_Data *s_data=NULL;            /* Stoch_Data pointer */
    Stoch_Model *s_model=NULL;          /* Stoch_Model pointer */
    usrCoreData *u_core=NULL;           /* usrCoreData pointer */
    double *solution;                   /* space for solution */
    long *osl_indx;                     /* space for osl indices */
    long rc, nf, ii, nblks, nmod=1;     /* local variables */
    long Stage, Scenario;               /* Stage#, Scenario# */
    long colbeg, colindex;              /* Column start, column index */

    /* filenames */
    char *splfile,*outfile;
    char myfile[80];
    long nchar,nsize;

    /* get application name from command line */
    if ( argc != 2 ) {
         printf(" Usage:  solv   , where "
             "{file.spl} is a valid SPL file.\n");
         return;
    }
    nchar = (long) strlen( argv[1] );
    nsize = (nchar+1)+4;
    splfile = (char *) malloc(sizeof(char)*nsize);
    if ( splfile==NULL ) {
         printf(" ERROR: solv is unable to store input file-name\n");
         return;
    }
    outfile = (char *) malloc(sizeof(char)*nsize);
    if ( outfile==NULL ) {
         printf(" ERROR: solv is unable to store output file-name\n");
         return;
    }
    strcpy( splfile,argv[1] );
    strcat( splfile,".spl" );
    sprintf(myfile,"%s",splfile);

    strcpy( outfile,argv[1] );
    strcat( outfile,".out" );

/* initialize OSL environment */
    rc = ekks_InitializeOSL();  if (rc>0) return;

/* initialize Stoch_Data with maximum MXSCN scenarios */
    ekks_InitializeApplication( &rc,&s_data,MXSCN );
     if ( rc>1 ) return;

/* Read SP-file */
    sprintf(splfile,"%s",myfile);
    ekks_GetStochasticData( &rc,s_data,splfile );
     if ( rc>1 ) return;

/* Declare stochastic model (full model) */
    ekks_DescribeModel( &rc,s_data,0,0,&s_model,&nblks,1 );
     if ( rc>1 ) return;

/* Solve -- using ekks_BendersLSolve */
    {
        long subproblems=10,iterations=100;
        long ipcut=2, inwmt=0, icrsh=3, ilpdc=0, iscal=0, iprsl=0,
                iprsltype=0, iprslunit=0, isalg=1;
        long large_bounds=0;
        long status;
        double optmaxmin=1.0, optgap=TARG_OPT_GAP;  /*  optmaxmin= 1.0: minimize
                                                        optmaxmin=-1.0: maximize */

        /* Solve */

        long bndrrtcod;
        ekks_BendersLSolve( &bndrrtcod, s_model, ipcut,
            subproblems, iterations, optmaxmin,
            inwmt, icrsh, ilpdc,
            iscal,
            iprsl, iprsltype, iprslunit,
            isalg, large_bounds,
            &status,&optgap );
        if ((rc = bndrrtcod) > 1) return;

    }

/* Get core data */
    ekks_GetCoreData(&rc, s_data, &u_core );
    if ( rc>1 ) return;

/* Get stage/scenario solution */
    Scenario = 3;
    Stage = 2;
/* Malloc space for solution and index arrays */
    solution = (double *) malloc( u_core->ncol*sizeof(double) );
    osl_indx = (long *) malloc( u_core->ncol*sizeof(long) );
    colbeg = u_core->colstart[Stage-1]-1;
    nf = u_core->colstart[Stage] - u_core->colstart[Stage-1];

/* Get the solution */
    ekks_GetSolution( &rc,s_model,Scenario,0,solution,osl_indx);
    if ( rc>1 ) return;

/* Print scenario/stage solution */
    printf("\n\nScenario %d has %d columns in stage %d\n\n",Scenario, nf, Stage);
    for (ii=0;iisortin2exc[ii+colbeg]-1;
     if ( abs(solution[colindex]) < 1.0e-10 ) solution[colindex] = 0.0;
     printf(" Column[%d] = %g\n",osl_indx[colindex], solution[colindex]);
    }

/* Free data structures */
    ekks_FreeStructure( s_model ); ekks_FreeStructure( s_data );
        return;
}

exsdst.c 

Solve and compute distribution of optimal values. 
/*
 *
 * exsdst.c:  Present distribution of optimal values.
 *
 *         Usage:  exsdst <file>
 *
 *         Reads problem from file <file>.spl and solves using
 *         ekks_BenderLsolve().  Then, as an example of user coding, computes
 *         the distribution of the optimal value.
 *         (compile w/ /DOSLMSDLL on INTEL platforms.
 */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "ekkstoch.h"

#define MXSCN 750               /* maximum scenarios */

#define TARG_OPT_GAP  1.0e-8   /* default target optimality gap */

void main( long argc, char *argv[] ) {

    Stoch_Data *s_data=NULL;             /* Stoch_Data pointer */
    Stoch_Model *s_model=NULL;           /* Stoch_Model pointer */
    usrCoreData *u_core=NULL;            /* usrCoreData pointer */
    double *solution=NULL;               /* space for solution */
    long rc, rtcod=0, nf, ii, nblks, nmod=1;        /* local variables */

     /* filenames */
     char *splfile;
     char myfile[80];
     long nsize,nchar;
     long len;

     /* get application name from command line */
     if ( argc != 2 ) {
          printf(" Usage:  exsdst <file>  , where "
              "{file.spl} is a valid SPL file.\n");
          return;
     }
     nchar = (long) strlen( argv[1] );
     nsize = (nchar+1)+3;
     splfile = (char *) malloc(sizeof(char)*nsize);
     if ( splfile==NULL ) {
          printf(" ERROR: exsdst is unable to store input file-name\n");
          return;
     }
         strcpy( splfile,argv[1] );
     strcat( splfile,".spl" );
     sprintf(myfile,"%s",splfile);


/*  Initialize OSL environments */
    rc = ekks_InitializeOSL();  if (rc > 0) return;

/* initialize Stoch_Data with maximum MXSCN scenarios */
    ekks_InitializeApplication( &rc,&s_data,MXSCN );
     if ( rc>1 ) return;

/* Read SP-file */
    sprintf(splfile,"%s",myfile);
    ekks_GetStochasticData( &rc,s_data,splfile );
     if ( rc>1 ) return;

/* Declare stochastic model (full model) */
    ekks_DescribeModel( &rc,s_data,0,0,&s_model,&nblks,1 );
     if ( rc>1 ) return;

/* Solve -- using ekks_BenderLSolve */
    {
        long subproblems=10,iterations=1;
        long ipcut=2, inwmt=0, icrsh=3, ilpdc=0, iscal=0, iprsl=0,
                iprsltype=0, iprslunit=0, isalg=1;
        long large_bounds=0;
        long status;
        double optmaxmin=1.0, optgap=TARG_OPT_GAP;   
                /* optmaxmin=  1.0: minimize 
                            = -1.0: maximize */


        /* Solve */

        long bndrrtcod;
        ekks_BenderLSolve( &bndrrtcod, s_model, ipcut,
            subproblems, iterations, optmaxmin,
            inwmt, icrsh, ilpdc,
            iscal,
            iprsl, iprsltype, iprslunit,
            isalg, large_bounds,
            &status,&optgap );
        if ((rc = bndrrtcod) > 1) return;

    }

/*
 * Example of User computation using scenario tree and OSL
 *
 */

     {
          usrStoch *stree = NULL;
          usrNodeData *ndata = NULL;
          usrCoreData *cdata = NULL;
          usrNode *pnode = NULL;
          double *node_obj = NULL;

          /* locals */
          long in, il, ibas, jj;
          double nodvalue, scenvalue, mean=0.0, variance=0.0;
          long ii, is, it;
          double *dnobj, *dcobj, *dsoln, dprob;

          /* get core data */
          ekks_GetCoreData( &rc,s_model,&cdata );   if ( rc>1 ) return;

          /* get stochdata information */
          ekks_GetScenarioTree( &rc, s_model, &stree ); if ( rc>1 ) return;

          /* malloc node objective value array */
          node_obj = (double *) malloc ( stree->nnode*sizeof(double));
          memset(node_obj,0,stree->nnode*sizeof(double));

          /* get OSL pointers to solution array */
          dsoln = (double *) malloc ( stree->node->ncol  * sizeof(double));
          if (dsoln == NULL) {
                      printf("\n Not enough memory available"); 
                          return;
                  }

          ekks_GetPointerToSolution(&rc, s_model, &nmod, &dsoln);
          if (rc > 1) return;

/*
 * Loop over nodes and compute objective value
 * once for each node.
 */
          for (ii=0;ii<stree->nnode;ii++) {
               pnode = stree->node + ii;

               /* get core objective */
               it = pnode->stag;
               dcobj = cdata->dobj + cdata->colstart[it-1]-1;

               /* node offset in OSL model */
               ibas = pnode->cbas;

               for(jj=0;jj<pnode->ncol;jj++){
                    node_obj[ii] += dsoln[ibas+jj] * dcobj[jj];
               }
               if ( ii > 0 ) {
                  /* get node objective */
                  ekks_GetNodeData( &rc, s_model, ii, &ndata );
                  dnobj = ndata->dobj;

                  if(dnobj!=NULL)
                        for(jj=0;jj<pnode->ncol;jj++){
                                node_obj[ii] += dsoln[ibas+jj] * dnobj[jj];
                  }
               }
          }


/*
 * For each leaf node,
 *      go up the tree adding in node_obj value along the way.
 */

          printf(" \n \t Distribution of Optimal Values \n \n");

          for(ii=0;ii<stree->nleaf;ii++) {

               scenvalue = 0.0;

               in = stree->leaf[ii];          /* name of leaf node */

               pnode = stree->node + in-1;     /* pointer to leaf info */
               is = pnode->scen;          /* leaf scenario */
               dprob = pnode->prob;          /* leaf probability */

               while (in>0) {

                    scenvalue += node_obj[in-1];

                    in = pnode->par;     /* parent pos'n in list */
                    if (in>0) pnode = stree->node + in - 1;
               }

               printf(" \t Scenario %d value %f probability %f \n",
                   is,scenvalue,dprob);

               mean += dprob*scenvalue;
               variance += dprob*scenvalue*scenvalue;

          }

          variance -= mean*mean;
          if (variance < 0.0) variance=0.0;

          printf(" \t \t Mean Value: \t \t %f \n",mean);
          printf(" \t \t Standard Deviation: \t %f \n",sqrt(variance));

    }

/* Free data structures */
    ekks_FreeStructure( s_model ); ekks_FreeStructure( s_data );

}

exsgbdt.c 

Assemble tree using ekks_*Tree() subroutines. 
/*  
 *
 * exsgbdt.c:  Generic driver using ekks_*Tree() subroutines.
 *
 *         Input distribution is sampled and processed by
 *         ekks_*Tree() subroutines.
 *
 *         Implements G.B. Dantzig's Aircraft Allocation model,
 *         in: G.B. Dantzig, Linear Programming and Extensions
 *         (Princeton U. Press, 1963) pp. 572-597.
 *
 *         See also: A.J. King, "Stochastic programming problems:
 *         examples from the literature",
 *          in: Yu. Ermoliev and R.J-B Wets, eds, Numerical Techniques for
 *         Stochastic Optimization (Springer-Verlag, 1988) pp. 543-567.
 *         (compile w/ /DOSLMSDLL on INTEL platforms.
 */

#include <stdlib.h>
#include "ekkstoch.h"

#define INF 1.0e31

#define N_SAMPLES  300      /* number of samples to be generated */

void main( long argc, char *argv[] ) {

    Stoch_Data *s_data;                     /* Stoch_Data pointer */

/* Model dimensions */
    long ncol=27, nrow=9, nels=44;

/* Sparse matrix data...organized by row */
    long mrow[]={ 1, 1, 1, 1, 1,
            2, 2, 2, 2,
            3, 3, 3,
            4, 4, 4, 4, 4,
            5, 5, 5, 5,
            6, 6, 6, 6, 6, 6,
            7, 7, 7, 7, 7,
            8, 8, 8, 8, 8, 8,
            9, 9, 9, 9, 9, 9 };
    long mcol[]={ 1, 2, 3, 4, 5,
            6, 7, 8, 9,
            10, 11, 12,
            13, 14, 15, 16, 17,
            1,        13, 18, 19,
            2, 6, 10, 14, 20, 21,
            3, 7,     15, 22, 23,
            4, 8, 11, 16, 24, 25,
            5, 9, 12, 17, 26, 27 };
    double dels[] = { 1.0, 1.0, 1.0, 1.0, 1.0,
                1.0, 1.0, 1.0, 1.0,
                1.0, 1.0, 1.0,
                1.0, 1.0, 1.0, 1.0, 1.0,
                16.0,              9.0, -1.0, 1.0,
                15.0, 10.0,  5.0, 11.0, -1.0, 1.0,
                28.0, 14.0,       22.0, -1.0, 1.0,
                23.0, 15.0,  7.0, 17.0, -1.0, 1.0,
                81.0, 57.0, 29.0, 55.0, -1.0, 1.0 };

/* Objective */
    double dobj[]={ 18.0, 21.0, 18.0, 16.0, 10.0, 15.0, 16.0, 14.0, 9.0,
              10.0,  9.0,  6.0, 17.0, 16.0, 17.0, 15.0, 10.0, 0.0,
              13.0,  0.0, 13.0,  0.0,  7.0,  0.0,  7.0,  0.0, 1.0 };
        double dqdg[]={ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
              0.0,  0.0,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
              13.0,  0.0, 13.0,  0.0,  7.0,  0.0,  7.0,  0.0, 1.0 };


/* Column bounds */
    double dclo[]={ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
              0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
              0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0 };
    double dcup[]={ INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,
              INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,
              INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF };

/* Row bounds */
    double drlo[]={ -INF, -INF, -INF, -INF,  0.0, 0.0, 0.0, 0.0, 0.0 };
    double drup[]={ 10.0, 19.0, 25.0, 15.0,  0.0, 0.0, 0.0, 0.0, 0.0 };

/* Stages */
    long nstg=2;
    long n_first_stg_rows=4;
    long n_first_stg_cols=17;
    long rstg[]={ 1,1,1,1,2,2,2,2,2 };
    long cstg[]={ 1,1,1,1,1,1,1,1,1,
            1,1,1,1,1,1,1,1,2,
            2,2,2,2,2,2,2,2,2 };

/* Stochastic data */
    long nindp=5;
    long nsamp[]={ 5, 2, 5, 5, 3 };
    double demand[]={ 200, 220, 250, 270, 300,
                50, 150,
                140, 160, 180, 200, 220,
                10, 50, 80, 100, 340,
                580, 600, 620 };
    double dprobs[]={ 0.2, 0.05, 0.35, 0.2, 0.2,
                0.3, 0.7,
                0.1, 0.2, 0.4, 0.2, 0.1,
                0.2, 0.2, 0.3, 0.2, 0.1,
                0.1, 0.8, 0.1 };

/* Tree structure */
    Tree *tree=NULL;

/* local variables */
    long rc,ns=1,ii,iii,jj,*indx,*incr;
    double dd,dp=1.0,dwgt;
    long irand,indx_beg,ev_indx;
    long nscen,ianc,istg;
    long *label;
    long evoffset;



/*  Initialize OSL environments */
    rc = ekks_InitializeOSL();
    if(rc > 0) return;

/* initialize Stoch_Data with ns scenarios */
    ekks_InitializeApplication( &rc,&s_data,N_SAMPLES );
        if ( rc>1 ) return;

    ekks_ProblemType( &rc,s_data,quadratic);
        if ( rc>1 ) return;

/* Create Core */
    ekks_CreateCore( &rc,s_data,nstg,rstg,cstg,nrow,ncol,nels,
          dobj,drlo,drup,dclo,dcup,mrow,mcol,dels,dqdg );


/* ...zero everything */
    memset( dobj,0,ncol*sizeof(double) );
    memset( drlo,0,nrow*sizeof(double) );
    memset( drup,0,nrow*sizeof(double) );
    memset( dclo,0,ncol*sizeof(double) );
    memset( dcup,0,ncol*sizeof(double) );
    nels = 0;

/* ...malloc event indexing array */
    indx = (long *) malloc( (1+nindp)*sizeof(long) );
    memset( indx,0,(1+nindp)*sizeof(long) );

/* ...malloc label array */
    label = (long *) malloc( (nstg-1)*sizeof(long) );

/* Generate Scenarios from sampling */

    for (ns=0;ns < N_SAMPLES;ns++) {

/*   ...for each indep rv, process random number to get index */
        indx_beg = 0;
        for (jj=0;jj<nindp;jj++) {
                ev_indx = 0;
                irand = rand();
                dd = (double)irand/(double)RAND_MAX;
                dp = dprobs[indx_beg+ev_indx];
                while ( irand > dp*RAND_MAX )
                        dp += dprobs[++ev_indx + indx_beg];
                indx[jj] = ev_indx;
                indx_beg += nsamp[jj];
        }

/*   ...process indep events into a unique "event label"
        (there should be one such label for each stage in the simulation
        -- this example has only one stage ) */

        label[0]= indx[0];
        for (jj=1;jj<nindp;jj++)
                label[0] = label[0]*nsamp[jj] + indx[jj];

        dp = 1.0e0;
        dp /= N_SAMPLES;

/*   ...initialize tree if necessary, or add event to tree */
        if (!tree)
                ekks_InitializeTree( &rc,&tree,label,nstg-1,N_SAMPLES,dp );
        else
                ekks_AddToTree( &rc,tree,label,nstg-1,dp);
        if (rc>0) return;

/*   ...get scenario branching information from tree */
        ekks_GetScenarioFromTree( &rc,tree,ns,&ianc,&istg,&dwgt,
                                                                &nscen,&label );
        if (rc>0) return;

        if ( istg == 1 )        /* always true in this example */
                ianc = 0;

        istg++;                 /* always need to do this */

/*   ...assign values to right-hand sides based on index */
        evoffset=0;
        for (jj=0;jj<nindp;jj++) {
                drlo[ n_first_stg_rows + jj ] = demand[ indx[jj]+evoffset ];
                drup[ n_first_stg_rows + jj ] = demand[ indx[jj]+evoffset ];
                evoffset += nsamp[jj];
        }

/* ...add sample */
        ekks_AddScenario( &rc,s_data,nscen,ianc,istg,dwgt,1,
                nrow,ncol,nels,dobj,drlo,drup,dclo,dcup,
                mrow,mcol,dels,ADD_TO_CORE_VALUES );
        if (rc>1) return;
    }


/* Save SP-file as "sgbd_smpl.spl" */
    ekks_PutStochasticData( &rc,s_data,"sgbd_smpl.spl" );

/* Free Stoch_Data structures */
    ekks_FreeStructure(s_data);
}

[ Top of Page | Previous (Reference Section) | Next (Bibliography) | Roadmap ]