\section{Introduction and the main result}
\label{sec:intro}

The setting of this paper is similar to that of Kelly et al. \cite{1998xxJORS_KelMauTan}. Consider a
network with $m$ directed link resources. Let $c(l)$ be the capacity of the link $l$. There are $n$
users and each has a single fixed path. Each user sends data along its associated path with the
first vertex of the path being the source of the user's data and the last vertex being its terminus.
Let $A$ be the $m \times n$ matrix with $A_{le}=1$ if the path $e$ uses link $l$ and $A_{le}=0$
otherwise. Let $[n]$ denote the set $\{1,2,\ldots,n\}$ of users and let $[m]$ denote the set
$\{1,2,\ldots,m\}$ of links. Let $w_e:\mathbb{R}_+ \rightarrow \mathbb{R},~e \in [n]$ be the
utility functions of the users. User $e$ derives a utility $w_e(x(e))$ when sending a flow of rate
$x(e)$. The functions $w_e, e \in [n]$ are assumed to be strictly concave and increasing. Let
$x=(x(e), e \in [n])$, $w= (w_e, e \in [n])$, and $c=(c(l), l \in [m])$. Let $\mathcal{A}=\{x|x \geq
0 ~\text{and} ~ Ax \leq c\}$. Throughout, we make the standing assumption that $\mathcal{A}$ has an
interior feasible point, i.e., there exists a point for which all inequalities are
strict. The system optimal operating point solves the problem:

\begin{align}
 \text{System}(w,A,c)&: \nonumber \\ \label{eqn:maximizeSumSeparableConcaveFunctions}
 \max_{x \in \mathcal{A}} ~ W(x) &:= \sum_{e=1}^{n}w_e(x(e)).
\end{align}
 The catches are that the network operator does not know the utility functions of the users, and the users know neither the rate choices of the other users nor the flow constraints on the network.

Kelly \cite{199701ETT_Kel} proposed the decomposition of the above problem into two subproblems, one
to be solved by each user, and the other to be solved by the network. Let $\lambda_e$ be the cost
per unit rate to user $e$ set by the network, and let $p_e$ be the price user $e$ is willing to pay.
The maximization problem solved by user $e$ is
\begin{align}
 \text{User}(w_e;\lambda_e):&\nonumber ~~~~~~~~~~~~~~~~~~~~~~~~~~\\ \label{eqn:userob}
 \max_{p_e:p_e \geq 0}~ &w_e \left(\frac{p_e}{\lambda_e}\right)-p_e .

\end{align}
If $p = (p_e, e \in [n])$ is known to the network, the network's optimization problem is
 \begin{align}
 \text{Network}(A,c;p):& \nonumber \\ \label{eqn:maximizeSumLOgFunctions}
 \max_{x\in \mathcal{A}}  ~& \sum_{e=1}^{n}p_e \log(x(e)).~~~~~~~~~~~
\end{align}
The solution to Network$(A,c;p)$ is well-known to satisfy the so-called proportional fairness
criterion: if $\mu_l, l \in [m]$ are the optimal dual price variables associated with the dual to
Network$(A,c;p)$, then
\begin{align}\label{eqn:propfairdual}
x(e)=\frac{p_e}{\sum_{l:l\in e}\mu_l}, e \in [n],
\end{align}
is the optimal solution to the network problem.
Kelly \cite{199701ETT_Kel} showed that there exist costs per unit rate $(\lambda^{\star}_e, e \in
[n])$, prices $(p^{\star}_e, e \in [n])$, and flows $(x^{\star}(e), e \in [n])$, satisfying
$p^{\star}_e = \lambda^{\star}_e \cdot x^{\star}(e)$ for $e \in [n]$ such that $p^{\star}_e$ solves
User$(w_e; \lambda^{\star}_e)$ for $e \in [n]$ and $(x^{\star}(e), e \in [n])$ solves
Network$(A,c;p^{\star})$; furthermore, $(x^{\star}(e), e \in [n])$ is the unique solution to
System$(w,A,c)$. The costs per unit rate satisfy $\lambda^{\star}_e=\sum_{l:l\in e}\mu_l^{\star}$
for some dual price variables.

In order to ensure operation at $x^{\star}$, taking the information asymmetry constraints into
account, Kelly et al. \cite{1998xxJORS_KelMauTan} proposed the following {\em fast user adaptation}
dynamics:

\begin{align}
 p_e(t)&= x(e,t) \cdot w_e^\prime(x(e,t)),~ e \in [n],\label{eqn:price}\\
 \frac{d}{dt}x(e,t)&=\kappa \cdot \left(p_e(t)-x(e,t) \cdot \sum_{l : l \in e} \mu_l(t)\right),
~e \in [n],\label{eqn:update}\\
\mu_l(t)&=\psi_l \left( \sum_{e:e \ni l} x(e,t)\right), \quad l \in [m] \label{eqn:cost},
\end{align}
where $\psi_l(y)$ is a penalty\footnote{Kelly et al. \cite{1998xxJORS_KelMauTan} suggest
two functions as examples, one of which is $\psi_l(y)=(y-c(l)+\varepsilon)^+/\varepsilon^2$ for
some $\varepsilon > 0$.} or cost per unit flow when the total flow in the link is $y$. It signifies
the level of congestion in that link. Thus $\mu_l(t)$ in (\ref{eqn:cost}) is the cost per unit flow
through link $l$, and may be interpreted as a dual variable of the network problem. The optimal
dual variables for Network$(A,c;p)$ are such that the net cost of user $e$ flow matches the price
$p_e$ paid by that user; see (\ref{eqn:propfairdual}). The network, adapts the flow
$x(e,t)$ using an additive-increase multiplicative-decrease scheme as in (\ref{eqn:update});
the network therefore attempts to equalize, albeit slowly, the instantaneous
net cost of user $e$ flow, $x(e,t) \cdot \sum_{l : l \in e} \mu_l(t)$, to the instantaneous price
paid by that user, $p_e(t)$. On the other hand, if we differentiate (\ref{eqn:userob}) with respect
to $p_e$ and use the relation $p_e(t)=\lambda_e(t) \cdot x(e,t)$, we get that $p_e(t)$ in
(\ref{eqn:price}) maximizes User$(w_e; \lambda_e(t))$. So the users adapt {\em instantaneously} (in
comparison to the network's slower speed of adaptation) to the congestion signal. Kelly et al.
\cite{1998xxJORS_KelMauTan} provided a Lyapunov function for the dynamical system defined by
(\ref{eqn:price})-(\ref{eqn:cost}). The stable equilibrium point of the dynamical system maximizes a
relaxation of the system problem, as determined by the choice of $\psi_l(\cdot)$  in
(\ref{eqn:cost}).

Despite the popularity of this approach, there are two issues we would like to highlight.
\begin{itemize}
\item $x(t)$ may not remain feasible at all times $t$.
\item $x(t)$ converges to the optimal value of a relaxation of the system problem.
\end{itemize}
The first of these issues was highlighted in Johansson et al. \cite{2006xx_JohSolJoh}. The dynamics
(\ref{eqn:price})-(\ref{eqn:cost}) cannot then be used in systems where feasibility has to be
ensured at all times. The second of these issues is often circumvented via iterative algorithms
where the Lagrange
multipliers or penalty functions are also adapted (see for example, Arrow and Hurwicz
\cite{1958xxEJES_ArrHur}, Low and Lapsley
\cite{1999xx_LowLap}, Chiang et al. \cite{2007xxIEEE_ChiLowCalDoy}, Palomar and Chiang
\cite{2006xxJSAC_PalChi}). Such approaches either assume knowledge of the utility functions at the
network end or may encounter infeasible iterates, or both.

The literature on network utility maximization is vast and we will not be able to do justice to all
of
them. However, there are two works, Hochbaum \cite{1994xxMOR_Hoc} and La and Anantharam
\cite{2002xxTON_LaAna}, that are very relevant to our contribution which we bring to the reader's
attention. A greedy algorithm proposed by Hochbaum \cite{1994xxMOR_Hoc} can be adapted to solve the
system problem with iterates remaining feasible at all times and without full knowledge of the
utility functions at the network side. Though the algorithm circumvents both the problems
highlighted above, it works only when the  set of feasible flows forms the ``independent set of a
polymatroid". This is the case when the network has, for example, a single source and multiple sinks
or when the network has multiple sources but a single sink.

Mo and Walrand \cite{2000xxTON_MoWal} proposed a window-based rate control mechanism that converges
to the solution to Network$(A,c;p)$ for a fixed $p$. The window update rule of
\cite{2000xxTON_MoWal} uses only delay information provided to the user (propagation and
round-trip delays). La and Anantharam \cite{2002xxTON_LaAna}
proposed two algorithms that solve the system problem using the decomposition of Kelly et al.
\cite{1998xxJORS_KelMauTan}. The first algorithm incorporates the solution to the user problem into
the window update rule of \cite{2000xxTON_MoWal}. The second algorithm of La and Anantharam
\cite{2002xxTON_LaAna} explicitly finds the solution to the user problem and the network problem in
each iteration. Although their simulations showed the convergence of the algorithm for general
networks, a mathematical proof was given only for the case of a network with a single link. Their
algorithm additionally imposes more stringent conditions on the utility functions than those assumed
in this paper.

In this paper, we propose a discrete-time algorithm (see Algorithm \ref{alg:mainalgorithm} below)
that remains feasible at all times and converges to the global maximum of System$(w,A,c)$. The
corresponding continuous time dynamics also shares the same properties. Our algorithm can be
considered as a generalization of the algorithm of La and Anantharam \cite{2002xxTON_LaAna} to 1)
more
general networks, and 2) a larger class of objective functions.

