\documentclass[11pt]{amsart}
\setlength{\baselineskip}{1cm}
\renewcommand{\baselinestretch}{1.4}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{0cm}
\setlength{\textheight}{24.2cm}
\setlength{\textwidth}{16.2cm}
\setlength{\topmargin}{.0cm}
\pagestyle{empty}
\thispagestyle{empty}

\usepackage[dvips]{graphics}
\usepackage{amsmath,amsthm}
\usepackage{amsfonts}

%Russian font defs
\font\cyrs=wncyr8
\font\cyr=wncyr10 at 12pt
\font\cyrb=wncyr10 at 16pt

\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{ax}[thm]{Axiom}
\newtheorem{res}[thm]{Result}
\newtheorem{conv}[thm]{Convention}
\newtheorem{defn}[thm]{Definition}
\newtheorem{notation}[thm]{Notation}
\newtheorem{example}[thm]{Example}

\theoremstyle{remark}
\newtheorem{rem}{Remark}[section]

%\numberwithin{equation}{section}

\newcommand{\thmref}[1]{Theorem~\ref{#1}}
\newcommand{\secref}[1]{\S\ref{#1}}
\newcommand{\lemref}[1]{Lemma~\ref{#1}}

\newtheorem{Assumption}{Assumption}
\newcommand{\eg}{{\em e.g. }}
\newcommand{\etal}{{\sl et al.}}
\newcommand{\etc}{{\sl etc.}}
\newcommand{\cf}{{\sl cf }}
\newcommand{\ie}{{\em i.e. }}
\hyphenation{con-straints}
\def\N{\hbox{I \kern-.5em N}}
\def\R{\hbox{I \kern-.53em R}}


\renewcommand{\thetable}{\Roman{table}}
 
\title{Relaxation Algorithms in Finding Nash Equilibria\thanks{Research supported by VUW GSBGM and IGC.}}
 
\author{Steffan Berridge\thanks{Institute of
Statistics and Operations Research, Victoria University of Wellington.
Email Steffan.Berridge@vuw.ac.nz} \and Jacek B. Krawczyk\thanks{Commerce \& 
Administration Faculty,
Victoria of Wellington, 
PO Box 600, Wellington, New Zealand; fax +64-4-4712200. Email: 
Jacek.Krawczyk@vuw.ac.nz.}}



\begin{document} 
\begin{abstract} {\large
Relaxation algorithms provide a powerful method of finding noncooperative
equilibria in general synchronous games.  Through use of the 
Nikaido-Isoda function, the Nash solution to a broad category of constrained,
multiplayer, non-zerosum games can easily be found.  We provide solutions to 
some simple games using this procedure and extend ourselves to more difficult 
games involving coupled constraints and multiple discrete time periods using a program 
developed in Matlab.}
\end{abstract}
\maketitle
\begin{center}
\vspace{.5cm}
Working Paper \\
{\em Version 2.3} \\
\end{center}

\vspace{.25cm}
\begin{small}

\begin{itemize}
\item[{\em JEL}] Keywords: C63 (Computational Techniques),
C72 (Noncooperative Games), C87 (Economic Software), E62 (Taxation),
Q25 (Water; Air). 

\item[{\em AMS}] Keywords: 90-08 (Computational Methods), 
90A14 (Equilibrium: general theory),
90A30 (Environmental economics), 90A56 (Special types of equilibria), 90A58
(Models of real-world systems), 90A70 (Macroeconomic policy-making, taxation),
90D05 ($2$-person games), 90D06 ($n$-person games), 90D10 (Noncooperative games),
90D50 (Discrete-time games), 90D80 (Applications of game theory). 

\item[{\em Authors}'] Keywords: Computational economics; Nash normalised 
equilibrium, coupled constraints, Nikaido-Isoda function, 
open-loop Nash equilibrium. \end{itemize} 
\end{small}

\pagestyle{empty} 
\thispagestyle{empty} 
\begin{center}
\vspace{.5cm}
Research supported by VUW GSBGM and IGC. \\
\vspace{1cm}
{\em Correspondence should be addressed to:} ~Jacek B. Krawczyk.
\end{center}
\hrulefill
\vspace{.5cm}

\begin{small}
\noindent About the authors:\,{\sc Steffan Berridge.} {Institute of 
Statistics and Operations Research ({\em Masters' Student}), 
Victoria University of Wellington.
Email: \mbox{\tt Steffan Berridge@vuw.ac.nz}} \,;\\
\underline{http://www.vuw.ac.nz/~berridge} \\
{\sc Jacek B. Krawczyk.} {Commerce \&
Administration Faculty,
Victoria University of Wellington,
PO Box 600, Wellington, New Zealand; fax +64-4-4712200. Email:
\mbox{\tt Jacek.Krawczyk@vuw.ac.nz}} \,;\\
\underline{http://www.vuw.ac.nz/~jacek} 

\end{small}
\newpage

\tableofcontents
\newpage
 
\thispagestyle{empty} 
 
\bigskip 
\begin{center} 
\begin{Large} 
Relaxation Algorithms in finding Nash Equilibria
\end{Large} \end{center} 
 
\pagestyle{plain} 


\section{Introduction}
Game theory problems are renown for being difficult to solve, especially
many player, non-zero sum and dynamic games.  Computation of Nash equilibria 
using the Nikaido-Isoda function \cite{pjm1} and the relaxation algorithm 
\cite{bas1}, \cite{UR} seems to be an 
attractive possibility in that the most advanced technique required is 
minimisation of a multivariate function; a well studied topic.  
The aim of this paper is to use the relaxation algorithm to find Nash 
equilibria in games of varying complexity. 

The motivating paper \cite{UR} sets the foundation for the usage of the
relaxation algorithm to find Nash normalised equilibria, and promises
its success for a certain class of games. The contribution of the present paper
is in applyinging this idea to some specific game theory problems
and reporting on equilibrium computation 
experiments. The experiments below were conducted using customarily developed
software in the Matlab programming environment.

In the latter stages of the paper, we are able to solve coupled 
constraints games, which force
players to obey constraints based on their collective actions.  Obeying such
constraints relies on cooperation, however taxes can be enforced to ensure
that players meet these constraints (see \cite{HK}.)  
We solve a coupled constraints game
in a static and then a dynamic setting using the open loop information pattern.
As analytical solutions to such games hardly ever exist, the strength 
of this paper is in studying a new {\em computational-economics}\,
method suitable for finding equilibria in intertemporal (``open-loop'')
conflict situations.  


The organisation of the paper is as follows.  Sections \ref{defcon}
and \ref{relalg} have a tutorial character and 
provide an introduction to some elementary concepts. 
In Section \ref{simex}, we  consider some simple examples which provide an insight 
into how the algorithm works.  Section \ref{dynam}  applies the
relaxation algorithm to a dynamic game. Concluding remarks and an Appendix, which describes
the software we developed, are at the end of the paper.

All definitions, theorems and such are numbered communally and consecutively in 
each section.


\section{Definitions and concepts} \label{defcon}

\subsection{Conventions}
\begin{conv}
 In this paper we will be considering two different levels of vector; each
 player in a game will have an {\sf action}, and all players together 
 will have a {\sf collective action}.  To avoid notational confusion 
 between the two and offer the reader some comfort, we will adopt the 
 convention that a player's action will be in normal type and have a subscript
 (\eg $x_{i}\in{\R}^{m_i}$), and the collective action will be in boldface 
(\eg ${\bf x}\in{\R}^{m_1}\times\cdots\times{\R}^{m_n}$).
 In this notation, ${\bf x}=(x_1,\ldots,x_n)$.  
\end{conv}

\begin{conv}
 Political correctness is a theme of this decade, and singular
pronouns are often required to denote a non gender specific person.  We adopt
the convention that the words ``he'' and ``she'' have the common meaning 
``he or she'' where appropriate.
\end{conv}

\subsection{Nash equilibrium and the Nikaido-Isoda function} %\noindent 
This paper concerns game theory, 
we therefore present some basic definitions in this area.

\begin{defn}
 A {\sf game in normal form} is a three-tuple $\{{\mathcal N},
(X_{i})_{i\in {\mathcal N}},(\phi_{i})_{i\in {\mathcal N}}\}$ 
 where 
 ${\mathcal N}$ is the set of players, ${\mathcal N}=\{1,2,\ldots,n\}$, $X_{i}$ 
is the action space of player $i$ and
 $\phi_{i}$ is the payoff function of player $i$, $\phi_{i}:X_{1}\times
 X_{2}\times\ldots\times X_{n}\rightarrow {\R}$.
\end{defn}

\begin{notation}
 Let ${\bf x}=(x_{1},\ldots ,x_{n})$ and ${\bf y}=(y_{1},\ldots ,y_{n})$ 
 be elements of the
 collective action space $X_{1}\times\ldots\times X_{n}$.  Denote the element
 $(x_{1},\ldots,x_{i-1},y_{i},x_{i+1},\ldots,x_{n})$ by $(y_{i}\vert {\bf x})$.
\end{notation}

\begin{defn}
  Let $X\subset{\R}^{m_{1}}\times\ldots\times{\R}^{m_{n}}=
  {\R}^{m}$ be the collective action space, 
  and the function $\phi_{i}:X\rightarrow{\R}$ be
  the payoff function for player $i$, $i\in N$.  Then the point 
  ${\bf x^{\ast}}=(x_{1}^{\ast},
  \ldots,x_{n}^{\ast})$ is called the {\sf Nash equilibrium point} if, 
  for each $i$,
  \begin{equation}
    \phi_{i}({\bf x^{\ast}})=\max_{x_{i}\in{\R}^{m_{i}}}\phi_i(x_i|{\bf x^*}).
  \end{equation}
  We will use the notation 
  \begin{eqnarray*}
    {\bf x}^* = \arg{\rm equil}_{{\bf x}\in{\R}^m}\left\{(X_i)_{i\in {\mathcal N}},(\phi_{i})_{i\in {\mathcal N}}\right\}
  \end{eqnarray*}
  to denote the equilibrium point of
  the game where player $i$ has action space and payoff functions $X_i$
  and $\phi_i$ respectively.
\end{defn}

We now introduce the Nikaido-Isoda function \cite{UR}.  This function is not so
elementary in the realms of game theory, and bears some thinking about.

\begin{defn}
 Let $\phi_{i}$ be the payoff function for player $i$.  Then the {\sf Nikaido-Isoda
 function} \nolinebreak{$\Psi:(X_1\times\cdots\times X_n)\times(X_1\times\cdots
\times X_n)\rightarrow{\R}$} is defined as
 \begin{equation}
  \Psi({\bf x},{\bf y})=\sum_{i=1}^{n}[\phi_{i}(y_{i}\vert {\bf x})-\phi_{i}({\bf x})].
 \end{equation}
\end{defn}

\begin{res}
 \label{gr0}\cite{UR}
 \begin{equation}
  \Psi({\bf x},{\bf x})\equiv 0\hspace{20mm}{\bf x}\in X.
 \end{equation}
\end{res}

Let us take a brief look at the significance of this function.
Each summand can be thought of as the improvement in payoff a player will
receive by changing his action from $x_{i}$ to $y_{i}$ while all other
players continue to play according to ${\bf x}$.
The function thus represents the sum of these improvements in payoff.
Note that the {\em maximum} value this function can take, for a given ${\bf x}$,
is always nonnegative,
owing to Result \ref{gr0} above.  Also the function is everywhere
nonpositive
when either ${\bf x}$ or ${\bf y}$ is a Nash equilibrium point, since in an 
equilibrium
situation no player can make a unilateral improvement to their payoff, and
so each summand in this case can be at most zero.

From here, we reach the conclusion that when the Nikaido-Isoda
function cannot be made (significantly) positive for a given ${\bf y}$, 
we have (approximately) reached the Nash equilibrium point.  We use this
observation in constructing a termination condition for our algorithm; that
is, we choose an $\varepsilon$ such that, when $\max_{{\bf y}\in{\R}^m}
\Psi({\bf x}^s,{\bf y})<\varepsilon$, we have achieved the equilibrium to a
sufficient degree of precision.

\begin{defn}
 \label{norm_nash}
 An element ${\bf x^{\ast}}\in X$ is referred to as a {\sf Nash normalised 
equilibrium
 point} if 
 \begin{equation}
  \max_{{\bf y}\in X}\Psi({\bf x^{\ast}},{\bf y})=0.
 \end{equation}
\end{defn}

Notice that Rosen \cite{rosen65} also introduces the notion of a
normalised equilibrium when examining equilibrium point uniqueness
in concave n-person games.
However, he defines it in a different
context which is one of a game with coupled constraints.  His
definition will be presented when we deal with coupled constraints
games in Section \ref{morex}.

The two following lemmas are presented in \cite{aubin}:
\begin{lem}
 A Nash normalised equilibrium point is also a Nash equilibrium point.
\end{lem}

\begin{lem}
 A Nash equilibrium point is a Nash normalised equilibrium point if the
collective action space $X$ satisfies
 \begin{equation}
  X=X_{1}\times\ldots\times X_{n}\hspace{20mm}X_{i}\subset{\R}^{m_{i}}.
 \end{equation}
\end{lem}

An algorithm which uses the Nikaido-Isoda function to find the
Nash normalised equilibrium will be developed in the next section.
Here we note that at each point of the 
algorithm we wish to move towards a point which is
an ``improvement'' on the one that we are at.  To this end, let us put 
forward the following definition.

\begin{defn}
  The {\sf optimum response function} at point ${\bf x}$ is
  \begin{equation}
    Z({\bf x})=\arg\max_{{\bf y}\in X}\Psi({\bf x},{\bf y}).
  \end{equation}
\end{defn}

In brief terms, this function returns the set of players' actions whereby 
they all attempt to unilaterally maximise their payoffs.

It is interesting to note that we can consider the algorithm as either 
performing a static optimisation or as calculating successive actions
in a convergence to equilibrium in a real time process.
We have been considering the case where
all payoffs are known to us, in which case we can directly find the Nash
equilibrium using the relaxation algorithm.  However, if we only had
access to one player's payoff function and all players' past actions, 
then at each stage in the real time process we choose the optimum response 
for that player, assuming that the other players will play as they had
in the previous period.  In this way, convergence to the Nash normalised
equilibrium will occur as $t\rightarrow\infty$.

We now introduce some more technical definitions to be used in the main
theorem.

\begin{defn}
A function of one argument $f({\bf x})$ is {\sf weakly convex} if, $\forall{\bf x},{\bf y}\in X$
\begin{eqnarray}
 \alpha f({\bf x})+(1-\alpha)f({\bf y})\geq f(\alpha{\bf x}+(1-\alpha){\bf y})+
 \alpha(1-\alpha)r({\bf x},{\bf y}) \\
 \nonumber 0\leq \alpha\leq 1,\;\;{\rm and}\;\;
 \frac{r({\bf x},{\bf y})}{\Vert {\bf x}-{\bf y}\Vert}\rightarrow 0\;\; {\rm as}\;\;
 \Vert {\bf x}-{\bf y}\Vert\rightarrow 0\;\; \forall {\bf x}\in X
\end{eqnarray}
\end{defn}

\begin{defn}
A function of one argument $f({\bf x})$ is {\sf weakly concave} if, $\forall{\bf x},{\bf y}\in X$
\begin{eqnarray}
 \alpha f({\bf x})+(1-\alpha)f({\bf y})\leq f(\alpha{\bf x}+(1-\alpha){\bf y})+
 \alpha(1-\alpha)\mu({\bf x},{\bf y}) \\
 \nonumber 0\leq \alpha\leq 1,\;\;{\rm and}\;\;
 \frac{\mu({\bf x},{\bf y})}{\Vert {\bf x}-{\bf y}\Vert}\rightarrow 0\;\; {\rm as}\;\;
 \Vert {\bf x}-{\bf y}\Vert\rightarrow 0\;\; \forall {\bf x}\in X
\end{eqnarray}
\end{defn}

\begin{defn}
A function of two vector arguments, $f({\bf x},{\bf y})$ is referred to as {\sf weakly convex-concave}
\cite{UR} if it satisfies weak convexity with respect to its first argument and weak concavity
with respect to its second.  

That is, for fixed ${\bf z}\in X$,
\begin{eqnarray}
 \alpha f({\bf x},{\bf z})+(1-\alpha)f({\bf y},{\bf z})\geq f(\alpha{\bf x}+(1-\alpha){\bf y},{\bf z})+
 \alpha(1-\alpha)r({\bf x},{\bf y};{\bf z}) \\
 \nonumber {\bf x},{\bf y}\in X,\;\;0\leq \alpha\leq 1,\;\;{\rm and}\;\;
 \frac{r({\bf x},{\bf y};{\bf z})}{\Vert {\bf x}-{\bf y}\Vert}\rightarrow 0\;\; {\rm as}\;\;
 \Vert {\bf x}-{\bf y}\Vert\rightarrow 0\;\; \forall {\bf x}\in X
\end{eqnarray}
and
\begin{eqnarray}
 \alpha f({\bf z},{\bf x})+(1-\alpha) f({\bf z},{\bf y})\leq f({\bf z},\alpha {\bf x}+(1-\alpha){\bf y})+
 \alpha(1-\alpha)\mu({\bf x},{\bf y};{\bf z}) \\
 \nonumber {\bf x},{\bf y}\in X,\;\;0\leq \alpha\leq 1,\;\;{\rm and}\;\;
 \frac{\mu({\bf x},{\bf y};{\bf z})}{\Vert {\bf x}-{\bf y}\Vert}\rightarrow 0\;\;{\rm as}\;\;
 \Vert {\bf x}-{\bf y}\Vert\rightarrow 0\;\; \forall {\bf x}\in X
\end{eqnarray}
$r({\bf x},{\bf y};{\bf z})$ and $\mu({\bf x},{\bf y};{\bf z})$ are referred to as the residual terms.
\end{defn}

The notions of weak convexity and concavity are weakenings of strict convexity and concavity.  The residual
terms, to be chosen at will, ensure that there are many concave functions which are weakly convex and many
convex functions which are weakly concave.

In fact, the family of weakly convex-concave functions includes the family of smooth
functions\footnote{Recall that a function is smooth iff its derivatives of all orders
are continuous.} \cite{UR}, and so this condition is a minimal restriction for practical purposes.

We now present an elementary example to illustrate the above ideas.

\noindent \underline{Example}: 
{\em A convex function which is weakly concave}

\noindent Consider the function 
$f:{\R}\rightarrow {\R}$ defined by
\begin{eqnarray*}
 f(x)=x^2.
\end{eqnarray*}
This function is convex in all senses.  To show that it is weakly concave, we must find an $r(x,y)$ such that, for all $x,y\in{\R}$ and $\alpha\in[0,1]$,
\begin{eqnarray*}
 \alpha f(x)+(1-\alpha)f(y)\leq f(\alpha x+(1-\alpha)y)+\alpha(1-\alpha)r(x,y)
\end{eqnarray*}
That is if and only if
\begin{eqnarray*}
 & & \alpha x^2+(1-\alpha)y^2\leq (\alpha x+(1-\alpha)y)^2+\alpha(1-\alpha)r(x,y) \\
 & \Leftrightarrow \hspace{20mm} & \alpha x^2+(1-\alpha)y^2\leq \alpha^2x^2+(1-\alpha)^2y^2+2\alpha(1-\alpha)xy+\alpha(1-\alpha)r(x,y) \\
 & \Leftrightarrow \hspace{20mm} & \alpha(1-\alpha)x^2+\alpha(1-\alpha)y^2-2\alpha(1-\alpha)xy\leq \alpha(1-\alpha)r(x,y) \\
 & \Leftrightarrow \hspace{20mm} & (x-y)^2\leq r(x,y) \\
\end{eqnarray*}
So it is sufficient to select $r(x,y)=(x-y)^2$.  Also
\begin{eqnarray*}
 \frac{r(x,y)}{\|x-y\|} = \frac{(x-y)^2}{|x-y|} =  |x-y| \rightarrow & 0 \hspace{20mm}\mbox{as $|x-y|\rightarrow 0$.}
\end{eqnarray*}
So we conclude that $f(x)=x^2$ is weakly concave.  (However, as this is a one variable function, it cannot be weakly convex-concave.) ~~~~~ \vspace{1.5cm}$\diamond$

The function $r({\bf x},{\bf y};{\bf z})$ was introduced with the concept of weak convex-concavity.
In the case of $\Psi({\bf x},{\bf y})$ being a twice continuously differentiable function with
respect to both arguments on $X\times X$, the residual terms satisfy
\begin{equation}
r({\bf x},{\bf y};{\bf y})=\frac{1}{2}\langle A({\bf x},{\bf x})({\bf x}-{\bf y}),{\bf x}-{\bf y}\rangle+o_{1}(\Vert {\bf x}-{\bf y}\Vert^{2})
\end{equation}
and
\begin{equation}
\mu({\bf y},{\bf x};{\bf x})=\frac{1}{2}\langle B({\bf x},{\bf x})({\bf x}-{\bf y}),{\bf x}-{\bf y}\rangle+o_{2}(\Vert {\bf x}-{\bf y}\Vert^{2})
\end{equation}
where $A({\bf x},{\bf x})=\Psi_{{\bf x}{\bf x}}({\bf x},{\bf y})|_{{\bf y}={\bf x}}$ is the Hessian of the Nikaido Isoda function with
respect to the first argument and $B({\bf x},{\bf x})=\Psi_{{\bf y}{\bf y}}({\bf x},{\bf y})|_{{\bf y}={\bf x}}$ is the Hessian of the Nikaido-Isoda function with respect to the second argument, both evaluated at ${\bf y}={\bf x}$.

To prove the inequality of Condition (5) of the following Theorem \ref{condit} under the assumption that $\Psi({\bf x},{\bf y})$
is twice continuously differentiable, it suffices 
to show that $Q({\bf x},{\bf x})=A({\bf x},{\bf x})-B({\bf x},{\bf x})$ is strictly positive.

\begin{defn}
The set $X$ is convex if, for every two points ${\bf x}$ and ${\bf y}$ in $X$, the point 
\newline \mbox{$\alpha{\bf x}+(1-\alpha){\bf y}$} 
is also in $X$ for $0<\alpha<1$.
\end{defn}

\begin{defn}
The set $X$ is compact if every Cauchy sequence has a limit in $X$.
In a finite dimensional space this is equivalent to $X$ being closed and bounded.  
\end{defn}



\section{The relaxation algorithm} \label{relalg}

\subsection{Statement of the algorithm}%\noindent 
Suppose we wish to find the Nash equilibrium for some game
and we have some initial estimate of it, say ${\bf x_{0}}$.
The relaxation algorithm we use is
\begin{eqnarray}
 {\bf x}^{s+1}=(1-\alpha_{s}){\bf x}^{s}+\alpha_{s}Z({\bf x}^{s})\hspace{20mm}&0<\alpha_{s}\leq 1\\
                                               \nonumber &s=0,1,2,\ldots 
\end{eqnarray}

The iterate at step $s+1$ is constructed by a weighted average of the
improvement point $Z({\bf x}^s)$ and the current point ${\bf x}^s$.
This averaging ensures convergence of the algorithm under certain
conditions, as stated in the following Theorem \ref{condit}.

By taking sufficiently many iterations of the algorithm, it is our aim to 
determine the Nash 
equilibrium ${\bf x^{\ast}}$ with arbitrary precision.  Indeed Theorem 
\ref{condit} in the following subsection provides conditions for the 
convergence of this relaxation algorithm.

\subsection{Conditions for existence of a Nash equilibrium and convergence of
relaxation algorithm} %\noindent 
The following theorem states the conditions of convergence for the 
relaxation algorithm.  The conditions may look rather restrictive, but in fact a large
class of games satisfy them.

\begin{thm}\label{condit}\cite{UR}
There exists a unique Nash equilibrium point to which the algorithm converges if
:
\begin{enumerate}
 \item $X$ is a convex compact subset of ${\R}^{m}$,
 \item the Nikaido-Isoda function $\Psi:\;X\times X\rightarrow{\R}$ 
  is a weakly convex-concave
  function and $\Psi({\bf x},{\bf x})=0\;\;{\bf x}\in X$,
 \item the optimum response function $Z({\bf x})$ is single valued and continuous on
  $X$,
 \item the residual term $r({\bf x},{\bf y};{\bf z})$ is uniformly continuous 
  on $X$ wrt ${\bf z}$,
 \item the residual terms satisfy
  \begin{equation}
   r({\bf x},{\bf y};{\bf y})-\mu({\bf y},{\bf x};{\bf x})\geq\beta(\Vert {\bf x}-{\bf y} \Vert)\;\;\;\;{\bf x},{\bf y}\in X
  \end{equation}
  where $\beta$ is a strictly increasing function,
 \item the relaxation parameters $\alpha_{s}$ satisfy
  \begin{enumerate}
   \item $\alpha_{s}>0$,
   \item $\sum_{s=0}^{\infty}\alpha_{s}=\infty$,
   \item $\alpha_{s}\rightarrow 0 \;\;as \;\;s \rightarrow\infty$.
  \end{enumerate}
\end{enumerate}
\end{thm}

\subsection{Step size optimisation} %\noindent
In order for the algorithm to converge, we may choose any sequence
$(\alpha_s)$ satisfying the final condition of Theorem \ref{condit}.  However,
it is of computational importance to attempt to optimise the convergence rate.

Suitable step-sizes may be obtained by trial and error, and we have found that
using a constant step of $\alpha_s\equiv 0,5$ leads to a quick convergence in
most of our experiments. In this case, we can think of the $\alpha_s$ as being
constant until our convergence conditions are reached, and thereafter decaying
with factors $\dfrac{1}{2},\dfrac{1}{3},\dfrac{1}{4},\ldots$.

We suggest here a method for step-size optimisation.

\begin{defn}
Suppose that we reach ${\bf x}^s$ at iteration $s$.  Then $\alpha_s^*$ is
{\sf one-step optimal} if it minimises the optimum response function at
${\bf x}^{s+1}$.  That is, if
\begin{equation}
  \alpha_s^* = \arg\min_{\alpha\in[0,1)}\left\{\max_{{\bf y}\in X}\Psi({\bf x}^{s+1}(\alpha),{\bf y})\right\}.
\end{equation}
(Recall that ${\bf x}^{s+1}$ depends on $\alpha$.)
\end{defn}
This method is intuitively appealing. Indeed, we are trying to force 
the maximum of the Nikaido-Isoda function that is an equilibrium condition.

Of course, we could also use two-step optimisation, or indeed
$n$-step.  However any number greater than one would be
computationally expensive.
We have found that one-step optimising the step sizes leads to fewer
iterations, but that each step takes longer to complete than when 
using constant step sizes.

When we refer to optimal step sizes throughout the paper, we 
generally mean one-step optimal, however in some cases where
optmisation is very simple, other more elementary
methods have been used.


\section{Some simple examples with static games} \label{simex}
In this section we will consider some simple games in order to gain an appreciation
for how the algorithm works in conjunction with the Nikaido-Isoda function.
\subsection{Simple two player games}
\subsubsection{A two player game with identical payoff functions}
\label{ex1}
 Consider a two player game where both players are maximisers, the action
space for each player is ${\R}$ and the payoff functions of players $1$ and $2$ are identical.
Consider the following simple payoff function on the region $X=\{(x_{1},x_{2}):
-10\leq x_{1},x_{2}\leq 10\}$.
\begin{equation}
\phi_{i}({\bf x})=-x_{1}^{2}-x_{2}^{2}\;\;\;\;\;\; i=1,2.
\end{equation}

In this case the Nikaido Isoda function is:
\begin{equation}
\Psi({\bf x},{\bf y})=x_{1}^{2}+x_{2}^{2}-y_{1}^{2}-y_{2}^{2}.
\end{equation}

Notice in this case that assumptions 1-5 of Theorem \ref{condit} are satisfied.

The Nash equilibrium for this game is clearly the maximum of the common
payoff function, that is ${\bf x}=(0,0)$.

The optimum response function $\displaystyle{Z({\bf x})=\arg\max_{{\bf y}\in X} \Psi({\bf x},{\bf y})}$
is determined
by simple calculus to be $(0,0)$ for any ${\bf x}$.  Hence if we were to choose
$\alpha_{0}$ to be $\varepsilon$ where $0<\varepsilon\ll 1$, it doesn't matter what 
values we assign to $\alpha_{1},\alpha_{2},\ldots$ since we can choose
$\varepsilon$ such that ${\bf x}^{1}$ satisfies our termination conditions, namely
that the Nikaido-Isoda function 
\begin{eqnarray*}
   {\bf x}^1    =      \alpha_{0}{\bf x}^{0}+(1-\alpha_{0})Z({\bf x}^{0})
                =      \varepsilon {\bf x}^{0}+(1-\varepsilon)(0,0)
                =       \varepsilon {\bf x}^{0}
                \simeq  (0,0).
\end{eqnarray*}

The convergence in the action space is displayed in Figure \ref{ex_1_fig}
for starting guess $(10,5)$.

\begin{figure}[ht]
  \begin{center}
    \leavevmode
    \scalebox{0.4}{\includegraphics{exampl_1.eps}}
    \caption{Convergence of the example in Section \ref{ex1} with starting guess (10,5)}
    \label{ex_1_fig}
  \end{center}
\end{figure}

\subsubsection{A slightly less trivial example}
\label{ex2}
Consider the situation of Section \ref{ex1}, but let the common payoff function now be
\begin{equation}
\phi_{i}({\bf x})=-\frac{(x_{1}+x_{2})^{2}}{4} - \frac{(x_{1}-x_{2})^{2}}{9}.
\end{equation}

The Nikaido-Isoda function in this case is

\begin{eqnarray}
\Psi({\bf x},{\bf y}) =&2\displaystyle{\left\{\frac{(x_{1}+x_{2})^{2}}{4}+\frac{(x_{1}-x_{2})^{2}}{9}\right\}}\\ 
\nonumber  &\displaystyle{-\left\{\frac{(y_{1}+x_{2})^{2}}{4}+\frac{(y_{1}-x_{2})^{2}}{9}\right\}}\\
\nonumber  &\displaystyle{-\left\{\frac{(x_{1}+y_{2})^{2}}{4}+\frac{(x_{1}-y_{2})^{2}}{9}\right\}}.
\end{eqnarray}
As in Section \ref{ex1}, the Nash equilibrium for this game is ${\bf x}=(0,0)$.

The optimum response function is calculated to be $Z({\bf x})=\displaystyle{-\frac{5}{13}}(x_{2},x_{1})$.
In this case it is also relatively simple to see how to optimise the $\alpha_s$; 
since both players have the same payoff function, the optimal $\alpha_{s}$ is the 
one which optimises the payoff function 
\begin{equation}
\phi_{i}({\bf x^{s+1}})=\phi_i\left(\alpha_{s} {\bf x^{s}}+(1-\alpha_{s})Z({\bf x^{s}})\right)\;\;\;\;\;
0<\alpha_{s}\leq 1.
\end{equation}

Our Matlab program (see Appendix) gives the results of Table \ref{ex2_conv}, with starting 
guess ${\bf x}=(10,5)$, and using optimised step sizes $\alpha_s$.

\begin{table}
\begin{center}
\begin{tabular}{|c|cc|}
\hline
Iteration (s) & ${\bf x_{s}}$ & $\alpha_{s}$\\
\hline
0 & (10,5)& 0.7309\\
1 & (1.2849,-1.4668) & 1.0000\\
2 & (0.5639,-0.4942) & 1.0000\\
3 & (0.1901,-0.2169) & 1.0000\\
4 & (0.0834,-0.0731) & 1.0000\\
5 & (0.0281,-0.0321) & 1.0000\\
6 & (0.0123,-0.0108) & 0.5001\\
7 & (0.0082,-0.0078) &       \\
\hline
\end{tabular}
\caption{Convergence of the example in Section \ref{ex2}}
\label{ex2_conv}
\end{center}
\end{table}

The equilibrium point is still ${\bf x}=(0,0)$, but the algorithm does not 
converge directly to it.  It rather comes in line with the equilibrium first, 
as seen in Figure \ref{ex_2_fig}.  We can now make a comparison between
the optimised and nonoptimised $\alpha_s$.  The first graph, with the
optimisation performed, shows a much quicker convergence.  In contrast
the second one, which has $\alpha_s\equiv 0.5$ shows a smoother but much
slower convergence.  The third shows step sizes of $\alpha_s\equiv 1$ (that
is, a non-relaxed algorithm.) We would clearly prefer to use optimised 
step sizes if this could be achieved easily.

Note that Theorem \ref{condit} gives no guarantee that the algorithm
will converge at all in the second or third cases, despite our success.

\begin{figure}[htp]
\begin{center}
\scalebox{0.25}{\includegraphics{exam_2_1.eps}}
\scalebox{0.25}{\includegraphics{exam_2_2.eps}}
\scalebox{0.25}{\includegraphics{exam_2_3.eps}}
\caption{Convergence of the example in Section \ref{ex2} with starting guess (10,5) using
 (i) optimised $\alpha$s, (ii) $\alpha\equiv 0.5$ and (iii) $\alpha\equiv 1$.
 The isolines are those of the identical payoff function}
\label{ex_2_fig}
\end{center}
\end{figure}


\subsubsection{The quantity setting duopoly } \label{example_qsd} 
\nopagebreak
In this model, the players' payoff functions 
are not necessarily the same.  Consider a
situation where two firms sell an identical product on the same market
\cite{ECON407}.  Each
firm will want to choose its production rates in such a way as to maximise
its respective profits.  Let $x_{i}$ be the production of firm $i$ and let
$\alpha$, $\lambda$ and $\rho$ be constants.  The market price is:
\begin{equation}
p({\bf x})=\alpha-\rho(x_{1}+x_{2})
\end{equation}
and the profit made by firm $i$, is:
\begin{eqnarray}
\nonumber \phi_{i}({\bf x}) &=& p({\bf x})x_{i}-\lambda x_{i}\\
                            &=& [\alpha-\lambda-\rho(x_{1}+x_{2})]x_{i} .
\end{eqnarray}
The Nikaido-Isoda function in this case is
\begin{eqnarray}
\Psi({\bf x},{\bf y})=&[\alpha-\lambda-\rho(y_{1}+x_{2})]y_{1}-[\alpha-\lambda-\rho(x_{1}+x_{2})]x_{1}\\
\nonumber &+[\alpha-\lambda-\rho(x_{1}+y_{2})]y_{2}-[\alpha-\lambda-\rho(x_{1}+x_{2})]x_{2}
\end{eqnarray}
leading to an optimum response function of
\begin{equation}
Z({\bf x})=\frac{\alpha-\lambda}{2\rho}\left(1,1\right)-\frac{1}{2}\left(x_{2},x_{1}\right).
\end{equation}

It is a classical result that the Nash 
equilibrium in this game is $\displaystyle{x_{i}^{N}=\frac{\alpha-\lambda}{3\rho}}$ 
with corresponding payoff 
$\displaystyle{\phi_{i}(x_{i}^{N})=\frac{(\alpha-\lambda)^{2}}{9\rho}}$
(see for example \cite{ECON407}).  To show
the convergence of the algorithm in this case let us assign values to the parameters.
Let $\alpha=20$, $\lambda=4$ and $\rho=1$, then $x^{N}=\displaystyle{\left(\frac{16}{3},\frac{16}{3}\right)}$.

The convergence of the algorithm for the quantity setting duopoly  
is displayed in Figure \ref{fig:qsd_conv}.
Note that $\alpha_s\equiv 0.5$ here.

\begin{figure}[ht]
  \begin{center}
     \leavevmode
     \scalebox{0.4}{\includegraphics{qsd.eps}}
    \caption{Convergence of the example in Section \ref{example_qsd}}
    \label{fig:qsd_conv}
  \end{center}
\end{figure}


\subsection{A static game with coupled constraints} %\noindent 
\label{morex}
We will now use our method to find the equilibrium of a game with 
coupled constraints.  This means that the players' action set is 
now a general convex set $X\subset{\R}^n$ rather than as
previously where we were confined to $X=X_1\times\cdots\times X_n\subset
{\R}^{m_1}\times\cdots\times{\R}^{m_n}$,
each separate $X_i$ being convex.


\subsubsection{Existence of the solution to a coupled constraints game}
\label{sec:rosen_exist}

The theorems and definitions to follow are presented and
proved in Rosen \cite{rosen65}, and
lead to the existence and uniqueness of an equilibrium in the
game we introduce in Section \ref{ex4}.

Consider the solution to a game where the $n$ players have 
payoff functions $(\phi_i)_{i=1\ldots n}$,
\begin{equation}
  \label{eqn:river_tosolve}
  \arg{\rm equil}_{{\bf x}\in X}\left\{\phi_1({\bf x}),\ldots,\phi_n({\bf x})\right\}
\end{equation}
such that $X$ is a convex subset of ${\R}^m$.  This 
shall be called a Rosen's coupled constraints game.
In the special case where $X=X_1\times\cdots\times X_n$,
the game is said to have uncoupled constraints.

\begin{thm}
An equilibrium point exists for every concave $n$-person game.
\end{thm}

So we know that if each player has a payoff function which is
concave with respect to his own action while other players'
actions remain fixed, then the game must have
at least one Nash equilibrium.

The conditions for uniqueness in the case of uncoupled
constraints rely on the concept of
diagonal strict concavity.  Put loosely, a game with
this property is one in which each player has more control
over his payoff than the other players have over it.

\begin{defn}
The function $f({\bf x},r)=\sum_{i=1}^nr_i\phi_i({\bf x})$ 
is called {\sf diagonally strictly concave} for \mbox{${\bf x}\in X$} 
and fixed $r\in {\R}_+^n$ if for every ${\bf x}^0,{\bf x}^1\in X$ we have
\begin{equation}
  \label{eqn:river_sdc}
  ({\bf x}^1-{\bf x}^0)'g({\bf x}^0,r)+({\bf x}^0-{\bf x}^1)'g({\bf x}^1,r)>0
\end{equation}

where $g({\bf x},r)$ is the pseudogradient of $f$

\begin{equation}
  \label{eqn:river_pg}
  g({\bf x},r) = \left[\begin{array}{c} r_1\nabla_1\phi_1({\bf x})\\ 
                 \vdots \\ r_n\nabla_n\phi_n({\bf x}) \end{array}\right].
\end{equation}
\end{defn}

\begin{thm}
In a game with uncoupled constraints,
if the joint payoff function $f({\bf x},r)=\sum_{i=1}^nr_i\phi_i({\bf x})$
is diagonally strictly concave for some $r>0$, then
there exists a unique Nash equilibrium.
\end{thm}

When the constraints are coupled, there are no such guarantees,
and we must define a special type of equilibrium.

We introduce the function $h:{\R}^m\rightarrow{\R}^n$
to describe the constraint set $R$,
where $m$ is the dimension of the collective action space and $n$
is the number of players.  Each component of $h$ is a concave function
of ${\bf x}$, and is defined such that
\begin{equation}
R=\{{\bf x} : h({\bf x})\geq 0\}.
\end{equation}
Denote the Lagrange multiplier vector for player $i$ by $u_i^0$, then
${\bf x}^0\in R$ is an equilibrium point if and only if it 
satisfies the Kuhn-Tucker conditions:
\begin{eqnarray}
  h({\bf x}^0) &\geq& 0 \\
  (u_i^0)^Th({\bf x}^0) &=& 0 \\
  \phi_i({\bf x}^0) &\geq& \phi_i(y_i|{\bf x}^0)+(u_i^0)^Th(y_i|{\bf x}^0).
\end{eqnarray}
\begin{defn}
The equilibrium point ${\bf x}^0$ is a {\sf Rosen normalised equilibrium point}
if, for some vector $r>0$ and constant $u^0\geq 0$,
\begin{equation}
  u_i^0=\dfrac{u^0}{r_i}
\end{equation}
for each $i$.
\end{defn}

\begin{thm}
There exists a normalised equilibrium point to a concave $n$-person game
for every $r>0$.
\end{thm}

\begin{thm}
Let $Q$ be a convex subset of ${\R}_+^n$.
If $f({\bf x},r)$ is diagonally strictly concave for every $r\in Q$,
then for each $r\in Q$ there is a unique normalised equilibrium point.
\end{thm}

This tells us that, if we can show that our game is diagonally strictly
concave, then we can find a unique normalised equilibrium.

We will not prove the applicability of Theorem \ref{condit} to the 
coupled constraints game.  We suggest this could be done if the
components of the Nikaido-Isoda function were the players' 
Lagrangians rather than their payoffs $\phi_i$.  Our numerical experiments 
do seem to indicate that the diagonal strict concavity, which guarantees a 
unique solution to this class of games, is a sufficient condition for the 
relaxation algorithm to find an equilibrium here.


\subsubsection{Formulation of a coupled constraints game}
\label{ex4}
{\em River Basin Pollution, \cite{HK}.}
In this game, we consider three players $j=1,2,3$ located along a 
river.  Each agent is engaged in an economic activity at 
a chosen level $x_j$, but
the players must meet environmental conditions set by a local authority.

Pollutants may be expelled into the river, where they disperse.  Two 
monitoring stations $\ell=1,2$ are located along the river, at which the 
local authority has set maximum pollutant concentration levels.

The revenue for player $j$ is
\begin{equation}
{\mathcal R}_j({\bf x}) = [d_1-d_2(x_1+x_2+x_3)]x_j
\end{equation}
with expenditure
\begin{equation}
{\mathcal F}_j({\bf x}) = (c_{1j}+c_{2j}x_j)x_j.
\end{equation}
Thus the net profit for player $j$ is
\begin{eqnarray}
\nonumber \phi_j({\bf x}) & = & {\mathcal R}_j({\bf x}) - {\mathcal F}_j({\bf x}) \\
\label{eqn:river_phi}
                          & = & [d_1-d_2(x_1+x_2+x_3)-c_{1j}-c_{2j}x_j]x_j.
\end{eqnarray}
The constraint imposed by the local authority at location $l$ is
\begin{equation}
\label{eqn:river_ql}
q_\ell({\bf x}) = \sum_{j=1}^3\alpha_{j \ell}e_jx_j \leq 100, ~~\ell = 1,2.
\end{equation}

The economic constants $d_1$ and $d_2$ determine the inverse demand
law and are set to 3.0 and 0.01 respectively.  
The values for constants $c_{1j}$ and $c_{2j}$ are
given values in Table \ref{ex3_consts}.

The $\alpha_{j \ell}$ are the decay-transportation coefficients from
player $j$ to location $l$, and $e_j$ is the emission coefficient of player 
$j$, also given in Table \ref{ex3_consts}.

\begin{table}
\begin{center}
\begin{tabular}{|c|ccccc|}
\hline
Player $j$ & $c_{1j}$ & $c_{2j}$ & $e_j$ & $\alpha_{j1}$ & $\alpha_{j2}$ \\
\hline
1 & 0.10 & 0.01 & 0.50 & 6.5 & 4.583 \\ 
2 & 0.12 & 0.05 & 0.25 & 5.0 & 6.250 \\
3 & 0.15 & 0.01 & 0.75 & 5.5 & 3.750 \\
\hline
\end{tabular}
\caption{Constants for the River Basin Pollution game}
\label{ex3_consts}
\end{center}
\end{table}


In this case, 
we can check conditions (\ref{eqn:river_sdc})  to show that problem 
(\ref{eqn:river_tosolve}) is a diagonally strictly concave one,
and so a unique equilibrium point exists.  
We can thus compute it for a given set of weights $r$ using our relaxation
algorithm, where the optimum response function $Z({\bf x}^s)$ is restricted 
to lie within the feasibility set $R$.  
Notice also that the region defined by Equation (\ref{eqn:river_ql}) is convex.


\subsubsection{Solution to the River Basin Pollution game} 
\label{RBPsol}
Here, we build the Nikaido-Isoda function.  According to Section 
\ref{sec:rosen_exist}, its components should be the players' 
corresponding Lagrangian payoffs.  However, we will formulate the 
Nikaido-Isoda function from the unconstrained game, letting the 
constraints be handled within the Matlab procedure.

To this end, the Nikaido-Isoda function is
\begin{eqnarray*}
\Psi({\bf x},{\bf(y)}) & = & \sum_{j=1}^3(\phi_j(y_i|{\bf x})-\phi_j({\bf x})) \\
                      & = & [d_1-d_2(y_1+x_2+x_3)-c_{11}-c_{21}y_1]y_1 \\
                      &   & + [d_1-d_2(x_1+y_2+x_3)-c_{12}-c_{22}y_2]y_2 \\
                      &   & + [d_1-d_2(x_1+x_2+y_3)-c_{13}-c_{23}y_3]y_3.
\end{eqnarray*}

We used a starting guess of ${\bf x}=(0,0,0)$ and $\alpha\equiv 0.5$ in 
our Matlab program.  The convergence is shown in Table \ref{fig:ex4.conv}.

\begin{table}
  \begin{center}
    \begin{tabular}{|c|ccc|c|}
      \hline
      Iteration (s) & $x_1^s$ & $x_2^s$ & $x_3^s$ & $\alpha_{s}$\\
      \hline
      0 & 0 & 0 & 0 & 0.5 \\
      1 & 9.68 & 8.59 & 1.90 & 0.5 \\
      2 & 14.85 & 12.62 & 2.655 & 0.5 \\
      3 & 17.65 & 14.49 & 2.913 & 0.5 \\
      4 & 19.18 & 15.35 & 2.961 & 0.5 \\
      5 & 20.03 & 15.73 & 2.934 & 0.5 \\
      \vdots & \vdots & \vdots & \vdots & \vdots \\
      10 & 21.07 & 16.03 & 2.762 & 0.5 \\
      \vdots & \vdots & \vdots & \vdots & \vdots \\
      20 & 21.14 & 16.03 & 2.728 & 0.5 \\
      \hline
    \end{tabular}
    \caption{Convergence in the river basin pollution game}
    \label{fig:ex4.conv}
  \end{center}
\end{table}

The convergence is displayed as a line in the 3D action space in Figure \ref{fig:ex4.3D}.  
This game was also solved in \cite{HK} using Rosen's algorithm, and was found to have
equilibrium ${\bf x}=(21.149,16.028,2.722)$, giving net profits ${\bf z}=
(48.42,26.92,6.60)$. The first constraint is active, \ie  ~
$q_1({\bf x})=100$; the second constraint is inactive ($q_2({\bf x})=81.17$).
Our solution of ${\bf x}=(21.14,16.03,2.728)$ is within ${\pm 0.01}$ of the 
solution achieved in \cite{HK}.

\begin{figure}[htp]
  \begin{center}
    \scalebox{0.4}{\includegraphics{exampl_4.eps}}
    \caption{Convergence in Section \ref{ex4} with starting guess (0,0,0), $\alpha\equiv 0.5$}
    \label{fig:ex4.3D}
  \end{center}
\end{figure}

\begin{figure}[htp]
  \begin{center}
    \scalebox{0.25}{\includegraphics{river_p1.eps}}
    \scalebox{0.25}{\includegraphics{river_p2.eps}}
    \scalebox{0.25}{\includegraphics{river_p3.eps}}
    \caption{Progression of payoffs of Figure \ref{fig:ex4.3D}}
    \label{ex_4_payoffs}
  \end{center}
\end{figure}
\clearpage
\subsubsection{Applying the Pigouvian taxes}
Now that the Nash normalised equilibrium has been found,
we can force the players to obey it by applying Pigouvian
taxes.  In this way we create a new, unconstrained
game whose equilibrium is the same as that of
the original.

For each active constraint, we place a tax 
on each player of the amount
\begin{equation}
  T_\ell({\bf x}) = \max(0,\lambda_\ell g_\ell({\bf x}))
\end{equation}
where $\lambda_\ell$ is the Lagrange multiplier
for constraint $\ell$ at the equilibrium
${\bf x}={\bf x}^*$.  

For our game, recall that $g_\ell=q_\ell$ is the
constraint at location $\ell$, $\ell=1,2$, and that
the only the first constraint is active.

Thus, in the river basin pollution game with
the parameters given,
the payoff function for player $j$ becomes
\begin{eqnarray}
\nonumber \phi_j^*({\bf x}) & = & R_j({\bf x}) - F_j({\bf x}) - \sum_\ell T_\ell({\bf x})\\
 & = & [d_1-d_2(x_1+x_2+x_3)-c_{1j}-c_{2j}x_j]x_j  
- \lambda_1\max\left(0,\sum_{j=1}^3\alpha_{j \ell}e_jx_j-100\right). 
\label{eqn:river_phi_l}
\end{eqnarray}

The coputations of  Section \ref{RBPsol} gave us the Lagrange multiplier value 
for the active constraint \, $\lambda_1=0.5774$.
Thus, applying the taxes as above leads to the
set of payoff functions $\phi_j^*$ described by Equation 
(\ref{eqn:river_phi_l}), which has the unconstrained
Nash equilibrium 
${\bf x}^{**}=\arg{\rm equil}_{{\bf x}\in\mathbb{R}^m}(\phi_1,\phi_2,\phi_3,q_1,q_2)={\bf x}^*$.
That is, the new (unconstrained) Nash equilibrium is equal
to the old (constrained) Nash normalised equilibrium.

We report that checking  numerically the 
validity of the above is difficult (using the Matlab {\tt constr} routine,
see Appendix) due to the derivative
discontinuities introduced by the $\max$ function.
However, cross-sectional graphs of the payoff
functions strongly vindicate the results (see Figure
\ref{ex_4_pigouvian}).
	
\begin{figure}[htp]
  \begin{center}
    \scalebox{0.25}{\includegraphics{river_l1.eps}}
    \scalebox{0.25}{\includegraphics{river_l2.eps}}
    \scalebox{0.25}{\includegraphics{river_l3.eps}}
    \caption{Payoff functions for players 1,2,3 with Pigouvian taxes applied}
    \label{ex_4_pigouvian}
  \end{center}
\end{figure}


\section{Finding equilibria in dynamic games} 
\label{dynam}

\subsection{Application to dynamic games}%\noindent 
The relaxation algorithm allows us to find solutions to a large number
of dynamic games with coupled constraints, so long as those games can 
be decomposed to a game where the present value of the payoff for each 
player is a function of the future collective actions and the initial
conditions of the game.  This is the case for deterministic games with 
the open-loop information pattern.

In particular, any finite horizon game where the coupled constraints 
and payoffs at time $t$ depend only on the players' states at that time,
and the states can be chosen by players through their actions, is
solvable by this method.

Let us illustrate this with an example where the game of Section
\ref{ex4} is played over a finite number of periods.

\subsection{Formulation of the game}
\label{subsec:rived}
{\em The River Basin Pollution game played over $T$ periods.}
Suppose that the agents of Section \ref{ex4} now repeat the
River Pollution game over $T$ discrete time periods.  They have an option to
invest in their products and capacities, which can  depreciate.

The values for all constants
will remain the same as in Section \ref{ex4} over the period
of the game.

{\em State equations. }
Let $x_i^{(t)}$ denote the fully utilised production capacity 
that player $i$ has at time $t$, governed by the state equation
\begin{equation}
  \label{eqn:rived_state}
  x_i^{(t+1)} = (1-\mu_i)x_i^{(t)} + u_i^{(t)}
\end{equation}
where $u_i^{(t)}$ is the investment that player $i$ makes at time
$t$, and $\mu_i$ is the depreciation that applies to player $i$,
and $x_i^{(0)}$ is fixed, $i=1,2,3$, $t=0,\ldots,T$.

Denote the collective action and state respectively by
\begin{eqnarray}
  {\bf u} &=& \left({\bf u}_1,{\bf u}_2,{\bf u}_3\right) \hspace{20mm}\mbox{and}\\
  {\bf x} &=& \left({\bf x}_1,{\bf x}_2,{\bf x}_3\right)
\end{eqnarray}
where ${\bf u}_i=(u_i^{(1)},\ldots,u_i^{(T)})$ is the vector
of all actions of player $i$ and ${\bf x}_i=(x_i^{(1)},\ldots,x_i^{(T)})$
is the collective-time state vector of player $i$.
Denote the collective action and state at time $t$ by
\begin{eqnarray}
  {\bf u}^{(t)} &=& \left(u_1^{(t)},u_2^{(t)},u_3^{(t)}\right) \hspace{20mm}\mbox{and}\\
  {\bf x}^{(t)} &=& \left(x_1^{(t)},x_2^{(t)},x_3^{(t)}\right)
\end{eqnarray}
respectively.

{\em Information structure. }
Suppose that the game is open-loop.  That is, the players cannot
observe the state or capacity of another player other than
at $t=0$.  This implies
that the game can have an open-loop Nash equilibrium\footnote{We consider this
solution concept valid for ``short'' horizon dynamic games, where the 
players are in the process of learning each others' behaviour. In that
case, the expression $\rho^Tkx_i^{(T)}$ in Equation (\ref{eqn:rived_payoff}) 
can represent the payoff in 
a game which is played from $T$ onwards, under a different information
pattern.}.

{\em Payoff function. }
The present value payoff for player $i$ is
\begin{equation}
  \label{eqn:rived_payoff}
  \Phi_i({\bf u}) = \sum_{t=0}^{T-1}\rho^t\left(\phi_i({\bf x}^{(t)})
                    - c_{3i}(u_i^{(t)})^2\right) 
+ \rho^Tk_ix_i^{(T)}\hspace{20mm}i=1,2,3
\end{equation}
where $\phi_i({\bf x}^{(t)})$ is defined through Equation 
(\ref{eqn:river_phi}), $c_{3i}(u_i^{(t)})^2$
is the investment cost at time $t$ and $\rho^Tk_ix_i^{(T)}$ 
represents the scrap value and time $T$ for player $i$,
$k_i\geq 0$, being the scrap value constant for player $i$.

{\em Constraints. }
The constraints of Equation (\ref{eqn:river_ql}) for $\ell=1,2$ are
to be satisfied by the ${\bf x}^{(t)}$ for all $t\in\{0,1,\ldots,T\}$.
The current production capacity $x_i^{(t)}$ must be nonnegative,
but the $u_i^{(t)}$ may be negative, in which case a withdrawal
of capital would be indicated.


\subsection{Solution to the game}%\noindent 
We want to find the normalised Nash equilibrium
\begin{equation}
  \arg{\rm equil}_{{\bf u}\in{\bf U}}\left\{\Phi_1,\Phi_2,\Phi_3\right\}
\label{dyneq}
\end{equation}
subject to the state equations and constraints given above.

To do this using the relaxation algorithm, it is necessary
to reformulate the payoff functions and constraints so that they contain
only constants and the actions.  That is, they need to be
free of the states for us to start.

The reformulation for the payoff functions is
\begin{eqnarray}
  \nonumber
  \Phi_i({\bf u}) &=& \sum_{t=0}^{T-1}\rho^t\left\{\phi_i\left( \left[ \begin{array}{c} (1-\mu_1) x_1^{(0)} + \sum_{s=1}^{t}(1-\mu_1)^{t-s} u_1^{(s)} \\ (1-\mu_2) x_2^{(0)} + \sum_{s=1}^{t}(1-\mu_2)^{t-s} u_2^{(s)}\\ (1-\mu_3) x_3^{(0)} + \sum_{s=1}^{t}(1-\mu_3)^{t-s} u_3^{(s)} \end{array} \right] \right)
                      - c_{3i}(u_i^{(t)})^2\right\} \\ 
                  & & + \rho^T\left\{ (1-\mu_i) x_i^{(0)} + \sum_{s=1}^{T}(1-\mu_i)^{t-s} u_i^{(s)}\right\}
                    \hspace{20mm}i=1,2,3
\end{eqnarray}
where the sums over $s$ are interpreted as zero if $t$ or $T$ are 
respectively zero.  The constraints become
\begin{equation}
  q_\ell({\bf u}^{(t)}) = \sum_{i=1}^3\alpha_{j \ell}e_j\left( (1-\mu_i) x_i^{(0)} + \sum_{s=1}^{t}(1-\mu_i)^{t-s} u_i^{(s)} \right) \leq 100
\end{equation}
for $\ell=1,2$ and $t \in \{0,1,\dots T-1\}$, where the sum 
over $s$ is again treated as zero if $t=0$.

In the rest of this section, we will be concerned with the players' first
period actions, \ie ~$u_1^{(0)}$. 
 The second period action depends on the parameters as follows:
\begin{equation}
  u_i^{(1)} = \dfrac{\rho k}{2c_{3i}}.
\end{equation}
We now solve the above game for some different parameter combinations.

{\em Two periods with no depreciation or scrap value. }
Let $T=2$, $\rho\in(0,1]$, $\mu_i\equiv 0$, $k_i\equiv 0$, 
$c_{3i}\equiv 1$ and let the
initial state vector be the Nash solution of the game over one period,
namely ${\bf x}^{(0)}=(21.149,16.030,2.722)$.

Given this information, and a starting guess of
${\bf u}=((0,0),(0,0),(0,0))$, which is the expected
outcome, the program output is the starting guess.

This is intuitively correct since we are just playing
an identical game twice, with a penalty but no reward
associated with investment, and thus expect the
Nash solution for the single period to be
repeated over both periods. The constraints' satisfaction
in each period is the same as in the static game.

{\em Two periods with depreciation but no discounting or scrap value.}\,
Let $T=2$, $\rho=1$, $\mu_i\equiv 0.1$, $k_i\equiv 0$, 
$c_{3i}\equiv 1$ and let the
initial state vector be the Nash solution of the game over one period,
namely ${\bf x}^{(0)}=(21.149,16.030,2.722)$.

The solution in this case, with the same starting guess as before,
 is \[{\bf u}=((0.9577,0),(0.4305,0),(1.1782,0)). \]
That is, the players adjust the states to ${\bf x}^{(1)}
=(19.992,14.858,3.628)$
(in period 1), and make no further investment.
The payoffs  corresponding to the above actions are $(93.79,52.77,14.02)$.
The fact that the state does not revert to the initial state is caused 
by the presence of investment costs.  
The constraints for period 1 are inactive ($q_1=-1.49$,
$q_2=-20.77$.)

{\em Two periods with scrap value and discounting but no depreciation. }
Let $T=2$, $\mu_i\equiv 1$, $c_{3i}\equiv 1$, $k=1$ and 
${\bf x}^{(0)}=(21.149,16.030,2.722)$ as before,
but now let $\rho$ be less than one.

With starting guess ${\bf u}=((0,0),(0,0),(0,0))$,
$\rho$ was varied with the resulting payoff evolution 
shown in Figure \ref{fig:payoff_vary_rho}. As in one agent problems,
the larger $\rho$ the higher the payoffs. 
\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdrhp1-c.eps}}
    \scalebox{0.25}{\includegraphics{rivdrhp2-c.eps}}
    \scalebox{0.25}{\includegraphics{rivdrhp3-c.eps}}
    \caption{How payoff varies with $\rho$ for players 1,2,3 respectively}
    \label{fig:payoff_vary_rho}
  \end{center}
\end{figure}
The payoffs are achieved through actions shown in 
Table \ref{tab:rivedii} and resulting in the state evolution
as in Figure \ref{fig:state_vary_rho}. 
Figure \ref{fig:action_vary_rho} compares the first 
period actions of the players (\ie $u_i^{(1)}$, 
$i=1,2,3$) for different discounting factors.

\begin{table}[ht]
  \begin{center}
    \leavevmode
    \begin{tabular}{|l|ll|ll|ll|}
      \hline
      $\rho$  & \multicolumn{2}{c|}{${\bf u}_1$} & \multicolumn{2}{c|}{${\bf u}_2$} & \multicolumn{2}{c|}{${\bf u}_3$} \\
      \hline
      1      & 0.0203 & 0.500 & 0.2985 & 0.500 & -0.1064 & 0.500 \\
      0.99   & 0.0199 & 0.495 & 0.2927 & 0.495 & -0.1044 & 0.495 \\
      0.98   & 0.0195 & 0.490 & 0.2869 & 0.490 & -0.1023 & 0.490 \\
      0.97   & 0.0191 & 0.485 & 0.2812 & 0.485 & -0.1002 & 0.485 \\
      \hline
    \end{tabular}
    \caption{Nash equilibria for varying $\rho$ in Section \ref{subsec:rived}}
    \label{tab:rivedii}
  \end{center}
\end{table}

\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.3}{\includegraphics{rivdrhs1-c.eps}}
    \caption{Player 1 states for varying $\rho$}
    \label{fig:state_vary_rho}
  \end{center}
\end{figure}

For the static equilibrium of Section \ref{ex4} as well as for all dynamic
equilibria computed in this section, the first of the constraints ($\ell=1$)
has been active.  An increase in production by players 1 and 2 is
compensated by the production decrease of player 3 who is the heaviest
polluter and the highest cost producer (see Table \ref{ex3_consts}.)

\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdrha1-c.eps}}
    \scalebox{0.25}{\includegraphics{rivdrha2-c.eps}}
    \scalebox{0.25}{\includegraphics{rivdrha3-c.eps}}
    \caption{Player 1, 2 and 3 actions ($t=0$) for varying $\rho$}
    \label{fig:action_vary_rho}
  \end{center}
\end{figure}

{\em Two periods with barely active constraints. }
We look here at how the Nash equilibrium behaves when
a constraint is non-active  in the first period ($t=0$). 
 In this way we can observe the behaviour
of the players when constraints become active.

Redefine the constraints as
\begin{eqnarray}
  q_1({\bf u}^{(t)}) = \sum_{i=1}^3\alpha_{j1}e_j\left( (1-\mu_i) x_i^{(0)} + \sum_{s=1}^{t}(1-\mu_i)^{t-s} u_i^{(s)} \right) \leq 424 \\
  q_2({\bf u}^{(t)}) = \sum_{i=1}^3\alpha_{j2}e_j\left( (1-\mu_i) x_i^{(0)} + \sum_{s=1}^{t}(1-\mu_i)^{t-s} u_i^{(s)} \right) \leq 304,
\end{eqnarray}
then the constraints are nonactive when the game is played
for a single period.  However, when we play for two periods
and vary $\rho$, the constraints become active.

{\underline {Firstly}},
setting $\mu_j$ to zero for $j=1,2,3$, and letting
$\rho$ vary from 0.90 to 1.00 give the payoffs and actions 
of Figures \ref{fig:payoff_vary_rho_p}  and \ref{fig:action_vary_rho_p}
 respectively.
\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdrhp1-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdrhp2-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdrhp3-p.eps}}
    \caption{How payoff varies with $\rho$ for players 1,2,3 respectively for the new constraints}
    \label{fig:payoff_vary_rho_p}
  \end{center}
\end{figure} 

\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdrhc2-p.eps}}
    \caption{How constraint values vary with $\rho$ for the new constraints}
    \label{fig:constr2_vary_rho_p}
  \end{center}
\end{figure}

Again, the payoffs are high when the discount factor is large.  
Regarding the players' actions shown in Figure \ref{fig:action_vary_rho_p}, 
they are such that for a certain
value of $\rho$ (here, $\rho=.95$ see Figure \ref{fig:constr2_vary_rho_p})
the second constraint becomes binding. From this value of $\rho$
``onwards'' the least economic agent's actions (\ie\, Player 3's;
 see his parameteres in Table \ref{ex3_consts}) 
start diminishing. 
\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdrha1-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdrha2-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdrha3-p.eps}}
    \caption{Player 1, 2 and 3 actions ($t=0$) for varying $\rho$ with the new constraints}
    \label{fig:action_vary_rho_p}
  \end{center}
