
\chapter{Environments}\label{chap:environments}

% This chapter is an update to sections of the following papers:
%   ``Complementarity Problems in GAMS and the PATH Solver''
%       - GAMS models and documentation
%   ``Interfaces to PATH 3.0: Design, Implementation, and Usage''
%       - Basic AMPL, GAMS, NEOS, and MATLAB interfaces
% Rewrote the sections on AMPL, NEOS, and MATLAB to add examples and
% update the interfaces.

Environments are used by practitioners to easily specify large, complex models.
Modeling languages \cite{bisschop.meeraus:development,fourer.gay.ea:modeling},
such as AMPL \cite{fourer.gay.ea:ampl} and GAMS \cite{brooke.kendrick.ea:gams},
provide natural facilities for representing complementarity relationships
\cite{dirkse.ferris.ea:gams,ferris.fourer.ea:expressing}.  Additional 
environments, such as MATLAB \cite{matlab:users} and NEOS 
\cite{czyzyk.mesnier.ea:network-enabled,ferris.mesnier.ea:neos}, are also
available for expressing complementarity problems.

% In particular, mechanisms for solving complementarity problems using the 
% MATLAB and NEOS tools are now available.

This chapter looks at complementarity from a practitioner's perspective.  We 
first develop two example models illustrating complementarity relationships in
\secref{applications}.  These examples are then used in the subsequent survey 
of the AMPL, GAMS, MATLAB, and NEOS environments for expressing complementarity
problems and communicating them to solvers.  The discussion of each 
environment includes information on the input format, output generated, and 
unique features of the environment.  

\section{Example Applications}\label{sec:applications}

%%% Add references as needed!!!

%%% For the book, we will extend the applications section to cover
%%% three different case studies:
%%%    1.  An economics section as written
%%%    2.  An engineering section
%%%    3.  Financial section -- options pricing

