# Numerical Feedback Stabilization with Applications to Networks

Authors: Simone Göttlich; Peter Schillen
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

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


---

## Body

## 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 a2×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 Section2, 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 forx∈[0,l]:(1)∂tξ+Λ∂xξ+Bξ=0,ξ0,x=ξ0x,where B is a p×p matrix with constant entries, that is, B∈Rp×p, and λi≠0 for i∈1,…,p are the eigenvalues of the matrix Λ with Λ=diagλ1,…,λp. The eigenvalues are either positive or negative, so we order these values such that λi+≔λi>0 for i∈1,…,m and λi-≔λi<0 for i∈m+1,…,p and split the matrix Λ into Λ+=diagλ1+,…,λm+ and Λ-=diagλm+1-,…,λp-.For characteristic variablesξ∈Rp, we have ξ+∈Rm and ξ-∈Rp-m such that(2)ξ=ξ+ξ-. 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(3)ξ+t,0ξ-t,l=G+00G-︸≕Gξ+t,lξ-t,0,with G matrix being a p×p matrix, that is, G∈Rp×p, with blocks G+∈Rm×m and G-∈Rp-m×p-m. 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 controlsu+(t) and u-(t); see applications in Section 3. Then, the aim is to find suitable controls such that for any initial data ξ0(x) 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,x)=ξ0(x) for x∈(0,l). 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 C>0 such that, for every ξ0∈L2((0,l);Rp), the solution to the Cauchy problem satisfies (4)ξt,·L2≤Cexp-νtξ0L2,∀t∈R+.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].LetDp be the set of all real p×p diagonal matrices with strictly positive diagonal elements. We define ρ(G)≔inf{ΔGΔ-12:Δ∈Dp}. Then the following theorem holds.Theorem 1 (Theorem 2, [26]).
Ifρ(G)<1, there exists ε>0 such that, if B<ε, 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 function(5)Lt=∫0lξTPxξdx,with weighting matrix P(x)≔diagc1e-sign(λ1)μ1x,…,cpe-sign(λp)μpx∈Rp×p and positive constants μi>0, ci>0 for all i∈1,…,p, 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 G and the source term matrix B 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 ∂tξ+Λ∂xξ=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 ∂tξ=-Bξ.Therefore, we discretizeξ for any component i, with time step size Δt and space step size Δx, that is, ξi,jk for every discrete time k∈N>0, component i∈1,…,p, and cell index j∈-1,…,N. We consider the cells -1 and N as ghost cells to describe the boundary values. So Δx is the cell width of an equidistant spatial grid and N+2 the number of cells in the domain 0,l. The time step size Δt is then chosen according to the CFL condition; that is,(6)λΔtΔx≤1,λ≔maxi∈1,…,pλi.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 j∈0,…,N-1 and time k∈1,…,K-1:(7)ξ~i,jk=ξi,jk-ΔtΔxλiξi,jk-ξi,j-1k,fori∈1,…,m,(8)ξ~i,jk=ξi,jk-ΔtΔxλiξi,j+1k-ξi,jk,fori∈m+1,…,p,(9)ξi,jk+1=ξ~i,jk-ΔtBiξ~jk,fori∈1,…,p,(10)ξi,-1k=∑l=1mgi,lξl,N-1k+∑l=m+1pgi,lξl,0k,fori∈1,…,m,(11)ξi,Nk=∑l=1mgi,lξl,N-1k+∑l=m+1pgi,lξl,0k,fori∈m+1,…,p,(12)ξi,j0=1Δx∫xj-1/2xj+1/2ξ0,ixdx,fori∈1,…,pwith(13)ξ~j+k=ξ~1,jk,…,ξ~m,jkT,ξ~j-k=ξ~m+1,jk,…,ξ~p,jkT,ξ~jk=ξ~j,+k,ξ~j,-kT,where gi,l is the element in the ith row and lth column of the matrix G and Bi is the ith 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(14)Vk=∑j=0N-1∑i=1pξi,jk2ciexp-signλiμixi,jΔx,with ξi,jk evaluated at grid point j and time step k and positive real numbers ci>0 and μi>0.Now, we are able to present the discrete version of Theorem1 based on the numerical discretization (7)–(12) on the spatial domain [0,l] for a positive finite time horizon T=KΔt. The main result taken from [25] can be summarized as follows.Theorem 2.
The discretized system (7)–(12) is exponentially stable; that is,(15)∑j=0N-1∑i=1pξi,jk2Δx≤C~exp-νΔtk+1∑j=0N-1∑i=1pξi,j02Δx for k∈1,…,K and ν,C~>0, if there exist μi>0 and 0<ci≤1 such that
(C1) the matrices(16)G+TP0+Λ+G+-PN+Λ+,P-1-Λ--G-TPN-1-Λ-G-are negative definite, with Pj± defined as(17)Pj+=diagc1exp-μ1xi,j,…,cmexp-μmxi,j,Pj-=diagcm+1expμm+1xi,j,…,cpexpμpxi,j; (C2) the source term is restricted to (18)ξ~jkTΔtBTB-2Bξ~jk≤0,∀ξ~jk=ξ~1,jk,…,ξ~p,jkT.Additionally, the decay rate ν>0 can be explicitly given by(19)ν=mini1Δxλi1-ciexp-μiΔxand converges to miniλiμi for Δx→0 for all i.For the proof, the same idea is used as in Theorem1. The major difference is to show that the discrete time derivative of the discrete Lyapunov function (14) is negative; that is, Vk+1-Vk/Δt≤0, for the numerical discretization (7)–(12).In the next section we apply the results of Theorem2 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.


