Numerical Feedback Stabilization with Applications to Networks

The focus is on the numerical consideration of feedback boundary control problems for linear systems of conservation laws including source terms. We explain under which conditions the numerical discretization can be used to design feedback boundary values for network applications such as electric transmission lines or traffic flow systems. Several numerical examples illustrate the properties of the results for different types of networks

Many network applications are represented by onedimensional hyperbolic systems of balance laws [4,[26][27][28][29][30].For the analytical case it can be proven that feedback boundary values, designed under certain conditions, yield an exponential decay of a continuous Lyapunov function [3,9,31] also in the context of networks.
In this paper, we are mainly concerned with the numerical analysis of the stability of the steady-state of such networked systems.We explain how suitable boundary values must be designed such that a linear system of balance laws converges exponentially to a steady-state.Inspired by [24,25], we restrict ourselves to the numerical approximation of such systems and study the asymptotic behavior of that solution or to be more precise the conditions under which an exponential decay of the discrete solution to hyperbolic balance law can be attained.One way to influence a system of balance laws is through the boundary values.The interpretation of the boundary values in the case of hyperbolic systems is that they describe an outflow or an inflow of a quantity.One can measure the outflow, modify it, and let it flow back again as inflow.This kind of control is then called feedback boundary control.
Our considerations are done first for arbitrary linear systems of balance laws and are then extended to networks which are in fact coupled systems of boundary problems.For simulation purposes, we present two concrete examples.One example is the so-called telegrapher's equation, an equation described by a 2 × 2 hyperbolic balance law and normally used to model transmission lines [29].The other example is a model for road traffic, the so-called LWR-model [32].It is a scalar model describing traffic flow based on physical conservation laws.The question we address here is how a feedback control can be used that small errors in some measured data are dampened.
The outline of the article is as follows: in Section 2, we answer the question under which conditions an exponential decay of the numerical solution to the hyperbolic equations can be achieved.The network applications are given in Section 3. We design feedback boundary values for different networks and analyze the stability conditions.

Problem Description and Preliminary Work
Hyperbolic systems of balance laws are often used to describe physical applications such as gas, water, traffic, or power flow.Boundary conditions can be designed to ensure that the analytical solution for these systems converges towards a desired steady-state.To prove the convergence of the approximate solution to a desired steady-state, a technique based on Lyapunov functions [3,4,8,28] is used.The required boundary conditions are called feedback boundary conditions and the method is the boundary feedback stabilization of hyperbolic systems.This is a common approach for controlling systems of conservation laws.
In two recent works [24,25], the theoretical results mentioned above are transferred to the numerical discretization of systems of balance laws.It can be shown that the approximate solution also converges exponentially to a steady-state.This is done using a suitable discretized Lyapunov function.Furthermore, this approach allows for an explicit representation of the decay rates in contrast to many theoretical results.

Theoretical Results on Feedback Stabilization.
Our goal is the stability analysis of the numerical solution to systems of conservation laws with linear source terms.We assume that the hyperbolic systems under investigation are given in the so-called characteristic form.Typically, linear systems of balance laws can be expressed in a characteristic form in a straightforward way while nonlinear systems can be transformed into a characteristic form by linearization with respect to a steady-state solution; see [3,26] for examples and more details.
Let us consider the system of balance laws in characteristic form for  ∈ [0, ℓ]: where  is a  ×  matrix with constant entries, that is,  ∈ R × , and   ̸ = 0 for  ∈ {1, . . ., } are the eigenvalues of the matrix Λ with Λ = diag( 1 , . . .,   ).The eigenvalues are either positive or negative, so we order these values such that  +  fl   > 0 for  ∈ {1, . . ., } and  −  fl   < 0 for  ∈ { + 1, . . ., } and split the matrix So  ± are the components of  corresponding to the positive and negative eigenvalues of Λ.For the above system (1), linear feedback boundary conditions [7,9] can be prescribed by with  matrix being a  ×  matrix, that is,  ∈ R × , with blocks  + ∈ R × and  − ∈ R −×− .The name feedback originates from the fact that outputs of a system are fed back as inputs.Since the inputs and outputs of a system of balance laws are the boundaries, one speaks of boundary feedback control or boundary feedback stabilization, if the control is used to reach a steady-state (or equilibrium).Note that boundary control problems might depend on additional controls  + () and  − (); see applications in Section 3.Then, the aim is to find suitable controls such that for any initial data  0 () and boundary values (3) the solution of system (1) converges to the desired steady-state.
We are now concerned to analyze the exponential stability of system (1) under boundary conditions (3) with initial data (0, ) =  0 () for  ∈ (0, ℓ).From theory we know that the linear hyperbolic system (1) and ( 3) is exponentially stable around the equilibrium state  ≡ 0 if there exist ] > 0 and  > 0 such that, for every  0 ∈  2 ((0, ℓ); R  ), the solution to the Cauchy problem satisfies In the case of linear boundary conditions (3) the exponential stability of the equilibrium  ≡ 0 can be ensured by the following theoretical result [26,33].
Let D  be the set of all real  ×  diagonal matrices with strictly positive diagonal elements.We define () fl inf{‖ΔΔ −1 ‖ 2 : Δ ∈ D  }.Then the following theorem holds.
In [25] it is shown how this theoretical result can be reinterpreted from a numerical point of view.This means there exists a discrete version of Theorem 1 for the discretized system (1) and (3).The question is then under which additional conditions exponential stability of the discretized system can be ensured.In particular, corresponding conditions for the feedback matrix  and the source term matrix  need to be established.This question will be addressed in the next paragraph.