For a set of flows, abusing notation, write $p_e(x(e)):=x(e) \cdot w_e^{\prime}(x(e))$ as per
(\ref{eqn:price}), and set $p(x) = (p_e(x(e)), e \in [n])$. Write
\begin{align}
  \label{eqn:functionT}
  T(x) \coloneqq \arg \max_{y \in \mathcal{A}} \sum_{e=1}^{n}p_e(x(e)) \log(y(e))
\end{align}
for the solution to Network$(A,c;p(x))$. If $p_e(x(e))=0$ for some $e$, then the objective function
in (\ref{eqn:functionT}) is not strictly concave over $\mathcal{A}$. The optimization problem
(\ref{eqn:functionT}) may then have multiple solutions, and so $T(x)$ is to be viewed as a
set-valued mapping whose values are convex and compact subsets of $\mathcal{A}$. Define
$a_k\coloneqq \frac{1}{k+1}, k =0,1,2,\ldots$.

\begin{algorithm}\label{alg:mainalgorithm}
\begin{enumerate}
\item Initialize $x^{(0)} \in \mathcal{A}$  such that $x^{(0)}(e) > 0,~e \in [n]$. Initialize $k=0$.
\item User update:\label{step:step2}
 \begin{align}\label{eqn:userupdate}
  p_e^{(k)}=x^{(k)}(e) \cdot w_e^{\prime}(x^{(k)}(e)),~e \in [n].
 \end{align}
\item Network update:\label{step:step4}

Find a point $v \in T(x^{(k)})$ and set:
\begin{align}x^{(k+1)}= x^{(k)}+a_{k+1}(v-x^{(k)}).\label{eqn:stocapproxscheme}
\end{align}
\item Set $ k\leftarrow k+1$ and go to step \ref{step:step2}.
\end{enumerate}
\end{algorithm}

Our main result is the following theorem.
\vspace{2mm}
\begin{thm}
\label{thm:maintheorem}
Assume that $\mathcal{A}$ has an interior feasible point. The iterates $x^{(k)}$ of Algorithm
\ref{alg:mainalgorithm} converge to the optimal solution $x^{\star}$ to the system problem, i.e.,
$x^{(k)} \rightarrow x^{\star}$ as $k \rightarrow \infty$.
\end{thm}
\vspace{2mm}
Observe that, in Algorithm \ref{alg:mainalgorithm}, the users exhibit the same fast adaptation as in
the dynamics (\ref{eqn:price}). But in the network update, iterate $x^{(k+1)}$ is a convex
combination of $x^{(k)}$ and $v \in T(x^{(k)})$ which, by induction, remains in the feasible set for
all $k$. This resolves the feasibility issue in the dynamics of (\ref{eqn:price})-(\ref{eqn:cost}).
In the proof we will argue that the iterates track the differential inclusion
\begin{align}\label{eqn:mainalgdi}
\frac{d}{dt}x(t) \in T(x(t)) - x(t);
\end{align}
we will in fact see that the solution to this differential inclusion remains
feasible at all times.

Theorem \ref{thm:maintheorem} also states that the iterates converge to the global optimum of the
system problem. This resolves the issue that the dynamics (\ref{eqn:price})-(\ref{eqn:cost})
converge to the solution to a relaxed problem different from the original system problem.

The main technical issues to surmount are 1) the dynamics in (\ref{eqn:mainalgdi}) have
multiple fixed points and it is nontrivial to show convergence to the global optimum; 2) $T(x)$
is not necessarily a continuous function of $x$.

Algorithm \ref{alg:mainalgorithm} assumes the existence of a central entity that computes the
solution to Network$(A,c;p)$, i.e. $T(x)$. Section \ref{sec:complexity} describes algorithms to
compute $T(x)$ efficiently for some class of networks. The assumption regarding the existence of a
central entity is not crucial. An alternative is to use Mo and Walrand's algorithm
\cite{2000xxTON_MoWal} that finds $T(x)$ for a fixed $p$; that algorithm uses only the information
available at the user end. The $p$ can then be adapted (user updates) at a slower time scale. This
allows for the use of Algorithm \ref{alg:mainalgorithm} in large
scale networks in a distributed setting.

