DISCRETIZED FEEDBACK CONTROL FOR SYSTEMS OF LINEARIZED HYPERBOLIC BALANCE LAWS

. Physical systems such as water and gas networks are usually oper-ated in a state of equilibrium and feedback control is employed to damp small perturbations over time. We consider ﬂow problems on networks, described by hyperbolic balance laws, and analyze the stability of the linearized systems. Suﬃcient conditions for exponential stability in the continuous and discretized setting are presented. The analysis is extended to arbitrary Sobolev norms. Computational experiments illustrate the theoretical ﬁndings.

1. Introduction. Hyperbolic balance laws can be used to model flow dynamics on networks. Isothermal Euler and shallow water equations form 2×2 hyperbolic systems that describe the temporal and spatial evolution of gas and water flow. Boundary control of such systems is subject of current research [4,6].
An underlying tool for the study of these problems are Lyapunov functions that yield upper bounds for the deviation from steady states in suitable norms, e.g. for the Lebesgue norm L 2 and Sobolev norm H 2 . Exponential decay of a continuous Lyapunov function under so-called dissipative boundary conditions has been proven in [9,12,7,8]. In particular, analytical results have been presented in the case of gas flow [16,2,15,20] and water flow [5,14,17,23,11,19,21,18]. Comparisons to other stability concepts are presented in [10].
One may distinguish between uniform and space-varying steady states, conservation and balance laws, linear and nonlinear hyperbolic systems. A general theory for the stabilization of linear conservation laws is available [4,Sec. 3]. For nonlinear systems, however, results are still partial. A problem is posed by the fact that Lyapunov's indirect method [22] does not necessarily hold for hyperbolic systems. Furthermore, solutions of conservation laws exist in the classical sense only in finite time due to the occurence of shocks, even for smooth initial data [13]. Fortunately, if the linearized system is stable with respect to the Sobolev norm H 2 , one may obtain the stability of the nonlinear system with respect to the weaker L 2 -norm [8,4].
Most analytical results do not state explicitly decay rates of the Lyapunov function and the influence of the source term is assumed to be small or in intuitive terms, the considered balance laws are viewed as perturbations of conservation laws [9]. In practical applications, however, there may be a large influence of the source term.

Theoretical results and basic notations.
We briefly recall results from [4]. We consider physical systems that are described by a hyperbolic balance law ∂ ∂t with time and space variables (t, x) ∈ [0, ∞) × [0, L]. We assume a 2×2 system with strictly hyperbolic flux function f ∈ C 1 R 2 ; R 2 and source term S ∈ C 1 R 2 ; R 2 , but an extension to larger systems is possible. To simplify notation, we omit time and space variables and abbreviate system (1) by y t + f (y) x = −S(y), which is for smooth solutions equivalent to solving the quasilinear form where D y f (y) denotes the Jacobian of f (y). This system will be analyzed at a steady state.

Definition 2.1 (Steady State). A steady state (or equilibrium) is a time-invariant
possibly space-varying solutionȳ(x) of system (1). It satisfies the system of ordinary differential equations 0 = f (ȳ) x + S(ȳ).
Since we are interested in H 2 -stabilization, we assume for nowȳ ∈ H 2 (0, L). Sobolev's embedding theorem ensures a representantȳ ∈ C 1 (0, L). We first linearize the original system at steady state and then transform the corresponding linear system into Riemann invariants. The opposite approach -first transforming into Riemann invariants and then linearizing -has been discussed in [4]. Both approaches yield stability of the linearized system in Riemann coordinates. Assuming small and smooth perturbations ∆y(t, x) := y(t, x) −ȳ(x) of the steady state, the flux function and the source term is linearized by f (ȳ) + A∆y and S(ȳ) +S∆y, where A andS denote Jacobian at steady state. We denote by A x the matrix, where each entry of A(x) is differentiated with respect to x. Neglecting higher order derivatives, the linearized balance law (1) reads The matrices A(x),S(x) ∈ R 2×2 may depend on the space variable x ∈ [0, L]. For a given steady stateȳ we assume an initial perturbation ∆y(0, x) := ∆y 0 (x). Stabilization aims to damp the deviation ∆y.

