SELECTED MATHEMATICA TOOLS FOR MATHEMATICAL PROGRAMMING Derek Holmes Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, MI 48109-2117 Technical Report 93-26 September 1993

Selected Mathematica Tools for Mathematical Programming Derek Holmes Department of Industrial and Operations Engineering University of Michigan Ann Arbor, MI 48109-2117 September 1, 1993. Abstract: Mathematica is a widely used computer software system for computational mathematics which is easily programmable in a higher level language than either Fortran or C. The Mathematica language is well suited for algorithmic prototyping and model manipulation. This report describes several Mathematica functions which may be of use in linear and stochastic programming. These functions, organized into "packages," are divided into three categories: packages used for determinisitc linear programming, packages used for stochastic programming, and a package which may be used to manipulate LPs and SPs. The last package may be used to transform a general two-stage stochastic program into a multi-stage stochastic program of any size.

Contents 1 Introduction 2 2 General Determinisitic Programming 2 2.1 Problem manipulation in MPS format (mps.m)........................ 2 2.1.1 Reading an MPS file.................................. 3 2.1.2 Writing an MPS file..................................... 3 2.1.3 W riting a Quadratic MPS file................................ 3 2.1.4 Obtaining submatrices and subproblems......................... 4 2.1.5 Reading nonzero fragments................................. 4 2.2 Problem Plots (lp.m, lpeduc.m)................................. 4 2.2.1 Problem Visualization.................................... 4 2.3 Solving linear programs and their variants (Ip.m, lpeduc.m)...................... 5 2.3.1 Problem Manipulation Routines............................. 6 2.3.2 Linear Programming using the Primal Simplex Method................. 6 2.3.3 Parametric and Fractional Linear Programming..................... 7 2.3.4 Small scale Interior Point Algorithms.......................... 8 3 Stochastic Programming 9 3.1 Approximations.................................... 9 3.1.1 Upper Bounds....................................... 11 3.1.2 Lower Bounds........................................ 12 3.1.3 Refinement Procedures.................................... 12 3.1.4 Miscellaneous functions................................... 12 3.2 L-Shaped M ethod.......................................13 3.3 Miscellaneous Functions.................................. 14 4 Stochastic Program generation and manipulation 14 4.1 Reading in Data................................. 14 4.2 Describing stochasticity................................... 15 4.3 Outputting stochastic programs............................... 16 4.4 An exam ple............................................ 17 1

1 Introduction Mathematica [1], a computer package developed by Wolfram Research, is described as "a system for doing mathematics by computer." The package includes several very high level functions which are of interest to operations researchers, including linear and nonlinear programming, linear systems solutions, and matrix manipulation routines. Of more interest is Mathematica's ability to be programmed in a language flexible enough to allow succint and quick prototyping of any algorithm or concept. For example, complete set operations, string manipulation operations, and procedural programming constructs can be combined with numerical integration routines to read in an arbitrary linear program and solve its continuously distributed stochastic counterpart. This technical report briefly describes the design and use of several procedures which have been developed by the author for research in stochastic programming. These procedures are organized into sets of functions called "packages", which can be loaded into Mathemaiica and used as any other (built-in) function. (We assume that the reader has a working knowledge of Mathematica and is familiar with Mathematica terminology.) The packages are organized into three categories, corresponding to their scope. Packages in the first category are used for general deterministic programming, and include functions to read and write LP and QP MPS files (see, e.g. [2] or [3]), solve LPs using the simplex method and interior point methods, and solve fractional LPs and parametric LPs. Matrix visualization tools that utilize Mathematica's graphics capability are also included. Packages used for stochastic programming include routines for the Van Slyke and Wets L-shaped algorithm (Van-Slyke and Wets, [4]) and the Edmondson-Madansky and Jensen bounds (see, e.g [?]. Functions in the last category are used to generate arbitrary multistage stochastic linear programs (for an introduction, see, e.g. Ermoliev and Wets, Chapter 1 [5]) from determinisitic LP formulations. Note: The author makes no guarantees regarding the efficacy or stability of these functions. They are experimental only and are intended for individual modification to a given task. 2 General Determinisitic Programming This section includes Mathematica routines which are of interest to the general operations research practitioner. The routines contained in this section are designed for small problems, but can tackle mesium sized to large tasks given adequate hardware. (In the future, Matihematica can be linked with large scale optimizers to allow larger problems to be tackled.) 2.1 Problem manipulation in MPS format (mps.m) mps.m is a package which reads and writes linear programs in MPS format. Data in Mathematica are usually described as nested lists. Since the elements of the list are not restricted to be of a single data type, they can be used as complex data structures. The routines described in this package use a standard nested list (denoted here by (LP)) to describe the linear program minimize CT subject to Ai.Obi, i = 1,..., m z>0 where E {=,<,>} is{ c, A, b, rhstype, colnames, rovrames, problenmane }. Here, c is a (dense) list of objective elements, b is a list of the RHS elements and A can either be a dense or sparse representation of the constraint matrix. The lists corresponding to these representations are described below. Dense {{a1,..., aln},..., {aml...,amn}} Sparse {{il,ji,,j,},..., {ipjpai,j,}} where ait,ji..., aij, are the p nonzeros of the constraint matrix. We define rhstype to be a list of m elements, where rhstype[[i]] is in {0, 1,-1} when constraint i is an (equality, <, or ~) constraint. (The value correspond to the multiplier of a slack variable added to ensure

that each constraint is an equality.) rornames and colnames are lists of m + 1 and n strings describing row and column names. The first element of rovnames is always the name of the objective row. Finally, problemname is a string identifying the problem. Given this background, the functions performed by the routines in mps.m can be described. 2.1.1 Reading an MPS file An MPS file is a test file containing sections for describing rows, columns, RHSs, ranges, and bounds. A problem in theis format may be read in using (in mps.m) ReadMPS [stream, Options] Inputs stream An input stream for the MPS file, typically OpenRead[ile) Options _ Sparse->True Output matrix in Sparse format. Use False for dense format. Removelulls->False Remove rows and columns that contain no nonzeros. where the output of the function is a list of the form (LP), defined above. (Note:) Ranges and Bounds are incorporated into the constraint matrix within ReadMPS. New row and column names are formed from the source ranges and bounds. 2.1.2 Writing an MPS file An MPS file can be written using (in mps.m) WriteMPS[stream, c, A, b, Options] Inputs stream An output stream for the MPS file, typically OpenVrite[file] c A dense objective vector. A A sparse or dense constraint matrix. b A dense requirements vector. Options RHSformlist->{list} constraint form (see above). If missing, all constraints are assumed =. Rowlames->{list} (m + 1)-list of strings for the problem rows. Coluamnames->{list} (n)-list of strings for the problem columns. ProblemName->"" The problem name. Sparse->True Input matrix is in Sparse format. Sort->True Sort columns if constraint matrix is in Sparse form. The the output of the function is True if no problems occurred, False otherwise. If RowNames or Columnmmaes are not specified, then VriteMPS creates its own. 2.1.3 Writing a Quadratic MPS file The following function outputs only theiquadratic portion of an MPS file (for format guides, see, e.g. [3]): (in mps. a) WriteQuadraticMPS[stream, Q, Options] Inputs stream An output stream for the MPS file, typically OpenWrite file] Q The (dense) quadratic matrix, i.e. Q for min TQx + cTz. Options Columnlazues->{list} (n)-list of strings for the problem columns. ProblemName->"" The problem name. MulitplybyTwo->False Mulitplies all entries in Q by 2. where the output of the function is True if no problems occurred, False otherwise. 3