%%% Probably want to add details and basic proofs to the text...
%%% (Don't want the book to be all words with no substance.)

% This section is an update to sections of the following papers:
%   ``Complementarity Problems in GAMS and the PATH Solver''
%       - Transportation model development
%       - Walrasian equilibrium
%   ``Case Studies in Complementarity: Improving Model Formulation''
%       - Von Thunen Case Study
% Added the linear programming example (from GAMS/PATH User Guide)

Problems fitting into the complementarity framework occur in a wide variety 
of disciplines.  In this section, we introduce the nonlinear complementarity 
problem from an economics perspective using the optimality conditions for a 
transportation model.  The more general 
mixed complementarity problem is motivated using a Walrasian equilibrium
\cite{arrow.debreu:existence} model.

\subsection{Transportation Model}\label{sec:transportation}
The transportation model is a linear program where demand for a
single good must be satisfied by suppliers at a minimal transportation 
cost.  The underlying transportation network is given as a set 
$\arcs \subseteq \mathcal{S}\times\mathcal{D}$ of arcs from suppliers,
$\mathcal{S}$, to demand centers, $\mathcal{D}$.  The problems
variables are the amounts $x_{i,j}$ to be shipped over each arc
$(i,j)\in \arcs$.  A linear program to find an optimal shipment
schedule can be written mathematically as 
\begin{equation}
\label{eqn:transport}
\begin{array}{rc}
\min_{x \geq 0} & \sum_{(i,j) \in \arcs} c_{i,j} x_{i,j} \\
\text{subject to} & \sum_{j: (i,j) \in \arcs} x_{i,j} \leq s_i, \;
\forall i \in \mathcal{S} \\
& \sum_{i: (i,j) \in \arcs} x_{i,j} \geq d_j, \;
\forall j \in \mathcal{D}
\end{array}
\end{equation}
where $c_{i,j}$ is the unit shipment cost on the arc $(i,j)$ 
connecting supplier $i$ to demand center $j$, $s_i$ is the amount of
supply available at $i$, and $d_j$ is the demand at $j$.

The derivation of a complementarity problem for the transportation
model begins by associating with each constraint a multiplier,
alternatively termed a dual variable or shadow price.  These
multipliers represent the marginal price on changes to the 
corresponding constraint.  Labeling the prices on the supply
constraint $p^s$ and those on the demand constraint $p^d$, 
we then intuitively have at each supply node $i$
\[
0 \leq p^s_i, \;
s_i \geq \sum_{j: (i,j) \in \arcs} x_{i,j}.
\]
Consider the case when $s_i > \sum_{j: (i,j) \in \arcs} x_{i,j}$, that
is there is excess supply at $i$.  Then, in a competitive marketplace,
no rational person is willing to pay for more supply at node $i$; it
is already over-supplied.  Therefore, $p^s_i = 0$.  Alternatively,
when $s_i = \sum_{j: (i,j) \in \arcs} x_{i,j}$, that is node $i$
clears, we might be willing to pay for additional supply of the good.
Therefore, $p^s_i \geq 0$.  We write these two conditions succinctly
as:
\[
\begin{array}{rcll}
0 \leq p^s_i & \perp &
s_i \geq \sum_{j: (i,j) \in \arcs} x_{i,j}, & \forall i
\end{array}
\]
where the $\perp$ notation is understood to mean that at least one of 
the adjacent inequalities must be satisfied as an equality.  For 
example, either $0 = p^s_i$, the first case, or 
$s_i = \sum_{j: (i,j) \in \arcs} x_{i,j}$, the second case.  This
notation is used to represent the complementarity relationship and means
that $p^s_i$ and $s_i - \sum_{j: (i,j) \in \arcs} x_{i,j}$ are orthogonal.

Similarly, at each node $j$, the demand must be satisfied in any feasible 
solution, that is
\[
\sum_{i: (i,j) \in \arcs} x_{i,j} \geq d_j.
\]
Furthermore, the model assumes all prices are nonnegative, $0 \leq p^d_j$.
If there is too much of the commodity supplied,
$\sum_{i: (i,j) \in \arcs} x_{i,j} > d_j$, then, in a 
competitive marketplace, the price $p^d_j$ will be driven down to $0$.  
Summing these relationships gives the following complementarity condition:
\[
\begin{array}{rcll}
0 \leq p^d_j & \perp &
\sum_{i: (i,j) \in \arcs} x_{i,j} \geq d_j, & \forall j .
\end{array}
\]

The supply price at $i$ plus the transportation cost $c_{i,j}$ from
$i$ to $j$ must exceed the market price at $j$.  That is
$p^s_i + c_{i,j} \geq p^d_j$.  Otherwise, competitors would enter the
market and increase the available supply, driving down the market price.
This chain would repeat until the inequality is
satisfied.  Furthermore, if the cost of delivery strictly exceeds the
market price, that is $p^s_i + c_{i,j} > p^d_j$, then nothing is
shipped from $i$ to $j$ because doing so would incur a loss.  In this
case, $x_{i,j} = 0$.  Therefore, 
\[
\begin{array}{rcll}
0 \leq x_{i,j} & \perp &
p^s_i + c_{i,j} \geq p^d_j, & \forall (i,j) \in \arcs.
\end{array}
\]

Collecting all of the above relationships, we have the linear complementarity
problem
\begin{equation}
\label{eqn:transmcp}
\begin{array}{rcll}
0 \leq p^s_i & \perp &
s_i \geq \sum_{j: (i,j) \in \arcs} x_{i,j}, & \forall i \\
0 \leq p^d_i & \perp &
\sum_{i: (i,j) \in \arcs} x_{i,j} \geq d_j, & \forall j \\
0 \leq x_{i,j} & \perp &
p^s_i + c_{i,j} \geq p^d_j, & \forall (i,j) \in \arcs
\end{array}
\end{equation}
which is easily recognized as the complementary slackness conditions for the 
linear program \eqnref{transport}.  For linear programs the complementary 
slackness conditions are both necessary and sufficient for $x$ to be an 
optimal solution to \eqnref{transport}.  Furthermore, the conditions 
\eqnref{transmcp} are also the necessary and sufficient optimality conditions 
for a related problem in the variables $(p^s,p^d)$ 
\[
\begin{array}{rc}
\max_{p^s,p^d \geq 0} & \sum_{j \in \mathcal{D}} d_j p^d_j - \sum_{i
\in \mathcal{S}} s_i p^s_i \\
\text{subject to} & c_{i,j} \geq p^d_j - p^s_i, \; \forall (i,j) \in
\arcs
\end{array}
\]
termed the dual linear program, hence the nomenclature ``dual
variables''. 

Looking at these conditions a bit more closely we can gain further
insight into complementarity problems. A solution of \eqnref{transmcp}
tells the arcs used to transport goods. A priori we do not need to
specify the arcs to use, the solution itself indicates them.  This
property represents the key contribution of a complementarity problem
over a system of equations. If we know what arcs to send flow down, we
can just solve a simple system of linear equations.  However, the key
to the modeling power of complementarity is that it chooses which of
the inequalities in \eqnref{transmcp} to satisfy as equations.  In 
economics we can use this property to generate a model with different regimes 
and let the solution determine which ones are active.  Using the 
transportation model, we could specify two alternate routes from a supplier to
a demand center: one using a four-lane interstate highway and the other a 
two-lane state highway.  A regime shift could occur when the state highway 
becomes active because of congestion on the interstate highway, which slows 
traffic down and consequently drives up the transportation cost.  Congestion 
could be modeled, for example by making the cost of traveling along a road a 
function of the usage:
\[
\begin{array}{rcll}
0 \leq x_{i,j} & \perp &
p^s_i + c_{i,j}(x) \geq p^d_j, & \forall (i,j) \in \arcs.
\end{array}
\]
for some cost function, $c_{i,j}(\cdot)$.

% A regime shift could, for example, be a back stop technology like 
% windmills that become profitable (active) if a $CO_2$ tax is increased. 

While many interior point methods for linear programming exploit this
complementarity framework (so-called primal-dual methods
\cite{wright:primal-dual}), the real power of this modeling format is
the new problem instances it enables a modeler to create. We now show
some examples of how to extend the simple model \eqnref{transmcp} to
investigate other issues and facets of the problem.  

Demand in \eqnref{transmcp} is independent of the prices $p$.  Since
the prices $p$ are variables in the complementarity problem
\eqnref{transmcp}, we can easily replace the constant demand $d$ with
a function $d(p)$.  This function can be used to more accurately model
demand which typically is higher at a low price than at a high price.
Clearly, any function of $p$ can be added to the model.  For 
example, a linear demand function could be expressed using 
\[
\sum_{i: (i,j) \in \arcs} x_{i,j} \geq d_j(1 - p^d_j),\; \forall j.
\]
Note that the demand is rather strange if $p^d_j$ exceeds 1. Other
more reasonable examples for $d(p)$ can be derived.
For example, if we assume the demand for the commodity is isoelastic
\cite{rutherford:extensions}, then the demand for the commodity is
\[
     d_j(p) \define \frac{d_j}{(p^d_j)^{\alpha_j}}
\]
for fixed $d_j$ and $\alpha_j > 0$. Note that the resulting complementarity 
problem becomes nonlinear in the variables $p$ and that the function is
undefined when $p^d_j = 0$. 

Another feature that can be added to this model are tariffs or taxes.
In the case where a tax is applied at the supply point, the third
general inequality in \eqnref{transmcp} is replaced by 
\[
p^s_i ( 1 + t_i ) + c_{i,j} \geq p^d_j,\; \forall (i,j) \in \arcs.
\]
%The taxes can be made endogenous to the model by adding an additional
%complementarity condition.  For example, we could add the constraint 
%\[
%\begin{array}{rcll}
%0 \leq t_i & \perp & t_i p^s_i \sum_j x_{i,j} \geq T_i & \forall i
%\end{array}
%\]
%to require that the tax revenue will be at least $T_i$.  
Details about complementarity problems for more general transportation 
models can be found in 
\cite{dirkse.ferris:traffic,ferris.meeraus.ea:computing}. 

% Note that the matrix for the above extensions is copositive.  (The sum
% of a positive semidefinite matrix and a positive matrix is at least
% copositive.  (x^T(M + P)x = x^TMx + x^TPx \geq 0) for all $x \geq 0$.
%
% Can we show existence of a solution?
%
% What about a tax on the demand side of the equation?  To get copositivity
% we would need the constraint:
% \[
%   p^s_i + c_{i,j} \geq p^d_j ( 1 - t_j ),\; \forall (i,j) \in \arcs
% \]
% which seems reasonable, as the consumer pays p^d, but the producer
% only gets (1 - t)p^d with the rest going to taxes.

\textit{An important observation to make is that with either of these
modifications, the resulting complementarity problem is no longer the 
optimality conditions for a linear program.  In many cases, there is no 
optimization problem corresponding to the complementarity conditions.}

We now abstract from the particular example to describe more carefully
the complementarity problem in its mathematical form.  All the above
examples can be cast as nonlinear complementarity problems defined as
follows:
\begin{description}
\item[(Nonlinear Complementarity Problem \cite{cottle:nonlinear})]
  Let $F : \real^{n} \rightarrow \real^{n}$ be
  a continuously differentiable function.  The nonlinear complementarity
  problem is to find $z \in \real^n$ such that
  \[
	0 \leq z\; \perp\; F(z) \geq 0.
  \]
\end{description}
Recall that the $\perp$ sign signifies that one of the inequalities is
satisfied as an equality, so that componentwise $z_iF_i(z) = 0$.  We
frequently refer to this property as $z_i$ is ``complementary'' to
$F_i$.  A special case of the nonlinear complementarity problem that
has received much attention is when $F$ is a linear function, the
linear complementarity problem \cite{cottle.pang.ea:linear}.

\subsection{Walrasian Equilibrium}\label{sec:walrasian}
A Walrasian equilibrium \cite{arrow.debreu:existence} can also be 
formulated as 
a complementarity problem.  In this case, we want to find a  
price $p \in \real^m$ and an activity level $y \in \real^n$ such that
\begin{equation}
  \begin{array}{rclcl}
    0 \leq y & \perp & L(p) &\define & -A^T p \geq 0 \\
    0 \leq p & \perp & S(p,y) &\define & b + Ay - d(p) \geq 0. 
  \end{array}
  \label{eqn:walras}
\end{equation}
Here, $S(p,y)$ represents the excess supply function, and $L(p)$
represents the loss function.  The first complementarity
relationship chooses the activities $y_i$ to run (i.e. only those that do not
make a loss).  The second set of inequalities state that the price of
a commodity can only be positive if there is no excess supply.  These
conditions indeed correspond to the standard exposition of Walras'
law \cite{walras:elements} which states that supply equals demand if we 
assume all
prices $p$ will be positive at a solution. 
Formulations of equilibria as systems of equations do not allow the
model to choose the activities present, but typically make an a
priori assumption on this matter.

% Note we have copositivity of each of the jacobian matrices if
% d_p(p) is a decreasing function (for example).  

Many large scale models of this nature have been developed.  An
interested modeler could, for example, see how a large scale
complementarity problem was used to quantify the effects of the
Uruguay round of talks \cite{harrison.rutherford.ea:quantifying*2}.

In many modeling situations, a key tool for clarification is the use
of intermediate variables.  As an example, the modeler may wish to
define a variable corresponding to the demand function $d(p)$ in the
Walrasian equilibrium \eqnref{walras}.  Using the variables $d_p$ to 
store the value of the demand function referred to in the excess
supply equation, we have the following complementarity problem:
\begin{equation}
  \begin{array}{rclcl}
    0 \leq y & \perp & L(p) & \define & -A^T p \geq 0 \\
    0 \leq p & \perp & S(p,y) & \define & b + Ay - d_p \geq 0 \\
	     &	     &  & & d_p = d(p).
  \end{array}
  \label{eqn:walras2}
\end{equation}
The model now contains a mixture of equations and complementarity
constraints.  Since constructs like the above are prevalent in many
practical models, a more general definition of complementarity is needed.

A mixed complementarity problem is specified by three pieces of data,
namely the lower bounds $\ell$, the upper bounds $u$, and a function $F$.
\begin{description}
\item[(Mixed Complementarity Problem)]
Given lower and upper bounds
\begin{eqnarray*}
\ell \in \{ \real\cup\{-\infty\}\}^n\\
u \in \{ \real\cup\{+\infty\}\}^n
\end{eqnarray*}
and a continuously differentiable function $F:\real^n \rightarrow \real^n$,  
the mixed complementarity problem is to find $z \in \real^n$ such that 
\textit{precisely} one of the following holds for each 
$i \in \{1, \ldots, n\}$: 
  \[
\begin{array}{rcl}
F_i(z) = 0 & \text{ and } & \ell_i \leq z_i \leq u_i \\
F_i(z) > 0 & \text{ and } & z_i = \ell_i \\
F_i(z) < 0 & \text{ and } & z_i = u_i .
\end{array}
  \] 
\end{description}
These relationships define the general complementarity problem needed
for \eqnref{walras2}.  We will use the notation
\[
\begin{array}{lcr}
\ell \leq z \leq u & \perp & F(z)
\end{array}
\]
to denote the mixed complementarity relationship. Both the nonlinear
complementarity problem found in \eqnref{transmcp} and \eqnref{walras}
and nonlinear systems of equations are special cases of mixed
complementarity problems. Nonlinear complementarity problems have
$\ell = \{0\}^n$ and $u = \{\infty\}^n$ and nonlinear equations have $\ell =
\{-\infty\}^n$ and $u = \{\infty\}^n$.  A point misunderstood by many
experienced modelers is that \textit{the variable bounds determine the
relationships satisfied by the function} $F$.  Therefore, modifications 
to the lower and upper bounds on the variables can drastically change the
problem under consideration. 

An advantage of the mixed complementarity formulation described above
is the pairing between ``fixed'' variables (ones with equal upper and
lower bounds) and a component of $F$.  If a variable $z_i$ is fixed, 
then $F_i(z)$ is unrestricted since precisely one of the three
conditions in the MCP definition automatically holds when $z_i =
\ell_i = u_i$.  Thus if a variable is fixed, the paired equation can be
completely dropped from the model.  This convenient trick can
be used to remove particular constraints from a model.

In many cases, equations (which are complementary to free 
variables) can be written in any convenient way since they are internally 
``substituted out'' of the model.  In particular, for defining equations, 
such as those presented in \eqnref{walras2}, the choice appears to 
be arbitrary.  However, underlying model monotonicity is important.  For
example, if we have the linear program
\[
\begin{array}{rc}
\min_{x} &  c^T x \\
\mbox{subject to} & Ax = b,\; x \geq 0
\end{array}
\]
we can write the first order optimality conditions as either
\[
  \begin{array}{rcl}
0 \leq x & \perp & -A^T \mu + c \geq 0\\
\mu \mbox{ free} & \perp & Ax - b 
  \end{array}
\]
or, equivalently,
\[
  \begin{array}{rcl}
0 \leq x & \perp & -A^T \mu + c \geq 0\\
\mu \mbox{ free} & \perp & b - Ax 
  \end{array}
\]
because we have an equation.  The former is a linear (mixed) complementarity 
problem with a positive semidefinite matrix, while the latter is indefinite.  
If we need to add a diagonal perturbation to the problem for numerical 
reasons, the former system becomes positive definite, while the later 
remains indefinite and unlikely to be solvable. 

Modelers typically add bounds to their variables when
attempting to solve nonlinear problems in order to restrict the domain
of interest.  For example, many square nonlinear systems are
formulated as
\[ 
\begin{array}{c}
F(z) = 0 \\
\ell \leq z \leq u
\end{array}
\]
where typically, the bounds on $z$ are inactive at the solution.
This is \textit{not} a mixed complementarity problem, but is an
example of a ``constrained nonlinear system.''
It is important to note the distinction between complementarity
problems and constrained systems.  The complementarity problem uses
the bounds to infer relationships on the function $F$.  If a finite
bound is active at a solution, the corresponding component of $F$ is only
constrained to be nonnegative or nonpositive in the complementarity
problem, whereas in the constrained system $F$ must be zero.  Thus
there may be many solutions of the complementarity problem that 
do not satisfy $F(z) = 0$.  Only if $z^*$ is a solution of the
complementarity problem with $\ell < z^* < u$ is it guaranteed that
$F(z^*) = 0$.

Simple bounds on the variables are a convenient modeling tool that
translates into efficient mathematical programming tools.  For
example, specialized codes exist for the bound constrained
optimization problem
\[
\begin{array}{rc}
\min & f(x) \\
\text{subject to} & \ell \leq x \leq u.
\end{array}
\]
The first order optimality conditions of this problem are precisely
MCP($\nabla f(x),\ell,u$).  We can easily reason about this
condition in the one dimensional setting.  If we are at an
unconstrained minimizer, then $\nabla f(x) = 0$.  Otherwise, 
if $x$ is a local minimizer at its lower bound, then the function must be 
increasing as $x$ increases, so $\nabla f(x) \geq 0$.  Alternatively, if $x$ 
is a local minimizer at its upper
bound, then the function must be increasing as $x$ decreases, so that
$\nabla f(x) \leq 0$.  Thus,
\[
\begin{array}{lcr}
	\ell \leq x \leq u & \perp & \nabla f(x).
\end{array}
\]
The complementarity framework allows such models to be easily and 
efficiently processed.

Upper bounds can also be used to extend the utility of existing models.
For example, in the Walrasian equilibrium it may be necessary to
have an upper bound on the activity level $y$.  In this case, we
simply add an upper bound to $y$.  The model then becomes:
\begin{equation}
  \begin{array}{rclcl}
    0 \leq y \leq 10 & \perp & L(p) &\define & -A^T p \\
    0 \leq p & \perp & S(p,y) &\define & b + Ay - d(p)
  \end{array}
  \label{eqn:walras3}
\end{equation}
where we use the $\perp$ notation in the generalized sense of the
mixed complementarity problem.  \textit{Recall that the bounds on the
variables completely determine the relationship satisfied by the
equations.}  We can interpret the relationships that the above change 
generates as follows.  If $y_j = 0$, the loss function can be positive 
since we are not producing in the $j$th sector.  If $y_j$ is strictly between 
its bounds, then the loss function must be zero by complementarity; this is 
the competitive assumption.  However, if $y_j$ is at its upper bound, then 
the loss function can be negative. Of course, if the market does not 
allow free entry, some firms may operate at a profit (negative loss).

\begin{comment}
\subsection{Von Thunen Land Model}\label{sec:von}

The Von Thunen Land Model \cite{stevens:location} is cast in the
Arrow-Debreu framework \cite{arrow.debreu:existence} as an
equilibrium problem.  The basic model is a closed economy consisting
of three economic agents, a landowner, a worker, and a porter.  There
is a central market, around which concentric regions of land are
located.  Since the produced goods have to be delivered to the market,
this is an example of a spatial price equilibrium.  The key variables
of the model are the prices of commodities, land, labor, and
transport.  Given these prices, it is assumed that the agents demand
certain amounts of the commodities, which are supplied so as to
maximize profit in each sector.  Walras' law is then a consequence of
the assumed competitive paradigm, namely that supply will equal demand
in the equilibrium state.

We now describe the problems that the consumers and the producers
face.  We first look at consumption and derive a demand function for
each of the consumer agents in the economy.  Each of these agents has
a utility function that they wish to maximize subject to their
budgetary constraints.  As is typical in such problems, the utility
function, $u_a(d)$, is assumed to be Cobb-Douglas
\[
\begin{array}{rcl}
u_a(d) & \define & \prod_{c} d_c^{\alpha_{c,a}} \\
& & \sum_{c} \alpha_{c,a} = 1 \\
& & \alpha_{c,a} \geq 0
\end{array}
\]
where the $\alpha_{c,a}$ are given parameters dependent only on the
agent.  For each agent $a$, the variables $d_{c}$ represent quantities
of the desired commodities $c$.  In the Von Thunen model, the goods
are wheat, rice, corn and barley.  The agents endowments determine
their budgetary constraint as follows.  Given current market prices,
an agents wealth is the value of the initial endowment of goods at
those prices.  The agents problem is therefore
\[
\begin{array}{rc}
\max_{d \geq 0} & u_a(d) \\
\text{subject to} & \sum_c p_cd_c \leq \sum_c p_ce_{c,a}
\end{array}
\]
where $e_{c,a}$ is the endowment for agent $a$ of commodity $c$.
A closed form solution, corresponding to demand from agent $a$ 
for commodity $c$ is thus
\[
     d_{c,a}(p) \define \frac{\alpha_{c,a} \sum_cp_ce_{c,a}}{p_c}.
\]
Note that this assumes the prices of the commodities $p_c$ are positive.

The supply side of the economy is similar.  The worker earns a wage
$w_L$ for labor input. The land is distributed around the market
in rings with a rental rate $w_r$ associated with each ring $r$ of
land.  The area of land $a_r$ in each ring is an increasing function
of $r$. The model assumes that labor and land are substitutable via a
constant elasticities of substitution (CES) function.  

Consider the production $x_{c,r}$ of commodity $c$ in region $r$. In
order to maximize profit (or minimize costs), the labor $y_L$ and
land use $y_r$ solve 
\begin{equation}
\begin{array}{rc}
\min_{y_L,y_r \geq 0} & w_L y_L + w_r y_r \\
\text{subject to} & \phi_c y_L^{\beta_c} y_r^{1 - \beta_c} \geq x_{c,r}
\end{array}
\label{prodProblem}
\end{equation}
where $\phi_c$ is a given cost function scale parameter, and $\beta_c
\in [0,1]$ is the share parameter.
The technology constraint is precisely the CES function allowing a
suitable mix of labor and land use.
Again, a closed form solution can be calculated.  For example, the
demand for labor in order to produce $x_{c,r}$ of commodity $c$ 
in region $r$ is given by
\[
     x_{c,r} \frac{\beta_c 
       \left(\frac{w_L}{\beta_c}\right)^{\beta_c}
       \left(\frac{w_r}{1-\beta_c}\right)^{1-\beta_c}}{\phi_c w_L} .
\]
Considering all such demands, this clearly 
assumes the prices of inputs $w_L$, $w_r$ are positive.  A key point to
note is that input commodity (factor) demands to produce $x_{c,r}$ can be
determined by first solving (\ref{prodProblem}) for unit demand
$x_{c,r} \equiv 1$ and then multiplying these factor demands by the
actual amount desired. 
Let $\bar{y}_L$ and $\bar{y}_r$ denote the optimal solutions of
(\ref{prodProblem}) with $x_{c,r} \equiv 1$.
Using this fact, the {\em unit} production cost $\gamma_{c,r}$ 
for commodity $c$ in region $r$ can be calculated as follows:
\begin{eqnarray*}
\gamma_{c,r} & = & w_L \bar{y}_L + w_r \bar{y}_r \\
& = & w_L \frac{\beta_c 
       \left(\frac{w_L}{\beta_c}\right)^{\beta_c}
       \left(\frac{w_r}{1-\beta_c}\right)^{1-\beta_c}}{\phi_c w_L}
+ w_r \frac{(1-\beta_c) 
       \left(\frac{w_L}{\beta_c}\right)^{\beta_c}
       \left(\frac{w_r}{1-\beta_c}\right)^{1-\beta_c}}{\phi_c w_r} \\
& = & \frac{1}{\phi_c}\left(\frac{w_L}{\beta_c}\right)^{\beta_c}
                \left(\frac{w_r}{1-\beta_c}\right)^{1-\beta_c} .
\end{eqnarray*}

Transportation is provided by a porter, earning a wage $w_p$.  If we
denote the unit cost for transportation of commodity $c$ by $t_c$,
then unit transportation cost to market is
\[
T_{c,r}(w_p) \define t_c d_r w_p ,
\]
where $d_r$ is the distance of region $r$ to the market.
Spatial price equilibrium arises from the consideration:
\[ 0 \leq x_{c,r} \perp \gamma_{c,r}(w_L,w_r) + T_{c,r}(w_p)
  \geq p_c .
\]
This is intuitively clear; it states that commodity 
$c$ will be produced in region $r$ only if
the combined cost of production and transportation equals
the market price.

The above derivations assumed that the producers and consumers act
as price takers.
Walras' law is now invoked to determine the prices so that markets
clear.
The resulting complementarity problem is:
\begin{eqnarray}
& & \gamma_{c,r} =
                \frac{1}{\phi_c}\left(\frac{w_L}{\beta_c}\right)^{\beta_c}
                \left(\frac{w_r}{1-\beta_c}\right)^{1-\beta_c} 
\label{costdef} \\
0 \leq x_{c,r} 
& \perp & \gamma_{c,r} + T_{c,r}(w_p) \geq p_c 
\label{rational} \\
0 \leq w_L 
& \perp & e_L  \geq \sum_{r,c} x_{c,r}
\frac{\beta_c \gamma_{c,r} }{ w_L} 
\label{laborbnd} \\
0 \leq w_r 
& \perp & a_r \geq
                \sum_c \frac{x_{c,r} (1-\beta_c) \gamma_{c,r}}{ w_r } 
\label{areabnd} \\
0 \leq w_p 
& \perp & e_P \geq \sum_{r,c} t_c d_r x_{c,r} 
\label{transbnd} \\
0 \leq p_c 
& \perp & \sum_r x_{c,r} \geq
                   \frac{\alpha_{c,P} e_P w_p
                +  \alpha_{c,L} e_L w_L
                +  \alpha_{c,O} \sum_r w_r a_r}{p_c}
\label{mktclear}
\end{eqnarray}
Note that in (\ref{laborbnd}), (\ref{areabnd}) and  (\ref{transbnd}),
the amounts of labor, land and transport are
bounded from above, and hence the prices on these inputs are
determined as multipliers (or shadow prices) on the corresponding
constraints. 
The final relationship (\ref{mktclear})
in the above complementarity problem
corresponds to market clearance;
prices are nonnegative and can only be positive if supply equals
demand.
(Some modelers multiply the last inequality throughout by $p_c$.  This
removes problems where $p_c$ becomes zero, but can also
introduce spurious solutions.)

The Arrow-Debreu theory guarantees that the problem is homogeneous in
prices; $(x,\lambda w,\lambda p)$ is also a solution whenever
$(x,w,p)$ solves the above.
Typically this singularity in the model is removed by
fixing a numeraire, that is fixing a price (for example $w_L = 1$)
and dropping the
corresponding complementary relationship.

Unfortunately, in this formulation even after fixing a numeraire,
some of the variables $p$ and $w$
may go to zero, resulting in an ill-defined problem.  In the case of the
Von Thunen land model, the rental price of land $w_r$ 
decreases as the distance
to market increases, and for remote rings of land, it becomes zero.  A
standard modeling fix is to put artificial lower bounds on these variables.
Even with this fix, the problem typically remains very hard to solve.
More importantly, the homogeneity property of the prices used above to
fix a numeraire no longer holds,
and the corresponding complementary relationship (which was dropped
from the problem) may fail to be
satisfied.  It therefore matters which numeriare is fixed, and many
modelers run into difficulty since in many cases 
the solution found by a solver is
invalid for the originally posed model.

\subsubsection{Intervention Pricing}
The principal difficulty is the fact that the
rental prices on land go to zero as proximity to the market
decreases, and become zero for sufficiently remote rings.  Such a
property is unlikely to hold in a practical setting.
Typically, a landowner has a minimum rental price (for example, land in
fallow increases in value).  As outlined above, a fixed lower bound on
the rental price violates the well-established homogeneity property.
One suggestion from Tom Rutherford is to allow the
landowner to intervene and ``purchase-back'' land whenever the
rental cost gets smaller than a certain fraction of the labor wage. 

The new model adds a (homogeneous in price) constraint
\[
\begin{array}{lcr}
0 \leq i_r 
& \perp & w_r  \geq 0.0001 * w_L 
\end{array}
\]
and modifies (\ref{areabnd}) and (\ref{mktclear}) as follows:
\begin{eqnarray}
0 \leq w_r 
& \perp & a_r - i_r \geq
                \sum_c \frac{x_{c,r} (1-\beta_c) \gamma_{c,r}}{ w_r } 
\nonumber
\\
0 \leq p_c 
& \perp & \sum_r x_{c,r} \geq
                   \frac{\alpha_{c,P} e_P w_p
                +  \alpha_{c,L} e_L w_L
                +  \alpha_{c,O} \sum_r w_r (a_r - i_r)}{p_c} .
\label{mod-mclear}
\end{eqnarray}

Given the intervention purchase, we can now add a lower bound on $w_r$
to avoid division by zero errors.

\subsubsection{Nested Production and Maintenance}
Another observation that can be used to overcome the land price going
to zero is the fact that land typically requires some maintenance
labor input to keep it usable for crop growth.  Traditionally, in
economics, this is carried out by providing a nested CES function as
technology input to the model.  The idea is that 
commodity $c$ in region $r$ is made from labor and an intermediate
good.  The intermediate good is ``maintained land''.  Essentially, the
following production problem replaces (\ref{prodProblem}):
\[
\begin{array}{rc}
    \min_{y_M,y_L,y_r,g} & w_L (y_M + y_L) + w_r y_r \\
    \text{ subject to } 
    & y_r \geq (1 - \beta_c - \epsilon) g \\
    & y_M \geq \epsilon g \\
    & \phi_c y_L^{\beta_c} g^{1 - \beta_c} \geq 1,\\
    & y_M, y_L, y_r, g \geq 0 .
\end{array}
\]
Note that the variable $y_M$ represents ``maintenance labor'' and $g$
represents the amount of ``maintained land'' produced, an intermediate
good.  The process of generating maintained land uses a Leontieff
production function, namely
\[
       \min ( \lambda_r y_r, \lambda_M y_M ) \geq g .
\]
Here $\lambda_M = \frac{1}{\epsilon}$, $\epsilon$ small, corresponds
to small amounts of maintenance labor, while $\lambda_r =
\frac{1}{1-\beta_c-\epsilon}$ is chosen to calibrate the model
correctly.  A simple calculus exercise then generates appropriate
demand and cost expressions.  The resulting complementarity problem
comprises (\ref{rational}), (\ref{transbnd}), (\ref{mktclear}) and
\begin{eqnarray*}
& & \gamma_{c,r} =
                \frac{w_L^{\beta_c}}{\phi_c}
                \left(\frac{w_L\epsilon + w_r(1-\beta_c-\epsilon)}{1-\beta_c}\right)^{1-\beta_c} 
\\
0 \leq w_L 
& \perp & e_L  \geq \sum_{r,c} x_{c,r} \gamma_{c,r} 
\left(\frac{\beta_c}{w_L} + \frac{\epsilon (1 - \beta_c)}{w_L\epsilon + w_r(1-\beta_c-\epsilon)}\right)
\\
0 \leq w_r 
& \perp & a_r \geq
                \sum_c \frac{x_{c,r} \gamma_{c,r} (1-\beta_c)
                (1-\beta_c - \epsilon) }{w_L\epsilon +
                w_r(1-\beta_c-\epsilon)} 
\end{eqnarray*}
\end{comment}

\begin{comment}
\subsection{Variational Inequalities}\label{sec:variational}

Recall that a mixed complementarity problem is specified by three pieces 
of data, namely the lower bounds $\ell$, the upper bounds $u$, and a 
function $F$.
\begin{description}
\item[(MCP)]
Given lower bounds $\ell \in \{ \real\cup\{-\infty\}\}^n$, upper bounds
$u \in \{ \real\cup\{\infty\}\}^n$ and a function 
$F:\real^n \rightarrow \real^n$, find $z \in \real^n$ such that 
\textit{precisely} one of the following holds for 
each $i \in \{1, \ldots, n\}$: 
\[
   \begin{array}{rcl}
   F_i(z) = 0 & \text{ and } & \ell_i \leq z_i \leq u_i \\
   F_i(z) > 0 & \text{ and } & z_i = \ell_i \\
   F_i(z) < 0 & \text{ and } & z_i = u_i.
   \end{array}
\] 
\end{description}
Unfortunately, this representation of a complementarity problem is
not always conducive to analysis.

One alternative formulation typically found in the literature is as a
variational inequality defined by the function $F$ and a closed convex 
set $C$.
\begin{description}
\item[(VI)]
Given a closed convex set $C \subseteq \real^n$ and a function 
$F:\real^n \rightarrow \real^n$, find $z \in C$ such that
\[
    \ip{F(z)}{y - z} \geq 0 \; \mbox{ for all } y \in C
\]
\end{description}
Using this definition, we can see that MCP($F,\ell,u$) is equivalent to 
solving VI($F,[\ell,u]$) where, in this case, $C = [\ell,u]$ is the Cartesian 
product of closed (not necessarily bounded) intervals.  The variational
inequality can also be represented as a generalized equation
\cite{robinson:generalized}.  Find a $z \in C$ such that
\begin{eqnarray*}
0 \in F(z) + N_C(z)
\end{eqnarray*}
where by $N_C(z)$ the normal cone \cite{rockafellar:convex} to a closed 
convex set $C$ at the point $z$ defined as
\[
N_C(z) :=
\left\{
\begin{array}{cl} 
\{y \mid \ip{\bar{z} - z}{y} \le 0,\; \forall \bar{z} \in C\} &
        \mbox{if } z \in C \\
\emptyset & 
        \mbox{otherwise,}
\end{array}\right.
\]
These formulations will be used in parts of this thesis.

When $C = \{ z \mid Az \geq b \}$, we have a polyhedrally constrained
variational inequality.  By adding multipliers, we can write down a
box constrained problem:
\[
    0 \in \left[\begin{array}{l} F(x) - A^T\mu \\ Ax - b \end{array}\right]
          + \left[\begin{array}{l} N_[\ell,u](x) \\ N_{\real_+}(\mu)\end{array}\right]
\]
\end{comment}

\section{AMPL}\label{sec:ampl}

%%% Need work on slacks and squareness.

% Input format
% Non-matching constraints

AMPL \cite{fourer.gay.ea:ampl} is a general purpose language for
modeling mathematical programs including linear, nonlinear, and
integer programs.  The syntax for associating variables and functions
in complementarity relationships uses the \textbf{complements} keyword
\cite{ferris.fourer.ea:expressing} in constraint declarations.
An implementation of the transportation model from \secref{transportation} 
can be found in \figref{transmcp-ampl}; the actual data for the model is 
assumed to be given in the file \texttt{transmcp.dat}.  Example data is 
found in \figref{transmcp-ampl2}.

\begin{figure}
\begin{verbatim}
set S;                     # Suppliers
set D;                     # Markets
                 
param cap {S};             # Capacity of supplier
param dem {D};             # Demand at market
param tra {S,D};           # Transportation cost
                 
var x{S,D};                # Shipment quantity
var p_s{S};                # Price at plant
var p_d{D};                # Price at market

subj to 
supply {i in S}:           # Supply limit
    0 <= p_s[i] complements cap[i] >= (sum {j in D} x[i,j]);

subj to 
demand {j in D}:           # Demand satisfaction
    0 <= p_d[j] complements (sum {i in S} x[i,j]) >= dem[j];

subj to 
rational {i in S, j in D}: # Rational marketplace
    0 <= x[i,j] complements p_s[i] + tra[i,j] >= p_d[j];
\end{verbatim}
\caption{Transportation model in AMPL, \texttt{transmcp.mod}}
\label{fig:transmcp-ampl} 
\end{figure}

\begin{figure}
\begin{verbatim}
set S := seattle, san-diego;
set D := new-york, chicago, topeka;

param:     cap :=
seattle    325
san-diego  575;

param:     dem :=
new-york   325
chicago    300
topeka     275;

param tra: new-york  chicago   topeka :=
seattle    0.225     0.153     0.162
san-diego  0.225     0.162     0.126;
\end{verbatim}
\caption{Transportation data in AMPL, \texttt{transmcp.dat}}
\label{fig:transmcp-ampl2} 
\end{figure}

Each condition declared using the \textbf{complements} keyword must
specify exactly two inequalities.  In \figref{transmcp-ampl}, these 
correspond to one inequality on each variable and function
which precisely define the complementarity condition.  Any ordering of the 
inequalities can be used in the pairing.  For example, an alternative 
declaration would be to reorder the \texttt{supply} constraint as 
\single\begin{verbatim}
     subj to 
     supply {i in S}:           # Supply limit
         cap[i] >= (sum {j in D} x[i,j]) complements 0 <= p_s[i];
\end{verbatim}
\double
Recall that mixed complementarity conditions were defined by lower and upper 
bounds on the variables and an unrestricted function.  AMPL can be used to 
write such a representation by placing both inequalities on the variables.
\single\begin{verbatim}
     subj to supply {i in S}:
         0 <= p_s[i] <= Infinity complements 
             cap[i] - (sum {j in D} x[i,j]);
\end{verbatim}
\double
This formulation is particularly useful if a user wants to modify the bounds 
during an iterative process.  If we declare parameters, \texttt{p\_l} and 
\texttt{p\_u}, for the lower and upper bounds on the prices, then we can 
write the complementarity constraint as:
\single\begin{verbatim}
     subj to supply {i in S}:
         p_l[i] <= p_s[i] <= p_u[i] complements 
             cap[i] - (sum {j in D} x[i,j]);
\end{verbatim}
\double
The meaning of the constraint is completely determined by the bound
parameters which can be modified at run time.  Note that the function must be 
oriented correctly to preserve the meaning of the complementarity condition.  
Thus, the following declaration is incorrect:
\single\begin{verbatim}
     subj to supply {i in S}:
         p_l[i] <= p_s[i] <= p_u[i] complements 
             (sum {j in D} x[i,j]) - cap[i];
\end{verbatim}
\double

% Options
% Solver specification

To solve the \texttt{transmcp} problem with AMPL, the following
statements are issued from the AMPL prompt:
\single\begin{verbatim}
     ampl: model transmcp.mod;
     ampl: data transmcp.dat;
     ampl: solve;
\end{verbatim}
\double
The first two statements load the model and problem data respectively,
while the last passes an internal representation of the problem
to a solver.  The solver used can be changed with the statement:
\single\begin{verbatim}
     ampl: option solver path.ampl;
\end{verbatim}
\double
where \texttt{path.ampl} can be the name of any executable for solving 
complementarity problems using AMPL input.  Here, \texttt{path.ampl} refers
to the AMPL version of the PATH 4.x code.  To add solver specific 
options, we write the statement:
\single\begin{verbatim}
     ampl: option path_options "logfile foo.log";
\end{verbatim}
\double
The options list and name of the options string, \texttt{path\_options}, is 
specific to the algorithm.  For example, the SEMI code is called
using: 
\single\begin{verbatim}
     ampl: option solver semi.ampl;
     ampl: option semi_options "logfile foo.log";
\end{verbatim}
\double
PATH 4.x and SEMI use a common set 
of options for AMPL which are detailed in \tabref{options-ampl}.  
Note that an \texttt{optfile} can be set which contains a listing of 
algorithm specific options and the corresponding values to be set.
%These options are documented in \partref{software}.

\begin{table}
\centering
\caption{AMPL Specific Options}\label{tab:options-ampl}
\begin{tabular}{lll}
Option & Type & Explanation \\
\hline
logfile & string & Specify file where log will be written\\
optfile & string & Specify file containing solver options \\
statusfile & string & Specify file where detailed status information will be written \\
version & boolean & Controls printing of version information \\
sideineq & integer & Control the handling of side constraints \\
\end{tabular}
\end{table}

% Return codes and what they mean
% Multipliers

After the solve statement is encountered, AMPL will generate the model
and pass it to the algorithm for solution.  Assuming that PATH 4.x is used, 
the following output for the \texttt{transmcp} model will be displayed on 
exit:
\single\begin{verbatim}
     Path v4.4a: Solution found.
     15 iterations (1 for crash); 18 pivots.
     16 function, 16 gradient evaluations.
\end{verbatim}
\double
The message can be accessed from the \texttt{solve\_message} variable
within AMPL, and indicates the solution status, ``Solution
found'', and number of iterations taken.  The total number of
function and gradient evaluations performed are also displayed.
Furthermore, a result code is passed back to AMPL that can be checked
using the \texttt{solve\_result\_num} variable.  A complete listing of
the result codes and corresponding messages returned by PATH 4.x and SEMI
are listed in \tabref{codes-ampl}.  The result codes are
broken into the categories found in \tabref{intervals-ampl}
where the category string for the return code can be accessed with the 
\texttt{solve\_result} variable within AMPL.

\begin{table}
\centering
\caption{AMPL Return Codes}\label{tab:codes-ampl}
\begin{tabular}{ll}
Code & String \\
\hline
0 & Solution found \\
\hline
200 & Infeasible \\
201 & Inconsistent bounds \\
\hline
400 & Major iteration limit \\
401 & Cumulative iteration limit \\
402 & Time limit \\
403 & User interrupt \\
\hline
500 & Unexpected return code \\
501 & Solver error \\
503 & Domain error \\
504 & Non-square system \\
505 & Model contains side constraints \\
510 & Not enough progress \\
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{AMPL Return Code Meanings}\label{tab:intervals-ampl}
\begin{tabular}{lll}
Interval & String & Meaning\\
\hline
0  -- 99 & solved & The model has been solved \\
100--199 & solved? & A local solution has been found \\
200--299 & infeasible & The model is has no solution\\
300--399 & unbounded & The objective value is unbounded \\
400--499 & limit & A resource limit was reached by the code\\
500--599 & failure & The algorithm failed to find a solution \\
\end{tabular}
\end{table}

Once a solution has been found, we are returned to the AMPL prompt.
We can then investigate the values for the variables using the normal AMPL
syntax.   For example, to display the variable values for
transportation model, we can issue the command:
\single\begin{verbatim}
     ampl: display p_s, p_d, x;
     :             p_s       p_d      :=
     chicago        .      86.9571
     new-york       .      87.0291
     san-diego   86.8041      .
     seattle     86.8041      .
     topeka         .      86.9301
     ;

     x :=
     san-diego chicago      0
     san-diego new-york   300
     san-diego topeka     275
     seattle   chicago    300
     seattle   new-york    25
     seattle   topeka       0
     ;
\end{verbatim}
\double

%%% More careful description of the slacks needs to be present

% Associated with each variable and function is a slack value, that
% can be checked using the \textbf{slack} suffix.  For example, we can 
% display the values for the slack with the statement:
% \begin{verbatim}
%      ampl: display x.slack, supply.slack;
% \end{verbatim}

In cases where the algorithm fails to find a solution, it might be
useful to specify a \texttt{logfile} for the problem and investigate the
contents to gain information about the reasons for failure.
% Descriptions of the logfiles generated by the algorithms
% are presented in \partref{software}. 

The Walrasian equilibrium in \secref{walrasian} can be modeled
similarly as presented in \figref{walrasian-ampl}, where we have used
an equation to define the demand function, \texttt{d\_p}.
Equations count as both $\leq$ and $\geq$ inequalities.  Thus, the 
corresponding variable is unrestricted.  Note that the demand function
is undefined when the prices are zero.  If a starting point is not
supplied, AMPL uses the zero vector projected onto the simple variable
bounds.  Therefore, we must initialize the prices to positive
values, as done in the variable declaration of \figref{walrasian-ampl} in
order to start at a well-defined point.  Algorithms will typically
generate an error if a user supplies a starting point  
(either directly or indirectly) where the function is undefined.
Furthermore, \textit{a good starting point can help an algorithm to
rapidly find a solution, especially for nonlinear models.}
 
\begin{figure}
\begin{verbatim}
set I;
set J;

param A {I,J};
param b {I};
param c {I};

var p {I} := 1;   # Price of commodity
var y {J};        # Activity level
var d_p {I};      # Demand

subj to
S {i in I}:       # Excess supply function
    0 <= p[i] complements 
        b[i] + (sum {j in J} A[i,j]*y[j]) - d_p[i] >= 0;

subj to
L {j in J}:       # Loss function
    0 <= y[j] complements - (sum {i in I} p[i]*A[i,j]) >= 0;

subj to
D {i in I}:       # Demand function
    d_p[i]    complements 
        d_p[i] = c[i]*(sum {k in I} b[k]*p[k]) / p[i];
\end{verbatim}
\caption{Walrasian equilibrium as an MCP in AMPL, \texttt{walrasian.mod}}
\label{fig:walrasian-ampl}
\end{figure}

\begin{comment}
\begin{figure}
\begin{verbatim}
set I = '1', '2', '3';
set J = '1', '2';

param:     b      c :=
'1'        3      5
'2'        4      4
'3'        2      6;

param A:   '1'    '2' :=
'1'         7      -2
'2'         2      -6
'3'         1      -10;
\end{verbatim}
\caption{Walrasian equilibrium data in AMPL, \texttt{walrasian.dat}}
\label{fig:walrasian-ampl2}
\end{figure}
\end{comment}

In many cases, it is not significant to match free variables
explicitly to equations.  Therefore, when we have free variables
matched to equations, we \textit{can} remove the \textbf{complements}
from the declaration.  For example, we can write the demand function
in \figref{walrasian-ampl} as:
\single\begin{verbatim}
     subj to
     D {i in I}:                # Demand function
         d_p(i) = c[i]*(sum {k in I} b[k]*p[k]) / p[i];
\end{verbatim}
\double
The solver interface will match any unmatched, unrestricted variables
to these equality constraints.  The number of free variables and
unmatched equations must be the same in order to have a square
complementarity problem.  Note that this extension allows existing models 
consisting of a square system of nonlinear equations to be easily recast 
as a complementarity problem - the model definition is unchanged.

\textit{However, equations should be explicitly matched to free
variables in order to preserve problem structure, even though not
explicitly required.  Matrix properties, such as positive
semidefiniteness, could be lost with the arbitrary matching performed
by the solver interface, leading to poor performance (or even failure) by an 
algorithm.}  For example, when constructing the first order optimality 
conditions for a linear program, the matching of equations to multipliers 
should be done to preserve problem structure.

Upper bounds can be imposed on the activity levels in the Walrasian
equilibrium by declaring the loss function as
\single\begin{verbatim}
     subj to
     L {j in J};                # Loss function
         0 <= y(j) <= 10 complements -(sum {i in I} p[i]*A[i,j]);
\end{verbatim}
\double
The bounds on the variables completely determine the complementarity
problem.  Hence, for bounded variables, we do not know beforehand if the
corresponding function will be satisfied as an equation, less than 
inequality or greater than inequality, since this determination depends 
on the values of the solution variables.  Therefore, the function is
unrestricted.

% Pitfalls - presolve
%          - unmatched ``side constraints''

\subsection{Generalized Declaration}
The notation used in AMPL allows generalized complementarity conditions 
to be written between two functions.  For example we can declare the
following: 
\single\begin{verbatim}
     subj to 
     func {i in I}:
         0 <= (sum {j in J} A[i,j]*x[j]) - a[i] complements 
              (sum {j in I} B[i,j]*x[j]) - b[i] >= 0;
\end{verbatim}
\double
Internal to the solver interface, the condition is reformulated
into a standard mixed complementarity problem by adding variables and  
constraints.  One equivalent definition is:
\single\begin{verbatim}
     subj to 
     s_func {i in I}:
              x[i] complements y[i] = (sum {j in I} A[i,j]*x[j]) - a[i];

     subj to 
     func {i in I}:
         0 <= y[i] complements (sum {j in I} B[i,j]*x[j]) - b[i] >= 0;
\end{verbatim}
\double
Functions can also be doubly bounded as in the following definition.
\single\begin{verbatim}
     subj to 
     func {i in I}:
         x[i] complements
              l[i] <= (sum {j in I} A[i,j]*x[j]) - a[i] <= u[i];
\end{verbatim}
\double
A similar reformulation is performed by the interface:
\single\begin{verbatim}
     subj to
     s_func {i in I}:
         x[i] complements y[i] = (sum {j in I} A[i,j]*x[j]) - a[i];

     subj to 
     func {i in I}:
         l[i] <= y[i] <= u[i] complements x[i];
\end{verbatim}
\double

\subsection{Side Constraints}
AMPL provides extensions for complementarity problems with side
constraints.  For example, variables can be declared with bounds on
them.  However, the complementarity relationship is defined based upon
the bounds present in the constraint declaration.  For example, if
we define a complementarity problem as:
\single\begin{verbatim}
     set I;

     param A {I,I};
     param b {I};

     var x {I} >= 1, <= 10;

     subj to 
     func {i in I}:
         0 <= x[i] complements (sum {j in I} A[i,j]*x[j]) - a[i] >= 0;
\end{verbatim}
\double
then the $x$ variable will be constrained to be at least between 1
and 10.  These bounds are passed along as additional, ``side'', constraints
to the model.  Side constraints can be bounds on the variables
or other extraneous equations and inequalities.  These are dealt with in 
the solver interface by adding ``dummy'' variables to the problem.

However, complementarity problems with side constraints can be
difficult for solvers to process.  Therefore, the user will be warned
when side constraints are present in the model.  The behavior can be
modified by changing the value of the ``sideineq'' option which 
takes on values of:
\begin{enumerate}
\item[0] -- do not warn
\item[1] -- warn about side constraints
\item[2] -- exit with return code 503 when side constraints exist
\item[3] -- generate a run-time error when side constraints exist
\item[4] -- use alternative reformulation (without warnings)
\item[5] -- use alternative reformulation (with a warning)
\end{enumerate}
The default is to provide a warning message and reformulate the
problem so that the Jacobian matrix does not have empty columns.
The alternative reformulation adds variables to the problem, but
creates empty columns in the Jacobian matrix.

%%% Put in the reformulations here...

\textit{Since square complementarity problems are typically more amenable
to solution than complementarity problems with side constraints, we highly 
recommend that a user avoid generating a model containing side 
constraints.}

\subsection{Presolve}
AMPL does some preprocessing for all of the mathematical programs that
it generates.  For the above model, the declared bounds on $x$ are
tighter than the bounds in the complementarity condition.  Since
none of the bounds in the complementarity condition can ever be
satisfied, they are redundant and removed.  The complementarity
condition is then replaced by the single constraint:
\single\begin{verbatim}
     subj to 
     func {i in I}:
         x[i] complements (sum {j in I} A[i,j]*x[j]) - b[i] = 0;
\end{verbatim}
\double
The declared bounds, $1 \leq x \leq 10$, are added as side
constraints to the model.  If however, the bounds in the
complementarity conditions are tighter, the relevant bounds in the
variable declaration are dropped and not passed to the solver.

For complementarity problems, the AMPL presolve phase can create side
constraints not otherwise present.  Particularly, when the value for
a variable becomes known, the corresponding complementary function will 
be passed along as a side constraint.  Therefore, if a side constraint 
warning is encountered, and the user believes they have a square 
complementarity problem, the AMPL presolve should be turned off with
the command: 
\single\begin{verbatim}
     ampl: option presolve 0;
\end{verbatim}
\double
before the solve statement.

\textit{Note that the introduction of side constraints by the AMPL 
presolve can negatively impact an algorithm.  In such cases, turning 
the AMPL presolve off could lead to improved performance by the solver.}

\section{GAMS}\label{sec:gams}

% Input format
% Non-matching constraints

GAMS \cite{brooke.kendrick.ea:gams} declares complementarity conditions in the
model definition statement using the $.$ operator, a surrogate for the 
generalized $\perp$ notation used when defining a complementarity problem 
mathematically.  An implementation of the transportation model from 
\secref{transportation} can be found in \figref{transmcp-gams}; the actual 
data for the model is assumed to be given in the file \texttt{transmcp.inc}.  
Example data is found in \figref{transmcp-gams2}.

\begin{figure}
\begin{verbatim}
sets   i                Suppliers,
       j                Markets;

parameter
       cap(i)           Capacity of supplier,
       dem(j)           Demand at market,
       tra(i,j)         Transportation cost;

$include transmcp.inc

positive variables
       x(i,j)           Shipment quantity,
       p_s(i)           Price at plant,
       p_d(j)           Price at market;

equations
       supply(i)        Supply limit,
       demand(j)        Demand satisfaction,
       rational(i,j)    Rational marketplace; 
 
supply(i)..     cap(i) =g= sum(j, x(i,j));

demand(j)..     sum(i, x(i,j)) =g= dem(j);

rational(i,j).. p_s(i) + tra(i,j) =g= p_d(j);

model transport / rational.x, demand.p_d, supply.p_s /;

solve transport using mcp;
\end{verbatim}
\caption{Transportation model in GAMS, \texttt{transmcp.gms}}
\label{fig:transmcp-gams}
\end{figure}

\begin{figure}
\begin{verbatim}
set  i   / seattle, san-diego /,
     j   / new-york, chicago, topeka /;

parameter cap(i) /
     seattle    325,
     san-diego  575
/;

parameter dem(j) /
     new-york   325,
     chicago    300,
     topeka     275
/;

table tra(i,j)
           new-york   chicago    topeka
seattle     0.225      0.153     0.162
san-diego   0.225      0.162     0.126;
\end{verbatim}
\caption{Transportation data in GAMS, \texttt{transmcp.inc}}
\label{fig:transmcp-gams2}
\end{figure}

The model corresponds very closely to \eqnref{transmcp}.  In GAMS, the $\perp$
sign is replaced in the \texttt{model} statement with a ``.''.  The pairing of
variables and equations is performed precisely at this point in the GAMS code.
Therefore, in \figref{transmcp-gams}, the function defined by \texttt{rational}
is complementary to the variable \texttt{x}.  The standard GAMS statements on 
the variables can be used to inform a solver of the bounds, namely (for a 
declared variable $z(i)$):
\single\begin{verbatim}
     z.lo(i) = 0;
\end{verbatim}
\double
or 
\single\begin{verbatim}
     positive variable z;
\end{verbatim}
\double
Further information on the GAMS syntax can be found in 
\cite{rutherford:extensions}.  Note that positive variables are matched with 
\texttt{=g=} constraints in this model.  However, using \texttt{=l=} in an 
equivalent declaration for the \texttt{supply} constraint
\single\begin{verbatim}
     supply(i)..     sum(j, x(i,j)) =l= cap(i);
\end{verbatim}
\double
results in an illegal complementarity problem and GAMS gives a model 
generation error. 

A GAMS implementation of \eqnref{walras2} is given in \figref{walrasian-gams},
where we have introduced a defining variable, \texttt{d\_p(i)}, and 
equation, \texttt{D(i)}, for the demand function.  Note that free variables 
must be paired with equations.  Thus in \figref{walrasian-gams}, \texttt{d\_p}
is a free variable and its paired equation \texttt{demand} is an equality.  
Since \texttt{p} is nonnegative, its paired relationship \texttt{S} is a 
(greater-than) inequality.  

\begin{figure}
\begin{verbatim}
$include walras.inc

positive variables p(i), y(j); 
variables d_p(i);
equations S(i), L(j), D(i); 

S(i)..   b(i) + sum(j, A(i,j)*y(j)) - d_p(i) =g= 0 ;

L(j)..   -sum(i, p(i)*A(i,j)) =g= 0 ;

D(i)..   d_p(i) =e= c(i)*sum(k, b(k)*p(k)) / p(i) ;

model walras / S.p, L.y, D.d_p /;

p.l(i) = 1;
solve walras using mcp;
\end{verbatim}
\caption{Walrasian equilibrium as an MCP in GAMS, \texttt{walrasian.gms}}
\label{fig:walrasian-gams}
\end{figure}

The demand function in \figref{walrasian-gams} is undefined when the prices 
are zero.  Since GAMS uses the projection of zero onto the variable bounds 
as a default starting point, we must specify an alternate starting point
to avoid a GAMS execution error.  The statement:
\single\begin{verbatim}
     p.l(i) = 1;
\end{verbatim}
\double
is used to specify an initial point where the function is well-defined.  
\textit{In general, a starting point close to a solution will typically 
help the algorithm to quickly solve the model.}

\begin{comment}
\begin{figure}
\begin{verbatim}
set i /1*3/;
set j /1*2/;

parameter b(i) /
    1 3,
    2 4,
    3 5
/;

paramater c(i) /
    1 5,
    2 4,
    3 3
/;

table A(i,j)
     1    2
1    7   -2
2    2   -6
3    1   -10;
\end{verbatim}
\caption{Walrasian equilibrium data in GAMS, \texttt{walrasian.inc}}
\label{fig:walrasian-gams2}
\end{figure}
\end{comment}

A simplification is allowed to the model statement in
\figref{walrasian-gams}.  In many cases, it is not significant to match free 
variables explicitly to equations; we only require that there are the same 
number of free variables as equations.  Thus, in the example of 
\figref{walrasian-gams}, the model statement could be replaced by
\single\begin{verbatim}
     model walras / S.p, L.y, D /;
\end{verbatim}
\double
This extension allows existing GAMS models consisting of a square
system of nonlinear equations to be easily recast as a complementarity
problem - the model statement is unchanged.  GAMS generates a list of all 
variables appearing in the equations found in the model statement, performs 
explicitly defined pairings and then checks that the number of remaining 
equations equals the number of remaining free variables.  All variables that 
are not free and all inequalities must be explicitly matched.

\textit{Equations should be explicitly matched to free variables to preserve 
matrix properties, such as positive semidefiniteness, that can be lost when an 
arbitrary matching is performed.  The loss of structure can lead to failure or
poor performance by an algorithm.}  For example, when constructing the first 
order optimality conditions for a linear program, the matching of equations to
multipliers should be done to preserve the skew symmetric structure.

We can place an upper bound on the activity level in the Walrasian
equilibrium by simply adding an upper bound to $y$ 
\single\begin{verbatim}
     y.up(j) = 10;
\end{verbatim}
\double
and replacing the loss equation with the following definition:
\single\begin{verbatim}
     L(j)..   -sum(i, p(i)*A(i,j)) =e= 0;
\end{verbatim}
\double
For bounded variables, we do not know beforehand if the constraint will be 
satisfied as an equation, less than inequality or greater than inequality, 
since this determination depends on the values of the solution variables.  
We adopt the convention that all bounded variables are paired to equations.

Note that if a variable $z_i$ is fixed, then $F_i(z)$ is unrestricted since 
precisely one of the three conditions in the MCP definition automatically 
holds when $z_i = \ell_i = u_i$.  Thus if a variable is fixed in a GAMS model,
the paired equation is completely dropped from the model.  This convenient
modeling trick can be used to remove particular constraints from a model at
generation time.

Furthermore, we can scale particular equations and variables in GAMS using 
the \texttt{.scale} attribute.  For example, we can write the statements
\single\begin{verbatim}
     D.scale(i) = 5;
     d_p.scale(i) = 0.5;
\end{verbatim}
\double
to specify scalings for the \texttt{D} equation and \texttt{d\_p} variable.
To turn the scaling on for the model we set the \texttt{.scaleopt} model 
attribute.  
\single\begin{verbatim}
     walras.scaleopt = 1;
\end{verbatim}
\double     
In certain circumstances, scaling can help an algorithm to solve the problem 
specified.

% Options
% Solver specification

To solve the problem, the modeler executes the command:
\single\begin{verbatim}
     gams transmcp
\end{verbatim}
\double
where \texttt{transmcp} can be replaced by any filename containing a
GAMS model.  
Many command line options for GAMS exist; the reader is 
referred to \cite{brooke.kendrick.ea:gams} for further details as well as the 
on-line documentation at
\single\begin{verbatim}
     http://www.gams.com/
\end{verbatim}
\double
One useful command line option allows a user to change the solver by issuing 
the command:
\single\begin{verbatim}
     gams transmcp mcp=path
\end{verbatim}
\double
which sets the MCP solver to PATH 4.x.

The solver can be changed in many different ways, including the following:
\begin{enumerate}
\item Add the following line to the \texttt{transmcp.gms} file prior to the 
{\tt solve} statement
\single\begin{verbatim}
     option mcp = path;
\end{verbatim}
\double
where \texttt{path} can be any solver for complementarity problems.
In this example, PATH 4.x will be used instead of the default solver
provided.  Note: setting the solver in the GAMS code takes precedence 
to setting the solver with a command line option.
\item Rerun the \texttt{gamsinst} program from the GAMS system directory
and choose the desired algorithm as the default solver for MCP.
\end{enumerate}
The SEMI algorithm is used by changing the MCP solver to \texttt{semi}.

The default algorithm options should be sufficient for most models; the 
technique for changing these options are now described.  
To change the default options on the model \texttt{transport}, the modeler
is required to write a file \texttt{path.opt} for the PATH 4.x algorithm
in the working directory and either add a line
\single\begin{verbatim}
     transport.optfile = 1;
\end{verbatim}
\double
before the \texttt{solve} statement in the file \texttt{transmcp.gms},
or use the command-line option
\single\begin{verbatim}
     gams transmcp optfile=1
\end{verbatim}
\double
Unless the modeler has changed the WORKDIR parameter explicitly, the
working directory will be the directory containing the model file. 
The procedure is the same when using the SEMI algorithm.  However, the
options file will be called \texttt{semi.opt}.
% See \partref{software} for complete listings of the available options
% for each of the algorithms.

GAMS controls the total number iterations allowed via the
{\tt iterlim} option.  If more iterations are needed for a particular model
then either of the following lines should be added to the
{\tt transmcp.gms} file before the solve statement 
\single\begin{verbatim}
     option iterlim = 2000;
     transport.iterlim = 2000;
\end{verbatim}
\double
An iteration in the PATH 4.x and SEMI codes consists of finding a solution
to a linear system of equations $Ax = b$ for some matrix $A$ and right hand 
side $b$.

Similarly, if the solver runs out of memory, then the workspace allocated can 
be changed using
\single\begin{verbatim}
     transport.workspace = 20;
\end{verbatim}
\double
The above example would allocate 20MB of workspace for solving the model.

After a \texttt{solve} statement is encountered, control is handed over to 
the solver which creates a log providing information on what the solver is 
doing as time elapses.
% \partref{software} details the output generated in the log file 
% for typical executions of the available algorithms.
After the solver terminates, a listing file is generated containing the
problem solution.  We now describe the output in the listing file
specifically related to the complementarity problem.

\subsection{Listing File}\label{sec:listing}
% Return codes and what they mean
% Multipliers
% Objective function value

The listing file is the standard GAMS mechanism for reporting model
results.  This file contains information regarding the compilation
process, the form of the generated equations in the model, and a
report from the solver regarding the solution process.  

We now detail the last part of this output for the PATH 4.x solver, an example 
of which is given in \figref{listing-gams}.  We use ``...'' to indicate where
we have omitted continuing similar output.

\begin{figure}
\begin{verbatim}
               S O L V E      S U M M A R Y

     MODEL   transport
     TYPE    MCP
     SOLVER  PATH                FROM LINE  50

**** SOLVER STATUS     1 NORMAL COMPLETION
**** MODEL STATUS      1 OPTIMAL

 RESOURCE USAGE, LIMIT          0.145     1000.000
 ITERATION COUNT, LIMIT        18        10000
 EVALUATION ERRORS              0            0

 Work space allocated           --    0.06 Mb

---- EQU RATIONAL         

                      LOWER     LEVEL     UPPER    MARGINAL  

seattle  .new-york    -0.225    -0.225     +INF     50.000      
seattle  .chicago     -0.153    -0.153     +INF    300.000      
seattle  .topeka      -0.162    -0.126     +INF       .         

...

---- VAR X           shipment quantities in cases

                      LOWER     LEVEL     UPPER    MARGINAL  

seattle  .new-york      .       50.000     +INF       .         
seattle  .chicago       .      300.000     +INF       .         

...

**** REPORT SUMMARY :        0     NONOPT
                             0 INFEASIBLE
                             0  UNBOUNDED
                             0  REDEFINED
                             0     ERRORS
\end{verbatim}
\caption{Listing file solving \texttt{transmcp.gms} in GAMS}
\label{fig:listing-gams}
\end{figure}

After a summary line indicating the model name and type and the solver name, 
the listing file shows a solver status and a model status.  
\tabref{solstat-gams} and \tabref{modelstat-gams} display the relevant
codes that are returned under different circumstances.  A modeler can
access these codes within the \texttt{transmcp.gms} file using 
\texttt{transport.solstat} and
\texttt{transport.modelstat} respectively. 
An objective value is also reported, which can be accessed as 
\texttt{transport.objval}.  This objective measures the complementarity
error, which is zero at a solution to the problem and positive elsewhere.

\begin{table}
\centering
\caption{Solution Status Codes for GAMS}\label{tab:solstat-gams}
\begin{tabular}{llp{8cm}}
Code & String & Meaning\\
\hline
1 & Normal completion & Solver returned to GAMS without an error\\
2 & Iteration interrupt & Solver used too many iterations \\
3 & Resource interrupt & Solver took too much time \\
4 & Terminated by solver & Solver encountered difficulty and was unable to continue \\
8 & User interrupt & The user interrupted the solution process \\
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Model Status Codes for GAMS}\label{tab:modelstat-gams}
\begin{tabular}{lll}
Code & String & Meaning\\
\hline
1 & Optimal & Solver found a solution of the problem\\
6 & Intermediate infeasible & Solver failed to solve the problem\\
\end{tabular}
\end{table}

After the solver and model status is reported, a listing of the time and 
iterations used by the algorithm is given, along with a count of the number 
of evaluation errors encountered.  If the
number of evaluation errors is greater than zero, further information
can typically be found later in the listing file prefaced by
the ``****'' string.  Information provided by the solver is 
then displayed.

Next comes the solution listing starting with each of the equations in
the model.   For each equation passed to the solver, four
columns are reported, namely the lower bound, level, upper bound and
marginal.  GAMS moves all parts of a constraint involving variables to
the left hand side and accumulates the constants on the right hand
side.  The lower and upper bounds correspond to the constants that
GAMS generates.  For equations, these should be equal, whereas for
inequalities one of them should be infinite.  The level value 
of the equation (an evaluation of the left hand side of the constraint
at the current point) should
be between these bounds, otherwise the solution is infeasible and the
equation is marked as follows:
\single\begin{verbatim}
     seattle  .chicago     -0.153    -2.000     +INF    300.000 INFES      
\end{verbatim}
\double
The marginal column for an equation returns the level value of the variable 
that was matched to the equation.  If the modeler did not pair a particular
equation with a variable, the value returned is that of the variable that 
the solver interface paired with the equation.

For the variable listing, the lower, level, and upper columns indicate
the lower and upper bounds on the variables and the solution value.
The level value returned by PATH 4.x and SEMI will always be 
between these bounds.  The marginal column contains the value of the slack 
on the equation that was paired
to the variable.  If a variable appears in one of the constraints in
the model statement but is not explicitly paired to a constraint, the
slack reported here contains the internally matched constraint slack.
The definition of this slack is the minimum of \texttt{equ.l - equ.lower}
and \texttt{equ.upper - equ.l}, where \texttt{equ} is the paired equation.
 
Finally, a summary report is given that indicates how many errors were
found.  \figref{listing-gams} is a good case; when the model has
infeasibilities, these can be found by searching for the string
``INFES'' as described above.

\subsection{Redefined Equations}
Unfortunately, the above description is incomplete because some of the
equations may have the following form:
\single\begin{verbatim}
                 LOWER     LEVEL      UPPER    MARGINAL  
     new-york   325.000   350.000    325.000   0.225 REDEF
\end{verbatim}
\double
This notation should be construed as a warning from GAMS, as opposed to 
an error.  In principle, the
\texttt{REDEF} should only occur if the paired variable to the constraint
had a finite lower and upper bound and the variable is at one of those bounds.
In this case, at the solution of the
complementarity problem the ``equation (\texttt{=e=})'' may not be 
satisfied.  The
problem occurs because of a limitation in the GAMS syntax for
complementarity problems.  The GAMS equations are used to define the
function $F$.  The bounds on the function $F$ are derived from the
bounds on the associated variables.  Before solving the problem, for
finite bounded variables, we do not know if the associated function
will be positive, negative or zero at the solution.  Thus, we do not
know whether to define the equation as ``\texttt{=e=}'', ``\texttt{=l=}'' 
or ``\texttt{=g=}''.
GAMS therefore allows any of these, and informs the modeler via the
``REDEF'' label that internally GAMS has redefined the bounds so that
the solver processes the correct problem, but that the solution given
by the solver does not satisfy the original bounds.  
However, in practice, a \texttt{REDEF} can also occur when the
equation is defined using ``\texttt{=e=}'' and the variable has a single
finite bound.  This matching is allowed by GAMS, and as above, if the
variable is at its bound at a solution to the complementarity problem,
then the function $F$ may not satisfy the ``\texttt{=e=}'' relationship.

\begin{comment}
\tabref{gams-matching} provides the possible matchings between variables
and functions.

\begin{table}
\centering
\caption{Possible matchings between variables and functions in GAMS}\label{tab:gams-matchine}
\begin{tabular}{cc}
Variable Type & Equation Type \\
\hline
Two finite bounds & =g=, =e= , =l= \\
Finite lower bound, infinite upper bound & =g=, =e= \\
Finite upper bound, infinite lower bound & =l=, =e= \\
Two infinite bounds & =e=
\end{tabular}
\end{comment}

\textit{Note that
this is not an error, just a warning.  The solver has solved the
complementarity problem specified by this equation.}  GAMS gives this
report to ensure that the modeler understands that the complementarity
problem derives the relationships on the equations from the variable bounds,
not from the equation definition.

\begin{comment}
\subsection{C-Links}
There are two version of the algorithms available for GAMS.  For 
example, we have both PATH and PATHC.  The only difference between the two is
the link to GAMS with the remaining code identical.  The C-Links
\textit{do not support MPSGE models}, but enable the use of preprocessing 
and can be used to solve constrained systems of nonlinear equations.  The 
output for C-links is identical to the main distribution.

Many square nonlinear systems are formulated as
\[ F(z) = 0,\; \ell \leq z \leq u , \]
where typically, the bounds on $z$ are inactive at the solution.
This is {\em not\/} an MCP, but is an example of a 
``constrained nonlinear system'' (CNS).
It is important to note the distinction between MCP and CNS.
The MCP uses the bounds to infer relationships on the function
$F$.  If a finite
bound is active at a solution, the corresponding component of $F$ is only
constrained to be nonnegative or nonpositive in the MCP, whereas in CNS
it must be zero.  Thus there may be many solutions of MCP that
do not satisfy $F(z) = 0$.  Only if $z^*$ is a solution of MCP
with $\ell < z^* < u$ is it guaranteed that $F(z^*) = 0$.

Internally, the C-Links reformulate a constrained nonlinear system of 
equations to an equivalent complementarity problem.  The reformulation adds 
variables, $y$, with the resulting problem written as:
\[
\begin{array}{rcl}
\ell \leq x \leq u & \perp & -y \\
y \mbox{ free}& \perp & F(x).
\end{array}
\]
This is the MCP passed on to the algorithm.
\end{comment}

\subsection{Pitfalls}
The ordering of an equation is important in the
specification of an MCP.  Since the data of the MCP is the function
$F$ and the bounds $\ell$ and $u$, it is important for the modeler to
pass the function $F$ and not $-F$ to the solver.

\begin{figure}
\begin{verbatim}
variables x;
equations d_f;

x.lo = 0;
x.up = 2;

d_f.. 2*(x - 1) =e= 0;

model first / d_f.x /;
solve first using mcp;
\end{verbatim}
\caption{First order conditions as an MCP in GAMS, \texttt{first.gms}}
\label{fig:first}
\end{figure}

For example, if we have the optimization problem,
\[
\min_{x \in [0, 2]} (x - 1)^2
\]
then the first order optimality conditions are
\[
  \begin{array}{rcl}
0 \leq x \leq 2 & \perp & 2(x - 1)
  \end{array}
\]
which has a unique solution, $x = 1$.  \figref{first} provides correct
GAMS code for this problem.  However, if we accidentally write the valid 
equation
\single\begin{verbatim}
     d_f.. 0 =e= 2*(x - 1);
\end{verbatim}
\double
the problem given to the solver is
\[
  \begin{array}{rcl}
0 \leq x \leq 2 & \perp & -2(x - 1)
  \end{array}
\]
which has three solutions, $x = 0$, $x = 1$, and $x = 2$.  This problem
in fact formulates the stationary conditions for the non-convex quadratic 
problem,
\[
\max_{x \in [0, 2]} (x - 1)^2,
\]
not the problem we intended to solve.

\section{MATLAB}\label{sec:matlab}

MATLAB \cite{matlab:users} offers a convenient environment for exploration 
and visualization.  Both linear and nonlinear complementarity problems 
can be written in MATLAB and conveyed to the PATH 4.x and SEMI solvers
through provided m-functions.  

Furthermore, the MATLAB interfaces are unique in that polyhedrally 
constrained variational inequalities can be specified and solved.  
Mathematically, the polyhedrally constrained variational inequality 
is given a set $C \defined \{ z | Az \geq b\}$, and a function, $F$.  We 
then want to find a point, $z \in C$ where:
\[
     \ip{F(z)}{\bar{z} - z} \geq 0 \; \mbox{ for all } \bar{z} \in C.
\]
As mentioned in \secref{complementarity} when $C$ is the Cartesian product 
of closed (not necessarily bounded) intervals, the (box constrained) 
variational inequality is equivalent to the mixed complementarity problem.  
The general polyhedral constraints are handled by adding multipliers, 
resulting in a standard mixed complementarity problem.

\subsection{Linear Complementarity Problems}\label{sec:mat-lin}
The routines tailored for linear complementarity problems offer improved 
performance over the nonlinear counterparts, as they do not make any 
call-backs to MATLAB for function evaluations.  Furthermore, the preprocessor 
for complementarity problems developed in \chapref{preprocessing} is enabled 
when using the MATLAB interface for linear models.

The m-function for solving standard linear complementarity problems,
\[
\begin{array}{lcr}
0 \leq z & \perp & Mz + q \geq 0,
\end{array}
\]
is defined in \texttt{pathlcp.m} for PATH 4.x and \texttt{semilcp.m} for 
SEMI.  To solve a standard linear complementarity problem using PATH 4.x, 
we write the MATLAB command:
\single\begin{verbatim}
     >> z = pathlcp(M, q);
\end{verbatim}
\double
An implementation of the transportation model from \secref{transportation} 
can be found in \figref{transmcp-matlab}.  Note that the matrix $M$ passed
to the \texttt{pathlcp} routine will be automatically converted into a sparse 
matrix.
The conversion can be avoided if the matrix passed into the function is
already stored as a sparse matrix.

\begin{figure}
\begin{smallverb}
>> M = [
   0   0   0   0   0   0  -1   0   0   1   0;
   0   0   0   0   0   0   0  -1   0   1   0;
   0   0   0   0   0   0   0   0  -1   1   0;
   0   0   0   0   0   0  -1   0   0   0   1;
   0   0   0   0   0   0   0  -1   0   0   1;
   0   0   0   0   0   0   0   0  -1   0   1;
   1   0   0   1   0   0   0   0   0   0   0;
   0   1   0   0   1   0   0   0   0   0   0;
   0   0   1   0   0   1   0   0   0   0   0;
  -1  -1  -1   0   0   0   0   0   0   0   0;
   0   0   0  -1  -1  -1   0   0   0   0   0
];
>> q = [
   0.225;
   0.153;
   0.162;
   0.225;
   0.162;
   0.126;
-325.000;
-300.000;
-275.000;
 325.000;
 575.000
];
>> z = pathlcp(M, q)

z =

   25.0000
  300.0000
         0
  300.0000
         0
  275.0000
   87.0291
   86.9571
   86.9301
   86.8041
   86.8041
\end{smallverb}
\caption{Transportation problem in MATLAB}\label{fig:transmcp-matlab}
\end{figure}

Linear mixed complementarity problems can be specified by adding lower and 
upper bounds:
\single\begin{verbatim}
     >> z = pathlcp(M, q, l, u);
\end{verbatim}
\double
which solves the problem
\[
\begin{array}{lcr}
l \leq z \leq u & \perp & Mz + q.
\end{array}
\]
The numerical value used by the MATLAB interfaces for $\infty$ 
is $10^{20}$.

Many times, a good starting point can help an algorithm to quickly
find a solution.  Therefore, a starting point can be passed as an additional 
argument:
\single\begin{verbatim}
     >> z = pathlcp(M, q, l, u, z);
\end{verbatim}
\double
Note that a call of the form \texttt{pathlcp(M, q, [], [], z)} is
handled by using $\ell = \{0\}^n$ and $u = \{\infty\}^n$ as default 
values.  Therefore, 
this statement solves the standard linear complementarity problem using a
provided starting point.

Additional arguments to \texttt{pathlcp} are supplied for conveying
polyhedrally constrained variational inequalities.  Internally, multipliers 
on the constraints are added resulting in a linear mixed complementarity
problem.  The constraints $Az \geq b$ are passed to the solver with additional
arguments for $A$ and $b$:
\single\begin{verbatim}
     >> z = pathlcp(M, q, l, u, z, A, b);
\end{verbatim}
\double
The inequality type can be passed in an extra vector; a negative value 
indicates that we have a less than or equal to inequality, a positive value 
means greater than or equal, and zero means an equation.  The multipliers 
can be passed to and returned from the \texttt{pathlcp} function call.
\single\begin{verbatim}
     >> [z, mu] = pathlcp(M, q, l, u, z, A, b, t, mu);
\end{verbatim}
\double
The transportation model can be more succinctly posed as a
polyhedrally constrained variational inequality as given in
\figref{transmcp-matlab2}. 

\begin{figure}
\begin{smallverb}
>> A = [
   1   0   0   1   0   0;
   0   1   0   0   1   0;
   0   0   1   0   0   1;
  -1  -1  -1   0   0   0;
   0   0   0  -1  -1  -1
];
>> b = [
   325;
   300;
   275;
  -325;
  -575
];
>> q = [
    0.225;
    0.153;
    0.162;
    0.225;
    0.162;
    0.126
];
>> [z, mu] = pathlcp(sparse(6,6), q, [], [], [], A, b)

z =

   25.0000
  300.0000
         0
  300.0000
         0
  275.0000

mu =

   87.0291
   86.9571
   86.9301
   86.8041
   86.8041
\end{smallverb}
\caption{Transportation model as a variational inequality in MATLAB}\label{fig:transmcp-matlab2}
\end{figure}

The polyhedrally constrained variational inequality can be used to
easily write the optimality conditions for linear and convex quadratic 
programs.  For example, if we have the quadratic program:
\[
\begin{array}{rc}
\min_x & \frac{1}{2}x^TQx + c^Tx \\
\text{subject to} & Ax \geq b \\
& l \leq x \leq u
\end{array}
\]
with $Q$ a symmetric, positive semidefinite matrix, then the necessary
and sufficient optimality conditions can be formulated and passed with
the command:
\single\begin{verbatim}
     >> [x, mu] = pathlcp(Q, c, l, u, [], A, b)
\end{verbatim}
\double
For convenience, two m-functions for solving linear programs and
convex quadratic programs are provided:
\single\begin{verbatim}
     >> [x, mu] = lp_lcp(c, A, b, l, u, x, t);
     >> [x, mu] = qp_lcp(Q, c, A, b, l, u, x, t);
\end{verbatim}
\double
These functions pass along the above optimality conditions to the 
\texttt{pathlcp} function.

The default algorithm options should be sufficient for most models.  However,
solver specific options can be modified by writing a \texttt{path.opt} file 
for PATH 4.x or a \texttt{semi.opt} file for SEMI.
In cases where the algorithm fails to find a solution, the logfile
written in \texttt{logfile.tmp} can be investigated
to gain information about the reasons for failure.  Note that if the 
algorithm fails to solve the problem, a MATLAB error is
generated.  This behavior is part of the m-function, so a user could
modify the termination behavior if desired.

\subsection{Nonlinear Complementarity Problems}

Nonlinear mixed complementarity problems require that a user provide
an m-function that evaluates the nonlinear function and Jacobian at a 
specified point.  The calling syntax is:
\single\begin{verbatim}
     z = pathmcp(z, l, u, 'mcp_funjac');
\end{verbatim}
\double
which uses PATH 4.x where \texttt{mcp\_funjac} is the name of the nonlinear 
function evaluation routine.  The \texttt{z} vector must be explicitly 
provided in
order to determine the size of the complementarity problem.  For nonlinear 
models the algorithm may perform much better by choosing an appropriate
starting point close to a solution.  $\ell = \{0\}^n$ and $u = \{\infty\}^n$ 
are the 
default values used if the lower and upper bounds are not present.
That is, we solve a standard nonlinear complementarity problem.

The corresponding m-function \texttt{mcp\_funjac.m} contains the 
definition
\single\begin{verbatim}
     function [F, J, domerr] = mcp_funjac(z, jacflag) 
\end{verbatim}
\double
that computes the function, $F$, and, if \texttt{jacflag=1}, the 
sparse Jacobian $J$ at the point $z$.  The \texttt{domerr} variable returns 
the number 
of domain violations encountered during the evaluation.  A domain violation 
occurs for example when we have division by zero.  \textit{The Jacobian must 
be passed to the solver as a sparse matrix.}  As an example, a function for
the transportation model from \secref{transportation} with isoelastic 
demand is provided in \figref{matlab-isoelastic}.

\begin{figure}
\begin{verbatim}
function [F, J, domerr] = mcp_funjac(z, jacflag)

% initialize
z = z(:);
F = [];
J = [];
domerr = 0;

% obtain variable values
x = z(1:6);     % quantities
p_d = z(7:9);   % demand price
p_s = z(10:11); % supply price

% check for domain violation
if (p_d <= 0)
  domerr = 1;
  return;
end

cost   = [ 0.225; 0.153; 0.162; 0.225; 0.162; 0.126; ];
demand = [ -325; -300; -275; ] ./ (p_d.^0.5);
supply = [ 325; 575; ];

A = sparse( [
   1   0   0   1   0   0;
   0   1   0   0   1   0;
   0   0   1   0   0   1;
  -1  -1  -1   0   0   0;
   0   0   0  -1  -1  -1;
] );

F = [ (-A'*[p_d; p_s] + cost); A*x + [demand; supply] ];

if (jacflag)
  J = [ sparse(6,6) -A';
        A diag( [ -0.5 * [-325; -300; -275;] ./ (p_d.^1.5); 0; 0] ) ];
end

return
\end{verbatim}
\caption{Transportation model with isoelastic demand in MATLAB}\label{fig:matlab-isoelastic}
\end{figure}

Providing the derivative explicitly can be time consuming and prone
to error.  Therefore, we provide a convenience routine that uses the
\texttt{ADMAT} software \cite{coleman.verma:admat} to perform automatic
differentiation of a supplied m-function.  Assuming that \texttt{ADMAT}
has been installed according to the user guide, we can call 
PATH 4.x with the commands:
\single\begin{verbatim}
     >> mcp_set_adfunc('mcp_fun');
     >> z = pathmcp(z, l, u, 'mcp_adfunjac');
\end{verbatim}
\double
The first line of the code sets a global variable for the user defined 
function with syntax
\single\begin{verbatim}
     function F = mcp_fun(z, extra)
\end{verbatim}
\double
Note that the \texttt{extra} argument is needed by the \texttt{ADMAT}
software, but should not be used in the function evaluation code.
The second line of code tells the interface to use the provided 
\texttt{mcp\_adfunjac} function to perform the function and Jacobian
evaluations.  The implementation of \texttt{mcp\_adfunjac} is:
\single\begin{verbatim}
     function [F,J,domerr] = mcp_adfunjac(z,jacflag)

     global mcp_adfun;

     if jacflag == 1
       [F, J] = evalJ(mcp_adfun, z);
     else
       F = eval(mcp_adfun, z);
       J = [];
     end;

     domerr = 0;
     return;
\end{verbatim}
\double
The global variable \texttt{mcp\_adfun} is used to keep track of the
m-function that the user specified for their function.  Note that the
\texttt{evalJ} uses \texttt{ADMAT} to evaluate the function and Jacobian.
Currently, domain violations are not reported by the automatic 
differentiator.  Therefore, we always assume that no
domain violations occurred.

\begin{comment}
An implementation of the MATLAB function for the Walrasian equilibrium model 
from \secref{walrasian} can be found in \figref{walrasian-matlab}.  We have
used some global vectors to declare the problem parameters.  We can
run the model with the following:
\single\begin{verbatim}
     >> % Set up data
     >> mcp_set_adfunc('walras');
     >> x = ones(n+m+d);
     >> z = pathmcp(x,l,u,'mcp_adfunjac');
\end{verbatim}
\double
We remark the a domain violation occurs when the prices in the model
are zero.  Therefore, we need to specify a starting point.

\begin{figure}
\begin{verbatim}
function F = walras(z,Extra)

return;
\end{verbatim}
\caption{Walrasian equilibrium function in MATLAB, \texttt{walras.m}}\label{fig:walrasian-matlab}
\end{figure}
\end{comment}

The function and Jacobian at the returned solution can be obtained by adding
additional output parameters:
\single\begin{verbatim}
     [z, F, J] = pathmcp(z, l, u, 'mcp_funjac');
\end{verbatim}
\double
We also allow extensions for polyhedrally constrained variational inequalities 
with the syntax:
\single\begin{verbatim}
     [z, F, J, mu] = pathmcp(z, l, u, 'mcp_funjac', A, b, t, mu);
\end{verbatim}
\double
where \texttt{A}, \texttt{b}, and \texttt{t} are used to convey 
the polyhedral set as in \secref{mat-lin}.
These are internally reformulated by adding multipliers so that a mixed
complementarity problem is solved.  The values for the multipliers can be
obtained with the \texttt{mu} arguments.

Algorithm options can be changed using the same technique as presented
in the section on linear complementarity problems.  Specifically, a
\texttt{path.opt} or \texttt{semi.opt} file is written with option,
value pairs.  A log file is also written in \texttt{logfile.tmp} that
can be examined.  As with the linear interface, if the algorithm fails to 
solve the problem, a MATLAB error is generated.  This behavior can
be changed by modifying the m-function.

\subsection{Mex Functions}

Most users will call the PATH 4.x or SEMI codes through the provided 
m-functions
for linear and nonlinear complementarity problems.  However, a user could 
bypass these m-function and directly invoke the appropriate mex-function.
Essentially, a mex-function contains a compiled version of a code that can
be accessed within MATLAB. To call the mex-function explicitly, we write:
\single\begin{verbatim}
    status = lcppath(n, nnz, z, l, u, M, q);
\end{verbatim}
\double
for linear complementarity problems and
\single\begin{verbatim}
    [status, F, J] = mcppath(n, nnz, z, l, u, func);
\end{verbatim}
\double
for nonlinear models.  Note that \texttt{lcpsemi} and \texttt{mcpsemi} can 
be used
for the SEMI solver with the same arguments.  A status of $1$ is returned 
if the model is solved and $0$ is returned otherwise.  $z$ is used as the
initial starting point as well as the solution vector.
The $F$ and $J$ output arguments to \texttt{mcppath} are not necessary.
However, if present, they will contain the function and Jacobian 
evaluation at the $z$ vector returned by the algorithm.

The number of nonzeros passed into the mex-routine should be an
\textbf{over-estimate} of the exact number of nonzeros when using a
nonlinear function.  The reason for this is that a fixed amount of
space is initially allocated by the solver for storing the Jacobian.
Therefore, in situations where the number of nonzeros in the Jacobian
changes, we could run into serious problems.  For example, the
(sparse) Jacobian of the function $x^2$ at zero is empty.  A heuristic 
suggested by Jorge Mor\'{e} is employed in the \texttt{pathmcp} m-function that 
perturbs the initial starting point slightly and evaluates the Jacobian at 
this perturbed point to calculate the over-estimate of the number of nonzeros 
in the Jacobian.  This technique is not foolproof, but has worked well in 
practice.

\subsection{Pitfalls}
Interrupting the solver in the middle of a computation can lead to 
segmentation violations on subsequent calls.  Therefore, if the solver
has been interrupted, a user should issue the command
\single\begin{verbatim}
     >> clear mex;
\end{verbatim}
\double
before attempting to solve another problem.  This statement unloads the 
mex-functions; subsequent calls will reload them.

\section{NEOS}\label{sec:neos}

% Input format
% Options
% Solver specification
% Output...
% Return codes and what they mean
%
% CONDOR, automatic differentiation

NEOS \cite{czyzyk.mesnier.ea:network-enabled} enables users to submit 
optimization problems across the Internet for solution.  While complementarity
problems can be solved by sending models written in the AMPL and GAMS modeling
languages to the NEOS server, an additional interface is provided for users 
wanting to supply a function evaluation routine written in FORTRAN.  

%Essentially, a user provides the NEOS server with the function evaluation 
%routine.  The NEOS server then sends this routine to a remote machine
%where it is compiled, linked with appropriate libraries, and executed.  The 
%output from the execution is then returned to the user.  Derivatives for the 
%function are created using ADIFOR \cite{bischof.carle.ea:adifor}. 

The essential routine the user provides evaluates the function, $F$,
at a specified point.  The function header is as follows:
\single\begin{verbatim}
     subroutine fcn(n,x,F)
\end{verbatim}
\double
where $n$ is the number of variables in the problem, $x$ is a vector
of length $n$ containing the specified point, and $F$ is a vector of
length $n$ which returns the user defined function evaluation at the
point $x$.  For example, in the transportation model, we have 11
variables and the function would be implemented as provided in
\figref{transmcp-neos}.  By default, the complementarity problem uses
lower bounds of zero and upper bounds of infinity, a standard nonlinear
complementarity problem, and an initial starting point of zero.

\begin{figure}
\begin{verbatim}
       subroutine fcn(n,x,F)
       implicit none

       integer n
       double precision x(n), F(n)

C      Rational marketplace constraints
       F(1)  = x(10) - x(7) + 0.225d+0
       F(2)  = x(10) - x(8) + 0.153d+0
       F(3)  = x(10) - x(9) + 0.162d+0
       F(4)  = x(11) - x(7) + 0.225d+0
       F(5)  = x(11) - x(8) + 0.162d+0
       F(6)  = x(11) - x(9) + 0.126d+0

C      Capacity constraints
       F(7)  = x(1) + x(4) - 325.0d+0
       F(8)  = x(2) + x(5) - 300.0d+0
       F(9)  = x(3) + x(6) - 275.0d+0

C      Satisfaction of demand constraints
       F(10) = 325.0d+0 - x(1) - x(2) - x(3) 
       F(11) = 575.0d+0 - x(4) - x(5) - x(6)

       return
       end
\end{verbatim}
\caption{Transportation function for NEOS, \texttt{fcn.f}}\label{fig:transmcp-neos}
\end{figure}

Having written the function, we then submit the problem to the NEOS server.  
There are several different methods available for submission as documented on 
the NEOS homepage:
\single\begin{verbatim}
    http://www-neos.mcs.anl.gov/
\end{verbatim}
\double
We concentrate on e-mail submissions.  In order to solve
the transportation problem through e-mail, we send a message to
\texttt{neos@mcs.anl.gov} containing the following: 
\single\begin{verbatim}
    TYPE = CP
    SOLVER = PATH

    N = 11

    BEGIN.FCN
        subroutine fcn(n,x,F)
        implicit none

        integer n
        double precision x(n), F(n)

    C   Actual function definition here
      
        return
        end
    END.FCN

    END-SERVER-INPUT
\end{verbatim}
\double
The results will then be mailed back to the person sending the mail.

Once the mail has been sent to NEOS, a remote machine is contacted
which compiles the function, calculates the derivative with ADIFOR 
\cite{bischof.carle.ea:adifor}, links the code with appropriate libraries, 
and runs the executable.  The mail reported to the user appears as below
for the transportation model, where we have used ... to denote that some 
information has been deleted to conserve space.
\single\begin{verbatim}
     ...

     ***PATH Output***
     Path v4.4a: Neos Link
     11 row/cols, 24 non-zeros, 19.83% dense.

     Could not open options file: path.opt
     Using defaults.
     Path 4.4a (Sun Jun 18 04:24:05 2000)

     SOLUTION FOUND.

     Major Iterations. . . . 14
     Minor Iterations. . . . 18
     Restarts. . . . . . . . 0
     Crash Iterations. . . . 1
     Gradient Steps. . . . . 0
     Function Evaluations. . 16
     Gradient Evaluations. . 16
     Total Time. . . . . . . 0.060000
     Residual. . . . . . . . 2.980724e-10


     SOLUTION VECTORS:
     x(1   ): +2.5000000000e+01    F(1   ): -2.5010549187e-13
     x(2   ): +3.0000000000e+02    F(2   ): +1.7894019599e-13
     x(3   ): +0.0000000000e+00    F(3   ): +3.5999999999e-02
     x(4   ): +3.0000000000e+02    F(4   ): +2.3306356844e-13
     x(5   ): +0.0000000000e+00    F(5   ): +9.0000000007e-03
     x(6   ): +2.7500000000e+02    F(6   ): -1.1846079673e-13
     x(7   ): +8.7029079521e+01    F(7   ): -2.2907897801e-11
     x(8   ): +8.6957079521e+01    F(8   ): +2.0571633286e-10
     x(9   ): +8.6930079521e+01    F(9   ): -1.3722001313e-10
     x(10  ): +8.6804079521e+01    F(10  ): +9.1517904366e-11
     x(11  ): +8.6804079521e+01    F(11  ): -1.3710632629e-10

     ...
\end{verbatim}
\double
The last part of the listing contains the log file generate by the
algorithm. %, which is fully documented in \partref{software}.
The solution vectors list $x$ and $F(x)$ for the point returned by the
algorithm.  A complete listing of the return messages can be found in
\tabref{return-neos}.

\begin{table}
\centering
\caption{Return messages for NEOS}\label{tab:return-neos}
\begin{tabular}{ll}
Message & Meaning \\
\hline
INFEASIBLE & The model is infeasible \\
ITERATION LIMIT & The iteration limit was reached without finding a solution\\
OTHER ERROR & An error occurred within the algorithm\\
SOLUTION FOUND & A solution was found by the algorithm \\
TIME LIMIT & The time limit expired without finding a solution
\end{tabular}
\end{table}

Other keywords can be sent in the e-mail message, as given in
\tabref{keywords-neos}.  The keywords specifying multiple lines of
input (\texttt{FCN,INITPT,XBOUND,COMMENTS}) use \texttt{BEGIN.keyword}
and \texttt{END.keyword} to denote the start and end of the text block.
In particular, algorithm specific options can be set with the
\texttt{OPTIONS} keyword.  
%These are documented fully in
%\partref{software}.

\begin{table}
\centering
\caption{Keywords for NEOS e-mail
submissions}\label{tab:keywords-neos}
\begin{tabular}{lll}
Keyword & Type & Meaning \\
\hline
TYPE & model type & The type of model to solve \\
SOLVER & solver name & The name of the solver to use \\
N & integer & Number of variables in the problem \\
FCN & source code & Source code for the function evaluation routine \\
INITPT & text & Source code for the initial point routine \\
XBOUND & text & Source code for the bounds subroutine \\
OPTIONS & text & Algorithm specific options to be set \\
COMMENTS & text & Comments about the model
\end{tabular}
\end{table}

Optional routines for specifying a starting point and general bounds
can be provided by a user with specifications:
\single\begin{verbatim}
     subroutine xbound(n,l,u)
     subroutine initpt(n,x)
\end{verbatim}
\double
The user places the lower bounds in \texttt{l}, the upper bounds in
\texttt{u}, and the initial starting point in \texttt{x} respectively.
In all cases, \texttt{n} is the number of variables specified for
the problem.  On input to the \texttt{xbound} routine, the lower and
upper bounds are set to plus and minus infinity respectively.
An initial starting point of zero will be used.  If the routines are
not supplied, the default implementation ($\ell = \{0\}^n$, $u = \{\infty\}^n$,
and $x = \{0\}^n$) will be used.

\begin{comment}
The Walrasian equilibrium model can also be written in functional
form as given in \figref{walrasian-neos}.  This function is not
defined when the prices are zero, so we provide an initial starting
point with the initpt routine:
\single\begin{verbatim}
     subroutine initpt(n,x)
     implicit none
     
     integer i, n
     double precision x(n)

     do i = 1,n
        x(i) = 1.0
     enddo

     return
     end
\end{verbatim}
\double
Furthermore, upper bounds can be added by sending an xbound routine.
\single\begin{verbatim}
     subroutine xbound(n,l,u)
     implicit none
 
     integer i, n
     double precision l(n), u(n)

     do i = 1, n
        l(i) = 0
        u(i) = 10
     enddo
   
     return
     end
\end{verbatim}
\double
Initially, $\ell = \{-\infty\}^n$ and $u = \{\infty\}^n$ are passed into the
routine.  Therefore, only the elements with finite bounds need to
be specified.  The numerical value used for infinity can change
depending upon system parameters.

\begin{figure}
\begin{verbatim}
       subroutine fcn(n,x,f)
       implicit none

       integer n
       double precision x(n), f(n)

       return
\end{verbatim}
\caption{Walrasian equilibrium function for NEOS, \texttt{fcn.f}}\label{fig:walrasian-neos}
\end{figure}
\end{comment}

% The current specification allows CONDOR \cite{} to be used to solve long
% running jobs.  CONDOR uses the idle workstations at the University of
% Wisconsin to solve the problems.  

\subsection{Pitfalls}

% Pitfalls - jacobian elements

Evaluation errors in the function and Jacobian cannot currently be
reported to the algorithm.  Therefore, the results can be undefined
when there is a domain violation, such as division by zero, in the
function or Jacobian.

The automatic differentiator used, ADIFOR, will process
non-differentiable functions and even those containing
discontinuities.  However, the algorithm is likely to not work well in
these situations.  Therefore, we encourage the user to specify 
models utilizing smooth functions.

The Jacobian is passed to the solver as a sparse matrix.  
A fixed amount of space is initially allocated by the solver.
Therefore, in situations where the number of nonzeros in the Jacobian
changes, we could run into serious problems.  For example, the
(sparse) Jacobian of the function $x^2$ at zero is empty.  A heuristic 
suggested by Jorge Mor\'{e} is employed in the routines that perturbs the 
initial starting point slightly and evaluates the Jacobian at this perturbed 
point to calculate an over-estimate of the number of nonzeros in the Jacobian.
This technique is not foolproof, but has worked well in practice.

\section{Summary}
This chapter developed two example complementarity problems, a transportation 
model and a Walrasian equilibrium, and demonstrated how they are communicated 
to a solver in the AMPL, GAMS, MATLAB, and NEOS environments.  These 
environments are heavily used in practical applications because of the ease 
with which a model can be specified and conveyed to different solvers.
The solver interface to the PATH 4.x and SEMI codes used to support all of 
these environments is documented in the next chapter.

We note that an earlier version of \secref{applications} and \secref{gams} 
appeared in \cite{ferris.munson:complementarity}.  The original AMPL 
interface to PATH 
was developed by David Gay, while the original GAMS interface to PATH was 
developed by Steven Dirkse and Thomas Rutherford.  Both have been revised to 
support both PATH 4.x and SEMI.  The AMPL interface to PATH 4.x can use the 
preprocessor for complementarity problems developed in \chapref{preprocessing}.
An alternate version of the 
GAMS interface was written that enables the use of the preprocessor for 
complementarity problems, but does not support MPSGE models 
\cite{rutherford:applied}.  This new interface is available in GAMS by using 
the PATHC and SEMIC solvers.  (The `C' at the end of the solver name is used to 
indicate that support is provided for constrained systems of nonlinear 
equations, which are internally reformulated as mixed complementarity 
problems.)  Finally, we remark that the MATLAB and NEOS interfaces are new 
to this thesis.

% Other environments can also be used to specify complementarity problems.
% These include 
% NEOS \cite{czyzyk.mesnier.ea:network-enabled,ferris.mesnier.ea:neos} and 
% simplified stand-alone library interfaces \cite{ferris.munson:interfaces}.  
% We refer the reader to the cited documentation for further information on 
% these environments.