Riemann invariants.
The Jacobian A in the quasilinear form (2) has an eigenvalue decomposition with real eigenvalues [13]. Thus, there exists a diagonal matrixΛ := diag λ+ ,λ − , whose entries are eigenvalues of A and a matrix T , whose columns r ± are the corresponding eigenvectors, such that We introduce the assumption which holds in the case of isothermal Euler and shallow water equations as long as the flow remains subsonic and subcritical, respectively. Again, the eigenvaluesλ ± (x) and eigenvectors r ± (x) may be space-varying. We introduce the Riemann invariant such that for smooth solutions equation (2) is equivalent to ζ t +Λζ x = −Cζ. We prescribe linear feedback boundary conditions by a matrixḠ ∈ R 2×2 such that the initial boundary value problem (IBVP) is well-posed [4].
The idea is to reduce H d -to L 2 -stabilization by introducing the Riemann invariant R := ζ (0) , . . . , ζ (d) T , which contains all derivatives up to order d. We have Note that there is no straightforward extension to nonlinear systems, where additional compatibility conditions must be fulfilled to ensure well-posedness of the initial boundary value problem [4].
In Theorem 2.3, we derive an IBVP for the quantities R ∈ R 2(d+1) . Additional initial and boundary conditions have to be imposed for the quantities ζ (1) , . . . , ζ (d) . Whereas the choice R (j) (0, x) = ∂ j ∂x j ζ(0, x) for the j-th component is natural, we state boundary conditions implicitly by assuming that the source term is neglectable small near and inside pipe intersections. Hence, we consider the advection equations to obtain the inflow boundary conditions. This assumption is motivated by the fact that in typical applications the physical size of junctions, e.g. compressor stations, are very small compared to a pipeline [26]. Furthermore, it is used for proving existence and uniqueness of Riemann problems at junctions [2]. Later in Section 4, we will use this assumption also to obtain a consistent numerical scheme similarly to [25].
with Λ := diag Λ , . . . ,Λ and R ± := ζ (0),± , . . . , ζ (d),± T . For C,Λ sufficiently smooth, the source term is a triangular block matrix of the form where ( * ) denotes non-zero entries. In particular for d = 4, the source term reads . Under the assumption (A2), problem (S d ) is well-posed for the boundary conditions where G (k) i,j denotes the i, j-th entry of the matrix Proof. We prove statement (3) by induction. The first derivative satisfies Let the claim hold for j − 1 with 2 ≤ j ≤ d. Then, we obtain i ∈ R 2×2 . The claim follows by induction.
Since each entry R (j) ∈ R 2 consists of one wave R (j),+ with positive characteristic speedλ + and R (j),− with negative speed, they can be reordered to R ± . Assumption (A2) implies inductively for j = 1, . . . , d More generally, we consider IBVP on networks with n arcs. We extend the notation to a set of Riemann invariants R ∈ R m . We assume positive characteristic speeds λ (j) (x) > 0 for j = 1, . . . , p and λ (j) (x) < 0 for j = p + 1, . . . , m.
3. Stabilization in the continuous case. In this section, we present sufficient conditions for the exponential stability of the system (S d ). As basic tool, we use a Lyapunov function similar to [16]. But in extension we equip the left boundary for Riemann invariants with positive characteristic speeds by some strictly positive weight h (j) 0 > 0 and the right boundary by h (j) L > 0 for negative speeds, respectively. For now, we assume that these weights are given. In Section 5, we show how they can be determined to ensure exponential stability. We refer to these weights as ghostcell-weights. The reason for this naming will become clear in Section 4.
with strictly positive entries such that for x ∈ [0, L] fixed, the weighted inner product STEPHAN GERSTER AND MICHAEL HERTY and the Lyapunov function We show conditions that guarantee a possibly negative decay rate µ ∈ R, which may defer from the constantμ > 0. The basic idea of our proof is to ensure nonnegative eigenvalues of a matrix M, which we will define later in Lemma 3.2. Then, we deduce R, MR ≥ 0 for all R ∈ R m . Note that this is only possible if the matrix M is symmetric. Therefore, we define the weighted source terms The source term B(x) and the weighted source term B(x,μ) represent the same weighted quadratic form, i.e. .