\end{figure} 
\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.3}{\includegraphics{rivdrhs1-p.eps}}
    \caption{Player 1 states for varying $\rho$ with the new constraints}
    \label{fig:state_vary_rho_p}
  \end{center}
\end{figure} 

The evolution of the state is shown for
Player 1 in Figure \ref{fig:state_vary_rho_p}. We note that 
his capacities in period 1 surpass those of the ``static''
equilibrium, and grow in the discount factor. An apparently
big jump in capacities in period 2 is mainly due to 
the relatively large weight of the scrap value in 
the agent's objective function.


\underline{Secondly}, varying $\mu_j$ from
being 0 to 0.10 (for all players simultaneously) 
and keeping $\rho=0.96$,  we get
the results displayed in Figures \ref{fig:payoff.vary.mu.p},
\ref{fig:constr2.vary.mu.p} and   \ref{fig:action.vary.mu.p}.

\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdmup1-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdmup2-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdmup3-p.eps}}
    \caption{Player 1, 2 and 3 payoffs for varying $\mu_j$ with the new constraints}
    \label{fig:payoff.vary.mu.p}
  \end{center}
\end{figure}
 \begin{figure}[htp]
   \begin{center}
     \leavevmode
     \scalebox{0.25}{\includegraphics{rivdmuc2-p.eps}}
     \caption{How constraint values vary with $\mu_j$ for the new constraints}
     \label{fig:constr2.vary.mu.p}
   \end{center}
 \end{figure} 