Numerical Feedback Stabilization.
We use a splitting method [34] to discretize system (1) and (3).This means in a first step the advection part    + Λ   = 0 is solved using an upwind scheme either in positive or in negative direction.In a second step, an explicit Euler scheme is applied to resolve the source term equation    = −.
Therefore, we discretize  for any component , with time step size Δ and space step size Δ, that is,   , for every discrete time  ∈ N >0 , component  ∈ {1, . . ., }, and cell index  ∈ {−1, . . ., }.We consider the cells −1 and  as ghost cells to describe the boundary values.So Δ is the cell width of an equidistant spatial grid and  + 2 the number of cells in the domain [0, ℓ].The time step size Δ is then chosen according to the CFL condition; that is, Note that for stiff problems there might be an additional time step restriction due to the explicit Euler discretization of the source term equation.Hence, we end up with the following numerical scheme to solve system (1) and ( 3) for every cell index  ∈ {0, . . .,  − 1} and time  ∈ {1, . . .,  − 1}: with where  , is the element in the th row and th column of the matrix  and   is the th row of the source term matrix.As one can recognize, ( 7) and ( 9) are the numerical approximation to the advection and source term problem and ( 10) and (11) describe the boundary conditions.The initialization is given by ( 12).Since we are interested in boundary feedback control from a numerical point of view, we need to define a discrete Lyapunov function of (5).This is necessary to analyze the exponential stability of the discretized system under linear boundary conditions.Therefore we consider a discrete candidate Lyapunov function of the form with   , evaluated at grid point  and time step  and positive real numbers   > 0 and   > 0.
Now, we are able to present the discrete version of Theorem 1 based on the numerical discretization ( 7)-( 12) on the spatial domain [0, ℓ] for a positive finite time horizon  = Δ.The main result taken from [25] can be summarized as follows.
For the proof, the same idea is used as in Theorem 1.The major difference is to show that the discrete time derivative of the discrete Lyapunov function ( 14) is negative; that is, ( +1 −   )/Δ ≤ 0, for the numerical discretization ( 7)- (12).
In the next section we apply the results of Theorem 2 to networks [17,[26][27][28][29] represented by directed graphs.Networked problems, where the dynamics on edges is governed by systems of balance laws, are linked at intersections and can be therefore treated as coupled boundary value problems.Hence, the results from above can be extended in a rigorous way.

Application to Networks
A network is modeled by a directed graph and can be viewed as a spatial one-dimensional domain.To define hyperbolic balance laws on a network one has to introduce coupling conditions, which are a special form of boundary conditions at intersections.The associated feedback control problem is then to design control laws at the coupling interfaces to stabilize the system.
For application purposes, we use the results presented in Section 2 on two concrete network examples.One example is the so-called telegrapher's equation, an equation described by a 2 × 2 hyperbolic balance law and normally used to model transmission lines.The other example is a model for road traffic, the so-called LWR traffic flow network model.In this context, we will explain how a feedback control can be designed that small measurement errors are dampened.
In our numerical study, we discuss several scenarios and demonstrate the properties of Theorem 2 for small and larger networks.This means, in particular, to check the assumptions (C1) and (C2) for different network topologies.
Let us now start how the results from Section 2 fit into the context of networks.Definition 3. A network is described by a directed graph G = (V, E), where V is a set of nodes and E is a set of edges.We represent an edge  by an interval I ⊂ R + with uniform length ℓ.
To keep the notation simple, we assume that each edge (or interval I) can be mapped to the domain  ∈ [0, ℓ].
On each edge, we consider a system of balance laws (1).These systems are coupled by so-called coupling conditions prescribed as boundary values.For the coupling conditions we use the conservation of flow with inputs at the nodes.These inputs correspond to our control variables (cf.comments after (3)).So for a system of balance laws such as (1) the coupling has the following componentwise form: and here  ±  are external control inputs at fixed edges .So with this, the whole system in vector notation reads with One can describe the controls as so-called feedback controls.A feedback law designs the controls in such a way that they depend on the outflow at the boundaries.We prescribe linear feedback boundary conditions for (21) as with matrices G± ∈ R × , and get These are the full network equations including boundary conditions as coupling and as feedback.
In the following, we present two applications, in particular electrical transmission lines and vehicular traffic networks.For the first example, we describe how the network equations look like and which conditions must be ensured to obtain exponential stability while for the second example we explain how errors in the data can be dampened with the presented approach.

Telegrapher's Equation on Networks.
We are interested in a feedback control law for electrical transmission lines.Those settings are normally represented as networks, where on each edge the so-called telegrapher's equation is solved.The telegrapher's equations are a simplification of the Maxwell equation and include signal losses on transmission lines; see [29] for the modeling details.
Normally, the variables current (, ) and voltage (, ) are the conservation quantities that are described by the telegrapher's equation.They are expressed in terms of characteristic variables as (, ) =  + (, ) +  − (, ) and (, ) = √/( + (, ) −  − (, )).Due to the losses alongside the lines, we have an additional reaction term  depending on the positive constants of , , and  (resistance, inductance, and admittance) and  (capacitance).The equation written in characteristic variables on edge  is then for  ∈ [0, ℓ], For a directed graph with |E| edges we are able to rewrite (24) in a compact form provided that  ± is independent of V ∈ V.

Numerical Results
. For all computations we choose the length of all edges as one over two, that is, ℓ = 1/2, and also discretize the edges with an equidistant grid.We use the final time horizon  = 12.If some parameters change throughout the section, it will be stated.For simplicity, we skip the index  whenever the situation is clear and all parameters are fixed to the same value for all edges.
Example 5.As a first example we consider the network shown in Figure 1.
We intend to introduce the network framework introduced above step by step.On each edge , we have two quantities  +  and  −  .We first choose   =   = 0 and   =   = 1 for all .If   =   = 0, the telegrapher's equation reduces to the linear wave equation.Also note that  +  = 1 and  −  = −1.We control the inflow of  + 1 at node V 0 and the inflow of  − 2 at node V 2 depending on the outflow at V 2 and V 0 , respectively; that is,  + 1 (, 0) =  1  + 2 (, ℓ) and  − 2 (, ℓ) =  2  − 1 (, 0).The matrices  + and  − are composed of The corresponding eigenvalues of for   +  + 0  + −  +  , and for If we assume that  =   for all  and insert the values for the spatial discretization  1,0 = 0,  1, =  2,0 = 1/2,  2, = 1 on the strips [0, 1/2] and [1/2, 1], we get We remark that the determination of the eigenvalues is dependent on the network discretization, that is, the mapping of the edges to spatial intervals.For the discretization used here, we obtain two eigenvalues that are zero; that is,  1 ,  3 = 0.However, these eigenvalues will not influence the stability result since they are independent of the feedback parameters  1 ,  2 .So we have to ensure that the matrices   +  + 0  + −  +  and   −  − −1  − −  − −1 are negative definite for the feedback relevant eigenvalues  2 ,  4 .We claim  2 ,  4 < 0 which is the case if 2 )) .Since  > 0 it should also hold that 0 <  2  1 < 1 and 0 <  2 2 < 1.For the numerical computations we set as initial values  +  (0, ) = −1/2,  −  (0, ) = 1/2.We will vary the values for  1 and  2 and expect a rapidly dropping discrete Lyapunov function   (14) for small  1 and  2 and vice versa.This behavior can be observed in Table 1.Since the CFL constant is equal to one and also the grid size is very small, that is, Δ = 1/1600, the decay rate ] converges to  for |  | = 1 due to (19).Example 6.As a more complex network, we consider the cascade network in Figure 2 for a different number of layers .Each layer consists of seven edges and six nodes.The nodes are arranged such that three in a row are facing each other.A node is then connected to its opponent in the next row and the direct neighbor of the opponent.The number of edges is |E| = 7.The network has controls (inputs)  ±  at the outward nodes for  ∈ {1, 2, 3}.
By numbering the edges from left to right and from top to bottom, we have The matrices describing the network with  layers have the following form: ( Here,  ± ∈ R |E|×|E| and  1 ,  2 ,  ∈ R 7×7 .We set   =  for all edges .The parameters are   =   = 1,   =   = 0 for all edges.Then we have the following matrices: The eigenvalues of the first two matrices are in {0, −exp(±)} and the eigenvalues of the last two matrices are in With the same arguments used before, we restrict our investigations again to the feedback relevant eigenvalues.That means if the condition is satisfied, the assumption (C1) of Theorem 2 leads to negative definite matrices.Since  has to be greater than zero, (3/28)(max( 1 ,  2 )) −2 must be greater than one.So (max( 1 ,  2 )) 2 < 3/28.Therefore, we might choose  1 =  2 = 1/4.We expect that in our numerical experiments  goes to zero for increasing layers  → ∞; compare (41) and results in Table 2.This means, for a large number of layers, we can only expect a linear decreasing Lyapunov function   (cf.( 14)).
We also observe in Table 2 that the decay rate ( 19) is equal to  for our choice of parameters.For initial values we use  +  (0, ) = −2 and  −  (0, ) = 2 for all .

Traffic Flow on Road Networks.
As another application we use a model originated from traffic dynamics, the so-called LWR traffic network model [35][36][37][38][39]. From a mathematical point of view, we assume that on each edge of the network the car density is governed by the nonlinear conservation law where   describes the density of vehicles on each edge .
Here, the velocity of a car depends on the local density, that is, (  ), and the flux function   =   (  ) =   (  ) =   ( max (1 −   / max )) is a concave function, with   (0) = 0 and   ( max ) = 0, which is maximal at a value   .These values separate a free-flow state, that is,   ≤   from a trafficcongestion state, that is,   >   .If we only consider a freeflow state, the flux function is monotonically increasing.So we can rewrite (42) in the form of so-called kinetic wave equations, which is useful to define the coupling conditions.The coupling conditions in a road network have then the following componentwise form with distribution rates  + ẽ and  − ẽ and some additional input and output, either inflows  +  () from outside ramps or outflows  −  () from exits to the outside.But since we only consider a free-flow state, we do not need to couple back-traveling solutions and we can simplify the coupling conditions to with distribution rates  ẽ and inflows   ().Note, it has to hold that ∑ ẽ  ẽ = 1, in the sense of flux conservation.So considering a network G = (V, E) and the additional vector notation we get the whole LWR-networked system    +  ()    = 0,  (, 0) =  (, ℓ) +  () ,  (0, ) =  0 () .
(47) For a given state   of (42), which is perturbed with a perturbation   at some point in time and by linearization of the flux function   , we can rewrite the LWR-model (42) as a conservation law of the perturbation   as follows: A short calculation using (42) leads to Equation (49) will be now analyzed and associated with the stabilization framework introduced before.

Numerical Results.
We consider the LWR-model on the network in Figure 3. Here,  1 is just the inflow and  2 is a control input.Our goal is to analyze how we can control the inflow, such that small perturbations will vanish after some time.
As shown in the previous example, we can rewrite the model (42) in the form of kinetic wave equations, and couple the resulting conservation law as follows: (51) For given inputs  1 and  2 , there also exists a free-flow steadystate ( 1 ,  2 , . . .,  5 ).We can then define a feedback control law, such as

Table 1 :
The number of cells on each edge is  = 1600, with Δ = 1/1600.The CFL constant is equal to one; that is, Δ/Δ = 1.

Table 2 :
The number of cells on each edge is  = 200, with Δ = 1/200.The CFL constant is equal to one; that is, Δ/Δ = 1.