2.1.4 Obtaining submatrices and subproblems Parts of problems may be separated out by using (in mps.m) SelectSubLP[LP, BeginningCoords, EndCoords] Inputs LP A (Sparse) linear program in the format described for (LP) above. BeginningCoords (row,col) coordinates for the upper left hand corner of the subproblem. EndCoords (row,col) coordinates for the lower right hand corner of the subproblem. where the output of the function is the subproblem in the (LP) format and (in mps.m) Submatrix[A, BeginningCoords, EndCoords, Options] Inputs A A matrix in sparse (default) or dense format. BeginningCoords (row,col) coordinates for the upper left hand corner of the submatrix. EndCoords (row, col) coordinates for the lower right hand corner of the submatrix. Options Sparse->True Output matrix in Sparse format. Use False for dense format. where the output is the submatrix in the same form (Dense or Sparse) as the input matrix. 2.1.5 Reading nonzero fragments A common situation in stochastic programming occurs when a fragment of nonzeros in MPS format must be read and matched up with the column and row names found from a previously read LP. For example, stoch files frequently contain BLOCKS or SCENARIOS sections with submatrix realizations. The following function reads in these fragments. (in mps.m) onzeroFragment [stream, LP] Inputs stream An input stream with the nonzeros in a form similar to MPS's COLUMNS section. LP A linear program in (LP) format, described above. Output of the function is the submatrix corresponding to the nonzero fragment in sparse form. 2.2 Problem Plots (Ip.m, lpeduc.m) One of Mathematica's greatest strengths is its graphics capabilities. These capabilities can be used to visualize the nonzero structures of large problems, or to visualize a particular algorithm for educational purposes. We split these functions into problem visualization tools and algorithm visualization tools. 2.2.1 Problem Visualization The following functions use these capabilities to display the nonzero structure of a given problem. (in ups a) LPPlot LP, Options] Inputs LP A linear program in (LP) format, described in Section 1.1. Options Sparse->True Input matrix is in Sparse format. Use False for dense format. HighlistList->{} Output 3-D graph in which the (sparse) nonzeros in HighlistList are taller than the underlying LP. Output->False Function returns Graphics object if True, or shows graph and returns True. 4

*... * —----------- * * - Figure 1: CEPI matrix Figure 2: CEPI with first stage highlighted 5