The highest payoffs are (obviously) for no depreciation ($\mu=0$).
The impact of contraints on the equilibrium actions is noticeable
in this case too. For the depreciation parameter values close
to zero, a contraint becomes binding and the players have to diminish
their actions.

\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \scalebox{0.25}{\includegraphics{rivdmua1-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdmua2-p.eps}}
    \scalebox{0.25}{\includegraphics{rivdmua3-p.eps}}
    \caption{Player 1, 2 and 3 actions for varying $\mu_j$ with the new constraints}
    \label{fig:action.vary.mu.p}
  \end{center}
\end{figure}





% Moreover, the more important the future becomes ($\rho\uparrow$), the ``greedier'' the players become, and the more likely is an equilibrium with active constraints. In real life, this may mean that the players should be more strictly monitored when the value of $\rho$ becomes large.


\clearpage
\section{Conclusions} 
The Nikaido-Isoda function was introduced as a means of
making the problem of finding a Nash equilibrium one of
maximising a multivariate function.

This together with the relaxation methodology allows
us to find the Nash equilibrium of many non-zero sum 
multiplayer games in a much simpler fashion than has
previously been possible.  Software was developed to
implement these ideas and has been used successfully
to solve some nontrivial examples.

The games we solved included several elementary static 
games and some less elementary games including the dynamic
river basin pollution game based on \cite{HK}.  
The latter was originally solved
using Rosen's Algorithm.

