#### Abstract

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.

#### 1. Introduction

During the last few years, a huge amount of literature has emerged that deals with theoretical stability results for boundary control of conservation laws; see, for example, [1–15]. Further work also includes hyperbolic systems with source terms (so-called balance laws) for special applications such as gas dynamics [16–18] or water flow in open canals [19–23]. To the best of our knowledge only a few results are available, where these kinds of control problems are analyzed from a numerical point of view; see [24, 25].

Many network applications are represented by one-dimensional hyperbolic systems of balance laws [4, 26–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 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.

#### 2. 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.

##### 2.1. 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 :where is a matrix with constant entries, that is, , and for are the eigenvalues of the matrix with . The eigenvalues are either positive or negative, so we order these values such that for and for and split the matrix into and .

For characteristic variables , we have and such that 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 bywith matrix being a matrix, that is, , with blocks and . 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 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 for From theory we know that the linear hyperbolic system (1) and (3) is exponentially stable around the equilibrium state if there exist and such that, for every , the solution to the Cauchy problem satisfies

In the case of linear boundary conditions (3) the exponential stability of the equilibrium can be ensured by the following theoretical result [26, 33].

Let be the set of all real diagonal matrices with strictly positive diagonal elements. We define . Then the following theorem holds.

Theorem 1 (Theorem 2, [26]). *If , there exists such that, if , then the linear hyperbolic system (1) and (3) is exponentially stable.*

The key idea of the proof is to show that the time derivative of the Lyapunov functionwith weighting matrix and positive constants , for all , is negative definite along the solution of (1) and (3).

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.

##### 2.2. Numerical Feedback Stabilization

We use a splitting method [34] to discretize system (1) and (3). This means in a first step the advection part 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 , component , and cell index . We consider the cells and as ghost cells to describe the boundary values. So is the cell width of an equidistant spatial grid and the number of cells in the domain . 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 and time :withwhere 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 formwith evaluated at grid point and time step and positive real numbers and

Now, we are able to present the discrete version of Theorem 1 based on the numerical discretization (7)–(12) on the spatial domain for a positive finite time horizon . The main result taken from [25] can be summarized as follows.

Theorem 2. *The discretized system (7)–(12) is exponentially stable; that is, for and , if there exist and such that**(C1) the matricesare negative definite, with defined as (C2) the source term is restricted to Additionally, the decay rate can be explicitly given byand converges to for for all .*

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, , for the numerical discretization (7)–(12).

In the next section we apply the results of Theorem 2 to networks [17, 26–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.

#### 3. 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 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 , where is a set of nodes and is a set of edges. We represent an edge by an interval with uniform length .

To keep the notation simple, we assume that each edge (or interval ) can be mapped to the domain .

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 readswith

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 , and getThese 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.

##### 3.1. 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 thenfor , For a directed graph with edges we are able to rewrite (24) in a compact form provided that is independent of . We collect , and introduce the matricesand with The coupling condition is rewritten in terms of , where , , with . Then, the network equations can be written in a compact form such as (21) and (24).

*Remark 4. *Note that the time step size should be additionally restricted in such a way that the numerical scheme (1) and (3) is still valid and assumption (C2) in Theorem 2 is satisfied; that is, is negative definite for all edges . This is the case for (cf. [29] for the details of the proof).

###### 3.1.1. Numerical Results

For all computations we choose the length of all edges as one over two, that is, , and also discretize the edges with an equidistant grid. We use the final time horizon . 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 and for all . If , the telegrapher’s equation reduces to the linear wave equation. Also note that and . We control the inflow of at node and the inflow of at node depending on the outflow at and , respectively; that is, and . The matrices and are composed ofThe corresponding eigenvalues of and (cf. (C1) in Theorem 2) are for , and for . If we assume that for all and insert the values for the spatial discretization , , on the strips and , we getWe 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, . However, these eigenvalues will not influence the stability result since they are independent of the feedback parameters , . So we have to ensure that the matrices and are negative definite for the feedback relevant eigenvalues , We claim which is the case ifSince it should also hold that and .

For the numerical computations we set as initial values , . We will vary the values for and and expect a rapidly dropping discrete Lyapunov function (14) for small and 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, , the decay rate converges to for 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 . The network has controls (inputs) at the outward nodes for .

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: with Here, and . We set for all edges . The parameters are , for all edges. Then we have the following matrices: with So the eigenvalues of and of are determined by the eigenvalues of the matrices, The eigenvalues of the first two matrices are in 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 conditionis satisfied, the assumption (C1) of Theorem 2 leads to negative definite matrices. Since has to be greater than zero, must be greater than one. So . Therefore, we might choose . 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 and for all .

##### 3.2. Traffic Flow on Road Networks

As another application we use a model originated from traffic dynamics, the so-called LWR traffic network model [35–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 lawwhere 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 is a concave function, with and , which is maximal at a value . These values separate a free-flow state, that is, from a traffic-congestion state, that is, . If we only consider a free-flow 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 , in the sense of flux conservation. So considering a network and the additional vector notation we get the whole LWR-networked system

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.

###### 3.2.1. Numerical Results

We consider the LWR-model on the network in Figure 3. Here, is just the inflow and 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: For given inputs and , there also exists a free-flow steady-state . We can then define a feedback control law, such as with . The goal is to design a feedback law to ensure that small perturbations of the initial value condition are damped. For this purpose we analyze (49). To fulfill the assumption (C2) of Theorem 2, we have to ensure that Since is a negative constant for concave flux functions, it follows that has to be less than or equal to zero. For our numerical simulations, we use the flux function ; that is, and for all edges. Hence, if .

By linearizing the coupling conditions above, we get We set and . For the initial values we choose We can now compute the matrix We also set and for all edges . Then, the matrices and are So with the matrix is negative definite. For the simulation we set , , and . Also the perturbation is set to zero for and . The results of the computation are shown in Figures 4 and 5. One can observe that the perturbation is damped to zero, which illustrates the proven properties.

**(a)**

**(b)**

#### 4. Conclusion

We have introduced a feedback stabilization result for general network topologies and applied it to the telegrapher’s equation and the LWR traffic flow model. For the telegrapher’s equation we have used the boundary value coupling to set up the feedback stabilization in the case of interlinked networks. As a further example we have considered a scalar conservation law, namely, the LWR-model. We have seen that the proposed stabilization framework can be also used to dampen perturbations along edges to reach a desired steady-state.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the DFG project “Novel models and control for networked problems: from discrete event to continuous dynamics,” Grant no. GO 1920/4-1.