Here, the output depends on the Output option shown above. Two examples of the output of this function applied to a problem from the literature are shown below. The first is the problem without highlights. The second example shows the problem with the first period constraints highlighted. (Higlights don't work very well for large problems) Prototype: (in mps.m) MatPlot[c,A,b,Options] where the inputs are the elements of an LP and the options and output are identical to LPPlot. For educational purposes, it is frequently desirable to plot arbitrary two-dimensional polytopes, and show the progress of algorithms within these regions. The following function determines the Graphics objects to visulaize the polytope described by {z(Az < b}. (in lpeduc.m) PolytopeGraph CA,b] Inputs A,b The constraint matrix and RHS vector in Ax < b. and the output is a list of Graphics objects. To see the graph, type Show CPolytopeGraphCA,b, Axes->True. 2.3 Solving linear programs and their variants (Ip.m, lpeduc.m) Although Mathematica has several optimization functions, including one for Linear Programming, these are frequently inadequate for research purposes. The biggest problem with the LinearProgramming function is that it returns only a primal solution; no dual solution is available when (e.g.) a problem is unbounded. The first package, Ip.m is the author's own implementation of common simplex-based algorithms for linear programming and its variants. The algorithms are taken from (Bazaraa, Jarvis, and Sherali, [7]) and are given after some problem manipulation routines. 2.3.1 Problem Manipulation Routines A table of some useful problem manipulation routes is given below. Forming Problems AppendColumns [M1, M2,...] Returns a submatrix with matrices Mi stacked on top of each other. AppendRows [M1, M2,...] Returns the matrix (M1:, M,:,...). BlockMatrix CM] Returns a valid matrix of elements when M is a matrix formed from submatrices. ZeroMatrixEm) Returns a square zero matrix of m elements. Changing Problems MakeEqualityConstraintsE Returns an equality constrained c,A,b,rhstype, Sparse->True] problem formed from an LP in (LP) form. LTtoEQCc,A,bJ Returns an equality constrained problem formed from min{cTzlAx < b,x > 0}. Vector Tests VecGTCvec,no] returns True if vec[[Ci3 > no for all i. VecGECvec,no3 returns True if vecCCi]] > no for all i. VecLT vec,no] returns True if vec[[Ci] < no for all i. VecLE[vec,no3 returns True if vecCCi3] < no for all i. VecEQ vec,no] returns True if veCCi33 = no for all i. A useful tool for "debugging" linear programs and computer codes that manipulate them is to find the set of nonzeros that are not common to two linear programs. A function for finding these nonzeros is in mps.m) LPDiffC[LP1, LP2J Inputs LP1,LP2 Two (sparse) linear programs in (LP) format to be compared. whose output is an LP in (LP) format with nonzero elements that are in either LP1 or LP2, but not both.

2.3.2 Linear Programming using the Primal Simplex Method Several functions in this package fit together to solve an LP problem. Most users will call LP, which uses the primal simplex method. Feasibility is obtained using the Big-M method (M is user settable). The Big-M method implemented here adds an artificial variable to each constraint and attempts to force it out of the basis first. Also, LP prevents cycling with an implementation of Bland's Rule (discussed in, e.g. [7]). (A warning message will be printed when 40 zero steplength pivots are performed, and Bland's Rule will be invoked after 80 zero steplength pivots.) Inputs (in lp.m_) LPc,A,b, Options] c, A, b A (dense) linear program of the form min{crzxAz = b,z > 0}. Options IdentityColumns->{} Use an identity in A given by a list of column numbers as artificials. Plevel->O Print more detail about algorithmic progress. Maximum: 3. Math->False If True, use Mathematica's LinearProgramming function. BigM->100000 Penalty for infeasibilities. and the output is a nested list of the form { z, xstar, pistar }, where z is the optimal objective function value, and xstar and pistar are the optimal primal and dual solutions, respectively. LP sets up the problem by adding artificial columns and declaring these to be an intial basis. The routine then calls LPvithBasis, with prototype (in lp.m) LPithBasisCBasislist, c,A,b, Options] Inputs_ Basislist is a list of columns that define an invertible submatrix. c, A, b A (dense) linear program of the form min{cTzxAx = b, z > 0}. Options_ Artificials->{} List of artificial column numbers. Plevel->O Printing detail. Bland->False If True, start with Bland's Rule. whose output is a nested list of the form { z, xbasic, pistar, Basislist } where z and pistar are defined as above, xbasic is the optimal values of the basic variables contained in Basislist. 2.3.3 Parametric and Fractional Linear Programming A parametric linear program has the general form min{(c+ Ac')TzxAz = b,x > 0}, where A > 0. LPPpara solves problems of this sort by performing basis change analyses within the context of the simplex method. (in lp.nm) LPpara plevel,iform,c,c',A,b,Options] Inputs plevel is the level of output printed by the algorithm (0=none). iforn is 1 if the constraints are of the form Az > b, 0 if Ax = b. c, c', A, b describe the parametric program. Options Lilits->{lo,hi} Perform parametric anlaysis for A E { lo,hi }. The optimal value of a parametric linear program is a convex piecewise linear function. So, the solution from the algorithm is a list of p piecewise linear terms of the form z+Az' for some range of A E {A,, Ai}, i = 1,..., p. Specifically, the output is a nested list of the form { z, xsollist, pi, paralist } where xsollist = {x;z,., *} is a list of optimal primal solutions for the p piecewise linear terms and paralist is a list of p lists of the form {A, A, z, z'}. Here, the first two terms are the ranges of A for which a given basis is optimal and second two terms define the optimal objective value z + Az' over that range. A fractional linear program has the general form pTx +a minmq + qjx + /

subject to x E x{zlA < b, z > 0}. The method of Charnes and Cooper [8]) rewrites this nonlinear program as a linear program of the form min pTy + az subject to Ay < bz, nonnegativity constraints, and qTy + 0/ = 1. The optimal solution to the original program is z' = y*/z*. The following function solves this problem. (in lp.m) LPfracEplevel,iopt, numvec, numcons, denvec, dencons, A,b] Inputs plevel is the level of output printed by the algorithm (O=none). iopt is 1 if the denominator is always > 0, -1 otherwise. numvec,numcons describe the objective numerator. denvec,dencons describe the objective denominator. A,b describe the constraint set. The output is a list of the form {zstar, xstar} where zstar is the (scalar) optimal objective function value and xstar is an optimal solution. 2.3.4 Small scale Interior Point Algorithms Ipeduc.m contains several functions for experimenting and visualizing interior point algorithms. Two types of algorithms are implemented: the affine scaling algorithm and the projective scaling algorithm. Details of the development and notation used for the algorithms are due to (Birge [12]) and are included as an appendix to this technical report. The relevent routines are AffScaling and ProjScaling. AftScaling runs the affine scaling algorithm until the stopping criterion Ilskll2 < e is met, where sj is the kth search direction. (in lpeduc.m) AffScaling CAffScaling [A, b, c,x,gamma, eps, Options] Inputs c, A, b A (dense) linear program of the form min{c' zlAz = b,z > 0}. x is an intial interior starting point z E {zxAz < b}. gamma is the fraction of the largest stepsize to take. eps is the stopping criteria mentioned above. Options Outputs->{} List of output options (see below). Print->{} List of printing options (see below). EllipseRanges->{lo,hi} Range of iterations for which ellipses are printed (if enabled). MatrixUpdate->0. 8 Update projection matrix only if relative change exeeeds this number. EllipseStepSize->True Use ellipse step sizes if True, 7 of max if False A list of the output options which can be included in the Outputs and Print options is given below: Option Meaning Solution Final Solution and objective Ellipses Ellipses at each iteration EllipseGraphics Ellipse graphics at each iteration Iterates Sequence of Iterates (points) Bounds Lower/Upper Bounds sNorms Norms of Search Directions StepSizes Actual StepSizes SearchDirections Search Directions Duals Dual iterates All All of the above The projective scaling implementation is similar. For details, refer to (Birge, [12]). 8

