Thomas F Rutherford
email: rutherford@colorado.edu

MILES: A Mixed Inequality and nonLinear Equation Solver

Keywords: R n 6R n, R, u 0 R n Find: z, w, v

MILES is a solver for nonlinear complementarity problems and nonlinear systems of equations. This solver can be accessed indirectly through GAMS/MPSGE or GAMS/MCP. This paper documents the solution algorithm, user options, and program output. The purpose of the paper is to provide users of GAMS/MPSGE and GAMS/MCP an overview of how the MCP solver works so that they can use the program effectively.

Introduction

MILES is a Fortran program for solving nonlinear complementarity problems and nonlinear systems of equations. The solution procedure is a generalized Newton method with a backtracking line search. This code is based on an algorithm investigated by [START_REF] Mathiesen | Computation of economic equilibria by a sequence of linear complementarity problems[END_REF] who proposed a modeling format and sequential method for solving economic equilibrium models.

The method is closely related to algorithms proposed by [START_REF] Robinson | A quadratically-convergent algorithm for general nonlinear programming problems[END_REF], Hogan (1977), [START_REF] Eaves | A locally quadratically convergent algorithm for computing stationary points[END_REF] and [START_REF] Josephy | Newton's method for generalized equations[END_REF]. In this implementation, subproblems are solved as linear complementarity problems (LCPs), using an extension of Lemke's almost-complementary pivoting scheme in which upper and lower bounds are represented implicitly. The linear solver employs the basis factorization package LUSOL, developed by [START_REF] Gill | Maintaining LU factors of a general sparse matrix[END_REF].

The class of problems for which MILES may be applied are referred to as "generalized" or "mixed" complementarity problems, which is defined as follows:

-2 -When R=-4 and u=4 MCP reduces to a nonlinear system of equations.

When R=0 and u=+4, the MCP is a nonlinear complementarity problem. Finite dimensional variational inequalities are also MCP. MCP includes inequality-constrained linear, quadratic and nonlinear programs as special cases, although for these problems standard optimization methods may be preferred. MCP models which are not optimization problems encompass a large class of interesting mathematical programs. Specific examples of MCP formulations are not provided here. See Rutherford (1992a) for MCP formulations arising in economics.

Other examples are provided by [START_REF] Harker | Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications[END_REF] and Dirkse (1993).

There are two ways in which a problem may be presented to MILES:

(i) MILES may be used to solve computable general equilibrium models generated by MPSGE as a GAMS subsystem. In the MPSGE language, a model-builder specifies classes of nonlinear functions using a specialized tabular input format embedded within a GAMS program. Using benchmark quantities and prices, MPSGE automatically calibrates function coefficients and generates nonlinear equations and Jacobian matrices. Large, complicated systems of nonlinear equations may be implemented and analyzed very easily using this interface to MILES. An introduction to general equilibrium modeling with GAMS/MPSGE is provided by Rutherford (1992a).

