\section{Introduction}
Distribution networks are systems composed of several inventories where goods/materials/products are stored, and transportation lines along which the goods are transported. 
The control problem in distribution networks is to assign the flows to the transportation lines such that an external supply of goods, arriving at some of the inventories, is distributed to inventories in such a way that a requested demand can always be satisfied. 
The static problem of distributing a constant flow to transportation lines according to some cost functions has been at the heart of the mathematical theory called \emph{network optimization} \cite{Ford1956}, \cite{Bertsekas1998}, \cite{Rockafellar1998}.
However, real networks must react dynamically on changes in the external conditions, such as time-varying supply and demand.
If the networks operate on a fast time scale 
continuous feedback controllers are required to adapt the flow instantaneously to external changes.
Beyond that, feedback controllers also promise to provide desirable robustness properties to the network.
In addition to achieving an optimal routing, those controllers should at best also distribute the stored goods equally to all inventories.

Feedback control of inventory systems has recently found significant attention.
The control of distribution networks with unknown supply and demand has been investigated in \cite{Bauso2006}, using a combination of robust optimization and control theory. A control strategy is derived that ensures robust feasibility as well as a min-max type optimality.
Robust control strategies for distribution networks have been considered, e.g., in \cite{Blanchini2001}, \cite{Baric2012}.
The aspect of balancing in distribution networks has been studied in \cite{Danielson2013} using a one-step model predictive control scheme.
In  \cite{vanderSchaft2012},  \cite{Wei2013}, 
it was shown that the balancing problem can be formulated in a Hamiltonian framework.

Network flow and balancing problems are intimately related to consensus or output agreement problems. In \cite{Burger2013a}, \cite{Burger2013} we showed that networks of passive systems with constant external disturbances solve asymptotically certain network flow problems. Furthermore, the phenomena of output agreement and clustering are mathematically related to optimal network flow problems  \cite{Burger2011}.
The connection between output agreement and network flow control has been further studied from an internal model control perspective in \cite{DePersis2013}, \cite{Burger2013b}, \cite{Burger2013c}.
Considering demands and supplies that are generated by an external dynamical system (exo-system), distributed dynamic controllers were proposed that ensure an exact routing as well as a balancing of the inventory levels. In \cite{Burger2013c}, also controllers were proposed that route the flow asymptotically optimal through the network, provided the underlying cost functions are quadratic. 

We address here the problem of balancing and optimal routing  in time-varying distribution networks with convex routing costs and capacity constraints on the transportation lines. 
The time-varying optimal flow problem is formulated, assuming that the supply/demand is generated by an exo-system. 
It is shown that for a certain class of convex cost functions, namely \emph{Legendre-type} functions, the time-varying optimal flow can be characterized in terms of the time-varying dual solutions. 
In particular, an explicit differential equation for the dual solutions is derived. The construction of this differential equation requires to find the analytic solution to a \emph{generalized Dirichlet problem} \cite{Biggs1997}.
Based on this characterization of the time-varying dual variables, a feedback controller is proposed that achieves a balancing of the storage levels and routes the flow asymptotically optimal through the network.

The remainder of the paper is organized as follows. 
The inventory system model and the relevant assumptions on the supply are introduced in Section \ref{sec.OptimalDistribution}.
The time-varying optimal flow problem is presented in Section \ref{sec.OptRouting}, where also a differential equation for the optimal dual variables is derived.
A feedback controller, ensuring balancing and asymptotic optimal routing is proposed in Section \ref{sec.Feedback}. 
The design procedure is illustrated on an example system in Section \ref{sec.Example}. \\