(in lpeduc. m) ProjScaling [ProjScaling CA,b,c,x,gamma, eps,UB, Options] Inputs c, A, b A (dense) linear program of the form min{c" xAx = b, x > 0}. x is an intial interior starting point x E {xzAx < b}. gamma is the fraction of the largest stepsize to take. eps is the stopping criteria mentioned above. UB is an initial upper bound on zTx*. Options Outputs->{ List of output options (see below). Print->{} List of printing options (see below). EllipseRanges->{lo,hi} Range of iterations for which ellipses are printed (if enabled). MatrixUpdate->0.8 Update projection matrix only if relative change exeeeds this number. PolynomialStepSize->alpha Use polnomial stepsize with a given a. A list of the output options which can be included in the Outputs and Print options is given below: Option Meaning Solution Final Solution and objective Ellipses Ellipses at each iteration xInnerEllipseGraphics Inner Ellipse graphics (in x-space) xOuterEllipseGraphics Outer Ellipse graphics (in x-space) zInnerEllipseGraphics Inner Ellipse graphics (in z-space) zOuterEllipseGraphics Outer Ellipse graphics (in z-space) Iterates Sequence of Iterates (points) Bounds Lower/Upper Bounds slorms Norms of Search Directions StepSizes Actual StepSizes SearchDirections Search Directions UBs Sequence of Upper Bounds Duals Dual iterates All All of the above A particularly common linear program solved in stochastic programming is an LP of the form Q(z) = min(qTyIWy = h - T, y > 0}. This program may be solved by calling (in lp.m) Q[x,h,Y,T,q] where the inputs are the (dense) vectors and matrices described above. 3 Stochastic Programming This section describes a few routines which may be of interest to stochastic programmers. Mathematica is especially useful for approximation algorithms, since it has powerful list manipulation capabilities. These capabilities are most useful for extremal approximations such as the Edmundson-Madansky upper bound and the Jensen lower bound. (See, e.g. [11]) Currently, two main packages have been developed. The first is for approximation, and the second implements the Van Slyke and Wets' L-Shaped algorithm for solving two stage linear programs. 3.1 Approximations One of the more difficult tasks in stochastic porgramming is calculating the expected value of the second stage solution. This calculation requires a multidimensional integration of an implicitly defined function. The most common method for approaching this problem is to use approximations to the second stage functional. An extensive set of approximations have been developed (see, e.g. Birge and Wets [131). The most common procedure is to find a set of discrete realizations which closely approximates the underlying distribution. Doing so allows the integral to be written as a finite sum of solutions to second stage optimization problems. 9

The two most common discretization procedures are based on Jensen's inequality and the EdmundsonMadansky inequality. (again for detail, see [13]). The problem we wish to approximate is a two-stage stochastic linear program with recourse, which can be written as min cxT + Q(z) s. t. Ax=b,z >O where Q(z) = E[Q(x, )] and Q(xz,) = min{qTylWy = - Tz, y > O}. The Jensen inquality finds a lower bound on Q(z) using the inquality Q(x) = E[Q(x,f)] > p'Q(x, E[I'l]), 1=1 where (S1,..., S') is a partition of ~'s support (_) and p1 = P[f E Sl]. The Edmundson-Madansky upper bound implemented here assumes that the support S is of the form [a,, b] x,..., x, [aN, bN], and that the components of ( are independent. The bound is W ext =1 - Refinements are taken by partitioning the support i into several subrectangles (S,... S^) and calculating the conditional bounds over each partition. The partitioning method implemented follows that of Frauendorfer and Kall [11], which (1) splits the cell with the greatest error between upper and lower bounds, (2) considers all 2m+l points to find the point and direction which maximizes the minimum linearization error. (The functions implemented here also assume the underlying continuous distribution is multiariate uniform.) Although both bounds require a description of the support, keeping more data can greatly speed up the calculation of the bound. Specifically, keeping the LP solution at each extremal point of the support and the mean solution eliminates the need to recalculate them each time a new bound is calculated. Each of these sets of data must be kept for each partition. The following lists are the "data structures" which hold this data: Name Contents Format PARTN Extreme points defining partition { pt 1,.., pt 2-n } SUPPTS Set of partitions { PARTN(1),.., PARTN(K) } MEANS Set of means of a parition { {mean 1},.., {ean k} } PROBS Marginal probabilities of a parition { pl,..,pk } SUPPORT Complete description of a support { SUPPTS, MEANS, PROBS } SOLN LP Solution to a point (z,', r*) { z, xstar, pistar } PTSOLS List of solutions of extrema of a partition. { SOLN(1),..,SOLN(2") } SOLVEDSUPPORT SUPPORT with extrema solutions. { SUPPORT, { PTSOLS(1),,.., PTSOLS(K) } } TIMES Cumulative timing data { Tine, Time, Time } (LB times, UB times, Partitioning times) TIMEDSUP SOLVEDSUPPORT with timing data { SOLVEDSUPPORT, TIMES } These lists serve as input to the main partitioning functions CP and CompletePartition: 10

(in approx.m) CP[iteration, timedsup, x, hconst, W,T,q, totalarea, streams] (in approx.m) CompletePartition[iteration, solvedsupport, x, hconst, W,T,q, totalarea, cumt, streams] Inputs___ iteration The iteration or refinement number. timedsup A list in TIMEDSUP form (shown above). x The current primal iterate. hconst (hl,...,hk) in the second stage RHS (hl,...,hk, hA+l,..., h). W,T,q Describe the second stage linear program. totalarea A divisor for normalizing the total initial probability to 1. streams A set of 3 streams for outputting (all dual solutions, not used, not used). Note: Production versions will eliminate streams. solvedsupport A list in SOLVEDSUPPORT form (shown above). cumt A list in TIMES form (shown above). Both calculate an upper and lower bound based on the refined support given and refine the support for another iteration. The output is a list of the form TIMEDSUP. These functions call several other functions to find the bounds and the partition. These functions are described in the following sections. 3.1.1 Upper Bounds The main function that calculates the upper bound is UBEM, which has the prototype (in approx.m) UBEM[itertation, x, support, hconst, W,T,q,probsup, opts] Inputs x, hconst, W,T,q are defined as above. iteration The iteration or refinement number. support A list in SUPPTS form (shown above). probsup A list in PROBS form (shown above). Options Streams->{} A list of three streams for outputting dual information (see below). LBs->{} A list of Lower Bounds for each cell to determine the cell to split. Optpi->{} A list of the optimal dual mulitpliers at x for comparison purposes. PrintLevel->l is the level of debugging detail. The Optpi and Streams options are used to determine the distance of the approximated dual solutions from the actual dual solutions. Solutions to each new point are output to Streams [[]]. Expected duals within each parition and their 2-norm distances to the optimal (marginal) r* are printed to Streams C[2]3, while the overall expected dual estimate and its distance to r* are printed on Streams I[3]]. UBEM returns a compund list { UB, newpis, UBlist, errorlist, celltosplit}, where the overall upper bound and new duals are the first two items. The third item is a list of upper bounds for each partition, and the last two items are the differences UBi - LBi in each cell (which is nonempty only if LBs is specified), and the cell which maximizes this error. To save computation, UBEM develops the upper bound from a list of solutions which are generated using the function: (in approx.m) SolvePointList[solnlist, support, x, hconst, W,T,q] Inputs ___ x, hconst, V,T,q are defined as above. so1nlist A list of solutions of the form PTSOLS. support A list in SUPPTS form (shown above). The output is a new solnlist in PTSOLS form. 11

3.1.2 Lower Bounds The main function that calculates the Jensen lower bound is LB, which has the prototype (in approx.m) LB[itertn, hbar, hconst, probsup, x, W,T,q, opts] Inputs iteration, x, hconst, W,T,q are defined as above. hbar A list of expected RHSs in each partition. probsup A list in PROBS form (shown above). Options Streams->{} A list of three streams for outputting dual information (see below). Optpi->{} A list of the optimal dual mulitpliers at x for comparison purposes. PrintLevel->1 is the level of debugging detail. Again, the Optpi and Streams options print the lower bounds and differences between the approximated dual solution and the optimal dual solution. LB returns a compound list of the form {LB, LBlist}, which contains the overall lower bound and the partition-specific lower bounds. 3.1.3 Refinement Procedures As mentioned earlier, the single partitioning strategy implemented is due to Fraundorfer and Kall, who assume a rectangular support. Several functions determine the refinement to be taken, all of which are called by (in approx.m) FindFKRefinement[Sup, sollist, partition,pivec,hbar, __________probsup, totalarea, opts] Inputs_____...... Sup A list in SUPPTS form (shown above). sollist A list of solutions of the form PTSOLS. partition The partition number to split. pivec A list of duals for the cell to be partitioned. (Use TrimPiCubout C2]],hconst] where ubout is from UBEM). hbar A list of expected RHSs in each partition. probsup A list in PROBS form (shown above). totalarea A divisor for normalizing the total initial probability to 1. Options PrintLevel->1 is the level of debugging detail. The returned list is of the form SOLVEDSUPPORT, with the new partitions replacing the old partition in the list. FindFKRefiineaent calls the function MaptoExtrea[Cpointlist, extremalist], which maps a list of references to extreme points to their actual values. 3.1.4 Miscellaneous functions A few other functions have been developed to assist researchers to experiment with approximations. We review their prototypes here. * (in approx.a) Findsupportdata[support, option, totalareal returns the mean value and probability of partitions of a rectangular support, assuming multivariate uniform distributions, support and totalarea are defined above, and option is 1 for the first call to the function,where and 2 if not. * (in approx.m) FindPartitionCSup, chi, hconst, options] returns the partition number in Sup (of the form SUPP defined above) which contains the right-handside (Chi, hconst). If the RHS is not in Sup, the function returns zero. The only option for this function is PrintLevel. 12

