A x ru, and bl
x bu. Make sure that the declaration of the size of your DSPACE array matches the value of the NWORDS parameter in the call to EKKDSCA. Note that the value of NWORDS can be less than the number of entries in DSPACE, but the application will use only the first NWORDS of the DSPACE array. This technique may be used to have multiple applications within the same workspace.
The library does not include a subroutine that gives execution times. However, some limited timing information can be obtained by setting the real control variable Rprintcpu. The default value for Rprintcpu is zero, and while Rprintcpu is zero, message number EKK0197I will not be printed. However, if you reset Rprintcpu to some positive number, this new value will be interpreted as a limiting execution time interval for library routines, beyond which timing information will be output by EKKPRTS, the standard Library output module. If Rprintcpu has a positive value, then message number EKK0197I will be generated after the return from each library subroutine whose execution time in seconds exceeds the value of Rprintcpu.
Message EKK0197I gives the total execution time, the name of the subroutine that ran for more than Rprintcpu seconds, and the execution time for this routine. The times are only given to two significant digits, so setting Rprintcpu smaller than 0.01 might be of doubtful utility. To obtain execution times, FORTRAN and C callable system timing routines are used on the various operating systems for which the library is available. These routines are not all alike, and thus there may be some variation in timing output under different operating systems. In addition, resetting Rprintcpu in the midst of a computation may adversely affect reported execution times.
Be careful about when and where you set control variables. Setting them prior to making a call to EKKDSCA is allowed and will create a template for each call to EKKDSCM within that application. This is helpful if, for example, you have multiple maximization models and do not want to set Rmaxmin in each of those models. Note that setting control variables after a call to EKKDSCA, but before a call to EKKDSCM is fruitless; the variables are reset when EKKDSCM is called. Also note that it is possible to let EKKMPS or EKKLMDL call EKKDSCM for you (with a default of 1 block), but the same care with control variables is needed. That is, you must set them either prior to calling EKKDSCA or after calling one of: EKKDSCM, EKKLMDL, and EKKMPS.
Make sure that the name in the NAME record of the EKKQMPS input file matches the name of the problem as specified by the character control variable Cname. If you have read in the linear part of a quadratic programming problem with EKKMPS and you have not subsequently changed the value of Cname, the problem name will be given by the character string in the NAME record of the MPS file from which EKKMPS has read the linear part of the problem data. In addition to Cname, EKKMPS sets several other character control variables. For more information, see the notes following the description of EKKCSET.
EKKMPS will set bounds on columns to defaults, if column bounds are not specified in your MPS file. The defaults for integer variables are 0 and 1. The default bounds for continuous variables are 0 and 1031. Refer to "Default Bounds" for more information.
Be careful when using DSPACE and MSPACE with array indices obtained using EKKNGET. MSPACE is an INTEGER*4 array which is equivalenced to DSPACE. Note that for indices into REAL*8 arrays, the index will apply to DSPACE. For example, reduced costs are REAL*8, so Ncolrcosts is a FORTRAN index into DSPACE. However, for INTEGER*4 arrays, the index is into MSPACE. Thus, Nrowstat is a FORTRAN index into MSPACE, since the status arrays are INTEGER*4.
To suppress the EKKwxyz message numbers, do the following: First, use EKKIGET/EKKISET to set integer control variable number 25, Imsgpos to 1. This moves the messages to the right side of the output line. Then, call EKKMSET and set NONUM=1 to turn off the message numbers.
Severe error messages are sometimes generated because the size of the DSPACE declared is not large enough to run the problem. If you receive a message stating that additional DSPACE is needed or that DSPACE is being overwritten, try increasing the size of DSPACE. Also note that you can track DSPACE usage with EKKSMAP, or by accessing index control variables Nfirstfree and Nlastfree with EKKNGET. This can be done within your application program or within user exits.
When calling EKKQSLV with INIT=1, the algorithm will go through a decomposition scheme to get close to the optimal solution. Within this decomposition, EKKSSLV is called a number of times. Note that iteration counts may be confusing. There are actually three iteration counts: one for the number of major QP iterations, one for the number of EKKSSLV iterations within the decomposition, and one for the number of QP iterations after exiting from the decomposition. Note that changing integer control variable number 5, Imaxiter will only impact the last type of iterations.
This section discusses the dual variables generated by EKKSSLV and demonstrates how EKKNGET can be used to locate the dual vectors within the workspace. It also outlines the linear programming theory and explains the role that each component of these vectors plays in the dual linear program.
Every linear program (P) has an associated dual linear program (D) with the same optimal objective value as (P), provided that (P) has an optimal solution. The linear program addressed by modules of the Optimization Library can be expressed as follows:
minimize c T x
subject to: rl A x ru, and
bl x bu.
where:
The associated dual problem is:
maximize dlT rl - du T ru + ql T bl - qu T bu
subject to: dlT A - duT A + ql - qu = c, and
dl, du, ql, qu, 
0.
where:
The dual variables dl, du, ql, and qu can be used in post optimal analysis to calculate, for example, the increase in overall cost that would result from perturbing ru or rl.
All the dual variables are provided by EKKSSLV, but linear programming theory permits the vectors dl and du to be stored in one vector of length m and the vectors ql and qu to be stored in one vector of length n. The "complementarity conditions" of linear programming state that:
du(i) = 0 if Ax(i) < ru(i) (note the strict inequality) dl(i) = 0 if Ax(i) rl(i) qu(j) = 0 if x(j) < bu(j) ql(j) = 0 if x(j) bl(j),
Any optimal set of variables (x, du, dl, qu, ql) must satisfy these conditions. This means that for any optimal solution x, the corresponding optimal dual solution (du, dl, qu, ql) is such that for each i either du(i) = 0 or dl(i) = 0, and thus one dual vector d() can contain all the nonzero entries of both du() and dl(). Similarly, one vector q(), can contain all the nonzero entries of both qu() and ql()).
Given a vector d() that contains du() and dl(), a simple way to determine whether d(i)=du(i) or d(i)=dl(i) is to see if Ax(i)=ru(i). If so, then d(i)=du(i).
A program that evaluates the dual objective function by accessing the dual variables in this way is shown below. It is important to note that library modules negate some of the dual variables. Also note that the inner-products are calculated in such a way that equality constraints must be handled separately.
Being able to reference the dual solution from a FORTRAN program is one of the advantages of using the Optimization Library. "Sample FORTRAN Program EXDUAL" shows precisely how to access the dual variables.
Extra care should be taken using the procedure described above, which is considered advanced use of library modules.
Most mixed-integer programming problems involve constraining some or all of the variables in an LP to take on integer values. Problems of this type have many applications. For example, building materials and packing crates often must be purchased in whole number units, and a mixed-integer programming (MIP) problem must be solved to minimize their consumption. Also, "yes or no" decisions must often be made to decide, for example, whether a machine will be used to produce product A or product B. For the many users with problems like these, sample programs are provided with the download material, and some useful sample programs are listed below. Much of this discussion, however, deals with the more delicate problem of defining special ordered sets (SOSs) in a library application. Generally speaking, these are used to model non-convex objective functions and constraints, and the flexibility of EKKIMDL makes it possible to define SOSs as well as tune the EKKMSLV solution process.
For most users the best way to get started with EKKIMDL is to look at the example program EXIMDL3. In this example, certain variables from an LP are identified as integers, and the call to EKKIMDL puts a list of the variables and other needed information in the library workspace. The sample programs EXIMDL and EXIMDL2 give further examples showing the use of EKKIMDL to impose integrality on LP variables.
Special ordered sets are defined in "Special Ordered Set (SOS)" and in Appendix D. "Glossary". EKKIMDL can be used to define SOSs of type 1, 2, or 3. In addition, regular integer variables (i.e. those not in a SOS) are defined by using EKKIMDL and specifying type 4. If the lower bound for a variable is 0 and the upper bound is 1, then the variable is a 0-1, or decision, variable and may be used to model a "yes or no" decision. Problems with many 0-1 variables are usually best solved by EKKMSLV when preceded by a call to EKKMPRE. If an integer variable does not have the 0-1 bounds, the variable is known as a general integer variable, and it may take on any of the integer values between its bounds in an optimal solution. Integer variables that are not in a SOS have up- and down-pseudocosts associated with them. Pseudocosts are defined briefly in the glossary. Pseudocosts are always positive, and are an estimate of the degradation in the objective function that will result from imposing integrality on the corresponding variable.
The example program EXMSLV3 shows a technique involving special ordered sets of type 1 that may be used to minimize the polynomial (2 x2 y + xy2 + 25 x2 +100y2 + 20x - 200y) over the square region -2<=x<=2, -2<=y<=2. A minimizer (x,y) must lie on a grid that spans the box. The non-separability of the objective function is handled through a linking constraint in the program.
Special ordered sets of type 2 may be used for minimization along a piecewise linear curve Px(x). This is done in the sample FORTRAN program EXIMDL4 and in the sample C program EXIMDL6. The piecewise linear curve defined on the interval [0, 10] by: (0,2), (7/10,1), (1,5), (4,-1), (5,-3) and (10,4) (where the abscissas give x, and the ordinates give Px(x)) may be modeled by a special ordered set of type 2 with 6 members. The reference row entries for the members are (0, 7/10, 1, 4, 5, 10), and the objective function for the members is given by (2, 1, 5, -1, -3, 4). The reference row entries define an uneven grid, and a "value row", described in the code below, can be added that permits x to take on values between the grid points. This is in contrast to EXMSLV3, where the minimizer is constrained to lie on one of the grid points. EXIMDL4 and EXIMDL6 use two special ordered sets of type 2 to model two piecewise polynomials. The second SOS was added along with the constraint x+y<=4.5 so that the minimizer would lie between grid points.
EKKSOS should be called after model data has been loaded into the workspace, and before calling any of EKKSCAL, EKKPRSL, EKKMPRE, and EKKMSLV. EKKSOS attempts to add special ordered sets of integers. If such sets are found, additionial space may be required for the internal array in which integer information is stored. Therefore, when using EKKSOS, it is wise to allow extra space for this information. You can do this by setting the integer control variable Imaxintinfo (number 53) to the negative of an estimate of the extra space required. This should be done after the call to EKKDSCM, but before any calls to EKKMPS, EKKLMDL, EKKIMDL, or EKKSMDL. The necessary increment to the value of Imaxintinfo cannot exceed the number of nonzero coefficients of integer variables in the constraint matrix, but this may be a difficult datum to obtain, and it may still be a generous estimate. An even more generous estimate, which would be more easily obtained, would be the total number of nonzero coefficients in the matrix. However, before using this estimate, be aware that space set aside for holding integer information will not be used for another purpose.
Refer to the description of EKKSOS for more information on the calling sequence, and to the Sample FORTRAN Program EXSOS for an example of the use of EKKSOS with EKKMSLV.
The transformation in which the coefficients of the 0/1 variables are the powers of 2 is referred as the LOG2 transformation, and the transformation in which the coefficients of the 0/1 variables are the natural numbers is referred to as a SOS transformation. The choice of transformation is specified by a parameter in the calling sequence of EKKBMPR. The LOG2 transformation is recommended, because it results in fewer 0/1 variables. Our test results indicate that having a smaller tree to search is usually more advantageous than having special ordered sets on which to branch.
To use EKKBMPR, one must allocate additional workspace for the transformed 0/1 model. This should be done in an initial call to EKKDSCA that specifies (at least) one model more than the number of models being processed. EKKBMPR attempts to create a new copy of the transformed constraint matrix stored in vector-block format, overwriting the old matrix. If there is not enough space to create this copy, EKKBMPR returns immediately, with a non-zero return code.
Since it is the transformed problem that is solved, EKKBMPR should be called after the general integer model has been loaded into the workspace, and before any initial solution oriented processing has been done. In particular, the setting of control variables for EKKMSLV, and calls to other preprocessing routines, such as EKKPRSL, EKKSCAL, and EKKCRSH, should be done after EKKBMPR has made the transformation. We point out (the obvious) that all MIP user exit calls will then apply to the transformed model, and not to the original one.
Similarly, EKKBMPS should be called after the transformed problem has been solved by EKKMSLV (and postsolved by EKKPSSL, if EKKPRSL was called after EKKBMPR). After a normal return from EKKBMPS, both bounds of each of the original general integer variables will have been reset to the solution value for that variable.
Communication between EKKBMPR and EKKBMPS is via a file that is specified in the calling sequences of both routines by a FORTRAN unit number. Of course, a file accessed by EKKBMPS must be a file that was created by EKKBMPR, or by a user supplied substitute that mimics EKKBMPR. Since Optimization Library modules are written in C, one can not use the FORTRAN open statement to associate a named file with a unit number to be supplied to EKKBMPR or EKKBMPS. The library module EKKFOPN is provided to mimic the FORTRAN open statement without using the FORTRAN operating environment. If a named file has not been explicitly opened and associated it with unit by EKKFOPN before a call to EKKBMPR (EKKBMPS), then EKKBMPR (EKKBMPS) opens and uses the file named fort.xx, where "xx" is the value of unit. Each of EKKBMPR and EKKBMPS closes the file associated with unit before returning.
See the sample C program EXBMPR for an example of the use of EKKBMPR and EKKBMPS with EKKMSLV.
Designing your own branching strategy usually requires an in-depth understanding of the problem being solved. In particular, it usually requires knowledge of the variables that can have a large impact on the overall cost of an optimal solution. Alternatively, knowledge of the variables that can impact feasibility can lead to a good branching strategy.
Assuming that you know the details of the problem you're solving, this article can help you to write an EKKBRNU user exit. The example user exit given below is for a 0-1 problem for which the constraint matrix has a "lower triangular" form. The strategy is to branch first on variables at the left end of the constraint matrix (i.e. variables with low index numbers). This way the variables at the left are fixed early in the branch and bound process to achieve feasibility of the top constraints. They are then left unchanged for the rest of the optimization, since their values are fixed. This technique was used to solve in a few hours a pure 0-1 problem that had not been solved in three days using the default branching strategy.
Branch selection is done immediately after solving the relaxed LP at a node. During the selection process, EKKBRNU is called repeatedly to allow user code to override various choices made. The calls are made in the embedded fashion indicated in the section Mixed-Integer User Exit Subroutines.
Before beginning the branch selection, EKKBRNU is called with IREASON=1. This permits user code to initialize some "holder" variables. For instance, a branching strategy that involved searching for the non-integer variable with the most positive cost would require initializing some real control variable, Rmaxcost, to minus infinity. This particular strategy is outlined in the sample code below.
For each set of integer variables (4), EKKBRNU is called with IREASON=2. This allows the user code to save information about the type and priority of each set being processed. Priorities may be set at this time, but the effect of the change is not felt until the next round of branch selection. The user return code should be set to 0 here to switch on the calls to EKKBRNU with IREASON=3.
For each variable in the set indicated by the most recent call with IREASON=2, EKKBRNU is called with IREASON=3. This call allows user code to change pseudocosts or reference row entries associated with each variable. It also allows the user code to change the branch direction for each variable. The branch direction for a given variable indicates which branch will be evaluated first, when the nodes created by branching on this variable are evaluated. Probably the most important user coding associated with IREASON=3 is coding that resets a pointer to the best branching variable. The sample user exit below shows a number of possible strategies that might be used to select the best variable, and it also shows some important implementation details. First of all, the user code includes a check to make sure that it does not select a satisfied integer variable as the next branching variable. Such a choice would be rejected, and a program halt would occur. Second, the user code illustrates that the fifth element of MARRAY is an internal number corresponding to the integer variable under consideration. To map the fifth element of MARRAY to the index of the variable in the original problem formulation, use the index control variable Nintinfo as shown.
For each set of integer variables, EKKBRNU is called with IREASON=4 to permit the user code to change the choice of branching variable. If the choice of variable is to remain unchanged, but the branching direction needs to be changed, the variable number may be negated. Finally, there is a single call to EKKBRNU with IREASON=5. During this call, the user code should load the number of the set containing the next branching variable into the first element of MARRAY.
C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX C A BRANCH USER EXIT FOR LOWER TRIANGULAR CONSTRAINT MATRICES C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SUBROUTINE EKKBRNU(DSPACE,MSPACE,IREASON,MARRAY,DARRAY,JRTCOD) REAL*8 DSPACE(*),RMAXCOST REAL*8 DARRAY(5),RINFEAS,RMAXINFEAS,RMININFEAS INTEGER*4 MARRAY(6),ISET,RTCOD,JRTCOD,IBEST_INIT,IREASON,ISAVE INTEGER*4 MSPACE(*),I_INT_NUM,M2,IBEST SAVE IBEST,RMAXINFEAS,ISET,RMININFEAS,RMAXCOST SAVE IBEST_INIT,ISAVE INCLUDE (OSLI) INCLUDE (OSLR) INCLUDE (OSLN) IF (IREASON.EQ.1) THEN C This call to EKKBRNU permits the initialization of C various "holder" variables. CALL EKKRGET(RTCOD,DSPACE,OSLR,OSLRLN) CALL EKKNGET(RTCOD,DSPACE,OSLN,OSLNLN) C RMAXINFEAS = 0.0D0 C RMININFEAS = 9999.0D0 C RMAXCOST = -9999.0D0 IBEST_INIT = 9999999 IBEST = IBEST_INIT JRTCOD = 0 GOTO 6000 ENDIF IF (IREASON.EQ.2) THEN C This call permits the user to save information about the set. M2 = MARRAY(2) IF((M2.EQ.1).OR.(M2.EQ.2).OR.(M2.EQ.3)) THEN write (6,*) 'THIS EKKBRNU NOT FOR SOS' STOP ENDIF C Setting JRTCOD=0 causes EKKBRNU to be called with IREASON=3 JRTCOD = 0 GOTO 6000 ENDIF IF (IREASON.EQ.3) THEN C Get the number of the variable. The value in MARRAY(5) is the C internal integer number. I_INT_NUM = MSPACE(NINTINFO+4*(MARRAY(5)-1)) C Compute integer infeasibility for this integer. RINFEAS = DABS(DARRAY(1)-DNINT(DARRAY(1))) C Don't attempt to branch on satisfied integer variables. IF (RINFEAS.LT.(RTOLINT)) GOTO 6000 C Chose the variable with the largest integer infeasibility. C IF (RINFEAS.GT.RMAXINFEAS) THEN C Chose the variable with the smallest integer infeasibility. C IF (RINFEAS.LT.RMININFEAS) THEN C Chose the variable with the maximum cost. C IF (DSPACE(NOBJECTIVE+I_INT_NUM-1).GT.RMAXCOST) THEN C Chose the variable with the smallest index to utilize the lower C triangular form of the matrix. IF (I_INT_NUM.lt.IBEST) THEN C RMAXCOST = DSPACE(NOBJECTIVE+I_INT_NUM-1) C RMININFEAS = RINFEAS C RMAXINFEAS = RINFEAS IBEST = I_INT_NUM ISAVE = MARRAY(5) ISET = MARRAY(1) C write (6,*) ' selecting up ' MARRAY(6) = 1 C write (6,*) ' selecting down ' C MARRAY(6) = 0 ENDIF GOTO 6000 ENDIF IF (IREASON.EQ.4) THEN C Load the (negated) internal number of the selected branching C variable into MARRAY(5). The negation ensures that the C selected branching direction will be adopted. IF (IBEST.NE.IBEST_INIT) MARRAY(5) = -ISAVE IF (MARRAY(5).GT.0) THEN write (6,*) ' default branch: ', MARRAY(5) ELSE write (6,*) ' selected branch: ',-MARRAY(5) ENDIF GOTO 6000 ENDIF IF (IREASON.EQ.5) THEN C Load the (negated) set number for the selected branching C variable into MARRAY(1). The negation is optional. IF (IBEST.NE.IBEST_INIT) MARRAY(1) = -ISET GOTO 6000 ENDIF 6000 CONTINUE RETURN END
A wide variety of block structured problems problems can be solved with a call to EKKLPDC followed by a call to EKKSSLV. In effect, EKKLPDC provides a sophisticated crash procedure for the LP solver that is much easier user to code than a crash procedure involving multiple models. This section gives descriptions of problems that are candidates for use with EKKLPDC, and performance tuners that may provide important additions to application code. Code fragments are given to show techniques for specifying problem structure and for tuning EKKLPDC, and results for a few test problems are given. The subsections are:
Many large-scale LP problems have constraint matrices that have an identifiable block structure. The blocks of the matrix may be connected by coupling rows, coupling columns, or neither.
An example of the first type might be a traffic planning problem in which the morning commute time for an entire city is to be minimized. The blocks of the constraint matrix would correspond to the network of accessible traffic routes for all the suburbs. In addition to the blocks, there would be a set of coupling constraints that impose, for example, capacity constraints that guaranteed that the sum of the traffic from all suburbs does not exceed the capacities of certain key routes. Problems with this block-angular structure are referred to here as Dantzig-Wolfe problems. The blocks are composed of disjoint sets of rows and columns, while the coupling rows may have non-zeros in the columns of multiple blocks. Multi-commodity flow problems often have Dantzig-Wolfe structure, but the number of coupling rows is comparable to or even exceeds the total number of rows in the separable blocks. Such problems require special consideration.
An example of the second type of block structured problems could be an inventory problem in which the distribution of goods from several supply points to many demand points must be optimized for each month, but demand for one month may be satisfied, to some extent, by supply from a previous month. A linear inventory cost would be incurred for each unit of the good that was produced in one month and distributed in the next month. Such a problem gives rise to a constraint matrix with blocks of constraints corresponding to the distribution network for each of the months. The inventory of goods that remain in stock at the end of a month may be modeled by columns with non-zeros in rows of two adjacent blocks. The complete constraint matrix has a staircase structure made up of blocks with coupling columns that "connect" two adjacent blocks.
The third type of block structured problem would be a completely decomposable problem, for which the constraint matrix has a disjoint block structure in which neither coupling rows nor columns "connect" adjacent blocks. These completely separable models are sometimes generated by modeling languages that provide substantial ease of use, but do not produce the most succinct of possible models. The optimization problems represented by such models can be completely solved by sequentially solving the smaller, uncoupled subproblems represented by the individual blocks.
To make good use of the EKKLPDC module, you should know that your problem has a constraint matrix with one of these common underlying block structures, and be able to identify the type of structure to the subroutine. The type of structure to be identified is communicated to EKKLPDC in the third parameter (type) of its five parameter calling sequence.
Besides the usual (rtcod, and dspace), the parameters in the calling sequence of EKKLPD are: type, strategy, and nblocks. The value of the type parameter identifies the type of structure to be found as follows, 1 => staircase, 2 => Dantzig-Wolfe (only one iteration), 3 => Dantzig-Wolfe (multiple iterations, to optimality or to the limit imposed by the value of integer control variable number 5, Imajorits). In theory, the constraint matrix must have either staircase or Dantzig-Wolfe structure in order for a call to EKKLPDC to be appropriate. However, in fact, one can (in theory, at least) always separate a sparse matrix into a staircase structure by identifying enough "coupling rows." At the other extreme, completely separable problems are also identified to EKKLPDC as staircase problems.
The strategy parameter may take values 0 - 3. How these are to be interpreted depends on the value of the type parameter. For an explanation of the twelve possible combinations of values of type and strategy, see the documentation of the EKKLPDC module.
If the value of nblocks is 0, then EKKLPDC will attempt to partition the matrix into as many blocks as possible. If the value of nblocks is some positive number, then EKKLPDC will combine small blocks (if necessary) so that there are no more than this number of blocks in the decomposition. Calling EKKLPDC with a non-negative value of nblocks invokes code that attempts to find blocks in the constraint matrix, beginning with coupling rows (or columns), and then to decompose the problem accordingly. If the value of type is 1, then EKKLPDC will search for coupling rows, and if not, then EKKLPDC will search for coupling columns. EKKLPDC will search for coupling rows (or columns) first among the lowest numbered rows (or columns), then among the highest numbered rows (or columns), and finally among the interior rows (or columns). The identification of blocks is accomplished by entering block numbers in the row status region for each row, identified by index control variable number 5 (Nrowstat).
For staircase problems, the process of identifying matrix blocks can be greatly facilitated by identifying the coupling rows in the row status region before calling EKKLPDC. That this has been done is communicated to EKKLPDC by assigning the value -1 to the nblocks parameter, and then setting the value of integer control variable number 71, Ilpdcflag, to the maximum number of blocks into which the matrix is to be decomposed. What must be done is to mark coupling rows by a value of 0 in their row status vector entry. If block numbers can be confidently assigned to some rows, these block numbers should be put in the corresponding row status vector entries. Rows that can not be confidently assigned a block number should have -1 put in their row status vector entries. See the code fragment presented below, for solving a completely decoupled problem.
For a simple, nonseparable staircase problem, the simple code fragment
CALL EKKLPDC(IRTCOD,DSPACE,1,1,0) CALL EKKSSLV(IRTCOD,DSPACE,1,0)
For a Dantzig-Wolfe problem, the code fragment
CALL EKKLPDC(IRTCOD,DSPACE,2,2,0) CALL EKKSSLV(IRTCOD,DSPACE,1,0)
may be used to solve the problem. The value of the type parameter in the calling sequence of EKKLPDC indicates that the problem is a Dantzig-Wolfe problem. The value of the strategy parameter indicates that the right hand side of the coupling rows will be partitioned and divided among the subproblems in an attempt to achieve feasibility in those rows. The value of the nblocks parameter indicates that the problem should be decomposed into as many blocks as possible. As with staircase problems, if the value of nblocks is some positive number, then EKKLPDC will combine small blocks (if necessary) so that there are no more than this many blocks in the decomposition.
EKKLPDC may have trouble identifying problem structure, especially for staircase problems. For the user that has a good understanding of the structure of his problem, there is an unambiguous method for specifying structure. This is done by loading block numbers into the status region. For staircase problems, this may look like:
C DATA CHARR /'01', '02', '03', '04', '05', '06', '07', & '08', '09', '10', '11', '12'/ C CALL EKKDSCA(RTCOD,DSPACE,MAXSPC,1) C CALL EKKMPS(RTCOD,DSPACE,8,2,0) C CALL EKKNGET(RTCOD,DSPACE,OSLN,OSLNLN) CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) C DO 133 III=1,INUMROWS QNAME = CSPACE(NROWNAMES+III-1) QNUM = QNAME(2:3) DO 134 II=1,12 IF (QNUM.EQ.CHARR(II)) MSPACE(NROWSTAT+III-1) = II 134 CONTINUE 133 CONTINUE C CALL EKKLPDC(RTCOD,DSPACE,1,1,-1) CALL EKKSSLV(RTCOD,DSPACE,1,3) C
This code fragment assumes that the row names in the MPS file have block numbers stored in characters 2 and 3. For instance, a valid row name might be R0300056. The name indicates that row 56 belongs to block 3. The -1 in the nblocks parameter of EKKLPDC invokes a marking of the columns, and the coupling columns will be handled in the way indicated by the strategy parameter. If strategy=1, each subproblem is solved with the neighboring coupling columns free. If strategy=2, each subproblem is solved with the neighboring coupling columns fixed at lower bound. If strategy=3, EKKLPDC solves a sequence of embedded problems. The constraint system of the first problem is an end block. Successive problems are formed by adding a block to the previous constraint system and are solved using the solution to the previous problem as a warm start. So if strategy=3, EKKLPDC solves the problem, and the final call to EKKSSLV should serve only to verify optimality. Finally, if strategy=0 EKKLPDC simply marks the columns and exits.
Since the solution generated by EKKLPDC is unlikely to be dual feasible, it is probably best to chose the primal algorithm of EKKSSLV to solve the problem to optimality. If the solution is very nearly primal feasible, it might help to decrease the value of the real control variable Rpweight and/or increase the value of Rchangeweight. These control variables are equally important for Dantzig-Wolfe problems, and some results for adjusting Rpweight for Dantzig-Wolfe problems are given below.
The maximum number of iterations between factorizations in EKKSSLV is stored in the integer control variable number 3, Imaxfactor. For staircase and for Dantzig-Wolfe problems, it may be beneficial to increase Imaxfactor before calling EKKLPDC. The value of Imaxfactor can then be restored to its original value before making the final call to EKKSSLV, or it can be left unchanged. Some results are given below for a home-grown staircase problem with 12000 columns, 2400 rows and 6 periods. The default value for Imaxfactor is given at the left of the table. The best results seem to be achieved by adding 40 to the default value of Imaxfactor. (The value of Imaxfactor was reset to its default value before EKKSSLV was called).
|
Imaxfactor |
124 |
144 |
154 |
164 |
|
CPU TIME |
45 |
44 |
40 |
36 |
|
Imaxfactor |
174 |
184 |
194 |
204 |
|
CPU TIME |
40 |
37 |
40 |
43 |
Dantzig-Wolfe (or block angular) problems are handled by EKKLPDC by doing either one iteration of Dantzig-Wolfe decomposition or more than one iteration. If only one iteration is desired, the value of the type parameter should be set to 2. If the number of coupling rows is very small, this option is probably the best. EKKLPDC has code that tries to identify a block structure in the constraint matrix. For Dantzig-Wolfe problems, the coupling rows must be grouped at the beginning of the rows or at the end, otherwise the problem structure will not be detected. For the user that knows the block structure of his problem, the following code shows how to specify the structure unambiguously.
C DATA CHARR/'01','02','03','04','05','06'/ C CALL EKKDSCA(RTCOD,DSPACE,MAXSPC,1) C CALL EKKMPS(RTCOD,DSPACE,98,2,0) C CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN) CALL EKKNGET(RTCOD,DSPACE,OSLN,OSLNLN) C DO 100 I=1,INUMROWS RNAME = CSPACE(NROWNAMES+I-1) IF (RNAME(1:1).EQ.'B') THEN C This row belongs to a block. Check the name to get the C block number. CBLOCKNUM = RNAME(2:3) DO J=1,IBLOCK IF (CBLOCKNUM.EQ.CHARR(J)) MSPACE(NROWSTAT+I-1) = J ENDDO ELSE C All other rows will be assigned to block number 0 MSPACE(NROWSTAT+I-1) = 0 ENDIF 100 CONTINUE C CALL EKKLPDC(RTCOD,DSPACE,2,1,-1) CALL EKKSSLV(RTCOD,DSPACE,1,3)
The program above reads a problem with a Dantzig-Wolfe (block angular) structure from an MPS file. The names of the coupling rows begin with "C". The names of the non-coupling rows begin with 'B' and the block number for each row is found in characters 2 and 3 of the name. The name of each of the rows is examined. If the row is a coupling row, a 0 is put in the corresponding entry of the row status region. Otherwise, the block number is extracted from the name, and this number is put into the row status region. EKKLPDC is called with nblocks=-1, to indicate that a decomposition of the rows is given in the row status region. In EKKLPDC, columns are assigned to blocks using this information.
If multiple passes of Dantzig-Wolfe Decomposition are desired, EKKLPDC should be called with type set to 3. The number of iterations is then given by integer control variable number 5, Imajorits. The accuracy of the warm start generated by EKKLPDC can depend heavily on the number of iterations that are executed. Table 1 gives results from Kennington. To get the data, EKKLPDC was called with type set to 3 and strategy set to 1. EKKLPDC was then followed by the primal algorithm of EKKSSLV.
|
Iterations |
3 |
6 |
12 |
24 |
48 |
96 |
121 |
|
1 |
379.83 |
293.73 |
266.90 |
254.92 |
256.17 |
246.84 |
247.27 |
|
2 |
392.13 |
299.58 |
223.88 |
261.24 |
243.40 |
225.21 |
213.38 |
|
4 |
347.45 |
251.80 |
213.31 |
122.66 |
113.98 |
113.27 |
59.10 |
|
8 |
343.08 |
213.44 |
206.42 |
93.92 |
118.64 |
78.18 |
59.62 |
|
16 |
114.84 |
95.93 |
44.81 |
57.57 |
48.50 |
51.53 |
52.21 |
|
32 |
61.59 |
63.19 |
49.93 |
43.60 |
48.78 |
53.06 |
52.38 |
|
64 |
70.66 |
68.54 |
54.67 |
42.98 |
47.91 |
52.74 |
52.22 |
The number of iterations is given on the vertical axis. The number of blocks into which the problem is decomposed is given on the horizontal axis. The decomposition into 121 blocks is the "maximum" decomposition, selected by setting nblocks=0. The best time shown is 42.98 seconds, and this was achieved by doing 64 iterations of Dantzig-Wolfe decomposition, with the problem decomposed into 24 blocks. These results show that the overall CPU time depends on both the number of blocks in the decomposition, and on the number of iterations of Dantzig-Wolfe that are executed. Some understanding of the problem structure is probably needed to determine appropriate settings of these parameters for an unsolved problem. Some understanding of problem structure may be gained by calling EKKLPDC with strategy=0 to decompose the problem and exit.
If the solution generated by EKKLPDC is nearly optimal, overall run time can sometimes be improved by making the real control
variable Rpweight small before the call to EKKSSLV. The results below were obtained by solving OSA14 using EKKLPDC (with type set to 2, strategy set to 1, and with the number of blocks restricted to 6) and then setting Rpweight to 10n before calling EKKSSLV. The result for the default value of Rpweight (10-1) is given at the left end of Table 1.
|
Exponent(n) |
-1 |
-2 |
-3 |
-4 |
-5 |
-6 |
|
CPU TIME |
197 |
139 |
213 |
160 |
175 |
166 |
Since the overall CPU time corresponding to Rpweight=10-3 is larger than the CPU time for the default value of Rpweight, it's difficult to make a general statement about how to adjust Rpweight. It seems, however, that making Rpweight quite small (10-6) often makes an improvement. If EKKLPDC generates a primal feasible solution, then there is no need to adjust Rpweight.
The following code fragment shows how to indicate that all rows are known to be non-coupling rows belonging to unknown blocks. This information is used to identify a completely separable block structure that can be exploited by EKKLPDC. The value of the integer control variable number 71, Ilpdcflag, indicates the maximum number of subblocks into which the constraint matrix is to be decomposed. If the actual number of blocks found exceeds the value of Ilpdcflag, then blocks will be merged by EKKLPDC.
EQUIVALENCE (DSPACE(1), MSPACE(1))
CALL EKKDSCA(RTCOD,DSPACE,MAXSPC,1)
CALL EKKMPS(RTCOD,DSPACE,98,2,0)
CALL EKKNGET(RTCOD,DSPACE,OSLN,OSLNLN)
CALL EKKIGET(RTCOD,DSPACE,OSLI,OSLILN)
ILPDCFLAG = 30
DO I=1,INUMROWS
MSPACE(NROWSTAT+I-1) = -1
ENDO
CALL EKKISET(RTCOD,DSPACE,OSLI,OSLILN)
CALL EKKLPDC(RTCOD,DSPACE,2,1,-1)
CALL EKKSSLV(RTCOD,DSPACE,1,3)
Other Library subroutines can be used to enhance the performance of EKKLPDC. EKKSCAL or EKKPRSL with type set to 0 may be called before EKKLPDC. Also, EKKCRSH with type set to 3 may be used before calling EKKLPDC(RTCOD,DSPACE,3,3,n). This call to EKKLPDC uses the current solution in DSPACE to generate starting bases for the subproblems.
Since the behavior of EKKLPDC is always very problem dependent, it's difficult to make any general statements about the right way to select the calling parameters for EKKLPDC, adjust control variables, or call additional library subroutines. But the computational results show the potential impact of using the tuners, and they illustrate the importance of the flexibility of the library.
It is possible to define row and column names for a problem without having the data in MPS format. The sample FORTRAN program EXNAME shows how to do this. In brief, the tactic is as follows.
This same idea can be used to combine data from an MPS input file with that from another source. One could proceed very nearly as described above, except that the MPS data file would not be a dummy, and the offset for the second block need not be (0,0).
In this same vein, it is possible to read data from two MPS input files and to combine them into one model. It is easy to combine input from MPS files if one file represents an addendum to the other. For example, one file might contain additional constraints on the structural variables of a problem defined in another file. Ideally, these files would have the same columns corresponding to the same structural variables, and these columns would appear in the same order in the two files. Alternatively, one file might contain coefficients for new structural variables to be added to existing constraints defined in the other file. In this case, the two files would have coefficients to be included in the same rows, and the rows would be in the same order in the files. These easy alternatives allow you to add rows (or columns) to a problem without changing the number of columns (or rows) and without duplicating, eliminating, or rearranging rows or columns.
The following fragment of a FORTRAN program shows how to add to a matrix read from one MPS file the additional rows read from another. (In this example, the rows read from the first MPS file are added to those read from the second.) In principle, much more extensive melding of MPS files can be done to define a large problem. But this requires complicated searching on row and column names and extensive data rearranging. It might be easier just to write a utility program to combine multiple MPS files before calling one of the library I/O modules. Here is a fragment of a FORTRAN program to build a composite model from two MPS input files with the same columns in the same order.
C Bring in the needed INCLUDE files. INCLUDE ( OSLI ) INCLUDE ( OSLN ) C Describe work space and allow room for two models. CALL EKKDSCA( RTCOD, DSPACE, MAXSPC, 2 ) C Read model data from first MPS file (from unit 98). CALL EKKMPS( RTCOD, DSPACE, 98, 2, 0 ) C Get and save information for adding this input as the C second block. CALL EKKIGET( RTCOD, DSPACE, OSLI, OSLILN ) NROWS = Inumrows NCOLS = Inumcols NELEMS = Inumels CALL EKKNGET( RTCOD, DSPACE, OSLN, OSLNLN ) NR1 = Nblockrow NC1 = Nblockcol NE1 = Nblockelem NRL = Nrowlower NRU = Nrowupper NAMES = Nrownames C Describe model 2 as having two blocks. CALL EKKDSCM( RTCOD, DSPACE, 2, 2 ) C Allow space for the additional rows, "NROWS" of them, from C the first MPS file. (There are no additional columns.) CALL EKKIGET( RTCOD, DSPACE, OSLI, OSLILN ) Imaxrows = -NROWS CALL EKKISET( RTCOD, DSPACE, OSLI, OSLILN ) C Read in the second MPS file (from unit 97). CALL EKKMPS( RTCOD, DSPACE, 97, 2, 0 ) C Get the number of rows just read in to use as the row offset C in the call to EKKDSCB. CALL EKKIGET( RTCOD, DSPACE, OSLI, OSLILN ) C Add the new rows to the model as a block, using the C saved pointers. CALL EKKDSCB( RTCOD, DSPACE, 2, 0, MSPACE(NR1), & MSPACE(NC1), DSPACE(NE1), Inumrows, 0, NCOLS, NELEMS ) C Move the upper and lower bounds, and the row names for the C added rows. CALL EKKNGET( RTCOD, DSPACE, OSLN, OSLNLN ) Do 10 I = 0, NROWS-1 DSPACE(Nrowupper+Inumrows+I) = DSPACE(NRU+I) DSPACE(Nrowlower+Inumrows+I) = DSPACE(NRL+I) QSPACE(Nrownames+Inumrows+I) = QSPACE(NAMES+I) 10 Continue C Reset the number of rows in the model. (The number of C columns has not changed.) CALL EKKIGET( RTCOD, DSPACE, OSLI, OSLILN ) Inumrows = Inumrows + NROWS CALL EKKISET( RTCOD, DSPACE, OSLI, OSLILN ) C Make a clean copy of the new matrix. CALL EKKNWMT( RTCOD, DSPACE, 2 )
The Optimization Library includes an module, called EKKSMDL, that provides access to spreadsheet programs. EKKSMDL can read linear, quadratic and mixed integer programming problems formulated in Lotus 1-2-3 and stored in a worksheet file. The formulation of problems in Lotus 1-2-3 is easy, and it eliminates the coding needed to implement calls to EKKLMDL, EKKQMDL and EKKIMDL in an application program. An update feature of EKKSMDL puts the solution of a spreadsheet problem into the worksheet file for further processing and report generation.
The spreadsheet formulation technique uses simple expressions for the formulation of linear constraints, and also gives a simple way to represent integer data. Variables that must be 0 or 1 in an optimal solution need only be listed in spreadsheet cells. Variables that may take any integer value have a similar representation. For other types of integer problems, such as the optimization of a piecewise linear objective function, spreadsheets may be used to formulate any of the four special ordered sets described above.
The rest of this article gives an introduction to the Lotus 1-2-3 interface by giving an example problem. The pieces of an optimization problem, the variables, the constraints, and the objective function, are formulated using Lotus 1-2-3 cells, named ranges, and @Functions. In this example, cells A1, B1, C1, and D1 contain the values of the problem variables for the separable quadratic network given in the following table.
A Lotus 1-2-3 exampleminimize: 10a + a2 + 5b + 2b2 + 2c + 3c2 + d + 4d2 ,
subject to: |
a |
+ b |
= |
3 |
||
c |
+ d |
= |
4 |
|||
- a |
- c |
= |
- 2 |
|||
- b |
- d |
= |
- 5 |
|||
and |
a |
b |
c |
d |
|
|
Initial values of 0 are loaded into the four cells. The node constraints and lower bounds are loaded into cells A2 through D2 and A3 through D3. Objective function data is loaded into cells A4 through D4 and A5 through D5. With these initializations, the spreadsheet contains:
A1: 0 A2: @SUM(+A1+B1) = 3 A3: +A1=0 A4: 10 A5: 1 B1: 0 B2: @SUM(+C1+D1) = 4 B3: +B1=0 B4: 5 B5: 2 C1: 0 C2: @SUM(-A1-C1) = -2 C3: +C1=0 C4: 2 C5: 3 D1: 0 D2: @SUM(-B1-D1) = -5 D3: +D1=0 D4: 1 D5: 4
To use range names for mnemonics, the following names are created:
SOL: A1..D1, NODES: A2..D2, BNDS: A3..D3, LIN: A4..D4, QUAD: A5..D5.
Finally, the following objective function is loaded into cell A6, and the (single cell) range A6..A6 is given the name OBJ.
A6: @SUMPRODUCT(SOL,LIN) + @SUMPRODUCT(SOL,QUAD,SOL)
Before the call to EKKSMDL in an application program that solves this problem, the character control variables Cssolution, Csconstrts, and Csobjective are assigned "SOL", "NODES,BNDS" and "OBJ" respectively to indicate the location of these cells. Users of Lotus 1-2-3 for Windows may specify the cells in the Solver Definition Dialogue Box rather than coding the character control variable assignments.
Using EKKSMDL, a single application program can be used to read, solve and update a variety of problems in spreadsheet format. This simplifies the process of formulating and solving problems as well as giving spreadsheet users access to the power and function of the Optimization Library.
Footnote (4): Each integer variable is associated with at least one special ordered set of integer variables. The members of special ordered sets of types 1, 2, and 3 are known from the explicit definitions of these sets. The members of special ordered sets of type 4 (ordinary or general integer variables) are known from the definitions of these sets by use of EKKIMDL or EKKSMDL. When general integer variables that are not to be assigned to a special ordered set of type 1, 2, or 3 are encountered by EKKMPS, they are assigned to a new special ordered set of type 4.
[ Top of Page | Previous Page | Next Page | Table of Contents ]