## 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 forx∈[0,l]:(1)∂tξ+Λ∂xξ+Bξ=0,ξ0,x=ξ0x,where B is a p×p matrix with constant entries, that is, B∈Rp×p, and λi≠0 for i∈1,…,p are the eigenvalues of the matrix Λ with Λ=diagλ1,…,λp. The eigenvalues are either positive or negative, so we order these values such that λi+≔λi>0 for i∈1,…,m and λi-≔λi<0 for i∈m+1,…,p and split the matrix Λ into Λ+=diagλ1+,…,λm+ and Λ-=diagλm+1-,…,λp-.For characteristic variablesξ∈Rp, we have ξ+∈Rm and ξ-∈Rp-m such that(2)ξ=ξ+ξ-. 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(3)ξ+t,0ξ-t,l=G+00G-︸≕Gξ+t,lξ-t,0,with G matrix being a p×p matrix, that is, G∈Rp×p, with blocks G+∈Rm×m and G-∈Rp-m×p-m. 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 controlsu+(t) and u-(t); see applications in Section 3. Then, the aim is to find suitable controls such that for any initial data ξ0(x) 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,x)=ξ0(x) for x∈(0,l). 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 C>0 such that, for every ξ0∈L2((0,l);Rp), the solution to the Cauchy problem satisfies (4)ξt,·L2≤Cexp-νtξ0L2,∀t∈R+.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].LetDp be the set of all real p×p diagonal matrices with strictly positive diagonal elements. We define ρ(G)≔inf{ΔGΔ-12:Δ∈Dp}. Then the following theorem holds.Theorem 1 (Theorem 2, [26]).
Ifρ(G)<1, there exists ε>0 such that, if B<ε, 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 function(5)Lt=∫0lξTPxξdx,with weighting matrix P(x)≔diagc1e-sign(λ1)μ1x,…,cpe-sign(λp)μpx∈Rp×p and positive constants μi>0, ci>0 for all i∈1,…,p, 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 G and the source term matrix B 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 ∂tξ+Λ∂xξ=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 ∂tξ=-Bξ.Therefore, we discretizeξ for any component i, with time step size Δt and space step size Δx, that is, ξi,jk for every discrete time k∈N>0, component i∈1,…,p, and cell index j∈-1,…,N. We consider the cells -1 and N as ghost cells to describe the boundary values. So Δx is the cell width of an equidistant spatial grid and N+2 the number of cells in the domain 0,l. The time step size Δt is then chosen according to the CFL condition; that is,(6)λΔtΔx≤1,λ≔maxi∈1,…,pλi.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 j∈0,…,N-1 and time k∈1,…,K-1:(7)ξ~i,jk=ξi,jk-ΔtΔxλiξi,jk-ξi,j-1k,fori∈1,…,m,(8)ξ~i,jk=ξi,jk-ΔtΔxλiξi,j+1k-ξi,jk,fori∈m+1,…,p,(9)ξi,jk+1=ξ~i,jk-ΔtBiξ~jk,fori∈1,…,p,(10)ξi,-1k=∑l=1mgi,lξl,N-1k+∑l=m+1pgi,lξl,0k,fori∈1,…,m,(11)ξi,Nk=∑l=1mgi,lξl,N-1k+∑l=m+1pgi,lξl,0k,fori∈m+1,…,p,(12)ξi,j0=1Δx∫xj-1/2xj+1/2ξ0,ixdx,fori∈1,…,pwith(13)ξ~j+k=ξ~1,jk,…,ξ~m,jkT,ξ~j-k=ξ~m+1,jk,…,ξ~p,jkT,ξ~jk=ξ~j,+k,ξ~j,-kT,where gi,l is the element in the ith row and lth column of the matrix G and Bi is the ith 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(14)Vk=∑j=0N-1∑i=1pξi,jk2ciexp-signλiμixi,jΔx,with ξi,jk evaluated at grid point j and time step k and positive real numbers ci>0 and μi>0.Now, we are able to present the discrete version of Theorem1 based on the numerical discretization (7)–(12) on the spatial domain [0,l] for a positive finite time horizon T=KΔt. The main result taken from [25] can be summarized as follows.Theorem 2.
The discretized system (7)–(12) is exponentially stable; that is,(15)∑j=0N-1∑i=1pξi,jk2Δx≤C~exp-νΔtk+1∑j=0N-1∑i=1pξi,j02Δx for k∈1,…,K and ν,C~>0, if there exist μi>0 and 0<ci≤1 such that
(C1) the matrices(16)G+TP0+Λ+G+-PN+Λ+,P-1-Λ--G-TPN-1-Λ-G-are negative definite, with Pj± defined as(17)Pj+=diagc1exp-μ1xi,j,…,cmexp-μmxi,j,Pj-=diagcm+1expμm+1xi,j,…,cpexpμpxi,j; (C2) the source term is restricted to (18)ξ~jkTΔtBTB-2Bξ~jk≤0,∀ξ~jk=ξ~1,jk,…,ξ~p,jkT.Additionally, the decay rate ν>0 can be explicitly given by(19)ν=mini1Δxλi1-ciexp-μiΔxand converges to miniλiμi for Δx→0 for all i.For the proof, the same idea is used as in Theorem1. The major difference is to show that the discrete time derivative of the discrete Lyapunov function (14) is negative; that is, Vk+1-Vk/Δt≤0, for the numerical discretization (7)–(12).In the next section we apply the results of Theorem2 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 Section2 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 Theorem2 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 Section2 fit into the context of networks.Definition 3.
A network is described by a directed graphG=V,E, where V is a set of nodes and E is a set of edges. We represent an edge e by an interval I⊂R+ with uniform length l.To keep the notation simple, we assume that each edge (or intervalI) can be mapped to the domain x∈[0,l].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: (20)λe+ξt,0=∑dee~+λe~+ξt,l+us+t,λe-ξt,l=∑dee~-λe~-ξt,0-us-t, and here us± are external control inputs at fixed edges s. So with this, the whole system in vector notation reads(21)∂tξ+Λ∂xξ+Bξ=0ξ+t,0ξ-t,l=Λ+-1D+00Λ--1D-ξ+t,lξ-t,0+u+tu-t,with (22)D+=matrix with entriesdee~+λe~+,D-=matrix with entriesdee~-λe~-,u+=u1+,u2+,…,uE+T,u-=u1-,u2-,…,uE-T,ξ+=ξ1+,…,ξE+T,ξ-=ξ1-,…,ξE-T,ξ=ξ+,ξ-.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 (23)u+t=G~+ξt,l,u-t=G~-ξt,0, with matrices G~±∈Rn×n, and get(24)∂tξ+Λ∂xξ+Bξ=0,ξ+t,0ξ-t,l=Λ+-1D++Λ+G~+ξ+t,lΛ--1D-+Λ-G~-ξ-t,0=G+00G-ξ+t,lξ-t,0.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.