\textbf{Preliminaries:}
The notation we employ is standard. The set of (non negative) real numbers is denoted by $\mathbb{R}$ ($\mathbb{R}_{\geq}$). The distance of a point $q$ from a set $\mathcal{A}$ is defined as $\mbox{dist}_{\mathcal{A}} q = \mbox{inf}_{p \in \mathcal{A}} \|p - q\|$.
Given a matrix $B$ (vector $x$), the notation $[B]_{ik}$ ($[x]_i$) denotes its $ik$-th ($i$-th) entry.
The range-space and null-space of a matrix $B$ are denoted by $\mathcal{R}(B)$ and $\mathcal{N}(B)$, respectively.
A graph $\mathcal{G}=(V,E)$ is an object consisting of a finite set of nodes, $|V| = n$, and edges, $|E| = m$.
Throughout the paper we assume that a graph $\mathcal{G}$ is assigned an (arbitrary) orientation.
The incidence matrix $B \in \mathbb{R}^{n \times m}$ of the
graph $\mathcal{G}$, is a $\{0, \pm 1\}$ matrix with $[B]_{ik}$ having value
`+1' if node $i$ is the initial node of edge $k$, `-1' if it is the terminal
node, and `0' otherwise. 


For a convex function $\mathcal{P}: \mathbb{R}^{p} \rightarrow \mathbb{R}$, the \emph{effective domain} is $\mathrm{dom} \mathcal{P}  = \{ s: \mathcal{P}(s) < \infty\}$. 
The convex conjugate of $\mathcal{P}(s)$ is $\mathcal{P}^{\star}(r) = \sup_{s} \{r^{T}s - \mathcal{P}(s) \; : \; s \in \mathrm{dom} \mathcal{P} \; \}$.




\section{Time-varying distribution problems} \label{sec.OptimalDistribution}

The objective of this paper is to design distribution control laws for transportation networks with storage that achieve two objectives: (i) balancing of the storage levels, and (ii) an optimal routing of the flows trough the network.
Consider an inventory network with $n$ inventories and $m$ transportation lines, described by a graph $\mathcal{G}=(V,E)$, and let $B$ be the incidence matrix (with arbitrary orientation). The dynamics of the inventory system is given as
\begin{align}
\begin{split} \label{sys.Inventory1}
\dot{x}(t) &= B\lambda(t) + Pw(t),
\end{split}
\end{align}
where $x \in \mathbb{R}^{n}$ represents the storage level, $\lambda \in \mathbb{R}^{m}$ the flow along the lines, and $Pw$ an external in-/outflow of the inventories, i.e., the supply or demand. 
This basic model (or its discrete time version) is studied  in \cite{Bauso2006}, \cite{vanderSchaft2012}, \cite{Baric2012}, \cite{Danielson2013}, \cite{DePersis2013} and is a relevant model for supply chains (\cite{Alessandri2011}) or data networks (\cite{Moss1983}).


We assume throughout this paper that the supply and demand are varying over time and in particular that they as well as their derivatives are known, that is, all information about the evolution of the flow is known. In compliance with the control theoretic approach of this paper, we assume that the supply/demand is generated by a dynamical system
\begin{align}\label{sys.Inventory_w} 
\dot{w}(t) = s(w(t)), \quad w \in \mathbb{R}^{r},\end{align}
with initial condition $w(0) = w_{0}$.
\begin{assumption} \label{ass.Boundedness}
There exists a compact set $\mathcal{W}$ that is forward invariant under \eqref{sys.Inventory_w} and the initial conditions $w(0)$ are contained in this set $\mathcal{W}$.
\end{assumption}
Note that this assumption ensures the boundedness of the supply and demand vector.
Throughout this paper we restrict our attention to networks with balanced supply and demand.
\begin{assumption}\label{ass.BalancedFlow}
The supply/demand of the network is balanced, that is 
\begin{align*}
\1_{n}^{T}Pw(t) = 0 \;\mbox{{\rm for all }} t \geq 0.
\end{align*}
\end{assumption}
We assume that the flow on each line can flow in both directions, but is restricted to satisfy certain capacity bounds. That is, we assume that for each line $k \in \{1,\ldots,m\}$ there exists $\bar{\lambda}_{k} \in \mathbb{R}_{>0}$ and the flow must satisfy
\begin{align} \label{eqn.CapacityConstraints}
- \bar{\lambda}_{k} \leq \lambda_{k}(t) \leq \bar{\lambda}_{k} \quad \mbox{for\;all\;} t \geq 0.
\end{align}
If the flow on one line is unconstrained, we let $\bar{\lambda}_{k}$ be infinity. In the following, we will also write $- \bar{\lambda} \leq \lambda \leq \bar{\lambda}$, meaning that the inequalities hold for each component.


A \emph{cut} of the network is defined to be a signed edge set\footnote{We assume that the direction of the edges is chosen in accordance with the definition of the incidence matrix $B$.} $Q \subset E$ (with positive and negative parts denoted by $Q^{+}$ and $Q^{-}$) such that, for some node set $S \subset V$ 
\begin{align}
\begin{split} \label{eqn.InducedCut}
Q^{+} = \{ k \in E \; :\; k = (i,j)\; \mbox{with}\; i \in S, j \in V \setminus S\} \\
Q^{-} = \{ k \in E \; :\; k = (j,i)\; \mbox{with}\; i \in S, j \in V \setminus S\}.
\end{split}
\end{align}
We will denote the cut induced by a set of nodes $S \subset V$ with $Q(S)$.
The \emph{capacity} of a cut is
\begin{align*}
\bar{\lambda}_{Q}= \sum_{k \in Q} \bar{\lambda}_{k}.
\end{align*} 

To ensure feasibility of the problem, we impose the following condition.
\begin{assumption} \label{ass.Cut}
There exists $\epsilon > 0$ such that for any set of nodes $S \subseteq V$ it holds that
\begin{align*}
-\bar{\lambda}_{Q(S)} + \epsilon \leq \sum_{i \in S} [Pw(t)]_{i}  \leq \bar{\lambda}_{Q(S)} - \epsilon.
\end{align*}
\end{assumption}
The assumption ensures that the ``in/out-flow'' to any set of nodes is strictly smaller than the capacity of the corresponding induced cut.
The following result summarizes the conditions that ensure feasibility of the distribution problem.

 
\begin{proposition}[Feasibility] \label{prop.Feasibility}
Consider the inventory system \eqref{sys.Inventory1} with a time-varying supply/demand generated by the exo-system \eqref{sys.Inventory_w} and capacity constraints \eqref{eqn.CapacityConstraints}. Suppose Assumptions \ref{ass.Boundedness}, \ref{ass.BalancedFlow} and \ref{ass.Cut} hold. Then,
(i) for any flow vector $\lambda \in \mathbb{R}^{m}$ the average inventory level  $x^{\mathrm{avg}}(t) = \frac{1}{n} \1_{n}^{T}x(t)$ remains constant, that is $x^{\mathrm{avg}}(t) = x^{\mathrm{avg}}(0)$ for all $t \geq 0$; and 
(ii) there exists a continuous flow trajectory $\lambda : \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}^{m}$ that is strictly feasible, that means there is a vector valued function $\lambda : \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}^{m}$ such that for each edge $k$ the bounds \eqref{eqn.CapacityConstraints} hold with strict inequality and $B \lambda(t) + Pw(t) = 0$. 
\end{proposition} 
 
 
\begin{proof}
To see the first claim note that the average $x^{\mathrm{avg}}(t)$ evolves as 
$\dot{x}^{\mathrm{avg}} = \frac{1}{n}\1_{n}^{T}B\lambda+\frac{1}{n}\1_{n}^{T}Pw(t) = 0.$
To prove the second claim, we show first that at each fixed time $t$, there exists a feasible flow vector. To indicate that we consider a single time instant, we write in the following $w = w(t)$.
Assume without loss of generality that all transportation lines have a finite capacity, i.e., $\bar{\lambda}_{k} < \infty$.
Let $\Lambda = \mathrm{diag} \{ \bar{\lambda}_{1}, \ldots, \bar{\lambda}_{m}\}$.
Consider the linear optimization problem
\begin{align}
\begin{split}  \label{prob.LPPrimal}
\min_{\lambda} & \; \| \Lambda^{-1}\lambda\|_{\infty}\quad \mathrm{s.t.} \; B\lambda = - Pw.
\end{split}
\end{align}
The linear programming dual can be derived  as
\begin{align}
\begin{split} \label{prob.LPDual}
\max_{u} & \; (-Pw)^{T}u \quad \mathrm{s.t.} \;   \| \Lambda B^{T}u\|_{1} \leq 1.
\end{split}
\end{align}
 Note now that the extreme points of the polytope defined by $\| \Lambda B^{T}u\|_{1} \leq 1$ are multiples of the vectors whose elements are zeros and ones (see e.g., \cite{Strang1984}).
 %
  Thus, each extreme point must be of the form $u = c_{Q} e_{Q}$, where $Q$ is a cut induced by some node set $S \subseteq V$ according to \eqref{eqn.InducedCut}, $c_{Q} \in \mathbb{R}$ is a constant such that $|c_{Q}|= \bar{\lambda}_{Q}^{-1}$, and $e_{Q} \in \mathbb{R}^{n}$ is a vector with entries being $[e_{Q}]_{i} =1$ if node $i \in S$ and $[e_{Q}]_{i} =0$ if $i \in V \setminus S$.
 %
It follows now from Assumption \ref{ass.Cut} that \eqref{prob.LPDual} must always have an optimal value strictly smaller than one. From strong duality of linear programs follows that also \eqref{prob.LPPrimal} must have an optimal value strictly smaller than one. 
This shows that at each fixed time instant $t$, there exists a feasible flow vector $\lambda(t)$ that satisfies all constraints, i.e., $- \bar{\lambda} < \lambda(t) < \bar{\lambda}$. 

It remains to prove that there exists a continuous function $\lambda : \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}^{m}$ that represents a feasible flow trajectory.
 %
Consider now the set of feasible flows for a given supply/demand vector $w$, i.e.,
\[ \Delta(w) = \{ \lambda \in \mathbb{R}^{m} \; : \; B \lambda = - Pw, \, - \bar{\lambda} \leq \lambda \leq \bar{\lambda} \}.\]
The set $\Delta(w)$ is compact, convex and has a non-empty relative interior.
One can understand $\Delta$ as a \emph{set-valued map} from $\mathbb{R}^{r}$ to $\mathbb{R}^{m}$. 
As a set-valued map, $\Delta$ is \emph{lower semi-continuous}.\footnote{ The set valued function $\Delta$ is lower semi-continuous at a point $w \in \mathrm{dom}\Delta$ if and only if for any $\lambda \in \Delta(w)$ and any sequence $w_{k}$ converging to $w$, there exists a sequence $\lambda_{k} \in \Delta(w_{k})$ converging to $\lambda$. It is said to be lower semi-continuous if it is lower semi-continuous at every point $w \in \mathrm{dom} \Delta$. }
The existence of a continuous selection $\lambda(w)$ of $\Delta$ follows now from Micheal's Theorem (\cite[Theorem 9.1.2]{Aubin1990}).
Since additionally $w(t)$ is a continuous function of time (i.e., it respects the dynamics \eqref{sys.Inventory_w}), there exists a continuous function $\lambda(t)$  such that $\lambda(t) \in \Delta(w(t))$, i.e., there is a continuous feasible flow trajectory.
\end{proof} 
 
 
 %

 
 
\section{Optimal Time-Varying Routing} \label{sec.OptRouting}

In the following, we assume that the routing of flow along transportation lines of the network incorporates a cost. We seek a control strategy to route the time-varying flow in such a way through the network that the routing cost is minimal at each time instant. We characterize next a class of feed-forward controllers that achieve such an optimal routing. 


The flow cost on a transportation line is described by a convex $C^{\infty}$ cost function
 \[ \mathcal{P}_{k}(\lambda_{k}), \, k = 1,2,\ldots,m \]
  for the flow $\lambda_{k}$. 
 
We aim to characterize the steady state flows that are such that at each time instant the  following \emph{static optimal distribution problem} is solved
\begin{align}
\begin{split} \label{prob.OFP}
\lambda^{w} = \mbox{arg}\min \;& \sum_{k=1}^{m}\mathcal{P}_{k}(\lambda_{k}) \\
& 0 = B \lambda + Pw,
\end{split}
\end{align}
where $ w = w(t)$ is the supply/demand at time $t$. We will use in the following the short hand notation $\mathcal{P}(\lambda) = \sum_{k=1}^{m} \mathcal{P}_{k}(\lambda_{k})$. 

The convex cost functions are such that they incorporate the capacity constraints of the transportation lines, in the sense that $P_{k}(\lambda_{k}) = \infty$ for $|\lambda_{k}| > \bar{\lambda}_{k}$.
More precisely, we assume that the cost functions are of \emph{Legendre type} \cite[p. 258]{Rockafellar1997}.

\begin{assumption}\label{ass.Legendre}
The functions $\mathcal{P}_{k}$, $k \in \{1,\ldots,m\}$ are strictly convex and differentiable on the open set $\{ \lambda_{k} \; : \; - \bar{\lambda}_{k} < \lambda_{k}(t) < \bar{\lambda}_{k} \}$, and $\lim_{\ell \rightarrow \infty} |\nabla \mathcal{P}(\lambda^{[\ell]}_{k})| = + \infty$ for any sequence $\lambda_{k}^{[1]},\lambda_{k}^{[2]},\ldots$ such that $|\lambda_{k}^{[\ell]} | < \bar{\lambda}_{k}$ and $\lim_{\ell \rightarrow \infty} |\lambda_{k}^{[\ell]}| = \bar{\lambda}_{k}$.
\end{assumption}
With cost functions of this form, any finite solution to \eqref{prob.OFP}  represents a feasible flow. 
Note that if the assumptions of Proposition \ref{prop.Feasibility} hold, then the optimization problem \eqref{prob.OFP} has always a finite optimal solution.




We can now characterize the optimality conditions of \eqref{prob.OFP}. The Lagrangian function associated to \eqref{prob.OFP} with Lagrange multiplier $\zeta$ is 
\[ \mathcal{L}(\lambda,\zeta) = \mathcal{P}(\lambda) + \zeta^{T}(-B \lambda - P w).\]
Now,  $(\lambda^{w},\zeta^{w})$ is an optimal primal/dual solution pair to \eqref{prob.OFP} if the following nonlinear equations hold
\begin{align}
\begin{split} \label{eqn.OptimalityConditions}
\nabla \mathcal{P}(\lambda^{w}) - B^{T}\zeta^{w} &= 0 \\
B \lambda^{w} + Pw &= 0.
\end{split}
\end{align}
Note that the first condition simply implies $\nabla \mathcal{P}(\lambda^{w}) \in \mathcal{R}(B^{T})$ since $\zeta^{w}$ has no further constraints. 
We define the set of optimal routing/supply pairs as
\begin{align*}
\begin{split}
\Gamma = \{ (\lambda,w) \in \mathbb{R}^{m} \times \mathcal{W} \; : &\; \nabla \mathcal{P}(\lambda) \in \mathcal{R}(B^{T}), \\
&\; B \lambda + Pw = 0  \}.
\end{split}
\end{align*}
In particular, $(\lambda^{w},w) \in \Gamma$ if and only if $\lambda^{w}$ is an optimal solution to the static optimal distribution problem \eqref{prob.OFP} with the supply vector $w$.


\begin{proposition}[Feed-forward Controller] \label{eqn.FFfeasibility}
Suppose Assumptions  \ref{ass.Boundedness}, \ref{ass.BalancedFlow}, \ref{ass.Cut} and \ref{ass.Legendre} hold. Then, there exists functions $\Phi : \mathbb{R}^{n} \times \mathbb{R}^{r} \rightarrow \mathbb{R}^{n}$ and $\Psi : \mathbb{R}^{n} \rightarrow \mathbb{R}^{m}$, such that 
for any $w(t)$ being a solution to \eqref{sys.Inventory_w} there exists $\zeta^{w}(0)$ such that the solution $\zeta^{w}(t)$ to the dynamical system
\begin{align}
\begin{split} \label{eqn.FF}
\dot{\zeta} &= \Phi(\zeta,w),  \\
\lambda &= \Psi(\zeta)
\end{split}
\end{align}
originating at $\zeta^{w}(0)$ is such that $\lambda^{w}(t) = \Psi(\zeta^{w}(t))$ satisfies $(\lambda^{w}(t),w(t)) \in \Gamma$, where $w(t)$ is the solution to \eqref{sys.Inventory_w} originating at $w(0) = w_{0}$.
\end{proposition}


The proof of this result provides a method to construct the dynamical system \eqref{eqn.FF}. However, before we present the proof,  some additional discussion is required. \\

The main difficulty associated to the proposed characterization of the set $\Gamma$ relates to the constraint $\nabla \mathcal{P}(\lambda) \in \mathcal{R}(B^{T})$. To avoid this constraint, we express the optimality conditions in terms of the \emph{dual solutions} $\zeta$. 
Under the given assumption that the cost functions $\mathcal{P}_{k}$ are of Legendre type,  the gradients $\nabla \mathcal{P}_{k}$ are invertible on the open sets $\{\lambda_{k} \; : \; -\bar{\lambda}_{k} < \lambda_{k} < \bar{\lambda}_{k}\}$ \cite{Rockafellar1998}. Since the optimal solution must always be strictly feasible, the two optimality conditions \eqref{eqn.OptimalityConditions} can  be combined to
\begin{align} \label{eqn.DualCond} B \nabla \mathcal{P}^{-1}(B^{T}\zeta^{w}) + Pw = 0.\end{align}
Bearing this in mind, we define the set of optimal dual solutions as
\[
\Gamma_{D} = \{ (\zeta,w) \in \mathbb{R}^{n} \times \mathcal{W} : B \nabla \mathcal{P}^{-1}(B^{T}\zeta) + Pw = 0 \}.
\] 
We want to emphasize two properties of $\Gamma_{D}$. First, if $(\zeta^{w},w) \in \Gamma_{D}$ then $(\zeta^{w}+ c\1,w) \in \Gamma_{D}$ for any $c \in \mathbb{R}$. Second, $(\zeta^{w},w) \in \Gamma_{D}$ if any only if the routing strategy $\lambda^{w} = \nabla \mathcal{P}^{-1}(B^{T}\zeta^{w})$ satisfies $(\lambda^{w},w) \in \Gamma$. 
\begin{remark}
In the considered framework, the inverse of the gradient of $\mathcal{P}_{k}$ is the gradient of the convex conjugate $\mathcal{P}_{k}^{\star}$, i.e., $\nabla \mathcal{P}^{-1}(\cdot) = \nabla \mathcal{P}^{\star}(\cdot)$ \cite{Rockafellar1997}. 
For completeness, we want to point out that the dual problem to \eqref{prob.OFP} is
\begin{align}
\begin{split} \label{prob.OPP}
\min_{z, \zeta}  \; &\sum_{k=1}^{m} \mathcal{P}^{\star}_{k}(z_{k}) + w(t)^{T}P^{T}\zeta \\
& z = B^{T}\zeta,
\end{split}
\end{align}
where $\mathcal{P}^{\star}_{k}$ are the convex conjugates of $\mathcal{P}_{k}$. Any optimal solution $\zeta^{w}$ to \eqref{prob.OPP} is such that $(\zeta^{w},w(t)) \in \Gamma_{D}$.
\end{remark}
The main idea pursued here is to design the feed-forward controller for the \emph{dual variables}.
We need to design a dynamics for $\zeta(t)$ such that the constraint \eqref{eqn.DualCond} is satisfied for all times.
Suppose now \eqref{eqn.DualCond} is feasible at the initial time, that is
\begin{align} \label{eqn.vwDynamics0}
B \nabla \mathcal{P}^{-1}(B^{T}\zeta^{w}(0)) = -Pw(0). \;    
\end{align}
Then, $\dot{\zeta}^{w}(t)$ must be such that
\begin{align} \label{eqn.vwDynamics}
\left(B \mathrm{D} \mathcal{P}^{-1}(B^{T}\zeta^{w}(t)) B^{T}\right) \dot{\zeta}^{w} = - P\dot{w}(t),
\end{align}
where $\mathrm{D} \mathcal{P}^{-1}(\cdot) = \mathrm{diag} [\nabla^{2}\mathcal{P}_{1}^{-1}(\cdot), \ldots, \nabla^{2} \mathcal{P}_{m}^{-1}(\cdot)] $. Note that $\mathrm{D} \mathcal{P}^{-1}(\cdot)$ is the Hessian matrix of $\mathcal{P}^{\star}$. 

The equation \eqref{eqn.vwDynamics} defines an implicit differential equation for $\zeta^{w}$. For the design of the feed-forward controller \eqref{eqn.FF}, we need to transform it into an explicit differential equation.

Under the stated assumptions $\mathrm{D} \mathcal{P}^{-1}(\cdot)$ is a diagonal matrix, with all entries on the diagonal being positive. Thus, the matrix 
\begin{align}\label{eqn.WeightedLaplacian1}
 L(\zeta^{w}) := B \mathrm{D} \mathcal{P}^{-1}(B^{T}\zeta^{w}) B^{T}
 \end{align}
 is a weighted, state dependent \emph{Laplacian matrix} of the network, with the weights on the edges being the second derivatives of the cost function's convex conjugates. 
 %
The Laplacian matrix has always one eigenvalue at zero and is therefore not invertible. To find an explicit differential equation for $\zeta^{w}$, we make the following observation. 
 Since $\1_{n}^{T}Pw(t) = 0$ there exists always $\omega(t) \in \mathbb{R}^{m}$ such that
\[ Pw(t) = B \omega(t). \]
One can always define $\omega(t) = B^{T}(BB^{T})^{\dagger}Pw(t)$, where $(BB^{T})^{\dagger}$ is the \emph{generalized inverse} of the matrix $BB^{T}$.\footnote{By the properties of the generalized inverse, it is readily verified that $B\omega(t) = BB^{T}(BB^{T})^{\dagger}Pw(t) = (I- \frac{1}{n}\1_{n}\1_{n}^{T})Pw(t) =Pw(t)$. }
Now, the conditions \eqref{eqn.vwDynamics} can be written as
\begin{align} \label{eqn.GenDirichlet}
L(\zeta^{w}(t)) \dot{\zeta}^{w}(t) = - B\dot{\omega}(t).
\end{align}
Following \cite[Chap. 16]{Biggs1997}, we call the problem to find $\dot{\zeta}^{w}$ a \emph{generalized Dirichlet problem}. 

It is described in \cite[p. 659]{Biggs1997}, how an analytic solution for $\dot{\zeta}^{w}$ in \eqref{eqn.GenDirichlet}  can be constructed.
In particular, one can find analytically a matrix $X(\zeta^{w}) \in \mathbb{R}^{m \times n}$ such that
\begin{align} \label{eqn.TreeSol}
 \dot{\zeta}^{w}(t)  = -X(\zeta^{w}(t))^{T} \mathrm{D} \mathcal{P}(B^{T}\zeta^{w}(t))\dot{\omega}
\end{align}
solves \eqref{eqn.vwDynamics}. The algorithmic procedure for the construction of $X(\zeta^{w})$ is provided in \cite[p. 659]{Biggs1997} and is reviewed in the appendix. The construction of $X(\zeta^{w})$ has an explicit graph theoretic interpretation and the matrix has a characteristic structure.
The entries of $X(\zeta^{w})$ are fractions of polynomials, with each polynomial being the sum of products of different functions $\nabla^{2}\mathcal{P}_{k}^{-1}(\cdot)$.
Since, by assumption, all $\nabla^{2}\mathcal{P}_{k}^{-1}(\cdot)$ are strictly positive and continuously differentiable,  $X(\zeta^{w})$ is continuously differentiable.
Note that if the assumptions of Proposition \ref{prop.Feasibility} hold, the optimal flow problem \eqref{prob.OFP} has at each time instant a feasible solution, and thus admits a finite dual solution.
We are now ready to prove Proposition \ref{eqn.FFfeasibility}. \\


\begin{proof} \textit{(of Proposition \ref{eqn.FFfeasibility})}
We prove the proposition by constructing the corresponding feed-forward controller.  Consider the dynamical system
 \begin{align}
\begin{split} \label{sys.FFcontrollerRouting}
 \dot{\zeta} &=  -X(\zeta)^{T} \mathrm{D} \mathcal{P}(B^{T}\zeta)B^{T}(BB^{T})^{\dagger} Ps(w) \\
 \lambda &= \nabla \mathcal{P}^{-1}(B^{T}\zeta),
\end{split}
\end{align}
where $X(\zeta)$ is as constructed in \eqref{eqn.TreeSol} and $(BB^{T})^{\dagger}$ is the generalized inverse of $BB^{T}$.
If the initial condition $\zeta(0)$ is chosen as an optimal dual solution at time $t=0$, that is $(\zeta(0),w(0)) \in \Gamma_{D}$, then the solution $\zeta(t)$ to \eqref{sys.FFcontrollerRouting} originating at $\zeta(0)$ will be a solution to \eqref{eqn.vwDynamics}, and will therefore be an optimal dual solution for all times, i.e., $(\zeta(t),w(t)) \in \Gamma_{D}$ for all $t \geq 0$. 
Furthermore, if $\zeta(t)$ is an optimal dual solution at time $t$, then  $\lambda(t)= \nabla \mathcal{P}^{-1}(B^{T}\zeta(t))$ is an optimal primal solution, i.e., $(\lambda(t),w(t) )\in \Gamma$.
Thus, Proposition \ref{eqn.FFfeasibility} follows with $\Phi(\zeta,w) =  -X(\zeta)^{T} \mathrm{D} \mathcal{P}(B^{T}\zeta)B^{T}(BB^{T})^{\dagger} Ps(w)$, $\Psi(\zeta) = \nabla \mathcal{P}^{-1}(B^{T}\zeta)$ and $\zeta(0)$ such that $(\zeta(0),w(0)) \in \Gamma_{D}$.
\end{proof}


\section{Feedback Controller Design} \label{sec.Feedback}

We can use the proposed feed-forward controller now to design a dynamic feedback controller for the inventory system \eqref{sys.Inventory1},  that ensures convergence to the optimal routing without being initialized with an optimal solution. 
In addition to the optimal routing, the feedback controller should achieve a \emph{balancing} of the inventory levels of \eqref{sys.Inventory1}.
Network balancing is widely accepted as control objective in storage networks, as it promises to make best use of the available storage, \cite{Danielson2013}, \cite{DePersis2013}. 
\begin{definition}
The \emph{time-varying optimal distribution problem} is solvable for the system \eqref{sys.Inventory1}, if there exists a controller of the form
\begin{align} \label{sys.ControllerBasic}
\begin{split}
\dot{\zeta} &= \Phi(\zeta,w,x) \\
\lambda &= \Psi(\zeta,x)
\end{split}
\end{align} 
such that any solution $(x(t),w(t),\zeta(t))$ of the closed-loop system \eqref{sys.Inventory1}, \eqref{sys.Inventory_w}, \eqref{sys.ControllerBasic} is bounded and satisfies
(i) $\lim_{t \rightarrow \infty} B^{T}x(t) = 0$; and
(ii)  $\lim_{t\rightarrow \infty} \mathrm{dist}_{\Gamma}\bigl( \lambda(t),w(t) \bigr) = 0.$ 
\end{definition}

Note that the flow constraints are treated here as soft constraints, that must be satisfied only asymptotically. 
Furthermore, the proposed problem formulation allows to use the supply/demand vector $w$ directly in the controller. 



We will construct a controller \eqref{sys.ControllerBasic} solving the time-varying optimal distribution problem on the basis of the feed-forward controller designed in the previous section.
For the sake of a compact presentation, we define the matrix
\begin{align} \label{eqn.Xi}
 \Xi(\zeta):= X(\zeta)^{T} \mathrm{D} \mathcal{P}(B^{T}\zeta)B^{T}(BB^{T})^{\dagger}. 
 \end{align}
Furthermore, the projection matrix for the space orthogonal to the all ones vector is
\begin{align} \label{eqn.Projection}
Q = (I_{n} - \frac{1}{n}\1_{n}\1_{n}^{T}).
\end{align}
For designing the feedback controller, we need to relax the original assumption on the feasibility of the flows (i.e., Assumption \ref{ass.Cut}).
\begin{assumption} \label{ass.CutInfty}
For any set of nodes $S \subset V$ it holds that $\bar{\lambda}_{Q(S)} = \infty$.
\end{assumption}
Note that this assumption does not imply that all lines of the network have an infinite capacity, but rather that there exists a spanning tree in the network with infinite capacity.
 %
 The main result of this paper is a controller that routes the flow asymptotically optimal through the network.
\begin{theorem}[Feedback Controller]
Consider the inventory dynamics \eqref{sys.Inventory1} with the supply generated by the dynamical system \eqref{sys.Inventory_w}. 
Suppose the Assumptions \ref{ass.Boundedness}, \ref{ass.BalancedFlow}, \ref{ass.Legendre}, and \ref{ass.CutInfty} hold.
Then,  the dynamic controller
\begin{align}
\begin{split} \label{sys.PriceController}
\dot{\zeta} &= - Q\Xi(\zeta)\Bigl( Ps(w)  + x+  k \bigl(B \nabla \mathcal{P}^{-1}(B^{T}\zeta) + Pw \bigr) \Bigr) \\
\lambda & = \nabla \mathcal{P}^{-1}(B^{T}\zeta) - B^{T}x.
\end{split}
\end{align}
with $k > 0$ a scalar controller gain, solves the time-varying optimal distribution problem for any initial inventory levels $x(0)$ and any controller initialization $\zeta(0)$, that is, all trajectories remain bounded, and (i) $\lim_{t \rightarrow \infty} B^{T} x  = 0$ and (ii) $\lim_{t\rightarrow \infty} \mathrm{dist}_{\Gamma}\bigl( \lambda(t),w(t) \bigr) = 0.$
\end{theorem}

\begin{proof}
We recall that, due to the assumption of balanced supply and demand, it must hold that 
\begin{align} \label{eqn.xw}
x^{w}(t) = x_{0}^{w}\1_{n},\; \mbox{with}\; x_{0}^{w}= \frac{1}{n} \1_{n}^{T}x(0).
\end{align}
Note that the average inventory level is constant over time since $\1_{n}^{T}\dot{x} = 0$ and therefore
\begin{align*}
(x(t) - x^{w}(t)) \in \mathcal{R}(B).
\end{align*}
Furthermore, we have that the average of the controller variables $\frac{1}{n}\1_{n}^{T}\zeta(t)$ is constant since $\1_{n}^{T}\dot{\zeta} = 0$.
The optimal controller (i.e., dual) variables satisfy
\begin{align} \label{eqn.zetaw}
 B \nabla \mathcal{P}^{-1}(B^{T}\zeta^{w}(t) ) + Pw(t) = 0. 
 \end{align}

To prove the desired convergence result, consider a composed storage function of the form
\begin{align}
U(x-x^{w},\zeta,w) = V(x-x^{w}) + W(\zeta,w).
\end{align}
As the first component of the storage function, we consider
\begin{align}
V(x-x^{w}) = \frac{1}{2}\|x - x^{w}\|^{2}.
\end{align}
This function satisfies
\begin{align*}
\dot{V} =& (x - x^{w})^{T}(B \lambda + Pw) \\
=& -(x-x^{w})^{T}BB^{T}(x-x^{w})   \\
&+(x-x^{w})^{T}\left(B\nabla \mathcal{P}^{-1}(B^{T}\zeta) + Pw) \right),
\end{align*}
where we exploited the two relations \eqref{eqn.xw} and \eqref{eqn.zetaw}.

Now, as the second component of $U$ we consider the positive semidefinite function
\begin{align} \label{eqn.W}
W(\zeta,w) =  \frac{1}{2}\|B\nabla \mathcal{P}^{-1}(B^{T}\zeta) + Pw\|^{2}.
\end{align}
The function satisfies $W(\zeta,w) = 0$ if and only if $(\zeta,w) \in \Gamma_{D}$, that is if $\zeta = \zeta^{w}$.

In order to keep the notation compact, we introduce the vector valued function
\begin{align} \label{eqn.Z}
Z(\zeta,w) := B\nabla \mathcal{P}^{-1}(B^{T}\zeta) + Pw.
\end{align}

The directional derivative of $W$ satisfies
\begin{align*}
&\dot{W}=   Z(\zeta,w)^{T}\bigl(B \mathrm{D} \mathcal{P}^{-1}(B^{T}\zeta) B^{T}\dot{\zeta} + Ps(w)    \bigr) \\
&= Z(\zeta,w)^{T}\Bigl( - L(\zeta) \Xi(\zeta)\bigl( Ps(w) + x  + k Z(\zeta,w) \bigr)  + Ps(w) \Bigr),
\end{align*}
where we defined analogous to \eqref{eqn.WeightedLaplacian1} 
\[ L(\zeta) := \bigl(B \mathrm{D} \mathcal{P}^{-1}(B^{T}\zeta) B^{T} \bigr) \]
and used the fact that $L(\zeta)Q = L(\zeta).$

Recalling the definition of $\Xi(\zeta)$ in \eqref{eqn.Xi} and the definition of $X(\zeta)$ provided in the appendix, and using that $Ps(w) \in \mathcal{R}(B)$ it follows that  
\begin{align*}
 L(\zeta) \Xi(\zeta)Ps(w) &= L(\zeta)X(\zeta)^{T} \mathrm{D} \mathcal{P}(B^{T}\zeta)B^{T}(BB^{T})^{\dagger}Ps(w) \\
&= BB^{T}(BB^{T})^{\dagger}Ps(w) \\
&=(I - \frac{1}{n}\1_{n}\1_{n}^{T})Ps(w) \\
&= Ps(w).
\end{align*}

Futhermore, we have
\begin{align*}
 L(\zeta) \Xi(\zeta)x &= L(\zeta)X(\zeta) \mathrm{D} \mathcal{P}(B^{T}\zeta)B^{T}(BB^{T})^{\dagger}x \\
 &=L(\zeta)X(\zeta) \mathrm{D} \mathcal{P}(B^{T}\zeta)B^{T}(BB^{T})^{\dagger}(x -x^{w}) \\
 &=B B^{T}(BB^{T})^{\dagger}(x -x^{w}) \\
 &= (I - \frac{1}{n}\1_{n}\1_{n}^{T})(x -x^{w}) \\
 &= (x - x^{w}).
\end{align*}
To derive the previous conclusion, we exploited that $(BB^{T})^{\dagger}\1_{n} = 0$ (see \cite[Lemma 2]{Gutman2004}) and that $(x - x^{w}) \in \mathcal{R}(B)$.
Along the same line of argumentation, and by exploiting that $\1^{T}_{n}Z(\zeta,w) = 0$, we obtain
\[ L(\zeta) \Xi(\zeta)Z(\zeta,w) = Z(\zeta,w). \]
After these considerations it follows that
\begin{align*}
\dot{W}=& -Z(\zeta,w)^{T}(x - x^{w}) -k Z(\zeta,w)^{T} Z(\zeta,w).
\end{align*}
To complete the proof of convergence we consider now the composed (positive semidefinite) storage function
\begin{align*} 
U(x - x^{w},\zeta,w)& =  V(x-x^{w}) + W(\zeta,w) \\
&= \frac{1}{2}\|x - x^{w}\|^{2} + \frac{1}{2}\|Z(\zeta,w)\|^{2}.
\end{align*}
The derivative of this function satisfies
\begin{align}\label{eqn.Udot}
\begin{split}
\dot{U} =& -(x-x^{w})^{T}BB^{T}(x-x^{w})    + (x - x^{w})^{T}Z(\zeta,w) \\
&-Z(\zeta,w)^{T}(x - x^{w}) -k Z(\zeta,w)^{T}  Z(\zeta,w) \\
= & -(x-x^{w})^{T}BB^{T}(x-x^{w})   -k Z(\zeta,w)^{T}Z(\zeta,w).
\end{split}
\end{align}
Thus, the function $U$ is non-increasing under the closed-loop dynamics. 
Note that $\1^{T}(x-x^{w})= 0$ and therefore it holds that
\begin{align}
(x-x^{w})^{T}BB^{T}(x-x^{w}) \geq s_{2}\|x-x^{w}\|^{2}
\end{align}
where $s_{2} > 0$ is the smallest non-zero eigenvalue of the Laplacian matrix $BB^{T}$.\footnote{The scalar $s_{2}$ is known as the \emph{algebraic connectivity} of the network.} 

The derivative of $U(x-x^{w},\zeta,w)$ satisfies therefore
\begin{align}
\dot{U} \leq - \sigma U,
\end{align}
where $\sigma = \min \{s_{2},k\}$.
It follows now from the Comparison Lemma \cite[Lemma 3.4]{Khalil2002} that
\begin{align*}
U\bigr(x(t)-x^{w}(t),\zeta(t),w(t) \bigl)  \\
\leq \mathrm{e}^{-\sigma t}U\bigr(x(0)-x^{w}(0),\zeta(0),w(0) \bigl),
\end{align*}
and consequently
$\lim_{t \rightarrow \infty} U\bigr(x(t)-x^{w}(t),\zeta(t),w(t) \bigl)  = 0. $
This proves that $\lim_{t\rightarrow \infty} \| x(t)- x^{w}(t)\| = 0$, which implies $\lim_{t \rightarrow \infty} B^{T}x(t) = 0$ and $\lim_{t \rightarrow \infty} Z(\zeta(t),w(t)) = 0$. These conditions imply $\lim_{t \rightarrow \infty} \mathrm{dist}_{\Gamma}\bigl( \lambda(t),w(t) \bigr) = 0$. 

It remains to show boundedness of all trajectories. The uniform boundedness of $U\bigr( x(t)-x^{w},\zeta(t),w(t) \bigl)$ over time implies directly that $V = \|x(t)-x^{w}\|^{2}$ must remain bounded, and consequently $x(t)$ must remain bounded.
Now, the boundedness of $U$ implies boundedness of $W = \|B \nabla \mathcal{P}^{-1}(B^{T}\zeta) -Pw(t)\|^{2}$. By Assumption \ref{ass.Boundedness}, $w(t)$ is bounded and consequently also $\|B \nabla \mathcal{P}^{-1}(B^{T}\zeta(t))\|$ must remain bounded.
Thus, it holds that at each time instant there exists a vector $c(t)$ with finite norm such that
\begin{align*}
B \nabla \mathcal{P}^{-1}(B^{T}\zeta(t)) = c(t).
\end{align*}
It follows with the same argumentation as used previously to derive \eqref{eqn.DualCond}, that at each time instant the vector $\zeta(t)$ is an optimal dual solution to a network flow problem of the form
\begin{align*}
\min \; &\sum_{k=1}^{m} P_{k}(\tilde{\lambda}_{k}(t))\quad \mathrm{s.t.} \; B\tilde{\lambda}(t) = c(t). 
\end{align*}
Due to Assumption \ref{ass.CutInfty} the problem has a feasible solution, and due to strict convexity of the cost functions (Assumption \ref{ass.Legendre}) the optimal solution is unique.
From the optimality conditions follows that there exists at each time instant a unique $z(t)$ such that $\nabla \mathcal{P}(\tilde{\lambda}(t)) = z(t)$ and $z(t) = B^{T}\zeta(t)$. 
Since additionally at each time instant $\1^{T}\zeta(t) = \1^{T}\zeta(0)$, the vector $\zeta(t)$ is uniquely defined and must be bounded. This proves that all trajectories remain bounded.
\end{proof}

\begin{remark}
The controller  \eqref{sys.PriceController} has a characteristic structure. The matrix  $\Xi(\zeta)$, as defined in \eqref{eqn.Xi}, can be understood as the \emph{inverse Hessian}  of the dual problem \eqref{prob.OPP}, and the vector $Z(\zeta,w)$, as defined in \eqref{eqn.Z}, is the gradient of \eqref{prob.OPP}. 
This interpretation suggests that  \eqref{sys.PriceController} represents a continuous-time \emph{Newton algorithm} for the dual problem \eqref{prob.OPP}, augmented with a feed-forward component and a regulation term for the inventory levels.\footnote{The projection matrix $Q$ ensures only that the controller generates dual variables with a constant mean value. Note that the optimal dual variables are only defined up to a variation in $\mathrm{span} (\1)$. Thus, for any desired mean value, there exists an optimal dual solution having exactly this mean value.}
Furthermore, as the dual variables of the optimal flow problem \eqref{prob.OFP} are often given an economic interpretation as \emph{prices}, the controller \eqref{sys.PriceController} can be understood as an optimal pricing controller.
\end{remark}



\section{Design Example} \label{sec.Example}
We illustrate the proposed controller design on a small example network. We consider an inventory system with three inventories and three transportation lines as illustrated in Figure \ref{fig.Inventory}.

We assign an arbitrary orientation to the three connecting lines to obtain the incidence matrix
\begin{align*}
B = \begin{bmatrix} -1  &0 & 1 \\   1 & -1 & 0 \\  0  & 1  & -1 \end{bmatrix}.
\end{align*}

\begin{figure}
\begin{center}
\begin{tikzpicture}
[robot/.style={circle,draw=black,fill=white,minimum size=7pt, drop shadow={color=gray!80}, line width=1pt}]
\node[robot] (1) at (0,0) [label=above:$1$]{};
\node[robot] (2) at (1.5,1.8)[label=above:$2$] {};
\node[robot] (3) at (3,0 ) [label=above:$3$]{};
\draw[ line width=1pt](1)-- node[label=left:$\lambda_{1}(t)$]{}(2);
\draw[ line width=1pt](2)--node[label=right:$\lambda_{2}(t)$]{}(3);
\draw[  line width=1pt](3)--node[label=below:$\lambda_{3}(t)$]{}(1);
\draw[<->,  line width=1pt](-0.8,0)--node[label=below left:$w_{1}(t)$]{}(1);
\draw[<->,  line width=1pt](0.7,1.8)--node[label=above left:$w_{2}(t)$]{}(2);
\draw[<->,  line width=1pt](3.8,0)--node[label=below right:$w_{3}(t)$]{}(3);
\end{tikzpicture}
\caption{Illustration of the inventory system.}
\label{fig.Inventory}
\end{center}
\end{figure}

For the first two lines, we assume that the flows are unconstrained and cost functions are quadratic
\begin{align*}
P_{k}(\lambda_{k}) = \frac{q_{k}}{2}(\lambda_{k} - r_{k})^{2}, \; k \in \{1,2\}.
\end{align*}
For the simulations we choose $q_{1} = 0.2$, $q_{2} = 10$, $r_{1}=2$, and $r_{2} = 0$. The inverse of the gradients for these functions become simply
\begin{align*}
\nabla \mathcal{P}_{k}^{-1}(z_{k}) = q_{k}^{-1} z_{k} + r_{k}, \; k \in \{1,2\}.
\end{align*}

For the third transportation line, we assume that the the flow is constrained to $|\lambda_{e}(t)| \leq \bar{\lambda}_{3}$. We assign therefore to line three a strictly convex cost function that goes to infinity as $\lambda_{3}$ approaches $\pm \bar{\lambda}_{3}$. We choose here the cost function
\begin{align}
\mathcal{P}_{3}(\lambda_{3}) =  -c a \log \Bigl( \cos(\frac{\lambda_{3}}{a}) \Bigr) 
\end{align}
where $c > 0$ and $a = \frac{2}{\pi}\bar{\lambda}_{3}$.
For the simulation, we choose $c = 0.1$ and $\bar{\lambda}_{3} = 4$. 
Choosing this cost function, one obtains
\begin{align}
\nabla \mathcal{P}_{3}^{-1}(z) = a \tan^{-1}(\frac{z}{c}).
\end{align}
The function $\mathcal{P}_{3}(\lambda_{3})$ and the function $\nabla \mathcal{P}_{3}^{-1}(z_{3})$ are shown in Figure \ref{fig.dP3}.


\begin{figure}[b]
\begin{center}
   \begin{subfigure}{0.23\textwidth}
   \input{P3.tikz}
    \label{fig.P3}
   \end{subfigure}
   %
   \hspace{2em}
   %
   \begin{subfigure}{0.23\textwidth}
   \input{dP3inv.tikz}
   \label{fig.dP3inv}
   \end{subfigure}
\end{center}   
   \caption{The cost function $\mathcal{P}_{3}(\lambda_{3})$ and the corresponding inverse of the gradient $\nabla \mathcal{P}_{3}^{-1}(z_{3})$.}
   \label{fig.dP3}
  \end{figure}


To construct the controller, the matrix $X(\zeta)$ as defined in \eqref{eqn.TreeSol} is required. For convenience, we will use in the following the notation we write for convenience  $h_{k} = \nabla \mathcal{P}_{k}^{-1}(\zeta_{k}), k \in \{1,2,3\}$.  Following the procedure described in the appendix, one can determine $X$ as
\begin{align*}
X = \frac{1}{\bar{h}} \begin{bmatrix}
0 & h_{1}h_{3} + h_{1}h_{2} & h_{1}h_{2} \\
h_{2}h_{3} & 0 & h_{2}h_{3} + h_{1}h_{2} \\
h_{2}h_{3} + h_{1}h_{3} & h_{1} h_{3} & 0
\end{bmatrix}
\end{align*}
where $\bar{h} = h_{1}h_{2}+ h_{2}h_{3} + h_{1}h_{3}$ .

For the simulations, we assume that the flow is generated by linear dynamical systems and that it is consequently an harmonic signal. We assume furthermore, that the supply generated at one node is an instantaneous demand at another node, i.e.,
\begin{align*}
w_{i}(t) = \sum_{ k = 1}^{3} b_{ik} \bigl(\kappa_{k} \sin(\phi_{k}t + \rho_{k}) \bigr).
\end{align*}
For the simulation we have chosen $\kappa_{1}=2$, $\kappa_{2} = \kappa_{3}=4$, $\phi_{1}=2$, $\phi_{2} = 4$, $\phi_{3} = 8$, and $\rho_{1} = 0$, $\rho_{2} = 2$, $\rho_{3} = 3.14$.
We consider a scenario, where the supply/demand evolves first continuously, but chances instantaneously at a certain time instant. That is, at time $t=3$, we change the phase shift of the reference signal to 
$\rho_{1} = 4$, $\rho_{2} = 6$, $\rho_{3} = 2$.
The simulation results are shown in Figure \ref{fig.Sim}. The controller balances all inventory levels and routes the flows such that they all approach the optimal flows. Note that the flow $\lambda_{3}$ remains always within its capacity bounds. It seems remarkable, that the transportation line $2$, which has a fairly high cost, is mainly used when the flow in line 3 is close to the capacity bounds.
At time $t=3$ the reference signal is changed instantaneously as described above, and the controller needs to readjust. However, the feedback controller can tolerate this change and regulates the flows to the new optimal strategy. The simulation study illustrates the robustness of the proposed feedback controller.



\begin{figure}[t]
\begin{center}
\begin{subfigure}{0.23\textwidth}
\input{x.tikz}
\end{subfigure}
\hspace{2em}
\begin{subfigure}{0.23\textwidth}
\input{Flow1.tikz}
\end{subfigure}

\begin{subfigure}{0.23\textwidth}
\input{Flow2.tikz}
\end{subfigure}
\hspace{2em}
\begin{subfigure}{0.23\textwidth}
\input{Flow3.tikz}
\end{subfigure}
\end{center}
\caption{Inventory levels (top left) and flows on the lines $\lambda_{k}$ (solid) together with the optimal flows $\lambda_{k}^{w}$ (dashed).} \label{fig.Sim}
\end{figure}


\section{Conclusions}
We proposed a feedback controller for an optimal routing control in distribution networks with storage and time-varying supply. 
The controller is composed of a feed-forward component, generating the time-varying optimal dual variables, and a feedback component, regulating the system to the flows to the optimal trajectory. 
The construction of the feed-forward component of the controller utilizes the analytic solution to a generalized Dirichlet problem.
With certain choices of the cost functions, the proposed controller will also respect capacity constraints on the transportation lines. 
One can expect that the presented results can be generalized in various directions. A relevant extension will be a distributed implementation of the controller. 
Furthermore, future research will focus on a controller design that does note require a measurement of the external supply and demand. 