The developed software enabled us to find the Nash normalised
equilibrium of a two period dynamic game played under the
open loop information pattern.  In particular, we observed
how varying the discount factor and capital depreciation
influence the equilibrium.

We noted that if a one period game has a (static) Rosen
normalised equilibrium that makes a constraint active, this
constraint usually remains active in the dynamic game.
The players' investment behaviour is such that the heaviest
polluter should decrease production for the game to have an
equilibrium.  All players' payoffs decrease with a diminishing 
discount factor.

However, if a static game has an unconstrained static equilibrium,
the dynamic equilibrium generally implies an increased production by all
players and that (usually one of) the constraints becomes active.
Also, the more important the future becomes (\ie ~$\rho$ approaches 1) 
the ``greedier'' the players become and the more likely is an
equilibrium with active constraints.  In real life this may mean
that players should be more strictly monitored for larger values
of $\rho$.  We summarise that in dynamic games 
with coupled constraints \underline{and} 
a final payoff function, the 
equilibrium (\ref{dyneq}), with $\rho\rightarrow 1$, moves
{\em away}\, from the static, or ``repeated'', equilibrium
of Section \ref{RBPsol}.   

Notice that all dynamic equilibria are Rosen normalised.  In order
to force the players to invest according to these equilibria,
Pigouvian taxes can be imposed \cite{HK}.