### 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 currentI(t,x) and voltage U(t,x) are the conservation quantities that are described by the telegrapher’s equation. They are expressed in terms of characteristic variables as I(t,x)=ξ+(t,x)+ξ-(t,x) and U(t,x)=L/C(ξ+(t,x)-ξ-(t,x)). Due to the losses alongside the lines, we have an additional reaction term B depending on the positive constants of R, L, and Z (resistance, inductance, and admittance) and C (capacitance). The equation written in characteristic variables on edge e is then(25)∂tξe+Λe∂xξe+Beξe=0for x∈[0,l],(26)Be=be11be12be21be22≔12ReLe-1+ZeCe-1ReLe-1-ZeCe-1ReLe-1-ZeCe-1ReLe-1+ZeCe-1,Λe=λe+00λe-,withλe±=±LeCe-1. For a directed graph with E edges we are able to rewrite (24) in a compact form provided that D± is independent of v∈V. We collect ξ=(ξ+,ξ-)T, ξ±=ξ1±,ξ2±,…,ξE± and introduce the matrices(27)B=B1+,B2+,…,BE+,B1-,B2-,…,BE-∈R2E×2Eand Λ=diagΛ+,Λ-∈R2E×2E with (28)Be+=0,…,0︸e-1,be11,0,…,0︸E-1,be12,0,…,0︸E-e∈R2E,Be-=0,…,0︸e-1,be21,0,…,0︸E-1,be22,0,…,0︸E-e∈R2E,Λ±=diagλ1±,λ2±,…,λE±∈RE×E. The coupling condition is rewritten in terms of D±, where D±e,e~=dee~±λe~±, ∀e,e~∈E, with D±∈RE×E. Then, the network equations can be written in a compact form such as (21) and (24).Remark 4.
Note that the time step sizeΔt 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, ΔtBeTBe-2Be is negative definite for all edges e. This is the case for Δt≤2minCe/Ze,Le/Re (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,l=1/2, and also discretize the edges with an equidistant grid. We use the final time horizon T=12. If some parameters change throughout the section, it will be stated. For simplicity, we skip the index e 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 Figure1.
We intend to introduce the network framework introduced above step by step. On each edgee, we have two quantities ξe+ and ξe-. We first choose Ze=Re=0 and Le=Ce=1 for all e. If Ze=Re=0, the telegrapher’s equation reduces to the linear wave equation. Also note that λe+=1 and λe-=-1. We control the inflow of ξ1+ at node v0 and the inflow of ξ2- at node v2 depending on the outflow at v2 and v0, respectively; that is, ξ1+(t,0)=g1ξ2+(t,l) and ξ2-(t,l)=g2ξ1-(t,0). The matrices G+ and G- are composed of(29)G+=0g110,G-=01g20.The corresponding eigenvalues of G+TP0+G+-PN+ and G-TPN-1-G--P-1- (cf. (C1) in Theorem 2) are (30)κ1=-exp-μx2,0+x1,Nexpμx2,0-expμx1,N,κ2=-exp-μx1,0+x2,Nexpμx1,0-g12expμx2,N,for G+TP0+G+-PN+, and(31)κ3=expμx1,N-1-expμx2,-1,κ4=g22expμx2,N-1-expμx1,-1, for G-TPN-1-G--P-1-. If we assume that μ=μe for all e and insert the values for the spatial discretization x1,0=0, x1,N=x2,0=1/2, x2,N=1 on the strips [0,1/2] and [1/2,1], we get(32)κ1=0,κ2=-exp-μ1-g12expμ,κ3=0,κ4=g22expμ-1.
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 g1, g2. So we have to ensure that the matrices G+TP0+G+-PN+ and G-TPN-1-G--P-1- are negative definite for the feedback relevant eigenvalues κ2, κ4. We claim κ2,κ4<0 which is the case if(33)μ<minln1g12,ln1g22.Since μ>0 it should also hold that 0<g12<1 and 0<g22<1.
For the numerical computations we set as initial valuesξe+(0,x)=-1/2, ξe-(0,x)=1/2. We will vary the values for g1 and g2 and expect a rapidly dropping discrete Lyapunov function Vk (14) for small g1 and g2 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, Δx=1/1600, the decay rate ν converges to μ for |λe|=1 due to (19).Table 1
The number of cells on each edge isN=1600, with Δx=1/1600. The CFL constant is equal to one; that is, Δt/Δx=1.
g
1 and g2
V
T
/
V
0
μ
ν
0.1
1.00
E
-
24
4.61
E
+
00
4.60
E
+
00
0.25
3.55
E
-
15
2.77
E
+
00
2.77
E
+
00
0.5
5.96
E
-
08
1.39
E
+
00
1.39
E
+
00
0.75
1.00
E
-
03
5.75
E
-
01
5.75
E
-
01
0.9
7.98
E
-
02
2.10
E
-
01
2.10
E
-
01Figure 1
Linear network without branches.Example 6.
As a more complex network, we consider the cascade network in Figure2 for a different number of layers n. 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=7n. The network has controls (inputs) us± at the outward nodes for s∈{1,2,3}.
By numbering the edges from left to right and from top to bottom, we have(34)u1+=g12ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u2+=g13ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u3+=g12ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u1-=g22ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-,u2-=g23ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-,u3-=g22ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-. The matrices describing the network with n layers have the following form:(35)G+=0⋯0G1M⋱0⋮⋱⋱⋮0⋯M0︸n+1,G-=0MT⋯0⋮⋱⋱⋮0⋱MTG20⋯0︸n+1 with ι=1,2(36)Gι=gι2gι2gι3gι3gι3gι2gι2⋮⋮⋮⋮⋮⋮⋮gι2gι2gι3gι3gι3gι2gι2T,M=120120000120120000013013013001301301300130130130000012012000012012.Here, G±∈RE×E and G1,G2,M∈R7×7. We set μe=μ for all edges e. The parameters are Le=Ce=1, Ze=Re=0 for all edges. Then we have the following matrices:(37)P0+=diagP~0+,…,P~1/2n+,PN+=diagP~1/2+,…,P~1/2n+1+,P-1-=diagP~0-,…,P~1/2n-,PN-1-=diagP~1/2-,…,P~1/2n+1-, with(38)P~n±=diagexp∓μn,…,exp∓μn∈R7×7. So the eigenvalues of G+TP0+G+-PN+ and of G-TPN-1-G--P-1- are determined by the eigenvalues of the matrices,(39)MTP~n+M-P~n+,MP~n-MT-P~n-,G1TP~0+G1-P~1/2n+1+,G2TP~1/2n+1-G2-P~0-. The eigenvalues of the first two matrices are in 0,-exp±μn and the eigenvalues of the last two matrices are in (40)-1,-exp-μ12n+1,283g12-exp-μ12n+1,283g22expμ12n+1-1. With the same arguments used before, we restrict our investigations again to the feedback relevant eigenvalues. That means if the condition(41)μ<2n+1ln328maxg1,g2-2is satisfied, the assumption (C1) of Theorem 2 leads to negative definite matrices. Since μ has to be greater than zero, 3/28max(g1,g2)-2 must be greater than one. So max(g1,g2)2<3/28. Therefore, we might choose g1=g2=1/4. We expect that in our numerical experiments μ goes to zero for increasing layers n→∞; compare (41) and results in Table 2. This means, for a large number of layers, we can only expect a linear decreasing Lyapunov function Vk (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 ξe+(0,x)=-2 and ξe-(0,x)=2 for all e.Table 2
The number of cells on each edge isN=200, with Δx=1/200. The CFL constant is equal to one; that is, Δt/Δx=1.
j
V
T
/
V
0
μ
ν
1
9.13
E
-
23
5.39
E
-
01
5.38
E
-
01
2
1.13
E
-
16
3.59
E
-
01
3.59
E
-
01
4
1.90
E
-
09
2.16
E
-
01
2.16
E
-
01
8
2.31
E
-
05
1.12
E
-
01
1.12
E
-
01
16
3.50
E
-
03
6.34
E
-
02
6.34
E
-
02Figure 2
Cascade network withn layers and six controls us±, s=1,2,3.


### 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 law(42)∂tρe+∂xρeVρe=0∀e,where ρe describes the density of vehicles on each edge e.Here, the velocity of a car depends on the local density, that is,Vρe, and the flux function qe=fe(ρe)=ρeV(ρe)=ρe(Vmax(1-ρe/ρmax)) is a concave function, with fe(0)=0 and fe(ρmax)=0, which is maximal at a value ρc. These values separate a free-flow state, that is, ρe≤ρc from a traffic-congestion state, that is, ρe>ρc. 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,(43)∂tqe+∂qefe-1qe-1∂xqe=0, which is useful to define the coupling conditions. The coupling conditions in a road network have then the following componentwise form(44)qet,0=∑dee~+qe~t,l+us+t,qet,l=∑dee~-qe~t,0+us-t, with distribution rates dee~+ and dee~- and some additional input and output, either inflows us+(t) from outside ramps or outflows us-(t) 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(45)qet,0=∑de~qee~t,l+ust, with distribution rates dee~ and inflows us(t). Note, it has to hold that ∑e~‍dee~=1, in the sense of flux conservation. So considering a network G=(V,E) and the additional vector notation (46)D=matrix with entriesdee~,u=u1,u2,…,uET,q=q1,…,qET,Fq=diag∂q1f1-1q1-1,∂q2f2-1q2-1,…,∂qEfE-1qE-1, we get the whole LWR-networked system(47)∂tq+Fq∂xq=0,qt,0=Dqt,l+ut,q0,x=q0x.For a given stateρ¯e of (42), which is perturbed with a perturbation dρe at some point in time and by linearization of the flux function fe, we can rewrite the LWR-model (42) as a conservation law of the perturbation dρe as follows: (48)∂tρ¯e+dρe+∂xfρ¯e+dρe=0. A short calculation using (42) leads to(49)∂tdρe+f′ρ¯e∂xdρe+f′′ρ¯e∂xρ¯edρe=0.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 Figure3. Here, u1 is just the inflow and u2 is a control input. Our goal is to analyze how we can control the inflow, such that small perturbations will vanish after some time.Figure 3
Road network with|E|=5 edges and two controls u1, u2.As shown in the previous example, we can rewrite the model (42) in the form of kinetic wave equations,(50)∂tqe+∂qefe-1qe-1∂xqe=0, and couple the resulting conservation law as follows: (51)q1t,0=u1t,q2t,0=d12q1t,l,q3t,0=d23q2t,l+u2t,q4t,0=1-d23q2t,l+u2t,q5t,0=1-d12q1t,l. For given inputs u¯1 and u¯2, there also exists a free-flow steady-state q¯1,q¯2,…,q¯5. We can then define a feedback control law, such as (52)u1t=u¯1,u2t=u¯2+g1ξ3t,l with ξe=qe-q¯e. 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(53)Δtf′′ρ¯e∂xρ¯e2-2f′′ρ¯e∂xρ¯e<0. Since f′′ρ¯e is a negative constant for concave flux functions, it follows that ∂xρ¯e has to be less than or equal to zero. For our numerical simulations, we use the flux function f(ρe)=ρe1-ρe; that is, Vmax=1 and ρmax=1 for all edges. Hence,(54)4Δt∂xρ¯e2+4∂xρ¯e≤0, if -1/Δt<∂xρ¯e<0.By linearizing the coupling conditions above, we get(55)f1′ρ¯1t,0dρ1t,0=0,f2′ρ¯2t,0dρ2t,0=d12f1′ρ¯1t,ldρ1t,l,f3′ρ¯3t,0dρ3t,0=d23f2′ρ¯2t,ldρ2t,l+g1f3′ρ¯3t,ldρ3t,l,f4′ρ¯4t,0dρ4t,0=1-d23f2′ρ¯2t,ldρ2t,l+g1f3′ρ¯3t,ldρ3t,l,f5′ρ¯5t,0dρ5t,0=1-d12f1′ρ¯1t,ldρ1t,l. We set d12=0.8 and d23=0.5. For the initial values we choose(56)q10,x=0.2,q20,x=0.16,q30,x=0.08·1-x,q40,x=0.08·1-x,q50,x=0.04·1-x. We can now compute the matrix(57)G=000000.6000000.360.6·g10000.360.6·g1000.10000. We also set μe=μ=0.5 and l=1 for all edges e. Then, the matrices P0 and PN are (58)P0=diag1,exp-μ,exp-2μ,exp-2μ,exp-μ,PN=diagexp-μ,exp-2μ,exp-3μ,exp-3μ,exp-2μ. So with g1=3/4 the matrix GTP0G-PN is negative definite. For the simulation we set u¯1=0.2, u¯2=0, and Δx=0.01. Also the perturbation dρe is set to zero for e∈{1,2,4,5} and dρ3(0,x)=sin2πx/5. 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.Figure 4
(a) The decay of the initial perturbations over time at the end of edge 3 and edge 4 in log-scale. (b) The evolution of the perturbation in time and space.
(a)
(b)Figure 5
Evolution of the perturbation withΔx=0.01.


## 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 currentI(t,x) and voltage U(t,x) are the conservation quantities that are described by the telegrapher’s equation. They are expressed in terms of characteristic variables as I(t,x)=ξ+(t,x)+ξ-(t,x) and U(t,x)=L/C(ξ+(t,x)-ξ-(t,x)). Due to the losses alongside the lines, we have an additional reaction term B depending on the positive constants of R, L, and Z (resistance, inductance, and admittance) and C (capacitance). The equation written in characteristic variables on edge e is then(25)∂tξe+Λe∂xξe+Beξe=0for x∈[0,l],(26)Be=be11be12be21be22≔12ReLe-1+ZeCe-1ReLe-1-ZeCe-1ReLe-1-ZeCe-1ReLe-1+ZeCe-1,Λe=λe+00λe-,withλe±=±LeCe-1. For a directed graph with E edges we are able to rewrite (24) in a compact form provided that D± is independent of v∈V. We collect ξ=(ξ+,ξ-)T, ξ±=ξ1±,ξ2±,…,ξE± and introduce the matrices(27)B=B1+,B2+,…,BE+,B1-,B2-,…,BE-∈R2E×2Eand Λ=diagΛ+,Λ-∈R2E×2E with (28)Be+=0,…,0︸e-1,be11,0,…,0︸E-1,be12,0,…,0︸E-e∈R2E,Be-=0,…,0︸e-1,be21,0,…,0︸E-1,be22,0,…,0︸E-e∈R2E,Λ±=diagλ1±,λ2±,…,λE±∈RE×E. The coupling condition is rewritten in terms of D±, where D±e,e~=dee~±λe~±, ∀e,e~∈E, with D±∈RE×E. Then, the network equations can be written in a compact form such as (21) and (24).Remark 4.
Note that the time step sizeΔt 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, ΔtBeTBe-2Be is negative definite for all edges e. This is the case for Δt≤2minCe/Ze,Le/Re (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,l=1/2, and also discretize the edges with an equidistant grid. We use the final time horizon T=12. If some parameters change throughout the section, it will be stated. For simplicity, we skip the index e 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 Figure1.
We intend to introduce the network framework introduced above step by step. On each edgee, we have two quantities ξe+ and ξe-. We first choose Ze=Re=0 and Le=Ce=1 for all e. If Ze=Re=0, the telegrapher’s equation reduces to the linear wave equation. Also note that λe+=1 and λe-=-1. We control the inflow of ξ1+ at node v0 and the inflow of ξ2- at node v2 depending on the outflow at v2 and v0, respectively; that is, ξ1+(t,0)=g1ξ2+(t,l) and ξ2-(t,l)=g2ξ1-(t,0). The matrices G+ and G- are composed of(29)G+=0g110,G-=01g20.The corresponding eigenvalues of G+TP0+G+-PN+ and G-TPN-1-G--P-1- (cf. (C1) in Theorem 2) are (30)κ1=-exp-μx2,0+x1,Nexpμx2,0-expμx1,N,κ2=-exp-μx1,0+x2,Nexpμx1,0-g12expμx2,N,for G+TP0+G+-PN+, and(31)κ3=expμx1,N-1-expμx2,-1,κ4=g22expμx2,N-1-expμx1,-1, for G-TPN-1-G--P-1-. If we assume that μ=μe for all e and insert the values for the spatial discretization x1,0=0, x1,N=x2,0=1/2, x2,N=1 on the strips [0,1/2] and [1/2,1], we get(32)κ1=0,κ2=-exp-μ1-g12expμ,κ3=0,κ4=g22expμ-1.
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 g1, g2. So we have to ensure that the matrices G+TP0+G+-PN+ and G-TPN-1-G--P-1- are negative definite for the feedback relevant eigenvalues κ2, κ4. We claim κ2,κ4<0 which is the case if(33)μ<minln1g12,ln1g22.Since μ>0 it should also hold that 0<g12<1 and 0<g22<1.
For the numerical computations we set as initial valuesξe+(0,x)=-1/2, ξe-(0,x)=1/2. We will vary the values for g1 and g2 and expect a rapidly dropping discrete Lyapunov function Vk (14) for small g1 and g2 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, Δx=1/1600, the decay rate ν converges to μ for |λe|=1 due to (19).Table 1
The number of cells on each edge isN=1600, with Δx=1/1600. The CFL constant is equal to one; that is, Δt/Δx=1.
g
1 and g2
V
T
/
V
0
μ
ν
0.1
1.00
E
-
24
4.61
E
+
00
4.60
E
+
00
0.25
3.55
E
-
15
2.77
E
+
00
2.77
E
+
00
0.5
5.96
E
-
08
1.39
E
+
00
1.39
E
+
00
0.75
1.00
E
-
03
5.75
E
-
01
5.75
E
-
01
0.9
7.98
E
-
02
2.10
E
-
01
2.10
E
-
01Figure 1
Linear network without branches.Example 6.
As a more complex network, we consider the cascade network in Figure2 for a different number of layers n. 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=7n. The network has controls (inputs) us± at the outward nodes for s∈{1,2,3}.
By numbering the edges from left to right and from top to bottom, we have(34)u1+=g12ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u2+=g13ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u3+=g12ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u1-=g22ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-,u2-=g23ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-,u3-=g22ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-. The matrices describing the network with n layers have the following form:(35)G+=0⋯0G1M⋱0⋮⋱⋱⋮0⋯M0︸n+1,G-=0MT⋯0⋮⋱⋱⋮0⋱MTG20⋯0︸n+1 with ι=1,2(36)Gι=gι2gι2gι3gι3gι3gι2gι2⋮⋮⋮⋮⋮⋮⋮gι2gι2gι3gι3gι3gι2gι2T,M=120120000120120000013013013001301301300130130130000012012000012012.Here, G±∈RE×E and G1,G2,M∈R7×7. We set μe=μ for all edges e. The parameters are Le=Ce=1, Ze=Re=0 for all edges. Then we have the following matrices:(37)P0+=diagP~0+,…,P~1/2n+,PN+=diagP~1/2+,…,P~1/2n+1+,P-1-=diagP~0-,…,P~1/2n-,PN-1-=diagP~1/2-,…,P~1/2n+1-, with(38)P~n±=diagexp∓μn,…,exp∓μn∈R7×7. So the eigenvalues of G+TP0+G+-PN+ and of G-TPN-1-G--P-1- are determined by the eigenvalues of the matrices,(39)MTP~n+M-P~n+,MP~n-MT-P~n-,G1TP~0+G1-P~1/2n+1+,G2TP~1/2n+1-G2-P~0-. The eigenvalues of the first two matrices are in 0,-exp±μn and the eigenvalues of the last two matrices are in (40)-1,-exp-μ12n+1,283g12-exp-μ12n+1,283g22expμ12n+1-1. With the same arguments used before, we restrict our investigations again to the feedback relevant eigenvalues. That means if the condition(41)μ<2n+1ln328maxg1,g2-2is satisfied, the assumption (C1) of Theorem 2 leads to negative definite matrices. Since μ has to be greater than zero, 3/28max(g1,g2)-2 must be greater than one. So max(g1,g2)2<3/28. Therefore, we might choose g1=g2=1/4. We expect that in our numerical experiments μ goes to zero for increasing layers n→∞; compare (41) and results in Table 2. This means, for a large number of layers, we can only expect a linear decreasing Lyapunov function Vk (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 ξe+(0,x)=-2 and ξe-(0,x)=2 for all e.Table 2
The number of cells on each edge isN=200, with Δx=1/200. The CFL constant is equal to one; that is, Δt/Δx=1.
j
V
T
/
V
0
μ
ν
1
9.13
E
-
23
5.39
E
-
01
5.38
E
-
01
2
1.13
E
-
16
3.59
E
-
01
3.59
E
-
01
4
1.90
E
-
09
2.16
E
-
01
2.16
E
-
01
8
2.31
E
-
05
1.12
E
-
01
1.12
E
-
01
16
3.50
E
-
03
6.34
E
-
02
6.34
E
-
02Figure 2
Cascade network withn layers and six controls us±, s=1,2,3.


## 3.1.1. Numerical Results

For all computations we choose the length of all edges as one over two, that is,l=1/2, and also discretize the edges with an equidistant grid. We use the final time horizon T=12. If some parameters change throughout the section, it will be stated. For simplicity, we skip the index e 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 Figure1.
We intend to introduce the network framework introduced above step by step. On each edgee, we have two quantities ξe+ and ξe-. We first choose Ze=Re=0 and Le=Ce=1 for all e. If Ze=Re=0, the telegrapher’s equation reduces to the linear wave equation. Also note that λe+=1 and λe-=-1. We control the inflow of ξ1+ at node v0 and the inflow of ξ2- at node v2 depending on the outflow at v2 and v0, respectively; that is, ξ1+(t,0)=g1ξ2+(t,l) and ξ2-(t,l)=g2ξ1-(t,0). The matrices G+ and G- are composed of(29)G+=0g110,G-=01g20.The corresponding eigenvalues of G+TP0+G+-PN+ and G-TPN-1-G--P-1- (cf. (C1) in Theorem 2) are (30)κ1=-exp-μx2,0+x1,Nexpμx2,0-expμx1,N,κ2=-exp-μx1,0+x2,Nexpμx1,0-g12expμx2,N,for G+TP0+G+-PN+, and(31)κ3=expμx1,N-1-expμx2,-1,κ4=g22expμx2,N-1-expμx1,-1, for G-TPN-1-G--P-1-. If we assume that μ=μe for all e and insert the values for the spatial discretization x1,0=0, x1,N=x2,0=1/2, x2,N=1 on the strips [0,1/2] and [1/2,1], we get(32)κ1=0,κ2=-exp-μ1-g12expμ,κ3=0,κ4=g22expμ-1.
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 g1, g2. So we have to ensure that the matrices G+TP0+G+-PN+ and G-TPN-1-G--P-1- are negative definite for the feedback relevant eigenvalues κ2, κ4. We claim κ2,κ4<0 which is the case if(33)μ<minln1g12,ln1g22.Since μ>0 it should also hold that 0<g12<1 and 0<g22<1.
For the numerical computations we set as initial valuesξe+(0,x)=-1/2, ξe-(0,x)=1/2. We will vary the values for g1 and g2 and expect a rapidly dropping discrete Lyapunov function Vk (14) for small g1 and g2 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, Δx=1/1600, the decay rate ν converges to μ for |λe|=1 due to (19).Table 1
The number of cells on each edge isN=1600, with Δx=1/1600. The CFL constant is equal to one; that is, Δt/Δx=1.
g
1 and g2
V
T
/
V
0
μ
ν
0.1
1.00
E
-
24
4.61
E
+
00
4.60
E
+
00
0.25
3.55
E
-
15
2.77
E
+
00
2.77
E
+
00
0.5
5.96
E
-
08
1.39
E
+
00
1.39
E
+
00
0.75
1.00
E
-
03
5.75
E
-
01
5.75
E
-
01
0.9
7.98
E
-
02
2.10
E
-
01
2.10
E
-
01Figure 1
Linear network without branches.Example 6.
As a more complex network, we consider the cascade network in Figure2 for a different number of layers n. 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=7n. The network has controls (inputs) us± at the outward nodes for s∈{1,2,3}.
By numbering the edges from left to right and from top to bottom, we have(34)u1+=g12ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u2+=g13ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u3+=g12ξE-6++ξE-5++ξE-4++ξE-3++ξE-2++ξE-1++ξE+,u1-=g22ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-,u2-=g23ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-,u3-=g22ξ1-+ξ2-+ξ3-+ξ4-+ξ5-+ξ6-+ξ7-. The matrices describing the network with n layers have the following form:(35)G+=0⋯0G1M⋱0⋮⋱⋱⋮0⋯M0︸n+1,G-=0MT⋯0⋮⋱⋱⋮0⋱MTG20⋯0︸n+1 with ι=1,2(36)Gι=gι2gι2gι3gι3gι3gι2gι2⋮⋮⋮⋮⋮⋮⋮gι2gι2gι3gι3gι3gι2gι2T,M=120120000120120000013013013001301301300130130130000012012000012012.Here, G±∈RE×E and G1,G2,M∈R7×7. We set μe=μ for all edges e. The parameters are Le=Ce=1, Ze=Re=0 for all edges. Then we have the following matrices:(37)P0+=diagP~0+,…,P~1/2n+,PN+=diagP~1/2+,…,P~1/2n+1+,P-1-=diagP~0-,…,P~1/2n-,PN-1-=diagP~1/2-,…,P~1/2n+1-, with(38)P~n±=diagexp∓μn,…,exp∓μn∈R7×7. So the eigenvalues of G+TP0+G+-PN+ and of G-TPN-1-G--P-1- are determined by the eigenvalues of the matrices,(39)MTP~n+M-P~n+,MP~n-MT-P~n-,G1TP~0+G1-P~1/2n+1+,G2TP~1/2n+1-G2-P~0-. The eigenvalues of the first two matrices are in 0,-exp±μn and the eigenvalues of the last two matrices are in (40)-1,-exp-μ12n+1,283g12-exp-μ12n+1,283g22expμ12n+1-1. With the same arguments used before, we restrict our investigations again to the feedback relevant eigenvalues. That means if the condition(41)μ<2n+1ln328maxg1,g2-2is satisfied, the assumption (C1) of Theorem 2 leads to negative definite matrices. Since μ has to be greater than zero, 3/28max(g1,g2)-2 must be greater than one. So max(g1,g2)2<3/28. Therefore, we might choose g1=g2=1/4. We expect that in our numerical experiments μ goes to zero for increasing layers n→∞; compare (41) and results in Table 2. This means, for a large number of layers, we can only expect a linear decreasing Lyapunov function Vk (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 ξe+(0,x)=-2 and ξe-(0,x)=2 for all e.Table 2
The number of cells on each edge isN=200, with Δx=1/200. The CFL constant is equal to one; that is, Δt/Δx=1.
j
V
T
/
V
0
μ
ν
1
9.13
E
-
23
5.39
E
-
01
5.38
E
-
01
2
1.13
E
-
16
3.59
E
-
01
3.59
E
-
01
4
1.90
E
-
09
2.16
E
-
01
2.16
E
-
01
8
2.31
E
-
05
1.12
E
-
01
1.12
E
-
01
16
3.50
E
-
03
6.34
E
-
02
6.34
E
-
02Figure 2
Cascade network withn layers and six controls us±, s=1,2,3.


## 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 law(42)∂tρe+∂xρeVρe=0∀e,where ρe describes the density of vehicles on each edge e.Here, the velocity of a car depends on the local density, that is,Vρe, and the flux function qe=fe(ρe)=ρeV(ρe)=ρe(Vmax(1-ρe/ρmax)) is a concave function, with fe(0)=0 and fe(ρmax)=0, which is maximal at a value ρc. These values separate a free-flow state, that is, ρe≤ρc from a traffic-congestion state, that is, ρe>ρc. 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,(43)∂tqe+∂qefe-1qe-1∂xqe=0, which is useful to define the coupling conditions. The coupling conditions in a road network have then the following componentwise form(44)qet,0=∑dee~+qe~t,l+us+t,qet,l=∑dee~-qe~t,0+us-t, with distribution rates dee~+ and dee~- and some additional input and output, either inflows us+(t) from outside ramps or outflows us-(t) 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(45)qet,0=∑de~qee~t,l+ust, with distribution rates dee~ and inflows us(t). Note, it has to hold that ∑e~‍dee~=1, in the sense of flux conservation. So considering a network G=(V,E) and the additional vector notation (46)D=matrix with entriesdee~,u=u1,u2,…,uET,q=q1,…,qET,Fq=diag∂q1f1-1q1-1,∂q2f2-1q2-1,…,∂qEfE-1qE-1, we get the whole LWR-networked system(47)∂tq+Fq∂xq=0,qt,0=Dqt,l+ut,q0,x=q0x.For a given stateρ¯e of (42), which is perturbed with a perturbation dρe at some point in time and by linearization of the flux function fe, we can rewrite the LWR-model (42) as a conservation law of the perturbation dρe as follows: (48)∂tρ¯e+dρe+∂xfρ¯e+dρe=0. A short calculation using (42) leads to(49)∂tdρe+f′ρ¯e∂xdρe+f′′ρ¯e∂xρ¯edρe=0.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 Figure3. Here, u1 is just the inflow and u2 is a control input. Our goal is to analyze how we can control the inflow, such that small perturbations will vanish after some time.Figure 3
Road network with|E|=5 edges and two controls u1, u2.As shown in the previous example, we can rewrite the model (42) in the form of kinetic wave equations,(50)∂tqe+∂qefe-1qe-1∂xqe=0, and couple the resulting conservation law as follows: (51)q1t,0=u1t,q2t,0=d12q1t,l,q3t,0=d23q2t,l+u2t,q4t,0=1-d23q2t,l+u2t,q5t,0=1-d12q1t,l. For given inputs u¯1 and u¯2, there also exists a free-flow steady-state q¯1,q¯2,…,q¯5. We can then define a feedback control law, such as (52)u1t=u¯1,u2t=u¯2+g1ξ3t,l with ξe=qe-q¯e. 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(53)Δtf′′ρ¯e∂xρ¯e2-2f′′ρ¯e∂xρ¯e<0. Since f′′ρ¯e is a negative constant for concave flux functions, it follows that ∂xρ¯e has to be less than or equal to zero. For our numerical simulations, we use the flux function f(ρe)=ρe1-ρe; that is, Vmax=1 and ρmax=1 for all edges. Hence,(54)4Δt∂xρ¯e2+4∂xρ¯e≤0, if -1/Δt<∂xρ¯e<0.By linearizing the coupling conditions above, we get(55)f1′ρ¯1t,0dρ1t,0=0,f2′ρ¯2t,0dρ2t,0=d12f1′ρ¯1t,ldρ1t,l,f3′ρ¯3t,0dρ3t,0=d23f2′ρ¯2t,ldρ2t,l+g1f3′ρ¯3t,ldρ3t,l,f4′ρ¯4t,0dρ4t,0=1-d23f2′ρ¯2t,ldρ2t,l+g1f3′ρ¯3t,ldρ3t,l,f5′ρ¯5t,0dρ5t,0=1-d12f1′ρ¯1t,ldρ1t,l. We set d12=0.8 and d23=0.5. For the initial values we choose(56)q10,x=0.2,q20,x=0.16,q30,x=0.08·1-x,q40,x=0.08·1-x,q50,x=0.04·1-x. We can now compute the matrix(57)G=000000.6000000.360.6·g10000.360.6·g1000.10000. We also set μe=μ=0.5 and l=1 for all edges e. Then, the matrices P0 and PN are (58)P0=diag1,exp-μ,exp-2μ,exp-2μ,exp-μ,PN=diagexp-μ,exp-2μ,exp-3μ,exp-3μ,exp-2μ. So with g1=3/4 the matrix GTP0G-PN is negative definite. For the simulation we set u¯1=0.2, u¯2=0, and Δx=0.01. Also the perturbation dρe is set to zero for e∈{1,2,4,5} and dρ3(0,x)=sin2πx/5. 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.Figure 4
(a) The decay of the initial perturbations over time at the end of edge 3 and edge 4 in log-scale. (b) The evolution of the perturbation in time and space.
(a)
(b)Figure 5
Evolution of the perturbation withΔx=0.01.


## 3.2.1. Numerical Results

We consider the LWR-model on the network in Figure3. Here, u1 is just the inflow and u2 is a control input. Our goal is to analyze how we can control the inflow, such that small perturbations will vanish after some time.Figure 3
Road network with|E|=5 edges and two controls u1, u2.As shown in the previous example, we can rewrite the model (42) in the form of kinetic wave equations,(50)∂tqe+∂qefe-1qe-1∂xqe=0, and couple the resulting conservation law as follows: (51)q1t,0=u1t,q2t,0=d12q1t,l,q3t,0=d23q2t,l+u2t,q4t,0=1-d23q2t,l+u2t,q5t,0=1-d12q1t,l. For given inputs u¯1 and u¯2, there also exists a free-flow steady-state q¯1,q¯2,…,q¯5. We can then define a feedback control law, such as (52)u1t=u¯1,u2t=u¯2+g1ξ3t,l with ξe=qe-q¯e. 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(53)Δtf′′ρ¯e∂xρ¯e2-2f′′ρ¯e∂xρ¯e<0. Since f′′ρ¯e is a negative constant for concave flux functions, it follows that ∂xρ¯e has to be less than or equal to zero. For our numerical simulations, we use the flux function f(ρe)=ρe1-ρe; that is, Vmax=1 and ρmax=1 for all edges. Hence,(54)4Δt∂xρ¯e2+4∂xρ¯e≤0, if -1/Δt<∂xρ¯e<0.By linearizing the coupling conditions above, we get(55)f1′ρ¯1t,0dρ1t,0=0,f2′ρ¯2t,0dρ2t,0=d12f1′ρ¯1t,ldρ1t,l,f3′ρ¯3t,0dρ3t,0=d23f2′ρ¯2t,ldρ2t,l+g1f3′ρ¯3t,ldρ3t,l,f4′ρ¯4t,0dρ4t,0=1-d23f2′ρ¯2t,ldρ2t,l+g1f3′ρ¯3t,ldρ3t,l,f5′ρ¯5t,0dρ5t,0=1-d12f1′ρ¯1t,ldρ1t,l. We set d12=0.8 and d23=0.5. For the initial values we choose(56)q10,x=0.2,q20,x=0.16,q30,x=0.08·1-x,q40,x=0.08·1-x,q50,x=0.04·1-x. We can now compute the matrix(57)G=000000.6000000.360.6·g10000.360.6·g1000.10000. We also set μe=μ=0.5 and l=1 for all edges e. Then, the matrices P0 and PN are (58)P0=diag1,exp-μ,exp-2μ,exp-2μ,exp-μ,PN=diagexp-μ,exp-2μ,exp-3μ,exp-3μ,exp-2μ. So with g1=3/4 the matrix GTP0G-PN is negative definite. For the simulation we set u¯1=0.2, u¯2=0, and Δx=0.01. Also the perturbation dρe is set to zero for e∈{1,2,4,5} and dρ3(0,x)=sin2πx/5. 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.Figure 4
(a) The decay of the initial perturbations over time at the end of edge 3 and edge 4 in log-scale. (b) The evolution of the perturbation in time and space.
(a)
(b)Figure 5
Evolution of the perturbation withΔx=0.01.


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



---

*Source: 6896153-2017-05-22.xml*