* (in approx.m) FindEncSimplexCbasept,Basis,,b,options] returns a list of n+ 1 points which describe a simplex which encloses the region K = {zlAx > b, z > 0}. The function requires a starting extreme point of K called basept and a list of columns defining basespt's basis. The only option for this function is PrintLevel. * (in approx.m) FormSimplexCsontraints[Simplex, options] returns a feasible region of the form K = {xIAx > b} which describes the simplex Simplex (a list of n + 1 points). The returned compound list is of the form {A,b}. Again, the only option for this function is PrintLevel. 3.2 L-Shaped Method The most common specialized algorithm for solving two stage stochastic linear programs is the L-Shaped method (Van Slyke and Wets [4]), which is based on Bender's decomposition. The algorithm, which is also approximates the value of the recourse function (second-stage solution) by finding cuts, which are constraints on both the approximation and the decision variables. A full discussion of this algorithm is beyond the scope of this report, so we just describe the main functions and their subfunctions. There are two main functions which implement the two-stage cutting plane algorithm. The first solves discrete problems and the second uses the approximations described above to solve problems with uniformly distributed RHSs. The prototype for the first algorithm is: (in cpa.m) CutPlane[Supin,hconst,A,b,c,W,T,q, options] Inputs Supin is a list of realizations and their probabilities of the form {{hi,...,hk},{p,...,pk}} A,b,c First Stage data. W,T,q Deterministic Second stage data. hconst Deterministic RHS for second stage. ____ _(Deterministic rows come first in Wy = h - Tx). Options PrintLevel->i is the level of debugging detail. The function returns the optimal first stage decision z. Of course, the primary use of this function is to modify it to output other items of interest or try variants on the algorithm. The second implementation uses the Jensen and Edmunson-Madansky upper bounds in approx.m to solve a two stage stochastic linear program with uniformly distributed RHSs. The prototype for the function is (in cpa.n) CutPlanerApproxtaxits,Supin,hconst,A,b,c,W,T,q, options] Inputs maxits is the maximum number of complete iterations allowed. Supin is a support of the form SUPPTS. A,b,c First Stage data. V,T,q Deterministic Second stage data. hconst Deterministic RHS for second stage. (Deterministic rows come first in Wy = h - Tx). Options PrintLevelAlg->1 is the level of algorithmic debugging detail. PrintLevelBounds->l is the level of bounding debugging detail. CutPlanewApprox also returns the optimal first stage decision z. One of the most lucrative improvements suggested for the L-Shaped algorithm dur to (Birge and Louveax, [16]) is the inclusion of multiple cuts at each iteration. This modification has been implemented in cpa.m, and may be run using the functions: * (in cpa.m) HultiCutPlanewApproxCmaxits.Supi.n,hconst,A,b, c,W,T,q, options]