The relaxation methodology increases the ability of game theory 
to solve ever more complex problems, and does so using 
simple ideas.


\newpage


\appendix \label{appA}
\section{Notes on customarily designed Matlab software 
for solutions of games}
This section contains a brief description of the MATLAB 
program we designed to implement the relaxation algorithm.

\subsection{Program design}
The flow diagram is displayed in Figure \ref{fig:design}.

\begin{figure}[htp]
  \begin{center}
    \leavevmode
    \rotatebox{270}{\scalebox{0.5}{\includegraphics{design.eps}}}
    \caption{Program flow diagram}
    \label{fig:design}
  \end{center}
\end{figure}

\subsection{Description of the procedures}

\subsubsection{User interface {\tt interfac.m}}
This reads in required parameters which describe the game
to be solved, and feeds them to the algorithm.  It then outputs
the results to screen.

\subsubsection{Parameter input {\tt getparms.m, scriptfile}}
The input to the program can be either through a series
of prompted questions or through use of a MATLAB m-file 
script which takes no arguments and returns all required
parameters to the program.

The required parameters at the time of writing are
the names of the payoff function and constraint function,
a description of the game, a vector giving the dimension
of each player's action space, the lower and upper bounds
on the players' actions, the number of equality constraints,
the precision required, the step-size optimisation method
 and a starting guess.