Lemma 3.2. Define the matrix
and let d (j) (x,μ) denote the eigenvalues of the matrixB(x,μ). The smallest eigenvalue d min (μ) := min We define the guaranteed rate µ g and the estimated rate µ e as Proof. The spectral radius σ max satisfies the inequality σ max {A} ≤ A p for any p-norm and matrix A.
and specifically for the diagonal matrix with positive entries Let 1 ∈ R m×m denote the identity matrix. The eigenvalue decomposition of the has non-negative eigenvalues. Therefore, the congruent and symmetric matrix M(x,μ,μ) is positive semidefinite.
The following theorem states sufficient conditions that yield the possibly negative lower bounds µ e ≤ µ g ≤ µ for the decay rate µ of the Lyapunov function.

contains the given ghostcell-weights
Define the matrices and let d min (μ) denote the smallest eigenvalues of the matrixB(x,μ). Define the guaranteed rate µ g and the estimated rate µ e as Then, the Lyapunov function satisfies the inequality and in particular, for the estimated rates If the decay rate µ > 0 is positive for all t ∈ R + 0 .
Proof. The analysis for the advection part R t + ΛR x = 0 can be reduced to a single For positive characteristic 524 STEPHAN GERSTER AND MICHAEL HERTY speeds j = 1, . . . , p we obtain by using integration by parts It follows analogously for j = p + 1, . . . , m By summing up the equations (5) and (6), we obtain The boundary terms (7) and (8) read and are non-positive if the matrix H(μ) is negative semidefinite. The Lyapunov function for the reaction part R t = −BR reads We sum up the equations (9) and (10) to obtain with Lemma 3.2 the estimate The H d -stability of the system (S) follows from the norm equivalence (∼) and The guaranteed decay rate not only depends on the influence of the source term, but also on the weights of the Lyapunov function. There is a positive decay rate µ > 0 if the source term or the condition number cond √ W (μ) is sufficiently small.
A source term is called diagonally stable if the matrixB(x,μ) and the congruent matrix W (x)B(x,μ) are positive semidefinite. Lyapunov stability for diagonally stable source terms has been shown in [3]. For these particular source terms and for conservation laws, there is the positive decay rate µ :=μ > 0 guaranteed. Note that a not necessarily symmetric source term with positive eigenvalues does not imply the positive definiteness of the matrixB(x,μ).
The important statement of Theorem 3.3 are the bounds µ ≥ µ g ≥ µ e . There may be cases, where the bounds µ g , µ e ∈ R are negative. Then, a decaying Lyapunov function is not guaranteed. We give in Section 6.3 an illustrate example for this case.
4. Stabilization in the discretized case. The explicit bounds µ e ≤ µ g ≤ µ hold for general balance laws. More accurate results can only be achieved by considering the numerical values of the matrices M and H. In this section, we discuss the stabilization of the discretized system that results from upwind and downwind schemes. We will see that a naive discretization may cause instabilities, since the definiteness property of the matrices M and H may not be inherited to their discretized analogues.
Recall that the H d -stability of the systemS results from the L 2 -stability of holds, where λ (j) are eigenvalues of (S d ). Cell averages at t k are approximated by The advection part can be approximated by a left and right sided upwind-scheme and the reaction part by the explicit Euler method, i.e.
where the diagonal matrix Λ i := diag{λ Note that this discretization is in accordance with assumption (A2). The boundary conditions affect the upwind scheme only. The splitting method accounts for the source term in the interior of the interval, but not at the ghostcells. For general systems of Riemann invariants R k i ∈ R m with i = 1, . . . , N we assume positive characteristic speed λ and a straightforward discretized analogue in t k would be which is proposed in [1,25]. In [1], a blow up in the discretized derivatives occurs in contrast to the continuous case. In [25], results on discretized uniform steady states are extended to diagonally stable source terms. To circumvent this problem, we approximate the continuous derivative The one-sided difference quotients (12) and (13) explain the choices (14) and (15) in the next definition. Like in the continuous case, the discrete weights are set by the ghostcell-weights, which we will explain in Section 5. N +1 > 0 for j = p + 1, . . . , m be given. We define the diagonal matrix h (j) We introduce for all i = 1, . . . , N the weighted inner product and the discrete Lyapunov function From now on we assume the restriction ∆x ∈ (0, λmin /μ) with λ min := min which ensures strictly positive weights W i such that the inner product (16) and consequently, the Lyapunov function (17) are well-defined. We show a discretized analogue of Theorem 3.3. .

Define the matrix
and assume that for all k, i, there exists θ < ∞ such that Then, the discretized Lyapunov function satisfies the inequality Proof. The numerical scheme (11) reads in splitting form whereR k i is the discretization of the upwind-scheme. The analysis for the advection part is reduced to a single discretized Riemann invariant L k,(j) Then, the time derivative of the discretized Lyapunov function satisfies With similar calculations for j = p + 1, . . . , m we obtain altogether Analogously to the continuous case, boundary terms (19) and (20) are non-positive ifĤ is negative semidefinite. Assumption (18) yields Finally, we consider the full discretization R k+1 Estimates (21), (22) and the positive semidefiniteness of the matrix M i (µ,μ) imply The claim follows from We end this section by discussing the assumptions of Theorem 4.2. Assumption (18) is a restriction on the spatial discretization. For the given CFL-condition and given θ > 0, it is satisfied if the space discretization ∆x > 0 is small enough, since for sufficiently smooth solutions the left-hand sides scales as O(∆x). Since the ratio ∆t /∆x is fixed, the influence of the term θ ∆t ∆x W . Therefore, this term does not cause an additional restriction compared to Theorem 3.3 in the formal limit. Note that we have to account for the steady state R k i = 0 and hence, assumption (18) does not simplify to If the discretization is in the stability region of the explicit Euler method and if the decay rate is set as µ :=μ − θ ∆t ∆x W , then, the eigenvectors of the source term satisfy the desired inequality Wi ≥ 0. To sum up, the discretization in space and time makes the guaranteed rate smaller. But this is expected, since the stability region of the Euler scheme should be considered. Furthermore, we note that eigenvalues of B i do not coincide with those of B i , only in the sense As in the continuous case our analysis covers the pathological case µ ≤ 0.

Suitable boundary conditions for exponential decay.
So far, we have assumed a given constantμ > 0 such that the matrix H(μ) is negative semidefinite. A strictly positive decay rate exists for systems of linear conservation laws if the boundary conditions satisfy the dissipativity condition where D denotes the set of real m×m diagonal matrices with strictly positive entries [4,Sec. 3]. We state a generalization for balance laws with small source term. We extend this dissipativity condition to non-uniform linearized balance laws with possibly large source terms. The important point of this discussion is that one may use some D ∈ D with DGD −1 2 < 1 to specify the weights of the Lyapunov function by the ghostcell-weights H G := D 2 . To obtain the corresponding result in the discretized case we assume the restriction (A3), i.e. ∆x ∈ (0, λmin /μ) for λ min := min Proof. According to Definition 4.1, there is the relationĤ B (μ) = H G Π 2 (μ). We define D B (μ) := DΠ(μ) to obtain . The inverse Π −1 (μ) exists for ∆x < λmin /μ. Since the · 2 -norm of a matrix B satisfies B 2 2 = σ max B T B , the inverse is estimated by Then, the assumption σμ ,∆x (G, D) ≤ 1 implies that the matrixĤ(μ) is negative semidefinite: In practical applications, the decay rate and boundary conditions are prescribed and weights of a Lyapunov function can be determined by the minimization problem The set D may be restricted for technical reasons. Furthermore, the minimum must neither be unique nor actually exist as long as σμ ,∆x (G, D) ≤ 1 holds. Of course, the minimization problem must be solvable, which is a slightly stronger assumption in the discretized case. In the limit ∆x → 0, the assumption tends to the continuous case. To see this, we assume that D is the infimum, i.e. σ 2 (G) = DGD −1 2 . Then, there is convergence from above σμ ,∆x (·, D) σμ(·, D) σ 2 (·) for ∆x,μ → 0. To sum up, the decay of the Lyapunov function is split into two parts. Firstly, there is a possible increase in the interior of the interval (0, L) which is described in Theorem 3.3 and Theorem 4.2. This increase is smaller for a larger constantμ > 0. If this constant is too large, the dissipativity condition may be violated. This means that both parts are coupled. But we balance the increase in the interior by appropriate boundary conditions to obtain a preferable large decay rate µ.
6. Numerical results for linearized balance laws. First, we investigate if the rate of convergence of the proposed upwind, downwind discretization [24] is inherited to the larger system (S d ) of size 2(d + 1). Furthermore, the convergence of the discretized Lyapunov function itself with respect to spatial refinement is studied. Then, we give a practical application. We state the relationship between the boundary feedback law in terms of the physical variables and Riemann invariants. Finally, we give an example how source terms may cause instabilities. We illustrate the analysis for isothermal Euler and shallow water equations.

Isothermal Euler equations are defined as
where ρ(t, x) is the density of the gas, q(t, x) the mass flux in the pipe, f > 0 a friction factor, D > 0 the diameter of the pipe and a > 0 is the speed of sound. The steady state momentumq is constant and we assumeq ≥ 0. The steady state is determined by the initial value problem withρ 0 :=ρ(0) > 0. Linearization at steady state gives Assuming a subsonic flow with q /ρ < a the Jacobian A is diagonalizable with characteristic speedsλ ± =q/ρ ± a that satisfyλ − < 0 <λ + . According to [16,20], stationary states exist as smooth solutions on a finite space interval only, until a critical length is reached. There, a blow-up in the derivatives occurs. The density at steady state is decreasing, i.e.ρ x < 0 and consequently, eigenvalues are increasing asλ ± x = −qρ −2ρ x > 0. Characteristics do not cross and a stable solution exists on unbounded time domains as long as the length of the pipe is sufficiently small. Shallow water equations describe the conserved quantities height h and momentum q. Since we do not consider shocks, we use the equivalent formulation with respect to the velocity u := q /ρ.
where g > 0 denotes the gravitational constant, C describes friction and S B := −B x the slope of the possibly space-varying bottom topography B(x). To obtain a subcritical flow, we assume a Froude number Fr := |u| / √ gh < 1 such that the Jacobian of the flux function The momentum in steady stateq =hū is again constant. If the length of the pipe is too large, non-constant classical steady states may break down due to a blow up in the derivative [18].
For all simulations we use a unit time and space interval with CFL := 1 and initial perturbations in Riemann coordinates are ζ ± (0, x) := cos(2πx). Lyapunov functions are normalized such that L(0) = 1.
6.1. Spatial discretization. To investigate the spatial discretization, we consider isothermal Euler equations with the periodic boundary conditions G := 1. The upwind-scheme lets us expect a linear convergence, which is indeed observed in Figure 1 and Table 1, where the logarithm of the error is shown together with the empirical order of convergence (EOC). The reference solution R ref is calculated on a refined grid with ∆x = 2 −12 and is summed up on the corresponding coarser grid, which causes no additional error for finite volume methods. Indeed, we observe linear convergence also for the larger system (S d ) of size 2(d + 1). For the Lyapunov function, there is no theoretical linear convergence. Therefore, we show in Table 2 additionally an approximation of the L 2 -and L ∞ -error where the reference value L ref is again calculated on a refined grid with ∆x = 2 −12 and the coarser approximation is linearly interpolated on the refined grid. Integrals are approximated by the trapezoidal rule. We observe the rate of convergence of the upwind-scheme is inherited to corresponding Lyapunov functions. We refer the interested reader to [25], where this is only observed for CFL = 0.75.   Figure 2 shows the Lyapunov functions with respect to the L 2 , H 1 , . . . , H 4 norm. Although the pathological caseμ = 0 does not guarantee a stable system, we observe a very slight decay in the Lyapunov function probably due to friction effects and numerical viscosity. 6.2. Stabilization of gas networks. Since friction makes the pressure decrease along pipes, compressors are used to amplify pressure, as illustrated in Figure 3 for a pipeline with two compressors. Like in [16,15] compressors are modelled by assuming conservation of mass and by applying a compressor power u(t) ≥ 0, i.e. q in (t, L) = q out (t, 0) and u(t) = q out (t, 0) where the superindices (in) and (out) stand in the caseq ≥ 0 for the left-sided and right-sided pipe of the corresponding compressor. The specific heat ratio κ depends on the considered gas, we use κ = 0.5. Additional boundary conditions are needed to close the system. So we assume a given steady state at the left and right end. Compressor powers u 1 , u 2 are determined by the steady state and affect boundary conditions. An interesting question is how to set power to obtain dissipative boundary conditions. An analytical answer for the nonlinear system is     Figure 2. Lyapunov functions for G = 1 andμ = 0 given in [16] for a finite time horizon. In extension to [16], Figure 3 shows for the linearized system a stability domain for compressor power (u 1 , u 2 ) ∈ R + 0 × R + 0 ρμ ,∆x (G, D) ≤ 1 for some D that guarantees a dissipativity condition forμ := 1. We have used the Matlab function fmincon to approximate the minimum min  ρ out (t, 0) = ρ in (t, L) u q out (t, 0) + 1 1 /κ ≈ρ out (0) + c 1 (u),c 2 (u) ∆ρ in (t, L) ∆q out (t, 0) for some coefficientsc i (u) that depend on the energies u 1 and u 2 . Superindices (j) denote the j-th pipe and the j-th compressor. Recall that at the left and right end additional boundary conditions are needed to close the system, in our simulations ∆ρ (1) (t, 0) = 0 and ∆q (3) (t, L) = 0. The transform into Riemann invariants is given by the following theorem.  non-constant bottom slope. More specifically, Table 3 shows the decay rate, guaranteed in Theorem 3.3, with estimated smallest eigenvalue and the observed rate: • guaranteed rate: µ g :=μ + 2 d min (μ) • estimated rate: µ e :=μ − 2 cond √ W (μ) max We see that stability is only guaranteed if the perturbations are small enough, as for µ g < 0 exponential stability is not guaranteed by Theorem 3.3. The case p = 1, where stability is not guaranteed and which is gray shaded in Table 3, is shown in Figure 5. Furthermore, the estimated decay rate seems as a good lower bound and the observed rate is indeed larger than the guaranteed rate, i.e. µ e ≤ µ g ≤ µ o . This is because Theorem 3.3 takes only the smallest eigenvalue of the matrix M(x, µ,μ) into account, whereas positive eigenvalues make Lyapunov functions decrease further. 7. Summary. Exponential Lyapunov stability of linearized hyperbolic balance laws and explicit decay rates have been discussed by theoretical and numerical analysis. We have shown how boundary conditions and the weights of Lyapunov functions can be set to ensure exponential stability. Results are generalized to stability with respect to higher Sobolev norms and the influence of source terms has been investigated. A discretized Lyapunov function with an upwind and downwind discretization that preserves the exponential stability has been presented.