{'fraction_non_alphanumeric': 0.051154385188463526, 'fraction_numerical': 0.012060190307590174, 'mean_word_length': 4.1637783279375356, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal pricing control in distribution networks with time-varying supply and demand

This paper studies the problem of optimal flow control in dynamic inventory systems. A dynamic optimal distribution problem, including time-varying supply and demand, capacity constraints on the transportation lines, and convex flow cost functions of Legendre-type, is formalized and solved. The time-varying optimal flow is characterized in terms of the time-varying dual variables of a corresponding network optimization problem. A dynamic feedback controller is proposed that regulates the flows asymptotically to the optimal flows and achieves in addition a balancing of all storage levels.


I. 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 network optimization [1], [2], [3]. 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 [4], 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 [5], [6]. The aspect of balancing in distribution networks has been studied in [7] using a one-step model predictive control scheme. In [8], [9], 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 [10], [11] 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 [12]. The connection between output agreement and network flow control has been further studied from an internal model control perspective in [13], [14], [15]. 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 [15], also controllers were proposed that route the flow asymptotically optimal through the network, provided the underlying cost functions are quadratic.
M. Bürger  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 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 generalized Dirichlet problem [16]. 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 II. The time-varying optimal flow problem is presented in Section III, 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 IV. The design procedure is illustrated on an example system in Section V.

Preliminaries:
The notation we employ is standard. The set of (non negative) real numbers is denoted by R (R ≥ ). The distance of a point q from a set A is defined as dist A q = inf p∈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 R(B) and N (B), respectively. A graph 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 G is assigned an (arbitrary) orientation. The incidence matrix B ∈ R n×m of the graph G, is a {0, ±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 P :

II. TIME-VARYING DISTRIBUTION PROBLEMS
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 G = (V, E), and let B be the incidence matrix (with arbitrary orientation). The dynamics of the inventory system is given asẋ where x ∈ R n represents the storage level, λ ∈ R m the flow along the lines, and P w an external in-/outflow of the inventories, i.e., the supply or demand. This basic model (or its discrete time version) is studied in [4], [8], [6], [7], [13] and is a relevant model for supply chains ( [17]) or data networks ( [18]). 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ẇ with initial condition w(0) = w 0 .
Assumption 1 There exists a compact set W that is forward invariant under (2) and the initial conditions w(0) are contained in this set W.
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.

Assumption 2
The supply/demand of the network is balanced, that is 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 ∈ {1, . . . , m} there existsλ k ∈ R >0 and the flow must satisfy If the flow on one line is unconstrained, we letλ k be infinity. In the following, we will also write −λ ≤ λ ≤λ, meaning that the inequalities hold for each component.
A cut of the network is defined to be a signed edge set 1 Q ⊂ E (with positive and negative parts denoted by Q + and Q − ) such that, for some node set S ⊂ V We will denote the cut induced by a set of nodes S ⊂ V with Q(S). The capacity of a cut is To ensure feasibility of the problem, we impose the following condition.

Assumption 3 There exists ǫ > 0 such that for any set of nodes
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. (1) with a time-varying supply/demand generated by the exo-system (2) and capacity constraints (3). Suppose Assumptions 1, 2 and 3 hold. Then, (i) for any flow vector λ ∈ R m the average inventory level x avg (t) = 1 n 1 T n x(t) remains constant, that is x avg (t) = x avg (0) for all t ≥ 0; and (ii) there exists a continuous flow trajectory λ : R ≥0 → R m that is strictly feasible, that means there is a vector valued function λ : R ≥0 → R m such that for each edge k the bounds (3) hold with strict inequality and Bλ(t) + P w(t) = 0.

Proposition 1 (Feasibility) Consider the inventory system
Proof: To see the first claim note that the average x avg (t) evolves asẋ avg = 1 n 1 T n Bλ + 1 n 1 T n P w(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.,λ k < ∞. Let Λ = diag{λ 1 , . . . ,λ m }. Consider the linear optimization problem The linear programming dual can be derived as Note now that the extreme points of the polytope defined by ΛB T u 1 ≤ 1 are multiples of the vectors whose elements are zeros and ones (see e.g., [19]). 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 ⊆ V according to (4), c Q ∈ R is a constant such that |c Q | =λ −1 Q , and e Q ∈ R n is a vector with entries being [e Q ] i = 1 if node i ∈ S and [e Q ] i = 0 if i ∈ V \ S. It follows now from Assumption 3 that (6) must always have an optimal value strictly smaller than one. From strong duality of linear programs follows that also (5) must have an optimal value strictly smaller than one. This shows that at each fixed time instant t, there exists a feasible flow vector λ(t) that satisfies all constraints, i.e., −λ < λ(t) <λ.
It remains to prove that there exists a continuous function λ : R ≥0 → R m that represents a feasible flow trajectory. Consider now the set of feasible flows for a given supply/demand vector w, i.e., The set ∆(w) is compact, convex and has a non-empty relative interior. One can understand ∆ as a set-valued map from R r to R m . As a set-valued map, ∆ is lower semi-continuous. 2 The existence of a continuous selection λ(w) of ∆ follows now from Micheal's Theorem ([20, Theorem 9.1.2]). Since additionally w(t) is a continuous function of time (i.e., it respects the dynamics (2)), there exists a continuous function λ(t) such that λ(t) ∈ ∆(w(t)), i.e., there is a continuous feasible flow trajectory.

III. OPTIMAL TIME-VARYING ROUTING
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 ∞ cost function for the flow λ k . We aim to characterize the steady state flows that are such that at each time instant the following static optimal distribution problem is solved where w = w(t) is the supply/demand at time t. We will use in the following the short hand notation P(λ) = m k=1 P k (λ k ).
The convex cost functions are such that they incorporate the capacity constraints of the transportation lines, in the sense that P k (λ k ) = ∞ for |λ k | >λ k . More precisely, we assume that the cost functions are of Legendre type [21, p. 258].
We can now characterize the optimality conditions of (7). The Lagrangian function associated to (7) with Lagrange multiplier ζ is is an optimal primal/dual solution pair to (7) if the following nonlinear equations hold Note that the first condition simply implies ∇P(λ w ) ∈ R(B T ) since ζ w has no further constraints. We define the set of optimal routing/supply pairs as Bλ + P w = 0}.
In particular, (λ w , w) ∈ Γ if and only if λ w is an optimal solution to the static optimal distribution problem (7) with the supply vector w.
The proof of this result provides a method to construct the dynamical system (9). However, before we present the proof, some additional discussion is required.
The main difficulty associated to the proposed characterization of the set Γ relates to the constraint ∇P(λ) ∈ R(B T ). To avoid this constraint, we express the optimality conditions in terms of the dual solutions ζ. Under the given assumption that the cost functions P k are of Legendre type, the gradients ∇P k are invertible on the open sets {λ k : −λ k < λ k <λ k } [3]. Since the optimal solution must always be strictly feasible, the two optimality conditions (8) can be combined to Bearing this in mind, we define the set of optimal dual solutions as We want to emphasize two properties of Γ D . First, if (ζ w , w) ∈ Γ D then (ζ w + c1, w) ∈ Γ D for any c ∈ R. Second, (ζ w , w) ∈ Γ D if any only if the routing strategy λ w = ∇P −1 (B T ζ w ) satisfies (λ w , w) ∈ Γ.

Remark 1
In the considered framework, the inverse of the gradient of P k is the gradient of the convex conjugate P ⋆ k , i.e., ∇P −1 (·) = ∇P ⋆ (·) [21]. For completeness, we want to point out that the dual problem to (7) is where P ⋆ k are the convex conjugates of P k . Any optimal solution ζ w to (11) is such that (ζ w , w(t)) ∈ Γ D . The main idea pursued here is to design the feed-forward controller for the dual variables. We need to design a dynamics for ζ(t) such that the constraint (10) is satisfied for all times. Suppose now (10) is feasible at the initial time, that is Then,ζ w (t) must be such that where DP −1 (·) = diag[∇ 2 P −1 1 (·), . . . , ∇ 2 P −1 m (·)]. Note that DP −1 (·) is the Hessian matrix of P ⋆ . The equation (13) defines an implicit differential equation for ζ w . For the design of the feed-forward controller (9), we need to transform it into an explicit differential equation.
Under the stated assumptions DP −1 (·) is a diagonal matrix, with all entries on the diagonal being positive. Thus, the matrix is a weighted, state dependent 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 ζ w , we make the following observation. Since 1 T n P w(t) = 0 there exists always ω(t) ∈ R m such that P w(t) = Bω(t).
One can always define ω(t) = B T (BB T ) † P w(t), where (BB T ) † is the generalized inverse of the matrix BB T . 3 Now, the conditions (13) can be written as Following [16,Chap. 16], we call the problem to findζ w a generalized Dirichlet problem.
It is described in [16, p. 659], how an analytic solution forζ w in (15) can be constructed. In particular, one can find analytically a matrix X(ζ w ) ∈ R m×n such thaṫ solves (13). The algorithmic procedure for the construction of X(ζ w ) is provided in [16, p. 659] and is reviewed in the appendix. The construction of X(ζ w ) has an explicit graph theoretic interpretation and the matrix has a characteristic structure. The entries of X(ζ w ) are fractions of polynomials, with each polynomial being the sum of products of different functions ∇ 2 P −1 k (·). Since, by assumption, all ∇ 2 P −1 k (·) are strictly positive and continuously differentiable, X(ζ w ) is continuously differentiable. Note that if the assumptions of Proposition 1 hold, the optimal flow problem (7) has at each time instant a feasible solution, and thus admits a finite dual solution. We are now ready to prove Proposition 2.

IV. FEEDBACK CONTROLLER DESIGN
We can use the proposed feed-forward controller now to design a dynamic feedback controller for the inventory system (1), 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 balancing of the inventory levels of (1). Network balancing is widely accepted as control objective in storage networks, as it promises to make best use of the available storage, [7], [13].
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 (18) 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 Furthermore, the projection matrix for the space orthogonal to the all ones vector is For designing the feedback controller, we need to relax the original assumption on the feasibility of the flows (i.e., Assumption 3).

Assumption 5 For any set of nodes
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. (1) with the supply generated by the dynamical system (2). Suppose the Assumptions 1, 2, 4, and 5 hold. Then, the dynamic controlleṙ

Theorem 1 (Feedback Controller) Consider the inventory dynamics
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 ζ(0), that is, all trajectories remain bounded, and (i) lim t→∞ B T x = 0 and (ii) lim t→∞ dist Γ λ(t), w(t) = 0.
Proof: We recall that, due to the assumption of balanced supply and demand, it must hold that Note that the average inventory level is constant over time since 1 T nẋ = 0 and therefore (x(t) − x w (t)) ∈ R(B).
Furthermore, we have that the average of the controller variables 1 n 1 T n ζ(t) is constant since 1 T nζ = 0. The optimal controller (i.e., dual) variables satisfy To prove the desired convergence result, consider a composed storage function of the form As the first component of the storage function, we consider This function satisfiesV where we exploited the two relations (22) and (23). Now, as the second component of U we consider the positive semidefinite function The function satisfies W (ζ, w) = 0 if and only if (ζ, w) ∈ Γ D , that is if ζ = ζ w . In order to keep the notation compact, we introduce the vector valued function The directional derivative of W satisfieṡ where we defined analogous to (14) L(ζ) := BDP −1 (B T ζ)B T and used the fact that L(ζ)Q = L(ζ).
Recalling the definition of Ξ(ζ) in (19) and the definition of X(ζ) provided in the appendix, and using that P s(w) ∈ R(B) it follows that

After these considerations it follows thaṫ
To complete the proof of convergence we consider now the composed (positive semidefinite) storage function The derivative of this function satisfieṡ 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 where s 2 > 0 is the smallest non-zero eigenvalue of the Laplacian matrix BB T . 4 The derivative of U (x − x w , ζ, w) satisfies thereforė where σ = min{s 2 , k}. It follows now from the Comparison Lemma [23,Lemma 3.4] that and consequently lim t→∞ U x(t) − x w (t), ζ(t), w(t) = 0. This proves that lim t→∞ x(t) − x w (t) = 0, which implies lim t→∞ B T x(t) = 0 and lim t→∞ Z(ζ(t), w(t)) = 0. These conditions imply lim t→∞ dist Γ λ(t), w(t) = 0.
It remains to show boundedness of all trajectories. The uniform boundedness of U x(t) − x w , ζ(t), w(t) 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∇P −1 (B T ζ) − P w(t) 2 . By Assumption 1, w(t) is bounded and consequently also B∇P −1 (B T ζ(t)) must remain bounded. Thus, it holds that at each time instant there exists a vector c(t) with finite norm such that It follows with the same argumentation as used previously to derive (10), that at each time instant the vector ζ(t) is an optimal dual solution to a network flow problem of the form min m k=1 P k (λ k (t)) s.t. Bλ(t) = c(t).
Due to Assumption 5 the problem has a feasible solution, and due to strict convexity of the cost functions (Assumption 4) the optimal solution is unique. From the optimality conditions follows that there exists at each time instant a unique z(t) such that ∇P(λ(t)) = z(t) and z(t) = B T ζ(t). Since additionally at each time instant 1 T ζ(t) = 1 T ζ(0), the vector ζ(t) is uniquely defined and must be bounded. This proves that all trajectories remain bounded.

Remark 2
The controller (21) has a characteristic structure. The matrix Ξ(ζ), as defined in (19), can be understood as the inverse Hessian of the dual problem (11), and the vector Z(ζ, w), as defined in (27), is the gradient of (11). This interpretation suggests that (21) represents a continuous-time Newton algorithm for the dual problem (11), augmented with a feed-forward component and a regulation term for the inventory levels. 5 Furthermore, as the dual variables of the optimal flow problem (7) are often given an economic interpretation as prices, the controller (21) can be understood as an optimal pricing controller. Fig. 1: Illustration of the inventory system.

V. DESIGN 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 1.
We assign an arbitrary orientation to the three connecting lines to obtain the incidence matrix For the first two lines, we assume that the flows are unconstrained and cost functions are quadratic 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 . For the third transportation line, we assume that the the flow is constrained to |λ e (t)| ≤λ 3 . We assign therefore to line three a strictly convex cost function that goes to infinity as λ 3 approaches ±λ 3 . We choose here the cost function where c > 0 and a = 2 πλ 3 . For the simulation, we choose c = 0.1 andλ 3 = 4. Choosing this cost function, one obtains The function P 3 (λ 3 ) and the function ∇P −1 3 (z 3 ) are shown in Figure 2.  To construct the controller, the matrix X(ζ) as defined in (16) is required. For convenience, we will use in the following the notation we write for convenience h k = ∇P −1 k (ζ k ), k ∈ {1, 2, 3}. Following the procedure described in the appendix, one can determine X as 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., For the simulation we have chosen κ 1 = 2, κ 2 = κ 3 = 4, φ 1 = 2, φ 2 = 4, φ 3 = 8, and ρ 1 = 0, ρ 2 = 2, ρ 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 ρ 1 = 4, ρ 2 = 6, ρ 3 = 2. The simulation results are shown in Figure 3. The controller balances all inventory levels and routes the flows such that they all approach the optimal flows. Note that the flow λ 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.

VI. CONCLUSIONS
We proposed a feedback controller for an optimal routing control in distribution networks with storage and timevarying 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.

A. Analytic Solution to the Generalized Dirichlet Problem
We review here the analytic solution for generalized Dirichlet problems provided in [16]. Consider the problem of finding a vector z that solves an equation of the form where B ∈ R n×m is the incidence matrix of a graph G = (V, E) (with arbitrary orientation), H is a diagonal matrix with positive elements on the diagonal, and v is a given vector. It is known ( [16]) that a solution to (33) is The analytic construction of the matrix X is reported and justified in [16, p. 659], and recalled below. In what follows we say that the node i ∈ V is the positive end of edge k ∈ E if [B] k,i = 1 and the negative end if [B] k,i = −1.
We call a spanning tree of the graph G an acyclic, connected subgraph T ⊂ G, defined on the same node set, having an edge set T being exactly of size |T | = n − 1.
Consider now any edge l ∈ T . After removing l from T , the spanning tree will partition into two connected components. The connected component that contains the positive end of l will be denoted with T + l . For the set of nodes forming the connected component T + l , we will write U (T , l). Each edge k ∈ E has a weight h k , which is the kth diagonal element of the matrix H. The weight of a spanning tree is now defined as h(T ) = k∈T h k .
Furthermore, we define the sum of the weights of all spanning trees of G bȳ We can now construct the matrix X ∈ R m×n in (34) according to the following procedure. For each edge k ∈ E and each node i ∈ V define the corresponding entry of X as that is the sum of the weights of all spanning trees for which the node i is contained in the positive connected component.
The proposed method can naturally be used to solve the equation (13). There, the weighting matrix is H = DP −1 (B T ζ w (t)). The time-varying nature of the weights has no influence on the solution. In fact, since the solution is analytically constructed, we will obtain expressions for X and H that depend explicitly on the state.