* (in cpa.m) MultiCutPlane[Supin,hconst,A,b,c,W,T,q, options] The inputs, options, and outputs are the same as in CutPlane and CutPlanewApprox. 3.3 Miscellaneous Functions To date, only one function has been written which cannot be classified into the other categories. This function (and future ones like it) is in stopro. m and finds the size (in term of rows and columns) of a deterministic representation to a multistage stochastic program. The function returns a list {matrices, nodes, rows, cols} with the total number of matrices, nodes, rows, and columns in the determinsitic equivalent. The prototype and arguments are (in stopro.m) DetEquiv[scenlist, pds, Adims, WMdims, WEdims] Inputs scenlist is a list with the number of distinct realizations in each period. pds is the number of periods. Adim is a list with the dimensions of the first period matrix A. WMdims is a list with the dimensions of the middle period matrices. WEdims is a list with the dimensions of the last period matrix. 4 Stochastic Program generation and manipulation The last category of functions concerns the manipulation and generation of large multistage stochastic programs (MSLPs). Currently, there are two known programs for solving large MSLPs, MSLiP (Gassman, [14], [15]) and ND, (Donohue, [17]). MSLiP will accept stochastic programs in a standard format (Edwards, et. al. [18]), while ND uses a more customized format. The functions in makesto.m will take an arbitrary linear program in MPS format (see section 2), and, given a description of stochasticity and the coordinates necessary to split the problem into periods, generate the input files necessary for both programs. The presentation of these functions will be split into three parts: reading in data, forming stochastic descriptions, and writing output files. A full example will be presented in a fourth section. 4.1 Reading in Data The basic premise of the makesto package is that a multistage stochastic linear program can be split into periods, and that constraints on first period decisions, intermediate period decisions, and last period decisions may all be different. The constraints for all periods are generally contained in one MPS file (called a CORE file). Since decisions in each period cannot depend on future decisions, the CORE file must be lower block triangular. The functions in alkesto.m make the additional assumption that the problems are Markovian in nature. This implies that the CORE file will actually be block diagonal. Most practical problems will have CORE files that describe at most the first, intermediate, and final periods. (CORE file that describe more than that may be cut down for makesto.m's purposes. As a result, the only information beyond the CORE file needed to extend the stochastic program is the upper left-hand coordinate of the Repeating Period and the Ending Period. (These coordinates will henceforth be abbreviated as repc and endc.) In standard format, these coordinates are expresseed as row and column names, and are in a separate file called a TIME file. In some cases, constraints on all decisions in periods 2,..., T may all be the same. If this is the case, then the repeating period and the ending period have the same starting coordinates. Assuming that ReadMPS (described above) is used to read in the CORE file, repc and endc may be found by looking at the matrix with LPPlot, or by using FindCoords, which finds the coordinates of a named column and row, or by reading an existing TIME file. The following two functions may be useful in this regard. 14

(in makesto.m) FindCoordsELP, row, col] Inputs LP is a complete linear program in (LP) format. row is a string (8 characters or less) containing a row name. col is a string (8 characters or less) containing a column name. (in makesto.m) ReadTimeFileC[ilename, colnames, rownames] Inputs filename A TIME filename. colnames,rowames are lists of strings obtained from ReadMPS. The output of the first function is the coordinates of the desired element, and the output of the second function is a list of starting coordinates for each time period. Once these coordintes have been obtained, SplitCore can be used to separate out the CORE file into its constituent components. (in makesto.m) SplitCore[LP, RepeatPeriodStart, EndPeriodStart, options] Inputs ___ LP is a complete linear program in (LP) format. RepeatPeriodStart is the upper left hand corner of the repeating period constraints. EndPeriodStart is the upper left hand corner of the final period constraints. Options WriteLPform->True Writes each block as a full LP in (LP) format. AddNames->True includes row and column names in the output LPs. Extending stochastic programs may also require input of an existing description of stochasticity. In standard format, stochsticity is described using a number of formats in a STOCH file. Currently, only one format can be read using the functions in makesto.m. This format is the BLOCKS format, which is used to described problems with stochastic T (off-diagonal) matrices. The function ReadBlockFile, described below, will return a list of the form {{period, {{pl,T1},..., {pn,Tn})),... }. Here, the technology matrix in period period takes on one of n different realizations T., each with probability Pi. (in makesto.m) ReadBlockFile[ filename, nperiods, pdnamelist, LPlist] Inputs filename is the STOCH filename to be read. nperiods is the number of periods in the STOCH file. pdnamelist is a list of period names obtained from a TIME file. LPlist is a list of LPs, obtained directly by running _SplitCore with WriteLPfora->True. 4.2 Describing stochasticity Stochasticity in either the RHS or the technology matrix Ti can be described in several ways. To make these descriptions as succinct and as easy to use as possible, the following options have been developed: * IndependentRfS, which requires as input a list of the form {..., {t, i, {{bi(1),p( pi )},..., {b,(k), pi(k)}}},...}. Each RHS bi in period t is independently distributed with the takes on the value bi(k) with probability Pi(k). * SampledRHS which is a list of the form {seed, {t, p, n, b, e, iu},..} where seed is a random number seed, t is the period for the sample, p is the maringal probability for each item in the list, and n is the number of samples generated (each with probability p/n). Samples are generated using the formula bi = (1 + a)bi where a is unformly distributed between ti and ui. This option is best for generating huge problems. 15