\subsubsection{The payoff function}
This takes the collective action and the player number
and returns the corresponding payoff.

\subsubsection{The constraint function}
This takes the collective action vector and returns a
vector of numbers which must be nonpositive to satisfy
the constraints.  Some may be required to be equal
to zero if the number of equality constraints is
positive.

\subsubsection{The relaxation algorithm {\tt relax.m}}
The main procedure containing the control loop for the
algorithm.  This function takes the required parameters,
and iterates the algorithm until the termination
conditions have been reached.  It returns the
normalised Nash equilibrium together with the
points traversed in achieving the solution,
the individual payoffs at each step, the
step sizes used and the time taken.

\subsubsection{Next iterate {\tt nextx.m}}
This simple function takes the $x^s$, $Z(x^s)$ and $\alpha$
and returns $x^{s+1}$.

\subsubsection{Nikaido-Isoda function {\tt nik\_iso.m}}
Returns the value of the Nikaido-Isoda function at the
point $({\bf y},{\bf x})$.  The order of the arguments is the reverse
of that in the paper for technical reasons.

\subsubsection{Optimum response function {\tt z\_const.m}}
This returns the ${\bf y}$ which maximises the Nikaido-Isoda
function at the point ${\bf x}$.

\subsubsection{Step-size optimisation primitive method {\tt op\_alpha.m}}
This simple step-size optimisation only works when the players
have identical payoff functions.  It simply determines the
$\alpha$ which maximises the common payoff.