(ii) MILES may be accessed as a GAMS subsystem using 1 There is one other MCP solver available through GAMS: PATH (Ferris and [START_REF] Dirkse | Robust solution of mixed complementarity problems[END_REF].

-3 -variables and equations written in standard GAMS algebra and the syntax for "mixed complementarity problems" (MCP).

If more than one MCP solver is available 1 , the statement "OPTION MCP=MILES;" tells GAMS to use MILES as the MCP solution system. When problems are presented to MILES using the MCP format, the user specifies nonlinear functions using GAMS matrix algebra and the GAMS compiler automatically generates the Jacobian functions. An introduction to the GAMS/MCP modeling format is provided by [START_REF] Rutherford | Applied general equilibrium modeling using MPS/GE as a GAMS subsystem[END_REF].

The purpose of this document is to provide users of MILES with an overview of how the solver works so that they can use the program more effectively. Section 2 introduces the Newton algorithm. Section 3 describes the implementation of Lemke's algorithm which is used to solve linear subproblems. Section 4 defines switches and tolerances which may be specified using the options file. Section 5 interprets the run-time log file which is normally directed to the screen. Section 6 interprets the status file and the detailed iteration reports which may be generated. Section 7 lists and suggests remedies for abnormal termination conditions.

The Newton Algorithm

The iterative procedure applied within MILES to solve nonlinear complementarity problems is closely related to the classical Newton algorithm for nonlinear equations. This first part of this section reviews the classical procedure. A thorough introduction to these ideas is provided by [START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF]. For a practical perspective, see [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF].

Newton algorithms for nonlinear equations begin with a local (Taylor series) approximation of the system of nonlinear equations. For a point z in the neighborhood of , the system z of nonlinear functions is linearized:

LF(z) ' F( z) % LF( z)(z&
8 = 1. If , the step size 8 is adopted, 2F(z z % d k )2 # 2F(z k )2
otherwise 8 is multiplied by a positive factor ", " < 1, and the convergence test is reapplied. This procedure is repeated until either an improvement results or 8 < 8. When 8 = 0, a positive " and 8 correspond to user-specified tolerances DMPFAC and MINSTP, respectively.

-5 -

d d 2 F(z k % d k )2 < 0. F(z k )%LF(z k )d&w%v'0 R#d%z k #u, w$0, v$0 w T (d%z k &R) ' v T (u&d&z k ) ' 0 step is taken provided that: 2
Convergence theory for this algorithm is quite well developed.

See, for example, Ortega and Rheinbolt (1970) or [START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF]. The most attractive aspect of the Newton scheme with the backtracking line search is that in the neighborhood of a well-behaved fixed point, 8=1 is the optimal step length and the rate of convergence can be quadratic. If this method finds a solution, it does so very quickly.

The application of Newton methods to nonlinear complementarity problems involves a modification of the search direction. Here, d solves a linear complementarity problem (LCP) rather than a linear system of equations. For iteration k, d solves:

Conceptually, we are solving for d, but in practice MILES solves the linear problem in terms of the original variables z = z k + d: [START_REF] Kaneko | A linear complementarity problem with an n by 2n 'P'matrix[END_REF] provides some convergence theory for the linearized subproblem.

-6 -

F(z k )&LF(z k )z k % LF(z k ) z ' w & v R # z # u, w $ 0, v $ 0 w T (z&R)'0, v T (u&z)'0 After computing the solution z, MILES sets d k = z -z k .
The linear subproblem incorporates upper and lower bounds on any or all of the variables, assuring that the iterative sequence always remains within the bounds: ( R # z k # u). This can be helpful when, as is often the case, F() is undefined for some

z 0 R n .
Convergence of the Newton algorithm applied to MCP hinges on three questions: (i) Does the linearized problem always have a solution? (ii) If the linearized problem has a solution, does Lemke's algorithm find it? and (iii) Is it possible to show that the computed direction d k will provide an "improvement" in the solution? Only for a limited class of functions F() can all three questions be answered in the affirmative. For a much larger class of functions, the algorithm converges in practice but convergence is not "provable". 3 The answer to question (iii) depends on the choice of a norm by which an improvement is measured. The introduction of bounds and complementarity conditions makes the calculation of an error index is made more complicated. In MILES, the deviation associated with a candidate solution z, ,(z), is based on a

In the following x + = max(x, 0).

-7 -

w i ' F i (z) % , v i ' &F i (z) % . B i ' (z i &u i ) % % (R i &z i ) % , C i ' L i w i % U i v i .
measure of the extent to which z, w and v violate applicable upper and lower bounds and complementarity conditions.

Evaluating Convergence

Let * i L and * i U be indicator variables for whether z i is off its lower or upper bound. These are defined as: 4 L i ' min(1,(z i &R i ) % ), and U i ' min(1,(u i &z i ) % ).

Given z, MILES uses the value of F(z) to implicitly define the slack variables w and v:

There are two components to the error term associated with index i, one corresponding to z i 's violation of upper and lower bounds:

and another corresponding to violations of complementarity conditions:

Parameter p may be selected with input parameter NORM. The default value for p is +4.

-8 - (z) ' 2 B (z)% C (z)2 p Given: M 0 R nxn , q,R,u 0 R n Find: z,w,v 0 R n such that Mz % q ' w & v R # z # u, w $ 0, v $ 0 w T (z&R)'0, v T (u&z)'0
The error assigned to point z is then taken:

for a pre-specified value of p = 1, 2 or +4. 5

Lemke's Method with Implicit Bounds

A mixed linear complementarity problem has the form:

In the Newton subproblem at iteration k, the LCP data are given by q = F(z k ) -LF(z k ) z k and M = LF(z k ).

The Working Tableau

In MILES, the pivoting scheme for solving the linear problem works with a re-labeled linear system of the form:

B x B % N x N ' q where x B 0 R n , x N 0 R 2n , and the tableau [B | N] is a conformal "complementary permutation" of [-M | I | -I].
That is, every

x B i ' R i if x B i ' z i 0 if x B i ' w i or v i x B i ' u i if x B i ' z i 4 if x B i ' w i or v i
column i in B must either be the ith column of M, I or -I, while the corresponding columns i and i+n in N must be the two columns which were not selected for B.

To move from the problem defined in terms of z, w and v to the problem defined in terms of x B and x N , we assign upper and lower bounds for the x B variables as follows:

The values of the non-basic variables x i N and x i N +n are determined by the assignment of x i B :

x B i ' z i Y x N i 'w i '0 x N i%n 'v i '0 w i Y x N i 'z i 'R i x N i%n 'v i '0 v i Y x N i 'w i '0 x N i%n 'z i 'u i
In MILES, B 0 is chosen using the initially assigned values for z.

When z i # R i , then x i B = w i ; when z i $ u i , then x i B = v i ; otherwise x i B = z i .
7

The present version of the code simply sets B 0 = -I and x B = w when the user-specified basis is singular. A subsequent version of the code will incorporate the algorithm described by Anstreicher, Lee and [START_REF] Rutherford | Extensions of GAMS for variational and complementarity problems with applications in economic equilibrium analysis[END_REF] for coping with singularity.

-10 -

x B '(B 0 ) &1 (q & N x N )
In words: if z i is basic then both w i and v i equal zero. If z i is non-basic at its lower bound, then w i is possibly non-zero and v i is non-basic at zero. If z i is non-basic at its upper bound, then v i is possibly non-zero and w i is non-basic at zero. Conceptually, we could solve the LCP by evaluating 3 n linear systems of the form:

x B ' B &1 q & N x N
Lemke's pivoting algorithm provides a procedure for finding a solution by sequentially evaluating some (hopefully small) subset of these 3 n alternative linear systems.

Initialization

Let B 0 denote the initial basis matrix. 6 The initial values for basic variables are then:

If

, then the initial basis is feasible and the

x B # x B # x B complementarity problem is solved. 7 Otherwise, MILES introduces x B ' x B & h z 0
an artificial variable z 0 and an artificial column h. Basic variables are then expressed as follows:

where is the "transformed artificial column" (the h untransformed column is

). The coefficients of are h ' B 0 h h selected so that:

(i) The values of "feasible" basis variables are unaffected by

z 0 : ( ). x B i # x B i # x B i Y hi '0 (ii)
The "most infeasible" basic variable ( i=p) is driven to its upper or lower bound when z 0 = 1: hp '

x B p &x B p if x B p >x B p x B p &x B p if x B p <x B p (iii)
All other infeasible basic variables assume values between their upper and lower bounds when z 0 increases to 1:

x B i ' 1%x B i if x B i >&4, x B i '%4 x B i %x B i 2 if x B i >&4, x B i <%4 x B i &1 if x B i '&4, x B i <%4 -12 -

Exiting Variable Entering Variable

Change in Non-basic Values (i) z i at lower bound w i increases from 0 When z 0 enters the basis, it assumes a value of unity, and at this point (baring degeneracy), the subsequent pivot sequence is entirely determined. The entering variable in one iteration is determined by the exiting basic variable in the previous iteration. For example, if z i were in B 0 and introducing z 0 caused z i to move onto its lower bound, then the subsequent iteration introduces w i . Conversely, if w i were in B 0 and z 0 caused w i to fall to zero, the subsequent iteration increases z i from R i . Finally, if v i were in B 0 and z 0 's introduction caused v i to fall to zero, the subsequent iteration decreases z i from u i .

x i N = z i = R i (ii) z i at upper bound v i increases from 0 x i N +n = z i = u i (iii) w i at 0 z i increases from R i x i N = x i N +n = 0 (iv) v i at 0 z i decreases from u i x i N = x i N +n = 0
The full set of pivoting rules is displayed in Table 1. One difference between this algorithm and the original Lemke (type III) pivoting scheme (see [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF], Garcia and Zangwill If all structural variables are subject to finite upper and lower bounds, then no z i variables may be part of a homogeneous solution adjacent to a secondary ray. This does not imply, however, that secondary rays are impossible when all z i variables are bounded, as a ray may then be comprised of w i and v i variables.

-13 -(1981), or [START_REF] Cottle | The Linear Complementarity Problem[END_REF]) is that structural variables (z i 's) may enter or exit the basis at their upper bound values.

The algorithm, therefore, must distinguish between pivots in which the entering variable increases from a lower bound from those in which the entering variable decreases from an upper bound.

Another difference with the "usual" Lemke pivot procedure is that an entering structural variable may move from one bound to another. When this occurs, the subsequent pivot introduces the corresponding slack variable. For example, if z i is increased from R i to u i without driving a basic variable infeasible, then z i becomes non-basic at u i , and the subsequent pivot introduces v i . This type of pivot may be interpreted as z i entering and exiting the basis in a single step. 8 In theory it is convenient to ignore degeneracy, while in practice degeneracy is unavoidable. The present algorithm does not incorporate a lexicographic scheme to avoid cycling, but it does implement a ratio test procedure which assures that when there is more than one candidate, priority is given to the most stable pivot. The admissible set of pivots is determined on both an absolute pivot tolerance (ZTOLPV) and a relative pivot tolerance (ZTOLRP). No pivot with absolute value smaller than min(ZTOLPV, ZTOLRP 5V5) is considered, where 5V5 is the norm of the incoming column.

Termination on a Secondary Ray Lemke's algorithm terminates normally when the introduction of a new variable drives z 0 to zero. This is the desired outcome, but it does not always happen. The algorithm may be interrupted prematurely when no basic variable "blocks" an incoming variable, a condition known as "termination on a secondary ray". In anticipation of such outcomes, MILES maintains a record of the value of z 0 for successive iterations, and it saves basis information associated with the smallest observed value, z 0 * . (In Lemke's algorithm, the pivot sequence is determined without regard for the effect on z 0 , and the value of the artificial variable may follow an erratic (non-monotone) path from its initial value of one to the final value of zero.)

When MILES encounters a secondary ray, a restart procedure is invoked in which the set of basic variables associated with z 0 * are reinstalled. This basis (augmented with one column from the non-basic triplet to replace z 0 ) serves as B 0 , and the algorithm is restarted. In some cases this procedure permits the pivot sequence to continue smoothly to a solution, while in others cases may only lead to another secondary ray.

When invoking MILES from within GAMS it is possible to use one of several option file names. See the README documentation with GAMS 2.25 for details.

-15 -

The Options File

MILES accepts the same format options file regardless of how the system is being accessed, through GAMS/MPSGE or GAMS/MCP. The options file is a standard text file which is normally named MILES.OPT. 9 The following is a typical options file: is the density at which the Markowitz strategy should search only 1 column or (preferably) use a dense LU for all the remaining rows and columns.

Log File Output

The log file is intended for display on the screen in order to permit monitoring progress. Relatively little output is generated.

A sample iteration log is displayed in ). The next line reports the convergence tolerance.

The lines beginning 0 and 1 are the major iteration reports for those iterations. the number following the iteration number is the current deviation, and the third number is the Armijo step length. The name of the variable complementary to the equation with the largest associated deviation is reported in parenthesis at the end of the line.

Following the final iteration is a summary of iterations, refactorizations, amd final deviation. The final message reports the solution status. In this case, the model has been successfully processed ("Solved.").

Status File Output

The status file reports more details regarding the solution process than are provided in the log file. Typically, this file is written to disk and examined only if a problem arises. Within GAMS, the status file appears in the listing only following the GAMS statement "OPTION SYSOUT=ON;".

The level of output to the status file is determined by the options passed to the solver. In the default configuration, the status file receives all information written to the log file together a detailed listing of all switches and tolerances and a report of basis factorization statistics.

When output levels are increased from their default values using the options file, the status file can receive considerably more output to assist in debugging. Tables 3-6 present a status file generated with LEVOUT=3 (maximum), PIVLOG=T, and LCPECH=T.

The status file begins with the same header as the log file.

Thereafter is a complete echo-print of the user-supplied option file when one is provided. Following the core allocation report is a full echo-print of control parameters, switches and tolerance as specified for the current run. This report has two sections: $ROWS and $COLUMNS. The four columns of numbers in the $ROWS section are the constant vector (q), the current estimate of level values for the associated variables (z), and the lower and upper bounds vectors ( R and u).

The letters L and U which appear between the ROW and Z columns are used to identify variables which are non-basic at their lower -25 -and upper bounds, respectively. In this example, all upper bounds equal +4, so no variables are non-basic at their upper bound.

By convention, only variable (and not equation names) appear in the status file. An equation is identified by the corresponding variable. We therefore see in the $COLUMNS:

section of the matrix echo-print, the row names correspond to the names of z variables. The names assigned to variables z i , w i and v i are z-<name i>, w-<name i> and v-<name i>, as shown in the $COLUMNS section. The nonzeros for w-<> and v-<> variables are not shown because they are assumed to be -/+I.

The status file output continues on Table 5 where the first half of the table reports output from the matrix scaling procedure, and the second half reports the messages associated with initiation of Lemke's procedure.

The "lu6chk warning" is a LUSOL report. Thereafter are two factorization reports. Two factorizations are undertaken here because the first basis was singular, so the program install all the lower bound slacks in place of the matrix defined by the initial values, z.

Following the second factorization report, at the bottom of Table 5 is a summary of initial pivot. "Infeasible in 3 rows."

indicates that contains 3 nonzero elements. "Maximum h infeasibility" reports the largest amount by which a structural variable violates an upper or lower bound. "Artificial column with 3 elements." indicates that the vector contains 3 h ' B 0 h elements (note that in this case B 0 = -I because the initial basis was singular, hence the equivalence between the number of nonzeros in and h.). h Table 6 displays The convergence report for iteration 1 is no different from the report written to the log file, and following this is a second report of variable and function values. We see here that a solution has been obtained following a single subproblem. This is because the underlying problem is, in fact, linear.

Table 3

Status File with Debugging Output (page 1 of 4) -2.25000000E-01 0.00000000E+00 0.00000000E+00 1.00000000E+20 X_01_02 -1.53000004E-01 0.00000000E+00 0.00000000E+00 1.00000000E+20 X_01_03 -1.61999996E-01 0.00000000E+00 0.00000000E+00 1.00000000E+20 X_02_01 -2.25000000E-01 0.00000000E+00 0.00000000E+00 1.00000000E+20 X_02_02 -1.61999996E-01 0.00000000E+00 0.00000000E+00 1.00000000E+20 X_02_03 -1.25999998E-01 0.00000000E+00 0.00000000E+00 1.00000000E+20 W_01 -3.25000000E+02 0.00000000E+00 0.00000000E+00 1.00000000E+00 W_02 -5.75000000E+02 0.00000000E+00 0.00000000E+00 1.00000000E+00 P_01 3.25000000E+02 1.00000000E+00 0.00000000E+00 1.00000000E+20 P_02 3.00000000E+02 1.00000000E+00 0.00000000E+00 1.00000000E+20 P_03 2.75000000E+02 1.00000000E+00 0.00000000E+00 1.00000000E+20 ... 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 9 to replace column 20 ---Pivot element:

Iteration 0 values. ROW Z F -------- ------------ ------------ X_01_01 L 0.
$COLUMNS: Z-X_01_01 W_01 -1.00000000E+00 P_01 1.00000000E+00 Z-X_01_02 W_01 -1.00000000E+00 P_02 1.00000000E+00 Z-X_01_03 W_01 -1.00000000E+00 P_03 1.00000000E+00 Z-X_02_01 W_02 -1.00000000E+00 P_01 1.00000000E+00 Z-X_02_02 W_02 -1.00000000E+00 P_02 1.00000000E+00 Z-X_02_03 W_02 -1.00000000E+00 P_03 1.00000000E+00 Z-W_01 X_01_01 1.00000000E+00 X_01_02 1.00000000E+00 X_01_03 1.00000000E+00 Z-W_02 X_02_01 1.00000000E+00 X_02_02 1.00000000E+00 X_02_03 1.00000000E+00 Z-P_01 X_01_01 -1.00000000E+00 X_02_01 -1.00000000E+00 Z-P_02 X_01_02 -1.00000000E+00 X_02_02 -1.00000000E+00 Z-P_03 X_01_03 -1.00000000E+00 X_02_03 -1.00000000E+00 ... ... 0.00000000E+00
-3.250E+02 Within GAMS, insert the line "OPTION SYSOUT=ON;" prior to the solve statement and resubmit the program in order to pass the MILES solver status file through to the listing.

3.25E+02 1.00E+00 1 1.14E-13 1.00E+00 (W_02 ) =============================================================== Iteration 1 values. ROW Z F -------- ------------ ------------ X_01_01 2
-33 -

Termination Messages

Basis factorization error in INVERT.

An unexpected error code returned by LUSOL. This should normally not occur. Examine the status file for a message from LUSOL. 11 Failure to converge.

Two successive iterates are identical -the Newton search direction is not defined. This should normally not occur.

Inconsistent parameters ZTOLZ0, ZTOLZE.

ZTOLZ0 determines the smallest value loaded into the cover vector h, whereas ZTOLZE is the feasibility tolerance employed in the Harris pivot selection procedure. If ZTOLZ0 < -ZTOLZE, Lemke's algorithm cannot be executed because the initial basis is infeasible.

Insufficient space for linearization.

Available memory is inadequate for holding the nonzeros in the Jacobian. More memory needs to be allocated. On a PC, you probably will need to install more physical memory -if there is insufficient space for the Jacobi matrix, there is far too little memory for holding the LU factors of the same matrix.

Insufficient space to invert.

More memory needs to be allocated for basis factors.

  z) Solving the linear system LF(z) = 0 provides the Newton direction from which given by . z d ' &LF &1 F( z) Newton iteration k begins at point z k . First, the linear model formed at z k is solved to determine the associated "Newton direction", d k . Second, a line search in direction d k determines the scalar steplength 8 and the subsequent iterate: z k+1 = z k + 8 d k . An Armijo or "back-tracking" line search initially considers

  entries are of the form "<keyword> = <value>", where keywords have at most 6 characters. The following are recognized keywords which may appear in the options file, identified by keyword, type and default value, and grouped according to function:Termination control CONTOL r Default = 1.0E-6is the convergence tolerance. Whenever an iterate is encountered for which ,(z) < CONTOL, the algorithm terminates. This corresponds to the GAMS/MINOS parameter "Row tolerance". bound on the number of Newton iterations. This corresponds to the GAMS/MINOS parameter "Major iterations". bound on the cumulative number of Lemke iterations. When MILES is invoked from within a GAMS program (either with MPSGE or GAMS/MCP), <model>.ITERLIM can be used to set this value. This corresponds to the GAMS/MINOS parameter "Iterations limit".NORM r Default = 3 defines the vector norm to be used for evaluating ,(z).Acceptable values are 1, 2 or 3 which correspond to p = 1bound on the number of factorizations in any LCP. This avoids wasting lots of CPU if a subproblem proves difficult to solve. bound on the number of restarts which the linear subproblem solver will undertake before giving up. the maximum number of Lemke iterations between A number of tolerances are set on the basis of the machine precision, EPS. For Lahey Fortran running on 80486 processor, EPS = 2.2D-16.-17 -recalculation of the basis factors. This corresponds to the GAMS/MINOS parameter "Factorization frequency". with which the factorization is checked. The number refers to the number of basis replacements operations between refinements. This corresponds to the GAMS/MINOS parameter "Check frequency".Pivot SelectionPLINFYr Default = 1.D20 is the value assigned for "plus infinity" ("+INF" in GAMS notation).ZTOLPVr Default = 3.644E-11 (EPS**(2./3.)) 10 is the absolute pivot tolerance. This corresponds, roughly, to the GAMS/MINOS parameter "Pivot tolerance" as it applies for nonlinear problems. ZTOLRP r Default = 3.644E-11 (EPS**(2./3.)) is the relative pivot tolerance. This corresponds, roughly, to the GAMS/MINOS parameter "Pivot tolerance" as it applies for nonlinear problems. the subproblem solution to determine when any -18 -variable has exceeded an upper or lower bound. This corresponds to GAMS/MINOS parameter "Feasibility tolerance". Linearization and Data Control SCALE l Default = .TRUE. invokes row and column scaling of the LCP tableau in every iteration. This corresponds, roughly, to the GAMS/MINOS switch "scale all variables". ZTOLDA r Default = 1.483E-08 (EPS**(1/2)) sets a lower bound on the magnitude of nonzeros recognized in the linearization. All coefficients with absolute value less than ZTOLDA are ignored. to generate a printout of the LCP data after of debug output written to the log and status files. The lowest meaningful value is -1 and the highest is 3. This corresponds, roughly, to the GAMS/MINOS parameter "Print level" PIVLOG l Default =.FALSE. is a switch to generate a status file listing of the Lemke pivot sequence.LUSOL parameters (as with MINOS 5.4, except LUSIZE) estimate the number of LU nonzeros which will be stored, as a multiple of the number of nonzeros in the Jacobian matrix. the pivot row and column and the no. of rows and columns involved at each the maximum number of columns searched allowed in a Markowitz-type search for the next pivot element. For some of the factorization, the number of rows searched is maxrow = maxcol -1. which limits waste space in U. In lu1fac, the row or column lists are compressed if their length exceeds uspace times the length of either file after the last compression.

  the final section of the status file. At the top of page 6 is the Lemke iteration log. The columns are interpreted as follows: ITER is the iteration index beginning with 0, STATUS is a statistic representing the efficiency of the Lemke path. Formally, status is the ratio of the minimum number of pivots from B 0 to the current basis divided by the actual number of pivots. When the status is 1, Lemke's algorithm is performing virtually as efficiently as a direct factorization (apart from the overhead of basis factor updates.) Z% indicates the percentage of columns in the basis are "structural" (z's). Z0 indicates the value of the artificial variable. Notice that in this example, the artificial variable declines monotonically from its initial value of unity. ERROR is a column in which the factorization error is reported, when it is computed. For this run, ITCH=30 and hence no factorization errors are computed. INFEAS. is a column in which the magnitude of the infeasibility introduced by the artificial column (defined using the box-norm) is reported. (In MILES the cover vector h contains many different nonzero values, not just 1's; so there may be a large difference between the magnitude of the artificial variable and the magnitude of the induced infeasibility. PIVOTS reports the pivot magnitude in both absolute terms (the first number) and relative terms (the second number). The relative pivot size is the ratio of the pivot element to the norm of the incoming column. IN/OUT report the indices (not names) of the incoming and outgoing columns for every iteration. Notice that Lemke's iteration log concludes with variable z 0 exiting the basis.

Table 1

 1 Pivot Sequence Rules for Lemke's Algorithm with Implicit Bounds

	Pivoting Rules

Table 2

 2

	Initial deviation ........ 3.250E+02	P_01
	Convergence tolerance .... 1.000E-06
	0	3.25E+02	1.00E+00 (P_01	)
	1	1.14E-13	1.00E+00 (W_02	)
	Major iterations ........	1
	Lemke pivots ............	10
	Refactorizations ........	2
	Deviation ............... 1.137E-13
	Solved.
	Work space allocated	--	0.01 Mb
	Initial deviation ........ 5.750E+02	W_02
	Convergence tolerance .... 1.000E-06
	0	5.75E+02	1.00E+00 (W_02	)
	1	2.51E+01	1.00E+00 (P_01	)
	2	4.53E+00	1.00E+00 (P_01	)
	3	1.16E+00	1.00E+00 (P_01	)
	4	3.05E-01	1.00E+00 (P_01	)
	5	8.08E-02	1.00E+00 (P_01	)
	6	2.14E-02	1.00E+00 (P_01	)
	7	5.68E-03	1.00E+00 (P_01	)
	8	1.51E-03	1.00E+00 (P_01	)
	9	4.00E-04	1.00E+00 (P_01	)
	10	1.06E-04	1.00E+00 (P_01	)
	11	2.82E-05	1.00E+00 (P_01	)
	12	7.47E-06	1.00E+00 (P_01	)
	13	1.98E-06	1.00E+00 (P_01	)
	14	5.26E-07	1.00E+00 (P_01	)
	Major iterations ........	14
	Lemke pivots ............	23
	Refactorizations ........	15
	Deviation ............... 5.262E-07
	Solved.

. This output is from two cases solved in succession. This and subsequent output comes from program TRNSP.FOR which calls the MILES library directly. (When MILES is invoked from within GAMS, at most one case is processed at a time.)

Table 2

 2

	Sample Iteration Log

Table 4

 4

contiues the status file. The iteration-byiteration report of variable and function values is produced whenever LEVOUT $ 2. Table

4

also contains an LCP echo-print.

Table 4

 4 Status File with Debugging Output (page 2 of 4)

	-31 -

Table 5

 5 Status File with Debugging Output (page 3 of 4)

					LEMKE PIVOT STEPS
					=================
	ITER STATUS	Z%	Z0	ERROR	INFEAS. ----PIVOTS ----IN	OUT
	1	1.00	0 1.000		3.E+02 1.E+00 1 Z0	W	9
	2	1.00	9 1.000		1.E+00 1.E+00 2 Z	9 W	1
	3	1.00 18 0.997		9.E-01 9.E-01 1 Z	1 W	10
	4	1.00 27 0.997		1.E+00 1.E+00 1 Z	10 W	2
	5	1.00 36 0.996		9.E-01 4.E-01 1 Z	2 W	11
	6	1.00 45 0.996		1.E+00 1.E+00 1 Z	11 W	6
	7	1.00 55 0.479		2.E+00 1.E+00 1 Z	6 W	7
	8	1.00 64 0.479		1.E+00 1.E+00 1 Z	7 W	4
	9	1.00 73 0.000		1.E+00 1.E+00 2 Z	4 W	8
	10	1.00 73 0.000		1.E-03 2.E-03 1 V	8 Z0
	================================
	Convergence Report, Iteration 1
	===============================================================
	ITER	DEVIATION		STEP
	0
						-32 -

  FTRAN ...............: 0.0 sec. Update ..............: 0.1 sec.

		.50000E+01	-8.32667E-17
	X_01_02	3.00000E+02	-5.55112E-17
	X_01_03 L	0.00000E+00	3.60000E-02
	X_02_01	3.00000E+02	-8.32667E-17
	X_02_02 L	0.00000E+00	8.99999E-03
	X_02_03	2.75000E+02	2.77556E-17
	W_01	1.00000E+00	-1.13687E-13
	W_02	1.00000E+00	1.13687E-13
	P_01	1.22500E+00	0.00000E+00
	P_02	1.15300E+00	0.00000E+00
	P_03	1.12600E+00	0.00000E+00
	Major iterations ........	1
	Lemke pivots ............	10
	Refactorizations ........	2
	Deviation ............... 1.137E-13
	Solved.
	Total solution time .:	0.6 sec.
	Function & Jacobian..:	0.2 sec.
	LCP solution ........:	0.2 sec.
	Refactorizations ....:	0.1 sec.

Table 6

 6 Status File with Debugging Output (page 4 of 4)

Table 7

 7 Transport Model in GAMS/MCP (page 1 of 2)

The status file (for this case) concludes with an iteration summary identical to the log file report and a summary of how much CPU time was employed overall and within various subtasks.

(Don't be alarmed if the sum of the last five numbers does not add up to the first number -some cycles are not monitored precisely.)

Increase the value of LUSIZE in the options file, or assign a larger value to <model>.workspace if MILES is accessed through GAMS.

Iteration limit exceeded.

This can result from either exceeding the major (Newton) or minor (Lemke) iterations limit. When MILES is invoked from GAMS, the Lemke iteration limit can be set with the statement "<model>.iterlim = xx;" (the default value is 1000). The Newton iteration limit is 25 by default, and it can be modified only through the ITLIMT option.

Resource interrupt.

Elapsed CPU time exceeds options parameter RESLIM. To increase this limit, either add RESLIM = xxx in the options file or (if MILES is invoked from within GAMS), add a GAMS statement "<model>.RESLIM = xxx;".

Singular matrix encountered.

Lemke's algorithm has been interrupted due to a singularity arising in the basis factorization, either during a column replacement or during a refactorization. For some reason, a restart is not possible.

Termination on a secondary ray.

Lemke's algorithm terminated on a secondary ray. For some reason, a restart is not possible.

Unknown termination status.

The termination status flag has not been set, but the code POSITIVE VARIABLES W(I) shadow price at supply node i, P(J) shadow price at demand node j, X(I,J) shipment quantities in cases; EQUATIONS SUPPLY(I) supply limit at plant i, FXDEMAND(J) fixed demand at market j, PRDEMAND(J) price-responsive demand at market j, PROFIT(I,J) zero profit conditions;

PROFIT(I,J).. REPORT("EQUIL",I,J) = X.L(I,J); REPORT("EQUIL","Price",J) = P.L(J); REPORT("EQUIL",I,"Price") = W.L(I);

DISPLAY "BENCHMARK CALIBRATION", REPORT; * Compute a counter-factual equilibrium:

C("SEATTLE","CHICAGO") = 0.5 * C("SEATTLE","CHICAGO"); SOLVE FIXEDQTY USING MCP; REPORT("FIXED",I,J) = X.L(I,J); REPORT("FIXED","Price",J) = P.L(J); REPORT("FIXED",I,"Price") = W.L(I); * Replicate the fixed demand equilibrium: SOLVE EQUILQTY USING MCP; REPORT("EQUIL",I,J) = X.L(I,J); REPORT("EQUIL","Price",J) = P.L(J); REPORT("EQUIL",I,"Price") = W.L(I);

DISPLAY "Reduced Seattle-Chicago transport cost:", REPORT;