* Blocks which is a list of the form {I..., {t, {P1, T1 )},..., {pk, Tk}},...} where T, is a matrix realization in sparse form. This option is the same as the BLOCKS format in STOCH files. The lists described above are Maihemaiica lists. In most cases, BLOCKS lists will be generated using ReadBlockiile, SampledRHS lists will be generated by using ReadList or ReadMPS to read an existing RHS, and IndependentRHS lists will be generated manually, or by using the following function (in makesto.m) MakestochlistC period, element, mean, diff, n, options ] Inputs period,element The stochastic element. mean is a base value from which to generate realizations. diff is an increment value (see below) n is the number of realizations. Options Range->{lo,hi} is a range specification for values. Realizations are usually incrementally generated using the formula mean + i*diff, where i is incremented until n realizations have been formed. In other words, diff is successively added and subtracted from the highest and lowest realizations until n realizations have been formed. The result returned is this list of numbers. If Range has been specified, then mean and difi are ignored, and n relizations uniformly spaced between the ranges specified by Ranges are returned. 4.3 Outputting stochastic programs Once the split stochastic program has been loaded into memory and a suitable definition of the stochasticity has been created, the following functions will form an arbitarily sized stochastic program and output it in standard format. (in makesto.m) MSCoretoCore[ infile, outcorefile, outtimefile, repc, endc, _____________stages, options] Inputs inmile is an input CORE file. outcorefile is an output CORE file with stages stages. outtimefile is an output TIME file with stages stages. repc, endc are the repeating and ending block coordinates. stages are the number of stages to output. Options SplitCoreLP->{} An existing split core from SplitCore CoreLP->'} An existing core LP from ReadMPS. The result returned will be True unless an error occurs. The companion function to NDCoretoCore is WriteStochiile. (in makesto.m) WriteStochfileE outstochfile, stages, stochlist, options] Inputs outstochfile is an output STOCH file with stages stages. stages are the number of stages to output. stochlist is a list describing the stochasticity. Options StochlistForm->"IndependentRHS" Form of st9ohlist.