\subsubsection{Step-size optimisation general method {\tt op\_alph2.m}}
This step-size optimisation works by trying to minimise the
maximum of the Nikaido-Isoda function at the next step.

\subsubsection{Vector to string function {\tt vect2str.m}}
Takes a collective action and the dimensions of the players'
action spaces, and returns a string output of a nested vector 
which can be output to screen.  For example, in the dynamic 
river pollution game, where each player had two actions, the
solution was $[[0,0][0,0][0,0]]$ for $\rho=1$.


\newpage
\begin{thebibliography}{99} 
\begin{normalsize}
\bibitem{aubin}
{\sc J.P. Aubin}, 
1979,
``Mathematical Methods of Game and Economic Theory''
{\em Amsterdam, North Holland.} 
 
\bibitem{bas1}
{\sc T. Ba\c{s}ar} 
1987,
~``Relaxation techniques and asynchronous algorithms for online computation
of non-cooperative equilibria'',
{\em Journal of Economic Dynamics and Control}, Vol. 11, pp. 531-549.

\bibitem{HK}
{\sc Alain Haurie \& Jacek B. Krawczyk}, 
1997, ~``Optimal charges on river effluent from lumped 
and distributed sources'', 
{\em Environmental Modelling and Assessment}, vol. 2, no.3, pp 93-106.