The rest of the paper is organized as follows. In Section \ref{sec:convergence}, we prove Theorem
\ref{thm:maintheorem}. In Section \ref{sec:complexity}, we address the complexity of identifying the
proportionally fair solution point for the network problem. We provide an example of a network where
flows aggregate into a `main branch', reminiscent of traffic from the suburbs flowing into an
arterial highway leading to the downtown of a large city, for which the complexity to solve the
network problem is $\mathcal{O}(n)$. We also argue that this complexity is manageable
($\mathcal{O}(n \log n)$ plus computations for feasibility checks) in situations where the feasible
set is a polymatroid, for example, when all flows either originate or terminate at a single vertex.
We also demonstrate via simulations that the dynamics in (\ref{eqn:mainalgdi}) converge to the
equilibrium  at a faster rate than the dynamics of (\ref{eqn:price})-(\ref{eqn:cost}) for identical
speed parameters $\kappa$. In Section \ref{sec:conclusion}, we end the paper with some
concluding remarks.

\section{Convergence}
\label{sec:convergence}
The update equation in step \ref{step:step4} of Algorithm \ref{alg:mainalgorithm} is a standard
stochastic approximation scheme but without the stochasticity. A common method to analyze the
asymptotic behavior of such schemes is the dynamical systems approach based on the theory of
ordinary differential equations (ODE). But $T(x)$ being a set valued map necessitates the use of
differential inclusions.

The outline of the proof is as follows. We will first characterize the fixed points of the
mapping $T$. We will then argue that the system optimal point is one of the finitely many fixed
points of the mapping $T$. We will next show that the solution to the differential inclusion in
(\ref{eqn:mainalgdi}) models the asymptotic behavior of the iterates $x^{(k)}$.
Following this, We will show that every solution to the differential inclusion converges to
one of the fixed points of $T$ via Lyapunov theory. Finally, we will prove that the fixed point
to which the solution to the differential inclusion converges as $t \rightarrow \infty$ is the
system optimal point.

\subsection{Characterization of the fixed points of $T(x)$}
\begin{defn}
A point $x$ is a fixed point of the set valued map $T$ if $x \in T(x)$.
\end{defn}

Let $S \subset [n]$. Let $\mathcal{A}|_{S}$ be the subset of $\mathcal{A}$ whose points have support
contained within $S$.
Define a subproblem of the system problem as
\begin{align}
\text{Subsystem}(w,A,c,S):&~~~~~\nonumber \\
\max_{y \in \mathcal{A}|_{S}}~&~\sum_{e\in
S}w_e(y(e)).~~~~~~~\label{eqn:objsubproblem}~~~~~~~~~~~~~~~~~~
\end{align}
\vspace{.1mm}

\begin{lemma}\label{lem:finitefixedpoint}
Let $\overline{x}$ be a fixed point of the mapping $T$.  Let $S=\{e:\overline{x}(e) > 0\}$. Then
$\overline{x}$ is the unique optimal solution to the Subsystem$(w,A,c,S)$.
\end{lemma}
\begin{IEEEproof}
We have $\overline{x}(e)=0$ for all $e \in S^{c}$. If $\lim_{x(e) \downarrow
0}x(e)\cdot w_{e}^{\prime}(x(e))>0$ for some $e \in S^{c}$, then any element  $y \in T(\overline{x})$ has $y(e)>0$ which contradicts the fact that $\overline{x}$ is a fixed point. Hence $p_e(\overline{x}(e))=\overline{x}(e)\cdot w_{e}^{\prime}(\overline{x}(e))=0$ for all $e \in S^{c}$, and we may write

\begin{align}
\overline{x} \in T(\overline{x}) = \arg\max_{y \in \mathcal{A}} ~~\sum_{e\in S}p_e(\overline{x}(e)) \log(y(e)).\label{eqn:subnetprob}
\end{align}
Since $\overline{x} \in \mathcal{A}|_S \subset \mathcal{A}$, we also have that  $\overline{x}$ maximizes (\ref{eqn:subnetprob}) over $\mathcal{A}|_{S}$, i.e.,
\begin{align}
\overline{x} = \arg\max_{y \in \mathcal{A}|_S}~~\sum_{e\in S}p_e(\overline{x}(e))
\log(y(e)).\label{eqn:subnetprob2}
\end{align}
Let $\mu_l,l \in [m]$ and $\eta_e, e \in [n]$ be the optimal dual variables for the network
subproblem (\ref{eqn:subnetprob2}). The corresponding Karush-Kuhn-Tucker (KKT) conditions are
\begin{align}
&\frac{p_e(\overline{x}(e))}{\overline{x}(e)}=\sum_{l:l\in e}\mu_l - \eta_e,~e \in S, \label{eqn:lagderivsub}\\
 &\mu_l\cdot \left (\sum_{e:e \ni l}\overline{x}(e)-c(l)\right)=0,~l\in [m],\label{eqn:compslackmusub}\\
& \eta_e \cdot \overline{x}(e)=0,~e \in [n], \label{eqn:compslacketasub}\\
 &\eta_e \geq 0, ~e \in S, ~\mu_l \geq 0,~l\in [m] ~\text{and}~  \overline{x}\in \mathcal{A}|_{S}.\label{eqn:feaspriduasub}
\end{align}
Since $\frac{p_e(\overline{x}(e))}{\overline{x}(e)}=w_e^\prime(\overline{x}(e))$, it is easy to see
that equations (\ref{eqn:lagderivsub}-\ref{eqn:feaspriduasub}) are the KKT conditions of
Subsystem$(w,A,c,S)$ as
well. Since $\overline{x}/2$ is an interior feasible point of $\mathcal{A}|_S$, KKT conditions are
sufficient
for optimality in (\ref{eqn:objsubproblem}) and $\overline{x}$ is the optimal solution to
Subsystem$(w,A,c,S)$. Uniqueness follows from the strict concavity of $w_e, e \in S$.
\end{IEEEproof}

Is every solution to Subsystem$(w,A,c,S)$, $S \subset [n]$, a fixed point of the mapping
$T$? The possibility that $\lim_{x(e) \downarrow 0}x(e)\cdot w_{e}^{\prime}(x(e))>0$ for an $e
\in S^{c}$ and the first
step of the proof of Lemma \ref{lem:finitefixedpoint} says this is not always true. However, we can
assert the following.
\begin{lemma}
The global maximum of the system problem, $x^{\star}$, is a fixed point of the mapping $T$.
\end{lemma}
\begin{IEEEproof}
$x^{\star}$ solves the system problem. Let $S=\{e:x^{\star}(e) >0\}$. $S$
can be a proper subset of $[n]$. Let $\mu_l,l \in [m]$ and $\eta_e, e \in [n]$ be the optimal dual
variables of the system problem. We then have
\begin{align}
&w_e^{\prime}(x^{\star}(e))=\sum_{l:l\in e}\mu_l - \eta_e,~e \in S, \label{eqn:lagderivsys}\\
&w_e^{\prime}(0)=\sum_{l:l\in e}\mu_l - \eta_e,~e \in S^{c}, \label{eqn:lagderivsys2}\\
&\mu_l\cdot \left(\sum_{e:e\ni l}x^{\star}(e)-c(l)\right)=0,~l\in
[m],\label{eqn:compslackmusys}\\
 &\eta_e \cdot x^{\star}(e)=0, ~e \in [n], \label{eqn:compslacketasys}\\
&\eta_e \geq 0,~ e \in [n],~ \mu_l \geq 0,~l\in [m] ~\text{and}~ x^{\star} \in \mathcal{A}.
\label{eqn:feaspriduasys}
\end{align}
 Observe that $w_e^{\prime}(0)$ is finite for an $e \in S^{c}$; otherwise a small increase in
$x^{\star}(e)$ and a corresponding decrease in $x^{\star}(i)$ for a suitable $i \in S$ (which has
finite $w_i^{\prime}(x^{\star}(i))$)  will result in a feasible flow that has a larger objective
function value. Hence $p_e(x^{\star}(e))=0$ for $e \in S^{c}$. Since
$w_e^{\prime}(x^{\star}(e))=\frac{p_e(x^{\star}(e))}{x^{\star}(e)}$ for all $e \in S$, it follows
from (\ref{eqn:lagderivsys})-(\ref{eqn:feaspriduasys}) that $x^{\star},(\mu_l, l\in [m]),
(\tilde{\eta}_e=\eta_e, e \in S)$ and $(\tilde{\eta_e}=\eta_e+w_e^{\prime}(0), e \in S^{c})$
satisfy the KKT conditions of the problem (\ref{eqn:functionT}). Hence $x^{\star} \in T(x^{\star})$.
 \end{IEEEproof}
 \subsection{Need for the theory of differential inclusions}
We now describe the issues that make it necessary to use differential inclusions to study the
asymptotic behavior of $x^{(k)}$. $T(x)$ is the set of points that solve (\ref{eqn:functionT}). If
$p_e(x(e)) = 0$ for some $e$ at a point $x \in \mathcal{A}$, then the objective function in
(\ref{eqn:functionT}) is not strictly concave. Hence there can be multiple points that solve
(\ref{eqn:functionT}). A continuous selection from $T(x)$ allows the use of differential equations
to analyze the stochastic approximation scheme in (\ref{eqn:stocapproxscheme}). A natural question
that arises is whether there is such a continuous selection from $T(x)$. We give an example in the
Appendix showing that such a selection is not possible.

\subsection{Differential Inclusions: Preliminaries}
In this section, we define a differential inclusion and state relevant results from \cite{2005xxSIAM_BenHofSor} that are used to show the convergence of Algorithm \ref{alg:mainalgorithm}.
Let $F: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$ be a set valued map. Consider the following differential inclusion:
\begin{align}\label{eqn:DI}
\frac{dx}{dt}\in F(x).
\end{align}
A solution to the differential inclusion in (\ref{eqn:DI}) with initial condition $x_{0} \in
\mathbb{R}^{n}$ is an absolutely continuous function $x:\mathbb{R} \rightarrow \mathbb{R}^{n}$ that
satisfies (\ref{eqn:DI}) for almost every $t \in \mathbb{R}$. The following conditions are
sufficient for the existence of a solution to the differential inclusion (\ref{eqn:DI}):
\begin{enumerate}
\item \label{itm:property1}$F(x)$ is nonempty, convex and compact for each $x \in \mathbb{R}^{n}$.
\item \label{itm:property2}$F$ has a closed graph.
\item \label{itm:property3}For some $K>0$, for all $x \in \mathbb{R}^{n}$, $F$ satisfies the following condition
\begin{align}\label{eqn:Fboundedness}
\sup_{z \in F(x)} \vert\vert z\vert\vert \leq K (1+\vert\vert x\vert \vert).
\end{align}
\end{enumerate}
The stochastic approximation scheme in $\mathbb{R}^{n}$ is given as
\begin{align}
y^{k+1} \in y^{(k)}+a_{k+1}(F(y^{(k)}) +U^{(k+1)}),\label{eqn:stocapproxschemenoise}
\end{align}
where $a_k$ satisfy the usual conditions:
\begin{align}
\lim_{k \rightarrow \infty}a_k = 0, ~\sum_{k=1}^{\infty}a_k \rightarrow \infty,
\end{align}
and $U^{(k)} \in \mathbb{R}^{n}$ are deterministic or random perturbations.

Let $t(0)=0, t(k)=\sum_{i=1}^{k}a_i$. Let $r_{y}:\mathbb{R}_{+}\rightarrow \mathbb{R}^{n}$ be a
continuous piece-wise linear function formed by the interpolation of $y^{(k)}$ as in
\begin{align}
&r_{y}(t) = y^{(k)}+\frac{y^{(k+1)}-y^{(k)}}{t(k+1)-t(k)}\cdot (t-t(k)),\label{eqn:linintpol} \\
&~~~~~~~~~~~~~~~~~~~~~~~~~\forall ~t \in [t(k),t(k+1)).\nonumber
\end{align}
\begin{defn}
(A perturbed solution to (\ref{eqn:DI})). Let $U:\mathbb{R}_{+} \rightarrow \mathbb{R}^n$ be locally integrable function such that
\begin{align*}
\lim_{t\rightarrow \infty} \sup_{0 \leq v \leq T}\left\vert\left\vert\int_{t}^{v}U(s)ds\right\vert\right\vert = 0.
\end{align*}
Let $\delta:[0,\infty) \rightarrow \mathbb{R}$ be a function such that $\delta(t) \rightarrow 0$ as $t \rightarrow \infty$. Define
\begin{align}
F^{\epsilon}(x) = \{y \in \mathbb{R}^{m}: \exists z: \vert\vert z-x\vert\vert < \epsilon, d(y,F(z)) < \epsilon\},
\end{align}
where $d(y,F(z)) = \inf\{\vert\vert y-q\vert\vert:q \in F(z)\}$.
An absolutely continuous function $y:[0, \infty)\rightarrow \mathbb{R}^{n}$ is a perturbed solution to (\ref{eqn:DI}) if there exists $U:\mathbb{R}_+ \rightarrow \mathbb{R}^{n}$ and $\delta:\mathbb{R}_+ \rightarrow \mathbb{R}^{n}$ as above such that
\begin{align}
\frac{dy}{dt} - U(t) \in F^{\delta(t)}(y(t)).
\end{align}
for almost every $t \in \mathbb{R}_{+}$.
\end{defn}
The following lemma, taken from \cite{2005xxSIAM_BenHofSor}, gives conditions on $y^{(k)}$ and $U^{(k)}$ for $r_y$ to be a perturbed solution to (\ref{eqn:DI}).
\begin{lemma}\label{lem:conditionforps}
\cite[Prop. 1.3]{2005xxSIAM_BenHofSor} Suppose $y^{(k)}$ is bounded, i.e., $\sup_{k} \vert \vert y^{(k)} \vert \vert < M < \infty$, and for all $T > 0$,
\begin{align}
\lim_{s \rightarrow \infty} \sup \Biggl\{ &\left\vert\left\vert\sum_{k=s}^{i-1} a_{k+1}
U^{(k+1)}\right\vert\right\vert:i=s+1,s+2,\nonumber \\& ~~~~~~~~~~~~~~~~~~~~ \ldots, m(t(s)+T)
\Biggr \}=0, \label{eqn:conditionforps2}
\end{align}
where $m(t) = \sup\{m:t(m) \leq t\}$. Then $r_{y}(t)$ is a perturbed solution of the differential
inclusion (\ref{eqn:DI}).
\end{lemma}
\begin{defn}
A compact set $L$ is an internally chain transitive set if for any $x, y \in L$ and every $\epsilon
> 0 , T > 0$, there exists $l \in \mathbb{N}$, solutions $x_1,x_2,\ldots,x_l$ to (\ref{eqn:DI}) and
$t_i > T,~ \forall ~ i$, that satisfy the following.
\begin{enumerate}
\item $x_i(t) \in L$ for all $ 0 \leq t \leq t_i$ and for all $i \in [l]$,
\item $\vert\vert x_i(t_i)- x_{i+1}(0) \vert \vert \leq \epsilon $ for all $i \in [l-1]$,
\item $\vert\vert x_1(0)- x \vert \vert \leq \epsilon$ and $\vert\vert x_l(t_l)- y\vert \vert \leq \epsilon$.
\end{enumerate}
We shall call the sequence $(x_1,x_2,\ldots,x_l)$ as an $(\epsilon,T)$ chain in $L$ from $x$ to $y$.
\end{defn}
The following lemma, again taken from \cite{2005xxSIAM_BenHofSor}, characterizes the limit set of
a perturbed solution.
\begin{lemma}\cite[Thm. 3.6]{2005xxSIAM_BenHofSor}\label{lem:psisict}
Let $r$ be a perturbed solution to (\ref{eqn:DI}). Then the limit set of $r(\cdot)$
\begin{align*}
L(r)\coloneqq \bigcap_{t \geq 0}\{ r(s) : s \geq t\}
\end{align*}
is internally chain transitive.
\end{lemma}
\subsection{Convergence analysis}
We proceed to prove the convergence of $x^{(k)}$ to the optimal solution to the system problem. Observe that $T$ maps points in $\mathcal{A}$ to itself. For $x \in \mathbb{R}^{n}$, define
\begin{align} \label{eqn:Fdefinition}
F(x) \coloneqq T(P_{\mathcal{A}}(x))-x.
\end{align}
where $P_{\mathcal{A}}(x)$ is the projection of $x$ onto the set $\mathcal{A}$.
\begin{lemma}
For each $x \in \mathbb{R}^{n}$, $F(x)$ is nonempty, convex and compact. Furthermore, $F$ has the closed graph property and satisfies (\ref{eqn:Fboundedness}).
\end{lemma}
\begin{IEEEproof}
The objective function of the network problem is continuous and the constraint set $\mathcal{A}$ is
compact. The maximum exists due to the Weierstrass theorem. Also, the set of maximizers is closed
and convex. Thus $T(P_{\mathcal{A}}(x))$ is nonempty, convex and compact, and hence so is $F(x)$.

We next prove the closed graph property of $F$. A function has the closed graph property if it is
upper hemicontinuous. The objective function in (\ref{eqn:functionT}), $\sum_{e=1}^{n}p_e(x(e))\log
(y(e))$, is jointly continuous\footnote{Here we may take $\log(y(e))$ to be continuous at $y(e) =
0$ with $\log \, 0 \coloneqq -\infty$ because then $\lim_{x(e) \downarrow 0} \log(y(e))=\log\, 0$
and this sequential continuity is all that is needed to apply Berge's maximum theorem
\cite[p. 116]{1963xxTS_Ber}} in $x$ and $y$. Also, the constraint set of the network problem does
not vary with $x$. By Berge's maximum theorem \cite[p. 116]{1963xxTS_Ber}, $T(x)$ is upper
hemicontinuous. Since $P_{\mathcal{A}}(x)$, the projection onto the  convex set $\mathcal{A}$, is
continuous, the composition $T(P_{\mathcal{A}}(x))$ is upper hemicontinuous. Consequently, $F(x)$ is
upper hemicontinuous and hence has the closed graph property.

Finally,
\begin{align}
\sup_{z \in F(x)} \vert\vert z\vert\vert &= \sup_{z \in T(P_{\mathcal{A}}(x))} \vert\vert z-x \vert\vert \nonumber\\
&\leq \sup_{z \in \mathcal{A}}\vert\vert z-x \vert\vert\nonumber \\ &\leq \sup_{z \in \mathcal{A}} \vert\vert z \vert\vert + \vert\vert x \vert\vert \nonumber\\
& \leq K + \vert\vert x \vert\vert \nonumber \\ &\leq K (1 + \vert\vert x \vert\vert),\nonumber
\end{align}
where $K > \max\{\sup_{z \in \mathcal{A}}\vert\vert z \vert \vert,1\}$.
\end{IEEEproof}
\begin{lemma}\label{lem:rxisps}
Let $r_x(t)$ be obtained by the linear interpolation of $x^{(k)}$ as given in (\ref{eqn:linintpol}).
Then $r_x(t)$ is a perturbed solution to the differential inclusion (\ref{eqn:DI}) with $F$ defined
as in (\ref{eqn:Fdefinition}).
\end{lemma}
\begin{IEEEproof}
We first show that $x^{(k)} \in \mathcal{A}$. Observe that $x^{(0)} \in \mathcal{A}$. Assume
$x^{(k-1)} \in \mathcal{A}$. Since $T(x^{(k-1)}) \in \mathcal{A}$ and $x^{(k)}$ is a convex
combination of $x^{(k-1)}$ and $T(x^{(k-1)})$, we have $x^{(k)} \in \mathcal{A}$. It follows that
\begin{align*}
F(x^{k}) = T(P_{\mathcal{A}}(x^{(k)}))-x^{(k)} = T(x^{(k)})-x^{(k)} .
\end{align*}

We now see that the update equation in (\ref{eqn:stocapproxscheme}) is the same as the stochastic
approximation scheme in (\ref{eqn:stocapproxschemenoise}) with $U^{(k)}= 0$ for all $k$.
Observe that $x^{(k)}$ is bounded because $x^{(k)} \in \mathcal{A}$ for all $k$ and $\mathcal{A}$
is compact; since $U^{(k+1)}= 0$, the condition in (\ref{eqn:conditionforps2}) is trivially
satisfied.
Hence, by Lemma \ref{lem:conditionforps}, $r_x(t)$ is a  perturbed solution.
\end{IEEEproof}
\vspace{1mm}
We restrict our attention to solutions of (\ref{eqn:DI}) with initial condition $x(0) \in \mathcal{A}$. Since $T(x) \in \mathcal{A}$, $x(t)$ lies in $\mathcal{A}$ for all $t$.
Define \[\Phi_{t}(x_0)\coloneqq \{x(t):x ~\text{solves}~ (\ref{eqn:DI}),~x(0)=x_0\}.\]
\begin{defn}
Let $\Lambda$ be a subset of $\mathcal{A}$. Let $V:\mathcal{A} \rightarrow \mathbb{R}$ be a continuous function such that
\begin{itemize}
\item $V(y) < V(x),~ y \in \Phi_t(x), x \in \mathcal{A}\backslash \Lambda,$
\item $V(y) \leq V(x),~ y \in \Phi_t(x), x \in \Lambda.$
\end{itemize}
Then $V$ is called a Lyapunov function for $\Lambda$.
\end{defn}

Define $V:\mathcal{A} \rightarrow \mathbb{R}$ as
\begin{align}\label{eqn:objlyapunov}
V(x)\coloneqq \sum_{e=1}^{n} w_{e}(x^{\star}(e))-\sum_{e=1}^{n}w_{e}(x(e)).
\end{align}

\begin{lemma}\label{lem:Vlyapunov}
Let $\Lambda$ be the set of fixed points of $T$. The function $V$ in (\ref{eqn:objlyapunov}) is a
Lyapunov function for $\Lambda$.
\end{lemma}

\begin{IEEEproof}
Let $x \in \mathcal{A}$ and $v \in T(x)$. We have, from the definition of $T(x)$ in
(\ref{eqn:functionT}), that
\begin{align}\label{eqn:strictinequality}
\sum_{e=1}^{n}p_e(x(e)) \log (v(e)) \geq \sum_{e=1}^{n} p_e(x(e)) \log (x(e))
\end{align}
because $v \in T(x)$ maximizes the network problem.

If $p_e(x(e)) >0$ for all $e$, then the network problem has unique solution. Therefore, equality
holds in (\ref{eqn:strictinequality}) if and only if $v=x$, i.e, $x$ is a fixed point of the
mapping
$T$. Thus we have
\begin{align}
0 &\leq \sum_{e=1}^{n} p_e(x(e)) \log~\frac{v(e)}{x(e)}\nonumber \\
&\stackrel{(a)}{\leq} \sum_{e=1}^{n} p_e(x(e)) \left(\frac{v(e)}{x(e)}-1\right)\nonumber\\
&=\sum_{e=1}^{n}w_e^{\prime}(x(e))(v(e)-x(e))\nonumber\\
&=\nabla W(x)\bigcdot (v-x), \label{eqn:timederivativepositive}
\end{align}
where (a) uses the inequality $\log\,y\leq y-1$.

More generally, let $p_e(x(e))>0$ for $e \in S \subset [n]$ and $p_e(x(e))=0$ for $e \in S^{c}$; in
particular, $x(e)=0$ for $e \in S^c$. Define $\tilde{v}$ to be $\tilde{v}(e)=v(e)1_{S}(e), e \in
[n]$.

The value of the objective function in (\ref{eqn:functionT}) evaluated at $v$ and $\tilde{v}$ are equal. Hence $\tilde{v} \in T(x)$,
\begin{align}\label{eqn:strictinequality1}
\sum_{e\in S}p_{e}(x(e))\log(\tilde{v}(e))\geq \sum_{e\in S}p_{e}(x(e))\log(x(e)),
\end{align}
and $\tilde{v}$ must be the unique solution to the problem defined in (\ref{eqn:subnetprob2}).
Therefore,
(\ref{eqn:strictinequality1}) holds with equality if and only if $\tilde{v}=x$.
Following the steps leading to (\ref{eqn:timederivativepositive}), we have
\begin{align}\label{eqn:strictinequality2}
\sum_{e\in S}w_e^{\prime}(x(e))(\tilde{v}(e)-x(e)) \geq 0
\end{align} which is a strict inequality if $\tilde{v}\neq x$.
Since $v(e)-x(e) \geq 0$ for $e \in S^{c}$, this along with (\ref{eqn:strictinequality2}) yields
\begin{align}\label{eqn:strictinequality3}
\nabla W(x)\bigcdot(v-x) \geq 0;
\end{align}
since $w_e^{\prime}(x(e))=w_e^{\prime}(0) > 0$ for $e \in S^c$, equality holds in (\ref{eqn:strictinequality3}) if and only if $v=x$. Hence
\begin{align}\label{eqn:rateofchange}
\frac{dV(x(t))}{dt}= -\nabla W(x(t))\bigcdot (v-x(t)) \leq 0 , ~ \forall ~ v \in T(x(t)).
\end{align}
The inequality in (\ref{eqn:rateofchange}) holds with an equality if and only if $v=x$. Therefore $V$ is a Lyapunov function for $\Lambda$.
\end{IEEEproof}

\begin{lemma}\label{lem:ictsubsetlambda}
Let $\Lambda$ be the set of fixed points of $T$. Every internally chain transitive set for $F$ in
(\ref{eqn:Fdefinition}) is a singleton that is a subset of $\Lambda$.
\end{lemma}
\begin{IEEEproof}
By Lemma \ref{lem:finitefixedpoint}, there are at most finitely many fixed points of the mapping
$T$.
Hence the  cardinality of the set $\Lambda$ is finite and $V(\Lambda)$ has empty interior. Also, by
Lemma \ref{lem:Vlyapunov}, $V$ is a Lyapunov function for $\Lambda$. Proposition 3.27 of
\cite{2005xxSIAM_BenHofSor}  states that if $V$ is a Lyapunov function for $\Lambda$ and if
$V(\Lambda )$ has an empty interior, then every internally chain transitive set is a subset of
$\Lambda$.

Choose
$\epsilon$ small enough so that open balls of radius $\epsilon$ centered at each of the finite
points $\Lambda$ are disjoint. Fix $T \geq 0$. Since any $(\epsilon,T)$ chain involves remaining in
$\Lambda$ for all time and jumps of size at most $\epsilon$ to another point in $\Lambda$, by the
disjointedness of the $\epsilon$-balls covering $\Lambda$, there can be no $(\epsilon,T)$ chain
in $\Lambda$ joining two of its distinct points. It follows that the internally chain transitive
subsets of $\Lambda$ are singletons.
\end{IEEEproof}

\begin{lemma}
The iterates $x^{(k)}$ converges to a fixed point of the mapping $T$.
\end{lemma}
\begin{IEEEproof}
In Lemma \ref{lem:rxisps}, we showed that $r_x(t)$ is a perturbed solution to (\ref{eqn:DI}). By
Lemma \ref{lem:psisict}, the limit set of $r_x(t)$ is internally chain transitive. By Lemma
\ref{lem:ictsubsetlambda}, $L(r_x)$ is a singleton and $L(r_x) \subset \Lambda$. Let $\hat{x} \in
L(r_x)$. Since $\mathcal{A}$ is compact and $\hat{x}$ is the only limit point of the sequence
$x^{(k)}$, every
subsequence of $x^{(k)}$ has a further subsequence that converges to $\hat{x}$. Hence $x^{(k)}$
converges to $\hat{x}$.
\end{IEEEproof}
In the rest of this section, we show that the iterates converge to $x^{\star}$, the optimal
solution to the system problem.

Let the dual variables of the optimization problem Network$(A,c;p)$ be $\mu_l, l \in [m]$.
Kelly et al. \cite{1998xxJORS_KelMauTan} simplified the dual to this problem to be:
\begin{align}
&\text{Dual}(p,A,c): \label{eqn:dualT}\nonumber \\
&~~\min_{\mu_l \geq 0, l \in [m]} \left(\sum_{e=1}^{n}p_e\cdot \log\frac{1}{\sum_{l:l\in e}\mu_l} +
\sum_{l=1}^{m}\mu_lc(l) \right).
\end{align}
We now argue that the search for the optimal $\mu_l,l\in [m]$ may be restricted to a compact set.
\begin{lemma}\label{lem:dualcompact}
The optimization problem in (\ref{eqn:dualT}) with $p_e=p_e(x(e))$ is equivalent to the following
optimization problem. For any $x \in \mathcal{A}$,
\begin{align}\label{eqn:dualobjfun}
\max_{0\leq \mu_l \leq 2P/c(l)} \sum_{e=1}^{n}p_e(x(e))\cdot\log \left(\sum_{l:l\in e}\mu_l\right)-\sum_{l=1}^{m}\mu_lc(l),
\end{align} where
\begin{align*}
P\coloneqq \max_{x \in \mathcal{A}}\sum_{e=1}^{n} x(e)\cdot w_e^{\prime}(x(e)) < \infty. \nonumber
\end{align*}
\end{lemma}
\begin{IEEEproof}
Define $R(\mu)$ to be the objective function in (\ref{eqn:dualobjfun}). For any $\mu_l > 2P/c(l)$, by reducing $\mu_{l}$, we increase the objective function's value. To see this, it suffices to show that $\frac{\partial R(\mu)}{\partial \mu_l} < 0$ for any $\mu_l > 2P/c(l)$.
But this is easily checked as follows:
\begin{align}
\frac{\partial R(\mu)}{\partial \mu_l}&=\sum_{e:e\ni l}p_e(x(e))\frac{1}{\sum_{l^{\prime}:l^{\prime}\in e}\mu_{l^{\prime}}}-c(l) \nonumber\\ &\leq \sum_{e:e \ni l} p_e(x(e)) \frac{1}{\mu_l}-c(l) \nonumber \\ &\leq \frac{1}{\mu_l} \sum_{e=1}^{n}p_e(x(e)) -c(l) \nonumber \\ &\leq \frac{1}{\mu_l}\left[\max_{x \in \mathcal{A}}\sum_{e=1}^{n} p_e(x(e))\right]-c(l) \nonumber \\&= \frac{P}{\mu_l}-c(l) < 0,
\end{align}
where the last inequality follows if $\mu_l > 2P/c(l)$.
\end{IEEEproof}
\begin{lemma}\label{lem:convxstar}
Let $x^{(k)}$ converge to $\hat{x}$, a fixed point of the mapping $T$. Then $\hat{x}= x^{\star}$,
the optimal solution to the system problem.
\end{lemma}
\begin{IEEEproof}
Let $v^{(k+1)}$ solve problem (\ref{eqn:functionT}) with $p_e=p_e(x^{(k)}(e))$, and so
$v^{(k+1)} \in T(x^{(k)})$ satisfies the KKT conditions
\begin{align}
&\frac{p_e(x^{(k)}(e))}{v^{(k+1)}(e)} - \sum_{l:l \in e} \mu_{l}^{(k)} +\eta_e^{(k)}=0, ~ e \in [n],\label{eqn:lagderiv} \\
&\mu_l^{(k)}\cdot \left( \sum_{e:e \ni l}v^{(k+1)}(e)-c(l)\right) = 0,~ l \in [m], \label{eqn:compslackintermediate}\\
&\eta_e^{(k)}\cdot v^{(k+1)}(e) = 0, ~e \in [n], \label{eqn:compslackintermediate2}\\
&\eta_e^{(k)}\geq 0,~ e \in [n],~ \mu_l^{(k)} \geq 0 , ~ l \in
[m].\label{eqn:dualpositiveintermediate}
\end{align}
Let us first claim that $x^{(k)} > 0$ for all $k \geq 0$. This holds for $k=0$, the initial point,
in Algorithm \ref{alg:mainalgorithm}. If, for some $k$, $x^{(k)}>0$, then $p_e(x^{(k)}(e)) > 0$ for
all $e$ and so, $v^{(k+1)} > 0$ and consequently, $x^{(k+1)}$ being a convex combination of
$x^{(k)}$ and $v^{(k+1)}$ also satisfies $x^{(k+1)} >0$. The claim follows by induction.
Since $x^{(k)}(e) > 0$, we have $v^{(k+1)}(e) > 0$, and so, by (\ref{eqn:compslackintermediate2}),
$\eta_e^{(k)}=0, ~ e \in [n]$.
Thus (\ref{eqn:lagderiv}) simplifies to
\begin{align}\label{eqn:denmupositive}
\frac{p_e(x^{(k)}(e))}{v^{(k+1)}(e))}=\sum_{l:l \in e} \mu_l^{(k)}.
\end{align}
Since we also have $p_e(x^{(k)}(e)) > 0$, and $v^{(k+1)}(e) > 0$ in (\ref{eqn:denmupositive}), we
have $\sum_{l:l \in e} \mu_l^{(k)} >0$. Hence
\begin{align}
v^{(k+1)}(e) &= \frac{p_e(x^{(k)}(e))}{\sum_{l:l \in e} \mu_l^{(k)}}\nonumber\\
&=\frac{x^{(k)}(e)\cdot w_e^{\prime}(x^{(k)}(e))}{\sum_{l:l \in e} \mu_l^{(k)}}.\label{eqn:tautoinfty}
\end{align}

Suppose $\hat{x}(e)=0$. Observe that $x^{(k)}(e) > 0$ for all $k$. Hence $v^{(k+1)}(e) < x^{(k)}(e)$ infinitely often, which is the same as saying
\begin{align}
\frac{x^{(k)}(e)w_e^{\prime}(x^{(k)}(e))}{\sum_{l:l \in e} \mu_l^{(k)}} < x^{(k)}(e)
\end{align}
occurs infinitely often.
This implies
\begin{align}\label{eqn:strictdecrease}
w_e^{\prime}(x^{(k)}(e)) < \sum_{l:l \in e} \mu_l^{(k)}
\end{align}
infinitely often.
Consider the subsequence that satisfies (\ref{eqn:strictdecrease}). Henceforth, let $x^{(k)}$ denote
that subsequence. We now make the following observations. In Lemma \ref{lem:dualcompact}, we showed
 that $\mu_l^{(k)}$ takes values on a compact set, and so we can find a further subsequence such
that $\mu_l^{(k_\tau)}\rightarrow \mu_l^{\star}$ for some $\mu_l^{\star}$, but for all $l \in [m]$.
Since (\ref{eqn:strictdecrease}) holds for $k = k_\tau$, by letting $\tau \rightarrow \infty$, we have
\begin{align}
w_e^{\prime}(\hat{x}(e)) &\leq \sum_{l:l \in e} \mu_l^{\star}, ~ \forall ~ e~\text{such that}~ \hat{x}(e)=0 \label{eqn:lagderivativefinal}\\
w_e^{\prime}(\hat{x}(e)) &= \sum_{l:l \in e} \mu_l^{\star},~ \forall ~ e~\text{such that} ~\hat{x}(e)>0;
\end{align}
the latter follows by letting $\tau \rightarrow \infty$ in (\ref{eqn:tautoinfty}) and from
$x^{(k+1)}=x^{(k)}+a_k (v^{(k+1)}-x^{(k)})$.
Choose
\begin{align}
\eta^{\star}(e)& =0 ~\text{if} ~\hat{x}(e) > 0, \\
\eta^{\star}(e)&=\sum_{l:l \in e} \mu_l^{\star}-w_e^{\prime}(\hat{x}(e))~\text{if} ~\hat{x}(e) = 0.
\end{align}
Thus, from (\ref{eqn:lagderivativefinal}),
\begin{align}
\eta^{\star}(e)\cdot \hat{x}(e)=0~\forall ~ e \in [n] ~\text{and}~\eta^{\star}(e) \geq 0~ \forall ~
e \in [n].
\end{align}
Since (\ref{eqn:compslackintermediate}) and (\ref{eqn:dualpositiveintermediate}) are true for
indices $k_\tau$, taking limit as $\tau \rightarrow \infty$, we get
\begin{align}
\mu_l^{\star}\cdot \left( \sum_{e:e \ni l}\hat{x}(e)-c(l)\right)& = 0,~ l\in [m],\\
\mu_l^{\star} &\geq 0, ~ \forall ~l.\label{eqn:positivityfinal}
\end{align}
Equations (\ref{eqn:lagderivativefinal})-(\ref{eqn:positivityfinal}) are the KKT conditions for the system problem. Hence $\hat{x}=x^{\star}$, the optimal solution to the system problem.
\end{IEEEproof}

{\em Proof of Theorem \ref{thm:maintheorem}:}
Theorem \ref{thm:maintheorem} follows from Lemma \ref{lem:convxstar}. $\hfill \IEEEQEDclosed$

\section{Algorithmic Complexity and Speed of Convergence}
\label{sec:complexity}
In this section, we remark on the complexity of Algorithm \ref{alg:mainalgorithm}. Each iteration
of
the algorithm has 1) a user update which adapts the amount a user is willing to pay to the
network,  and 2) a network update which adapts the rates allocated to the users.

Since $w_e$ is known at the user end, $w_e^{\prime}$ is easy to obtain either numerically or
analytically. Hence the user update (\ref{eqn:userupdate}) can be implemented by each user in
$\mathcal{O}(1)$ steps.

The network update consists of solving the network problem (\ref{eqn:functionT}). Its complexity
depends on the network structure. We indicate the complexity of the network update for the
following simple networks: a polymatroidal network with a single source and multiple sinks or
multiple sources and a single sink; a flow aggregating network with the structure in Figure
\ref{fig:figure1}.

\begin{figure}[!ht]
\centering
\includegraphics[draft=false,scale=0.6]{network_structure.eps}
 \caption{Flow aggregating network.}
 \label{fig:figure1}
 \end{figure}

{\em Polymatroidal network}: Consider a network with a single source and $n$ sinks. The
source sends flows at a rate $x(1),x(2),\ldots,x(n)$ to the sinks. Megiddo \cite{1974xx_Meg} showed
that the set of feasible flows $(x(e), e \in [n])$ forms the independent set of a polymatroid.
Therefore, the network problem is a separable concave maximization over the independent set of a
polymatroid. The fastest known algorithm that solves this optimization problem is a scaling based
greedy algorithm proposed by Hochbaum \cite{1994xxMOR_Hoc}. The algorithm obtains an
``$\epsilon$-optimal" solution to the network problem in $\mathcal{O}\left(n (\log n +F) \log
\frac{B}{n\epsilon}\right)$ where $F$ is the complexity to check whether a certain increase in one
of the  components of $(x(e), e \in [n])$ would make the flow infeasible. $B$ is the total amount of
resource to be allocated and is $\mathcal{O}(n)$.

{\em Flow aggregating network}:
Let  $(x(e),e \in [n])$ denote the flow through the network in Figure \ref{fig:figure1}. The flow
constraints of the network are
\begin{align}
\begin{split}\label{eqn:lac}
x(e) & \geq 0, e \in [n],\\
x(1) &\leq \alpha(1),\\
x(1)+x(2) &\leq \alpha(1)+\alpha(2),\\
&~\,\vdots\\
x(1)+x(2)+\cdots+x(n) &\leq \alpha(1)+\alpha(2)+\cdots+\alpha(n).
\end{split}
\end{align}
where $\alpha(e) \geq 0$ for all $e$.
The constraints in (\ref{eqn:lac}) are referred to as linear ascending constraints. This problem
arises as the core optimization problem in several wireless communication problems (Padakandla and
Sundaresan \cite{200910TCOM_PadSun}, Viswanath and Anantharam \cite{200206TIT_VisAna}, Lagunas
et al. \cite{200405TSP_PalLagCio}, Sanguinetti et al. \cite{201205TSP_SanDam}) and operations
research problems (Clark and Scarf \cite{1960xxMS_ClaSca}, Wang \cite{201510OL_Ziz}). See
\cite{2016xx_AkhSun} for a survey and a discussion of several algorithms. The network problem is
the maximization of a so-called $d$-separable concave function over the linear ascending
constraints. Veinott Jr. \cite{1971xxMS_Vei} mapped this problem to the geometrical problem of
finding the ``concave cover" of the set of points $(\sum_{e=1}^{i}\alpha(e),\sum_{e=1}^{i}p_e),~ i
\in [n]$, in $\mathbb{R}^2$. The ``string algorithm" of Muckstadt and Sapra \cite{2010xxSPR_MucSap}
finds the concave cover of a set of points in $\mathbb{R}^{2}$ in $\mathcal{O}(n)$ steps. See
\cite{2016xx_AkhSun} for details.

{\em Simulations}:
We now discuss some simulation studies investigating the speed of convergence of our proposed
algorithm. In the plots, the solid curves correspond to the error plots of the system
\begin{align}
 \label{eqn:scaledmainalgdi} \frac{d}{dt}x(t)  \in \kappa \cdot \left ( T(x(t)) -x(t) \right ).
\end{align}
The dashed curves correspond to the error plots of the system (\ref{eqn:price})-(\ref{eqn:cost}).
The differential inclusion (\ref{eqn:scaledmainalgdi}) has scaling factor $\kappa$ when compared
with (\ref{eqn:mainalgdi}) and corresponds to a scaled version of Algorithm \ref{alg:mainalgorithm}.
The scaling is to enable comparison of (\ref{eqn:scaledmainalgdi}) with the system
(\ref{eqn:price})-(\ref{eqn:cost}) which already has the scaling factor  $\kappa$ in
(\ref{eqn:update}).

\begin{figure}[!ht]
\centering
\includegraphics[draft=false,scale=0.8]{figurescale.eps}
 \caption{(a) Error dynamics for (\ref{eqn:price})-(\ref{eqn:cost}) denoted (``KMT'') with
$\kappa = 1,2,10,50,100$; error dynamics for (\ref{eqn:scaledmainalgdi}) denoted ``Algorithm
\ref{alg:mainalgorithm}'' with $\kappa=1$ is also plotted for comparison. (b) The roles of
(\ref{eqn:price})-(\ref{eqn:cost}) and (\ref{eqn:scaledmainalgdi}) are swapped.}
 \label{fig:figure2}
 \end{figure}

 \begin{figure}[!ht]
\centering
\includegraphics[draft=false,scale=0.8]{figure11050100.eps}
 \caption{Comparison of error dynamics for the systems (\ref{eqn:price})-(\ref{eqn:cost}) (``KMT''
curve) and (\ref{eqn:scaledmainalgdi}) (``Algorithm \ref{alg:mainalgorithm}'' curve). Subplot (a)
has $\kappa=1$, (b) has $\kappa=10$, (c) has $\kappa=50$, and (d) has $\kappa=100$. Note that the
time axis is scaled differently in each of the subfigures.}
 \label{fig:figure3}
 \end{figure}

All figures are for the flow aggregating network with $n=m=10$, $\beta(e)=0.09 \times e, e \in
[n]$, and link capacities $c(l)=10 \times l, l \in [m]$. The initial point for both differential
inclusions is always the lexicographically maximal point\footnote{The lexicographically maximal
point is one where the minimum allocation (across users) is maximized among all feasible points;
further the second minimum is maximized among all points with equal minimum allocation, and so
on.}. This is the most natural starting point when the network does not know the users' utility
functions and considers all users to be equal. While we report the results only for this particular
$\beta$ and $c$, we have simulated several other settings, and the results are qualitatively the
same. We do not repeat them here for brevity.

Figure \ref{fig:figure2}(a) shows that as $\kappa$ scales up, the speed of convergence of the
system (\ref{eqn:price})-(\ref{eqn:cost}) increases. For comparison, we have included the solid
curve for (\ref{eqn:scaledmainalgdi}) with $\kappa=1$.

Figure \ref{fig:figure2}(b) shows that as $\kappa$ scales up, the rate of convergence of
(\ref{eqn:scaledmainalgdi}) also increases similarly. Again, for comparison, we have included the
dashed curve for (\ref{eqn:price})-(\ref{eqn:cost}) with $\kappa=1$.

These two subfigures show that convergence can be sped up similarly in the two systems,
(\ref{eqn:price})-(\ref{eqn:cost}) and (\ref{eqn:scaledmainalgdi}), by simply increasing $\kappa$.

Figures \ref{fig:figure3}(a)-\ref{fig:figure3}(d) compare (\ref{eqn:price})-(\ref{eqn:cost})
directly with (\ref{eqn:scaledmainalgdi}) for identical $\kappa$. The speeding up parameter
$\kappa$ equals 1, 10, 50 and 100 in Figures
\ref{fig:figure3}(a), \ref{fig:figure3}(b), \ref{fig:figure3}(c), and \ref{fig:figure3}(d)
respectively. These figures demonstrate that (\ref{eqn:scaledmainalgdi}) converges faster than
(\ref{eqn:price})-(\ref{eqn:cost}) as long as the $\kappa$ parameter is identical for the two
systems. We saw the same qualitative behavior across several randomly chosen problem parameters.

\section{Conclusion}
\label{sec:conclusion}
We considered a network utility maximization problem in a distributed
framework where the users do not know the network structure or utility functions of other users and
the network does not know the users' utility functions. We decomposed the system problem into user
subproblems and a network subproblem following the methodology of \cite{1998xxJORS_KelMauTan}.
Unlike the dual decomposition iterative methods of \cite{1958xxEJES_ArrHur},
\cite{1998xxJORS_KelMauTan}, \cite{1999xx_LowLap}, etc., the iterations proposed in Algorithm
\ref{alg:mainalgorithm} ensure feasibility at every step. The convergence of the algorithm was
shown
using the theory of differential inclusions. Efficient methods to solve the network problem for
some
special networks were also described. Finally, sample simulations show that, in several examples,
Algorithm \ref{alg:mainalgorithm}'s associated differential inclusion (\ref{eqn:mainalgdi})
converges faster to the system optimal point when compared with the iterates arising from the ODE
(\ref{eqn:price})-(\ref{eqn:cost}).


\section{Introduction}
We show that if complex analytic functions and the complex numbers themselves had not been discovered before the study of the two-dimensional Laplace equation, $\Delta h =  h_{xx} + h_{yy} = 0$, then they would have likely been discovered as a result of it. In this paper, a hypercomplex analytic function is an element of a class of differentiable functions, $C^1 (A)$, over a two-dimensional commutative unital algebra, $A$. Although, hypercomplex analytic functions have appeared in the literature, e.g. \cite{Ketchum},  \cite{Ward2}, we will explore the connections between these atypical theories of differentiable functions and certain partial differential equations. We define $\Gamma_{\alpha,\beta}$ as the partial differential operator,
\begin{center}
$\Gamma_{\alpha,\beta} = \partial_{x}^2 + \alpha \partial_{x}\partial_{y} + \beta \partial_{y}^2$ with $\alpha,\beta \in \mathbb R$.
\end{center}
We have selected $\Gamma_{\alpha,\beta}$ in this manner as it is a generalization of the harmonic operator $\Delta = \partial_x^2 + \partial_y^2$, in the special case, $\Gamma_{0,1}$. Now, we introduce and utilize these new classes of differentiable functions to answer the following important questions;

\begin{itemize}

\item[1.]	What kind of continuous change of variables preserve solutions to the problem,
\begin{equation} \label{Eqn1}
 \Gamma_{\alpha,\beta}h = 0
\end{equation}
for $h=h(x,y)$ and $(x,y) \in \Omega \subseteq \mathbb R^2$?

\item[2.]	What is a general form of solution to the generalized Bi-harmonic Equation,

\begin{equation} \label{Eqn2}
 \Gamma_{\alpha,\beta}^2 h = 0
 \end{equation}
 for $h=h(x,y)$ and $(x,y) \in \Omega \subseteq \mathbb R^2$ \label{Eqn2}

\end{itemize}

In fact, we have chosen to select these questions as they are naturally answered in the case of the harmonic operator via the use of complex analytic functions, thus, rendering the results obtained in the case $\Gamma_{\alpha,\beta}$ true generalizations.

We define the commutative algebra,
\begin{center}
$A_{ \lambda_1 , \lambda_2 } = \{ t + sj | (t,s) \in \mathbb R^2 \}$ where,
$ j^2 = - \lambda_2  - \lambda_1 j $ and $j \not\in \mathbb R$.

\end{center}

The numbers $z = t + sj \in A_{ \lambda_1 , \lambda_2 }$ will have a {\it real} part, $Re_{ \lambda_1 , \lambda_2 }(z) = t$, and an {\it imaginary } part $Im_{ \lambda_1 , \lambda_2 }(z) = s$.

The conjugate will also be defined in the analogous manner,
\begin{center}
$\bar{z}^{ \lambda_1 , \lambda_2 } = t - sj$.
\end{center}

Over these algebras we may define functions,
\begin{center}
$ F : A_{ \lambda_1 , \lambda_2 }  \rightarrow A_{ \lambda_1 , \lambda_2 } $, \\
$ F(x,y) = u(x,y) + jv(x,y), $
\\

\end{center}
where $u,v$ are functions from $\mathbb{R}^2$ into $\mathbb{R}$.
\\
Such functions will be called $ A_{ \lambda_1 , \lambda_2 }$-differentiable if they satisfy the generalized Cauchy-Riemann conditions,
\begin{center}
$ u_x = v_y + \lambda_1 v_x, $ \\
$u_y = - \lambda_2 v_x. $
\end{center}

See e.g. \cite{Ward1} for related work on generalized Cauchy-Riemann conditions.
If the reader replaces $\lambda_1 = 0 , \lambda_2 = 1$, they will obtain the system of complex numbers and complex analytic functions.

\subsection{Main Results}
We will now list the main results which serve as answers to our primary questions.
The following theorem answers question $1$.

\begin{Thm}  \label{thm1}
Let $\Omega$ be a domain in $\mathbb R^2$ and let $h: \Omega \rightarrow \mathbb R$ be a solution to,
\begin{center}
$ \Gamma_{\alpha,\beta}h = 0 $
\end{center}

for all $(x,y) \in \Omega$ and let $F = (u,v): \Omega \rightarrow \mathbb R^2 $. Then, the composition $H = h \circ F$ is also a solution to \ref{Eqn1} if, $F = u + jv $ defines an $ A_{ \frac{\alpha}{\beta},\frac{1}{\beta} }$-differentiable function on the region $\Omega$.
\end{Thm}

The next theorem answers question 2.
\begin{Thm} \label{thm2}
Let $f = u + jv$, and $g = k + jw$ be $ A_{ \frac{\alpha}{\beta},\frac{1}{\beta} }$-differentiable functions on a domain $\Omega $ in  $\mathbb R^2$. Then, $h(x,y) = \Im_{ \lambda_1 , \lambda_2 }( \bar{z}^{ \lambda_1 , \lambda_2 }f(z) + g(z) ),$  is a solution to the equation,
\begin{center}
$ \Gamma_{\alpha,\beta}^2 h  = 0 $
\end{center}
on the region, $\Omega$.

\end{Thm}

Again, if the reader considers the complex-harmonic case, $\lambda_1 = 0 , \lambda_2 = 1$, then they will see familiar theorems, \cite{Charles},\cite{{Rudin}},\cite{Kre}.
\\
The following Lemma is very useful and as well constitutes an informative generalization of complex analytic functions.
\begin{Lemma} \label{lemma1}
Let $F = u + j v$ be an $ A_{ \frac{\alpha}{\beta},\frac{1}{\beta} }$-differentiable function.
\\ Then,
\begin{center}
$\Gamma_{\alpha,\beta}u = \Gamma_{\alpha,\beta}v = 0$.
\end{center}
\end{Lemma}

\begin{proof} \ref{lemma1}
Consider the generalized Cauchy-Riemann condition.
\begin{center}
$ u_x = v_y + \ddfrac{ \alpha }{ \beta } v_x $  and $u_y = - \ddfrac{ 1 }{ \beta } v_x $
\end{center}
and so,
\begin{center}
$ v_{yy} + \ddfrac{ \alpha }{ \beta } v_{xy} = - \ddfrac{ 1 }{ \beta } v_{xx}$
\end{center}
gives,
\begin{center}
$  \Gamma_{\alpha,\beta}v = v_{xx}+ \alpha v_{xy} + \beta v_{yy} = 0$.
\end{center}
Then, similarly, we infer,  $\Gamma_{\alpha,\beta}u = 0$

\end{proof}

\section{Proofs of presented results}

\begin{proof}[ Proof \ref{thm1}]
Begin by calculating the derivatives and grouping terms.

$(u_{xx}+\alpha u_{xy} + \beta u_{yy} )h_u + (u_x^2+\alpha u_x u_y + \beta u_y^2)h_{uu} + (2u_xv_x+\alpha(u_xv_y+u_yv_x) +  2\beta u_y v_y )h_{uv} + (v_x^2 + \alpha v_x v_y + \beta v_y^2 )h_{vv} + (v_{xx} + \alpha v_{xy} + \beta v_{yy} )h_v  = **$
Then, recall by Lemma \ref{lemma1}, both $u$ and $v$ are solutions to the equation \ref{Eqn1}. Then we begin considering each factor of $h_{uu},h_{uv},h_{vv}$. We apply different variations of the Generalized Cauchy-Riemann conditions to see that the terms all simplify nicely.
\begin{itemize}

\item $u_x^2+\alpha u_x u_y + \beta u_y^2 = u_x(u_x+\alpha u_y)-u_y(-\beta u_y) = | J(F) | $

\item $ 2 u_x v_x+\alpha(u_x v_y+u_y v_x) +  2 \beta u_y v_y = \alpha \big ( \ddfrac{2}{ \alpha }(v_y - \alpha u_y )v_x + u_x v_y +  v_x u_y + \ddfrac{ 2 \beta }{ \alpha } ( \ddfrac{-1}{ \beta }v_x) v_y  \big ) = \alpha | J(F) |   $

\item $v_x^2+\alpha v_x v_y + \beta v_y^2 = \beta \big( - ( -\ddfrac{1}{\beta} v_x ) v_x +  ( \ddfrac{\alpha}{\beta}v_x + v_y )v_y \big ) = \beta |J(F) | $

\end{itemize}
Then so,

\begin{center}
$ ** = | J(F) |  h_{uu} + \alpha | J(F) | h_{uv} +\beta |J(F) | h_{vv} = |J(F) | \big (  h_{uu} + \alpha h_{uv} +\beta h_{vv} \big ) = 0 $
\end{center}

Thus completing the proof.
\end{proof}

\begin{proof}[ Proof \ref{thm2}]
Begin by calculating derivatives of $h(x,y) = \Im_{ \lambda_1 , \lambda_2 }( \bar{z}^{ \lambda_1 , \lambda_2 }f(z) + g(z) )$ and grouping terms.
If we denote $z = x+jy$ with $x,y \in \mathbb{R}$, then we first group terms by factors of $y$,
\begin{center}
$\ddfrac{ \alpha }{ \beta }  v^{(4,0)} -  u^{(4,0)} + 2  \ddfrac{ \alpha^2 }{ \beta }  v^{(3,1)}  - 2 \alpha  u^{(3,1)} + (\alpha^2 + 2 \beta)\ddfrac{ \alpha }{ \beta }  v^{(2,2)}  -  (\alpha^2 + 2 \beta) u^{(2,2)} + 2 \alpha^2 v^{(1,3)} -  2 \alpha \beta  u^{(1,3)} + \alpha \beta   v^{(0,4)}  - \beta^2  u^{(0,4)} = \ddfrac{ \alpha }{ \beta } \big (   v^{(4,0)} + 2  \alpha  v^{(3,1)} + (\alpha^2 + 2 \beta) v^{(2,2)}+ 2 \alpha \beta v^{(1,3)}+ \beta ^2  v^{(0,4)} \big )- \big (u^{(4,0)} + 2 \alpha  u^{(3,1)} +  (\alpha^2 + 2 \beta) u^{(2,2)} + 2 \alpha \beta  u^{(1,3)} + \beta^2  u^{(0,4)} \big ) = \ddfrac{ \alpha }{ \beta } \Gamma^2_{ \alpha , \beta } v -  \Gamma^2_{ \alpha , \beta } u = 0 $
\end{center}
In terms of $x$,
\begin{center}
$ v^{(4,0)} + 2 \alpha v^{(3,1)} + (\alpha^2 + 2 \beta )v^{(2,2)} + 2 \alpha \beta v^{(1,3)} + \beta^2 v^{(0,4)} = \Gamma^2_{ \alpha , \beta } v  = 0 $
\end{center}
Then, the function $w$ and its derivatives,
\begin{center}
$ w^{(4,0)} + 2 \alpha w^{(3,1)} + (\alpha^2 + 2 \beta )w^{(2,2)} + 2 \alpha \beta w^{(1,3)} + \beta^2 w^{(0,4)} = \Gamma^2_{ \alpha , \beta } v = 0 $
\end{center}
Finally, the functions $u,v$ and their derivatives with constant factors,
\begin{center}
$ 4v^{(3,0)} + 2 \ddfrac{ \alpha^2 }{ \beta }  v^{(3,0)} - 2 \alpha  u^{(3,0)} + 6 \alpha v^{(2,1)} + 2 \ddfrac{ \alpha^2 }{ \beta } (\alpha^2 + 2\beta) v^{(2,1)} - 2(\alpha^2 + 2\beta) u^{(2,1)} + 2(\alpha^2 + 2\beta) v^{(1,2)} + 6 \alpha^2 v^{(1,2)} - 6 \alpha \beta u^{(1,2)} + 2 \alpha \beta  v^{(0,3)} - 4 \beta^2 u^{(0,3)} + 4 \alpha \beta v^{(0,3)}.$

\medskip

\end{center}

Next we differentiate the Generalized Cauchy Riemann Conditions and substitute $u$,
\begin{center}
$ u^{(1,0)} = v^{(0,1)} +  \ddfrac{ \alpha }{ \beta } v^{(1,0)} $

$ u^{(0,1)} = - \ddfrac{ 1 }{ \beta } v^{(1,0)} $

\end{center}

\begin{center}
$ 4v^{(3,0)} + 2 \ddfrac{ \alpha^2 }{ \beta }  v^{(3,0)} - 2 \alpha  \big(  v^{(2,1)} + \ddfrac{ \alpha }{ \beta } v^{(3,0)}  \big )  + 6 \alpha v^{(2,1)}+ 2(\alpha^2 + 2\beta) \big( v^{(1,2)}  + \ddfrac{ \alpha }{ \beta } v^{(2,1)} \big ) - 2(\alpha^2 + 2\beta) u^{(2,1)} + 6 \alpha^2 v^{(1,2)} - 6 \alpha \beta \big ( v^{(0,3)} + \ddfrac{ \alpha }{ \beta }v^{(1,2)} \big )+ 2 \alpha \beta  v^{(0,3)} + 4 \beta^2 \big( \ddfrac{ 1 }{ \beta }  v^{(1,2)} \big ) + 4 \alpha \beta v^{(0,3)} = 4v^{(3,0)} + 4 \alpha v^{(2,1)} + 4 \beta  v^{(1,2)}  = 4 \partial_x  \big ( v^{(2,0)} +  \alpha v^{(1,1)} +  \beta  v^{(0,2)} \big ) = 4 \partial_x \Gamma_{ \alpha , \beta } v =  0$  .

\end{center}

This completes the proof.
\end{proof}

\section{Concluding Remarks}
We hope that this paper has helped to give perspective on the usefulness of the hyper-complex number system and the corresponding differentiable functions. As well, to put perspective on the place of the complex-analytic functions within the mathematical landscape relative to their, often equally significant, counterparts.