The multistage solver ND accepts problems using two files, a fort file, which contains diagonal constraint matrices Wt for each period t, and a scendata file, which contains the T matrices and RHSs. The scendata file does not require all nodes in a decision tree, but only those which can be repeated for a given stage. (in makesto.m) HSCoretoMD[ infile, outfortfile, outscenfile, repc, endc, stages, stochlist, options] Inputs infile is an input CORE file. outfortfile is an output fort file with stages stages. outscenfile is an output scendata file with stages stages. repc, endc are the repeating and ending block coordinates. stages are the number of stages to output. stochlist is a list describing the stochasticity. Options____ StochlistForm->"IndependentRHS" Form of stochlist. SplitCoreLP->{} An existing split core from SplitCore CoreLP->{} An existing core LP from ReadMPS. Again, the result returned is True if all files were sucessfully written, and False otherwise. 4.4 An example As a brief demonstration for how these routines fit together, the following notebook fragment will read in a stochastic program in standard format and output several versions of increasing size. Input files and partition <<d:\r\math\makesto.m <<d: \r\probs\ndprobs\agr7\stochs.m CORE=ReadMPS [str=OpenRead["e: core"], Sparse->True]; Close [str]; repc={16,21}; endc={35,41}; SPLITCORE=SplitCore(CORE,repc, endc]; Off [lumberForm::sigz]; Output 4 Stages NSCoretoCore "e: core", "e:p4. cor", "e:p4. ti",repc, endc,4, StochlistForm->"IndependentRHS", SplitCoreLP->SPLITCORE]; WriteStochfile " e: p4s2. sto", 4, s4s2, StochlistForl->"IndependentRHS"] WriteStochfile " e:p4s4. sto",4,s4s4, StochlistFor->" IndependentRHS"] WriteStochfile ["e:p4s8. sto", 4,s4s8, StochlistForu->"IndependentRHS"] WriteStochfile ["e:p4s16. sto", 4, s4s16, StochlistForm->"IndependentRHS"] WriteStochfile["e:p4s32. sto",4,s4s32,StochlistForm->"IndependentRHS"] The file stochs.a contains several lists describing the stochasticity to be modelled: (* --- 3 stage ---------------------------— _ —---- *s3s2:= List A { Makestochlirst 1,5,1700,200,2], Makestochlist [2,11,1700,200,2]}; (* --- 4 stage — --------------------------------- *) s4s2:= List ~C { Hakestochlist [1,5,1700,200,2], Makestochlist[2,5,1700,200,2], Makestochlist[3,11,1700,200,2 }; s48:= Likst h 3,111 { Makestochlist [1,5,1700,125,8, akestoc9ist 2,5, 1700,125,8],

Makestochlist 3, 11, 1700,125,8] }; s4s16:= List 06 { Makestochlist[1,5,1700,75,16], Makestochlist[2,5,1700,75,16], Makestochlist[3,11,1700,75,16 }; s4s32:= List 00 { Makestochlist[1,5,1700,50,32], Makestochlist[2,5,1700,50,32], Makestochlist[3, 11, 170,0,532]}; References [1] Wolfram, S., 1992. Mathematica, A System for doging mathematics by computer., Addison Wesley, New York, NY. [2] B.A. Murtaugh and M.A. Saunders, 1983. "MINOS 5.1 User's Guide," Technical Report SOL-83-20R, Systems Optimization Laboratory, Stanford University, Stanford, California. [3] International Business Machines Corp., 1991. "Optimization Subroutine Libary Guide and Reference, Release 2," document SC23-0519-02, International Business Machines Corp. [4] Van Slyke R., and R. Wets, 1969. " L-shaped linear programs with applications to optimal control and stochastic programming," SIAM Journal of Applied Mathematics, 17, pp. 638-663. [5] Ermoliev, Yu. and Wets, R. 1980. "Stochastic Programming, an Introduction" in Numerical Techniques for Stochastic Optimization, Y. Ermoliev and R. J.-B. Wets (eds.), Springer-Verlag, Berlin. [6] Birge, J., 1985. "Decomposition and Partitioning Methods for Stochastic Linear Programs, " Operations Research, 33, 989-1007. [7] Bazaraa, M. S., J. J. Jarvis, and H. D. Sherali, 1990. Linear Programming and Network Flows, Second Edition, Wiley and Sons, New York, NY. [8] Charnes, A. and W. W. Cooper, 1962. "Programming with Linear Fractionals," Naval Research Logistics Quarterly bf 9, pp. 181-186. [9] Bazaraa, M. S., and C. M. Shetty, 1979. Nonlinear Programming, Wiley and Sons, New York, NY. [10] Kall, P., A. Ruszczyhski, and K. Frauendorfer. 1980. "Approximation Techniques in Stochastic Programming" in Numerical Techniques for Stochastic Optimization, Y. Ermoliev and R. J.-B. Wets (eds.), Springer-Verlag, Berlin. [11] Kall, P., A. Ruszczynski, and K. Frauendorfer. 1980. "Approximation Techniques in Stochastic Programming" in Numerical Techniques for Stochastic Optimization, Y. Ermoliev and R. J.-B. Wets (eds.), Springer-Verlag, Berlin. [12] Birge, J. R.'Interior Point Notes," Class Materials, Industrial and Operations Engineering 510, Ann Arbor, MI 48109. [13] Birge, J. and R. Wets. 1986. "Designing approximation schemes for stochastic optimization problems, in particular for stochastic programming problems with recourse." Mathematical Programming Study, 27 54-102. [14] Gassman, H.I., 1990.' MSLip: A Computer Code for the Multistage Stochastic Linear Programming Problem, " Mathematical Prolmming 47, 407-423. [15] Gassman, H.I., 1990. " MSLip User's Guide, " Working Paper WP-90-4, School of Business Administration, Dalhousie University, Halifax, Canada. [16] Birge, J. R. and F. Louveuax, 1988. "A multicut algorithm for two-stage stochastic linear programs," European Journal of Operations Research 34, pp. 384-392. [17] Donohue, C., 1993. Private communication, techincal report under preparation, University of Michigan, Ann Arbor, MI 48109. [18] Edwards, J.. 1980. "A proposed Standard Input format for computer codes which solve stochastic programs wit recourse. " in Numerical Techniques for Stochastic Optimization, Y. Ermoliev and R. J.-B. Wets (eds.), Springer-Verlag, Berlin. 18