\bibitem{pjm1}
{\sc Nikaido Hukuhane and Isoda Kazuo},
1955,
~``Note on noncooperative convex games'', 
{\em Pacific Journal of Mathematics}, Vol. 5, Supp. 1, pp 807-815.

\bibitem{ECON407}
{\sc Jacek B. Krawczyk}, 
1995, 
~``ECON 407 - Economic Dynamics B,
Introduction to Dynamic Games with Application - Course Notes'', 
{\em Printing Unit, Works and Services, Victoria University of Wellington} 1774/95/5.

\bibitem{rosen65}
{\sc J. B. Rosen},
1965,
~``Existence and uniqueness of equilibrium points for concave
n-person games'',
{\em Econometrica}, Vol. 33, No. 3, pp 520-534.

\bibitem{uryasev90}
{\sc S. P. Uryas'ev},
1990,
-``{\cyr Adaptivnye Algoritmy Stohastiqesko{j0} Optimizacii i Teorii Igr }''
~(``Adaptive Algorithms for Stochastic Optimisation and Game Theory,'')
{\cyr Moskva <Nauka>}.

\bibitem{UR}
{\sc Stanislav Uryas'ev \& Reuven Y. Rubinstein}, 
1994, 
~``On Relaxation Algorithms in Computation of Noncooperative Equilibria'', 
{\em IEEE Transactions on Automatic Control}, Vol. 39, No. 6. pp. 1263-1267.

\end{normalsize}
\end{thebibliography}

\end{document}



\noindent {\em JEL}\, Keywords: C63 (Computational Techniques),
C72 (Noncooperative Games), C87 (Economic Software), E62 (Taxation),
Q25 (Water; Air). 

\noindent {\em AMS}\, Keywords: 90-08 (Computational Methods), 
90A14 (Equilibrium: general theory),
90A30 (Environmental economics), 90A56 (Special types of equilibria), 90A58
(Models of real-world systems), 90A70 (Macroeconomic policy-making, taxation),
90D05 ($2$-person games), 90D06 ($n$-person games), 90D10 (Noncooperative games),
90D50 (Discrete-time games), 90D80 (Applications of game theory). 

\noindent {\em Authors'}\, Keywords: Computational economics; Nash normalised 
equilibrium, coupled constraints, Nikaido-Isoda function, 
open-loop Nash equilibrium.