1 Regularity and Lyapunov stabilization of weak entropy solutions to scalar conservation laws Sébastien Blandin, Xavier Litrico, Maria Laura Delle Monache, Benedetto Piccoli and Alexandre Bayen Abstract— We consider the problem of Lyapunov boundary stabilization of the weak entropy solution to a scalar conservation law with strictly convex flux in one dimension of space, around a uniform equilibrium. We show that for a specific class of boundary conditions, the solution to the initial-boundary value problem for an initial condition with bounded variations can be approximated arbitrarily closely in the L1 norm by a piecewise smooth solution with finitely many discontinuities. The constructive method we present designs explicit boundary conditions in this class, which guarantee Lyapunov stability of the weak entropy solution to the initial-boundary value problem. The stabilization result accounts for the proper treatment of boundary conditions in the weak sense. We also design a greedy controller obtained by maximizing the instantaneous decrease rate of the Lyapunov function, and illustrate the limitations of such simple controllers. Finally, we design improved boundary feedback controller which guarantees Lyapunov asymptotic stability while accounting for proper treatment of weak boundary conditions. Controllers performance is illustrated on numerical benchmarks using the Godunov scheme. Index Terms— Lyapunov stabilization, control, conservation laws MSC: 93D05, 35L65 I. I NTRODUCTION A. Motivation The conservation principle is one of the most fundamental modeling principles for physical systems. Statements of conservation of mass, momentum, energy are at the center of modern classical physics. For distributed dynamical systems, this principle can be written in conservation law form with the use of partial differential equations (PDE). The problem of well-posedness of the partial differential equation is concerned with the existence, uniqueness, and continuous dependence of the solution to the problem data [17]. First existence results for scalar conservation laws in one dimension of space date back to [29]. For hyperbolic systems of conservation laws in one dimension of space, existence results were provided in [18] with the introduction of the random choice method. Existence and uniqueness in the scalar case for several dimensions of spaces were proven in [24], and to this date constitute the only general results known on well-posedness for several dimensions of space. Existence and uniqueness for n × n hyperbolic systems of conservation laws in one-dimension of space was shown only Sébastien Blandin is a Research Scientist, IBM Research, Singapore (email: [email protected]). Xavier Litrico is the Director of LyRE, R&D center of Lyonnaise des Eaux, Bordeaux, France ([email protected]). Maria Laura Delle Monache is Postdoctoral Researcher at the Department of Mathematical Sciences, Rutgers University, Camden ([email protected]). Benedetto Piccoli is the Joseph and Loretta Lopez Chair Professor of Mathematics, Department of Mathematical Sciences, and Director, Center for Computational and Integrative Biology, Rutgers University, Camden ([email protected]). Alexandre Bayen is the Chancellor’s Professor, Department of Electrical Engineering and Computer Sciences, Department of Civil and Environmental Engineering, University of California, Berkeley (e-mail: [email protected]). recently, see [11], [12] for systems with genuinely nonlinear or linearly degenerate characteristic families, and [5], [8] for the general case of systems of strictly hyperbolic conservation laws. The global well-posedness of solutions to hyperbolic systems of conservation laws in several dimensions of space is still largely open. The keystone of well-posedness results, and a standard argument for constructive existence proofs, is the consideration of a functional space in which small variations in the problem data, i.e. initial condition for the Cauchy problem, initial and boundary conditions for the initial-boundary value problem (IBVP), create only small variations in the tentative solution. The control problem is posed from a different perspective, and in different terms. In the case of boundary control [23], the control problem consists of an objective trajectory for the system in a given functional space, and the knowledge of an initial condition. The problem of control [13] or stabilization [22], [7], [31], [14] consists in the existence and design of a controller, i.e. boundary conditions in the case of boundary control or stabilization, for which the solution to the partial differential equation stays in a domain prescribed by the objective trajectory. In this article we propose to show Lyapunov stability of the solution to the initial-boundary value problem associated with a scalar conservation law under suitably designed boundary conditions. We consider a scalar conservation law with smooth strictly concave or convex flux in one dimension of space. This partial differential equation is called the Burgers equation [20] in the case of a quadratic convex flux, and in the case of a concave flux, corresponds to the Lighthill-Whitham-Richards (LWR) [28], [32] PDE in its various forms, used in particular for macroscopic traffic flow modeling. The problem is defined in the following section. B. Problem statement Consider the scalar conservation law in one dimension of space ∂t u + ∂x f (u) = 0 (1) . on the domain Ω = { (t, x)| t ≥ 0 and a ≤ x ≤ b}. The flux function f (·) is assumed to be smooth (infinitely differentiable) and strictly convex or strictly concave1 . The initial-boundary value problem (IBVP) for (1) in Ω with initial condition u0 : (a, b) 7→ R, and boundary conditions ua , ub : R+ 7→ R, reads ∂t u + ∂x f (u) = 0 (2) u(0, x) = u0 (x) (3) u(t, a) = ua (t), u(t, b) = ub (t). (4) The Lyapunov boundary stabilization problem can now be formulated. Definition 1: Given a stationary solution u∗ to the PDE (1), the Lyapunov (resp. asymptotic) boundary stabilization problem 1 this is equivalent to the condition of genuine nonlinearity of the characteristic field. 2 consists in the existence of boundary conditions ua , ub depending on initial condition u0 with bounded variations such that the following is true. The IBVP (2)-(3)-(4) is well-posed and its solution is (resp. asymptotically) stable in the sense of Lyapunov at u∗ . Note that Lyapunov (resp. asymptotic) stabilization consists in the existence of a positive definite function decreasing (resp. vanishing) in time along trajectories of the system. Stationary solutions to the PDE include constant solutions and solutions with a single stationary jump discontinuity, called shock. In this article we address the case of constant solutions. The well-posedness of the IBVP (2)-(3)-(4) is critical to the definition of the problem, since the design of arbitrary boundary conditions can make the problem ill-posed (see [36] for an illustration on the LWR equation in the case of traffic). This would lead to a discrepancy between the control implemented and its realized value in the system, and a divergence between the desired trajectory of the system and its real trajectory. In the case of traffic, it corresponds for instance to installing a green traffic light at the location of a traffic jam with the intended goal that cars in the jam adopt the corresponding free-flow speed. For well-posedness of the IBVP (2)-(3)-(4), the PDE (1) and the boundary conditions (4) must be understood in the weak sense. The weak formulation is presented in Section II-A. The remainder of the article is organized as follows. Section II defines the notations used later in the article and states useful lemmas. Section III contains the proof of well-posedness of the initial-boundary value problem with piecewise-smooth increasing datum with negative gradient concentrated at a finite number of locations. Preparatory derivations involving the Lyapunov function candidate can be found in Section IV. In Section V, we show stability of the Lyapunov function for weak entropy solutions to the scalar conservation law remaining in the special class with a finite number of shocks. In Section VI we design a controller that maximizes the instantaneous decrease of the Lyapunov function identified previously, but highlight several configurations in which asymptotic stability is not achieved by this controller. In Section VII we design a new controller which guarantees asymptotic stability, and guarantees the existence of solution to the IBVP in the special class of solutions with a finite number of shocks. Moreover we propose a nonlocal control to improve the convergence time. Numerical examples are proposed in Section VIII, and promising research avenues related to this work in Section IX. II. P RELIMINARIES It is well-known that jump discontinuities can arise in finite time in solutions to conservation laws [17]. Thus classical solutions do not exist in general, and it is necessary to consider a more general formulation of the conservation law. 1) Weak entropy solution to the Cauchy problem: The weak formulation of the conservation law is obtained by considering derivatives in the sense of distribution. Definition 2: A function u : [0, T ] × R 7→ R is a weak solution to the Cauchy problem (2)-(3) if for any continuously differentiable function φ with compact support contained in (−∞, T ) × R, Z TZ ∞ Z ∞ (u φt + f (u) φx ) dx dt + u0 (x) φ(0, x)dx = 0, −∞ −∞ and t 7→ u(t, ·) is continuous from [0, T ] into L1loc . σ ∆u = ∆f (u), (5) (6) where ∆u = ur − ul is the jump in u, with ur , respectively ul , the value of the right, respectively left, limit of u at the jump location. To isolate a unique weak solution to a Cauchy problem associated with the conservation law, an additional admissibility condition (see Section 4.5 of [15]) is required. Different conditions have been proposed in the literature. In the scalar case, one of the first admissibility conditions, due to Oleinik [29], states that for a shock joining a left state ul and a right state ur , the following inequality must be satisfied for all u between ul and ur : f (u) − f (ul ) f (u) − f (ur ) ≥σ≥ , u − ul u − ur (7) where σ is the Rankine-Hugoniot speed (6). Kruzkhov [24] showed that it was sufficient to satisfy the entropy inequality condition for a specific family of entropy-entropy flux pairs in the scalar case, yielding the Kruzkhov entropy condition. The Lax admissibility condition [25], which exhibits a convenient geometric interpretation, states that for a shock joining a left state ul and a right state ur , the following inequality must be satisfied: λ(ul ) ≥ σ ≥ λ(ur ) (8) where λ(u) is the characteristic speed of u (i.e. f 0 (u)), and σ is the Rankine-Hugoniot speed (6). For the case of systems one requires that condition (8) holds for a genuinely nonlinear or linearly degenerate i-th characteristic family with λ replaced by the i-th eigenvalue λi of the Jacobian matrix DF (u). In the scalar case, for a convex flux, these formulations have been proven to be equivalent (see Section 2.1 of [26]). The Lax admissibility condition allows the selection of a particular weak solution. Definition 3: A function u : [0, T ]×R 7→ R is the weak entropy solution to the Cauchy problem (2)-(3) if it is a weak solution (definition 2), that satisfies the Lax admissibility condition (8). The definition of weak conditions to the IBVP requires a corresponding statement of weak boundary conditions, presented in the following section. 2) Weak boundary conditions: The first statement of weak boundary conditions was introduced in [6] in the scalar case in multiple dimensions of space, with C 2 flux and C 2 initial and boundary datum, using a vanishing viscosity method. In one dimension, this formulation reads: max sgn (u(t, a) − ua (t)) (f (u(t, a)) − f (k)) = 0 (9) min sgn (u(t, b) − ub (t)) (f (u(t, b)) − f (k)) = 0 (10) k∈[α,β] A. Weak solutions to the initial-boundary value problem 0 Given that u is smooth around a jump discontinuity, integrating the weak formulation (5) yields the Rankine-Hugoniot relation [17] defining the speed σ of propagation of jump discontinuities k∈[γ,δ] for almost all t > 0, and where α = min(u(t, a), ua (t)), β = max(u(t, a), ua (t)), γ = min(u(t, b), ub (t)), δ = max(u(t, b), ub (t)), and sgn denotes the sign function. For the case of systems of conservation laws, the interested reader is referred to [9],[34]. In the scalar case, at a left boundary a, the corresponding statement of weak boundary conditions derived from the structure of the solution to the Riemann problem is the following. Definition 4: A function u : Ω 7→ R satisfies the boundary condition ua at a if for almost every time t (in the sense of Lebesgue measure), the solution to the Riemann problem centered at a with initial data ( ua (t) if x < a (11) u(t, a) if x > a 3 either does not contain any wave (when left and right initial states are the same), or contains waves with non positive speeds (a wave with zero speed is allowed). Notice that for weak solutions we consider the condition will hold for all times except a finite number. B. Well-posedness of the initial-boundary value problem In [6], a solution satisfying (5) in the scalar case is constructed using a vanishing viscosity method for the weak boundary conditions statement (9)-(10) and is shown to be the admissible solution according to Kruzkhov entropy condition [24]. More recently, an existence result for n × n systems using wavefront tracking was proposed in [1]. The standard Riemann semigroup (SRS) method, introduced in [10] for the Cauchy problem associated with a Temple system [37] of conservation laws, was extended to the IBVP in [2], [3], with the boundary conditions statement from definition 4. In [3], it is shown for a n × n system that if the SRS exist, its trajectories coincide with wavefront tracking solutions. Uniqueness and continuous dependence is obtained for the case of non-characteristic conditions, and uniqueness for the characteristic case. The SRS is constructed for 2 × 2 system in [2], and for the case of n × n system in [16]. The stability of the IBVP with two boundaries was established via vanishing viscosity for 2 × 2 systems in [35]. We state in the scalar case for a static boundary the main result of [2] for characteristic boundary conditions, obtained for 2 × 2 systems with continuous boundary (see theorem C of [2]). A general result for n × n systems was established in [4]. Theorem 1: [2] Let f be a smooth map such that the equation (1) is strictly hyperbolic with characteristic field linearly degenerate or genuinely nonlinear (i.e. f is linear, convex or concave). For every δ > 0 there exists L > 0 and a continuous semigroup S defined for data in L1 ∩ BV with total variation bounded by δ, such that • • • The map t 7→ u(t, ·) yields a weak solution to the IBVP (2)(3)-(4). For piecewise constant initial and boundary data, the trajectories of the semigroup coincide with the solution to the IBVP obtained by piecing together the standard solutions to the Riemann problems at the points of discontinuity of the initial condition and at the boundary. For initial data u00 , u000 , boundary data u0a , u00a , u0b , u00b in L1 ∩ BV with total variation bounded by δ, let u0 , u00 denote the corresponding trajectories of the semigroup S, and t0 , t00 > 0, then ku0 (t0 , ·)−u00 (t00 , ·)k1 ≤ L(|t0 −t00 |+ku00 −u000 k1 +ku0a −u00a k1 + ku0b − u00b k1 ). We also refer the interested reader to the work of Otto [30] for the case where the entropy solution does not have traces at the boundary, see also [33]. In the following section we use this result in the case of a left and a right boundary to show that we can restrict our Lyapunov analysis to the case of piecewise smooth data. III. A PPROXIMATION OF SOLUTION TO INITIAL - BOUNDARY VALUE PROBLEM BY PIECEWISE SMOOTH SOLUTION In this section, we present results on the approximability of the solution to an IBVP with initial condition in BV by the solution to an IBVP with piecewise smooth solution at all times. We show that the solution to the IBVP with BV data can be approximated arbitrarily closely in the L1 norm by the solution to an IBVP with piecewise smooth data. We define the required regularity class used throughout the article. Definition 5: We note PWS+ the class of piecewise smooth functions f : R 7→ R such that 0 • f is positive (where defined). 0 • f has only downward jumps (i.e. f seen as measure has only negative Dirac masses). We now state the approximability result. Theorem 2: Consider T > 0, a < b, and let u0 : (a, b) 7→ R, ua , ub : (0, T ) 7→ R be functions with bounded total variation. For every > 0, there exists u0 : (a, b) 7→ R in PWS+ and piecewise constant boundary data ua , ub : (0, T ) 7→ R such that kua (t) − ua (t)k1 ≤ , kub (t) − ub (t)k1 ≤ and the solution u to the IBVP for equation (1) and data (u0 , ua , ub ) and the solution u to the IBVP for equation (1) and data (u0 , ua , ub ) satisfy: ∀ 0 ≤ t ≤ T, ku(t, ·) − u (t, ·)k1 ≤ . Proof: In the compact domain [a, b], we can approximate the initial condition condition u0 arbitrarily closely in the L1 sense by a piecewise constant function u0 with a finite number of discontinuities and lower total variation. Also we can replace the upward jumps of u0 with smooth increasing functions without increasing the total variation. We can also approximate boundary data by piecewise constant functions of lower variation. Since u0 has only positive derivative, no new shock can form in the solution inside the domain. Moreover, a finite number of shocks will be introduced by the boundary conditions. Therefore the solution u will remain in the class PWS+ . Using the continuous dependence result of Theorem 1, the resulting trajectories u, u of the semigroup can be made arbitrarily close in the L1 norm by controlling the distance between the initial conditions. Under suitable boundary conditions, the solution to the IBVP with piecewise smooth data is piecewise smooth: Theorem 3: Let T, δ > 0, a < b, and let u0 : (a, b) 7→ R be in PWS+ and ua , ub : (0, T ) 7→ R be piecewise constant. Let u denote the solution to the IBVP (5)-(3)-(11). At all times 0 ≤ t ≤ T , u(t, ·) is piecewise smooth. In the next section, we present the Lyapunov stability analysis for functions in PWS+ . IV. LYAPUNOV ANALYSIS In this section, we propose a Lyapunov function and compute its derivative. In the following we call ũ = u − u∗ where u∗ is a constant, hence stationary, solution around which we want to stabilize the system, and u is the solution to the IBVP associated with the scalar conservation law (1). Following the results from section III, we assume that u is in PWS+ . A. Lyapunov function candidate We consider the classical Lyapunov function candidate [21], [23]: Z Z 1 b 2 1 b V (t) = ũ (x, t) dx = (u(x, t) − u∗ )2 dx, (12) 2 a 2 a where u is a weak solution to the scalar conservation law. From definition 3, we have t 7→ u(t, ·) continuous from [0, T ] to L1 , and the function V (·) is well-defined and continuous. We index the jump discontinuities of u(t, ·) in increasing order of their location at time t by i = 0, . . . , N (t), including for notational purposes the boundaries a, b, with x0 (t) = a and xN (t) = b. The Lyapunov function candidate can be rewritten as: N (t)−1 Z x i+1 (t) 1 X V (t) = ũ2 (t, x) dx. (13) 2 i=0 xi (t) 4 From Theorem 3, we know that for all integer i ∈ [0, N (t)), the function u(t, ·) is smooth in the domain (xi (t), xi+1 (t)), thus ∂t u(t, ·) exists and is continuous for t such that xi (t) < xi+1 (t). Since discontinuity trajectories are differentiable with speed given by the Rankine-Hugoniot relation (6), it follows that at any time t such that N (t) is constant in a neighborhood of t and the boundary trace is continuous, the function V (·) is differentiable. B. Differentiation of the Lyapunov function candidate In this section, we compute the derivative of the Lyapunov function candidate (12), at any time t such that N (t) is constant in a neighborhood and the boundary trace is continuous. Differentiating expression (13) yields: In equation (15) we gather the first four terms that depend on the boundary trace of the solution, and the last two terms that depend on the shock dynamics inside the domain. In the following section, we analyze the stability properties of the internal terms. C. Internal stability The last two terms of equation (15) correspond to jump discontinuity in the solution and are neither observable nor controllable from the boundaries. We now show that these terms have a stabilizing effect on the Lyapunov function candidate (12). Proposition 1: Given a constant solution u∗ to the scalar conservation law (1), we have the following inequality N (t)−1 Z x " # N (t)−1 i+1 (t) X ũ(t, xi −) + ũ(t, xi +) dV 1 X ∗ ∗ ∗ ∆i ũf (ũ + u ) − F (ũ + u ) − ∆i f (ũ + u ) ≤ 0 (t) = ∂t ũ2 dx 2 i=1 dt 2 i=0 xi (t) (16) N (t)−1 1 X dxi dxi+1 i.e. the jump discontinuity dynamics of the solution u to the + (t) − ũ2 (t, xi (t)+) (t) . ũ2 (t, xi+1 (t)−) 2 i=0 dt dt IBVP, contributes to the decrease of the Lyapunov function candi(14) date (12). Proof: In order to show that the term (16) is negative, we As detailed at the end of section IV-A, the term under the sum show that each term in the sum is negative. If we note ul , ur the is smooth, and we can write ∂t ũ2 = 2 ũ ∂t ũ. Since u satisfies the value of u on the left and on the right of the jump discontinuity, conservation law (1), we have ∂t ũ = −∂x f (ũ+u∗ ). The derivative respectively, and ũl , ũr the corresponding reduced variables, we of the Lyapunov function can be written as: want to show that N (t)−1 Z x X i+1 (t) dV (t) = − ũ ∂x f (ũ + u∗ ) dx dt x (t) i i=0 1 + 2 (ũr f (ũr + u∗ ) − F (ũr + u∗ )) N (t)−1 dxi+1 dxi ũ2 (t, xi+1 (t)−) (t) − ũ2 (t, xi (t)+) (t) . dt dt X i=0 Equivalently, in the original state variable u = ũ + u∗ , we have: By integrating by parts the sum terms, we obtain: N (t)−1 X dV x (t) (t) = − [ũ f (ũ + u∗ )]xii+1 (t) dt i=0 N (t)−1 Z x X i+1 (t) + f (ũ + u∗ ) ∂x ũ dx i=0 xi (t) N (t)−1 + 1 X 2 i=0 ũ2 (t, xi+1 (t)−) dxi+1 dxi (t) − ũ2 (t, xi (t)+) (t) . dt dt and if we note F (·) a primitive function of the flux function f (·) we have: dV (t) = ũ(t, a) f (ũ(t, a) + u∗ ) − ũ(t, b) f (ũ(t, b) + u∗ ) dt − F (ũ(t, a) + u∗ ) + F (ũ(t, b) + u∗ ) N (t)−1 X 1 dxi + ∆i (ũf (ũ + u∗ ) − F (ũ + u∗ )) − (t) ∆i ũ2 , 2 dt i=1 where ∆i is defined around the discontinuity xi (t) as in equation (6). Using the Rankine-Hugoniot relation, defined in equation (6), to write the speed of the jump discontinuity dxi (t)/dt as a function of the left and right jump values we obtain: dV (t) = ũ(t, a) f (ũ(t, a) + u∗ ) − ũ(t, b) f (ũ(t, b) + u∗ ) dt − F (ũ(t, a) + u∗ ) + F (ũ(t, b) + u∗ ) N (t)−1 + X ∆i (ũf (ũ + u∗ ) − F (ũ + u∗ )) i=1 N (t)−1 X ũ(t, xi −) + ũ(t, xi +) − ∆i f (ũ + u∗ ). 2 i=1 − (ũl f (ũl + u∗ ) − F (ũl + u∗ )) ũl + ũr − (f (ũr + u∗ ) − f (ũl + u∗ )) ≤ 0. 2 (15) [((ur − u∗ ) f (ur ) − F (ur )) − ((ul − u∗ ) f (ul ) − F (ul ))] ul + ur − 2 u∗ − (f (ur ) − f (ul )) ≤ 0. 2 This can be rewritten as: 1 (17) F (ul ) − F (ur ) + (ur − ul ) (f (ur ) + f (ul )) ≤ 0, 2 which can be obtained from the Oleinik condition(7), here by integration the left inequality of (7) between ul and ur . Thus any solution satisfying the Oleinik entropy condition benefits from stability of the jump discontinuity dynamics. V. W ELL - POSED BOUNDARY STABILITY In this section, we motivate and define the control space and propose a stabilizing controller. A. Control space Definition 6: Let us denote smin , smax the minimal and maximal speed of the waves composing the solution to the Riemann problem at the boundary. The control space at the left boundary is the set of pairs (ul , ur ) such that either no wave is generated by the Riemann problem (ul = ur ), or smin ≥ 0 and smax > 0. The control space at the right boundary is the set of pairs (ul , ur ) such that either no wave is generated by the Riemann problem (ul = ur ), or smin < 0 and smax ≤ 0. Proposition 2: Let m denote the minimizer of the strictly convex flux function f . The control spaces Ca , Cb at the left and right boundaries, respectively, are characterized as the set of pairs (ul , ur ) such that one of the following properties is satisfied 5 summarized in the following formula: . Ca = (ul , ur ) ul = ur ul ≥ m and u ≥ m and l s.t. ul = ur ul ≤ m and u ≥ m and l . Cb = (ul , ur ) s.t. ur ≥ m ur ≤ m and ur ≤ m ur ≤ m and f (ul ) > f (ur ), (18) f (ul ) < f (ur ). (19) Using the characterization of the control space introduced in this section, we show in the following section that the system is stabilizable. B. Lyapunov stabilization In this section, we prove that there exist boundary conditions in the control space (18)-(19) such that the candidate Lyapunov function (12) is strictly decreasing. Lemma 1: Let g : u 7→ (u − u∗ ) f (u) − F (u), with f a smooth strictly convex function, and F a primitive of f . Let m denote the minimum of f . The function g is smooth on the real line, and satisfies the following properties: g is strictly increasing in (−∞, min(m, u∗ )), strictly decreasing in (min(m, u∗ ), max(m, u∗ )), and strictly increasing in (max(m, u∗ ), +∞). • For u > v such that f (u) = f (v), we have g(u) > g(v). Proof: The fact that g is smooth results from the smoothness of f . The first property is obtained by computing the derivative g 0 (u) = (u − u∗ ) f 0 (u) of g, and noting that f is strictly convex with minimum at m. To prove the second property, let us consider u > v such that f (u) = f (v). The difference g(v) − g(u) reads g(v) − g(u) = F (u) − F (v) + (v − u) f (u) that is strictly negative by strict convexity of f . The function g is represented for the case of the Burgers flux function in figure 1 with the arbitrary choice of g(m) = 0. g(u) g(u) • 0 −1 u* 0 u 1 0 −1 0 u u* 1 Fig. 1. Representation of the variations of g: for a Burgers flux function f in the case u∗ < 0 (left) and in the case u∗ > 0 (right). The points u = m (m = 0 in this case) and u = u∗ (u∗ = ±0.5 in this case) are local extrema of g. Theorem 4: Let V (·) denote the candidate Lyapunov function (12) for the PDE (1). There exist boundary conditions ua (·), ub (·) in the control spaces (18) (19), respectively, such that the following holds. If the corresponding solution to the IBVP is in the class PWS+ then the function V (·) is strictly decreasing, thus the solution is stable in the sense of Lyapunov. Proof: The set of boundary values (ua (t), ub (t)) which guarantee stabilization can be computed using Lemma 1 and are [m, u(t, b)) × {u(t, b)} (m, +∞) × (−∞, m) s.t. g(ua (t)) < g(ub (t)) {u(t, a)} × {u∗ } {u(t, a)} × {u(t, b)} {u(t, a)} × {u∗ } if (u(t, a), u(t, b)) ∈ [m, +∞) × (m, +∞) if (u(t, a), u(t, b)) ∈ [m, +∞) × (−∞, m] if (u(t, a), u(t, b)) ∈ (−∞, m) × (−∞, m] if (u(t, a), u(t, b)) ∈ (−∞, m) × (m, +∞) and g(u(t, a)) < g(u(t, b)) if (u(t, a), u(t, b)) ∈ (−∞, m) × (m, +∞) and g(u(t, a)) ≥ g(u(t, b)) (20) We highlight that the boundary values (20) only provide stability, and not asymptotic stability in general. In the following section we instantiate a greedy controller which maximizes the instantaneous decrease of the Lyapunov function, and also illustrate that the greedy controller may fail to provide asymptotic stability. In Section VII we then design an improved controller which guarantees asymptotic stability and we show that the solution resulting from these boundary controls and an initial condition in the class PWS+ remains in PWS+ . VI. M AXIMIZING INSTANTANEOUS LYAPUNOV FUNCTION DECREASE RATE In this section, we characterize the values of the control, in the control space, that minimize the Lyapunov function derivative. A greedy boundary control (ua (t), ub (t)) maximizing the instantaneous decrease of the Lyapunov function, read as follows (for the case u∗ < m): {m} × {u(t, b)} {m} × {u∗ } ∗ {u(t, a)} × {u } if (u(t, a), u(t, b)) ∈ [m, +∞) × (m, +∞) if (u(t, a), u(t, b)) ∈ [m, +∞) × (−∞, m] if (u(t, a), u(t, b)) ∈ (−∞, m) × (−∞, m] {u(t, a)} × {u(t, b)} if (u(t, a), u(t, b)) ∈ (−∞, m) × (m, +∞) and g(u(t, b)) > g(u∗ ) {u(t, a)} × {u∗ } if (u(t, a), u(t, b)) ∈ (−∞, m) × (m, +∞) and g(u(t, b)) ≤ g(u∗ ) (21) While the greedy controller (21) maximizes the instantaneous decrease of the Lyapunov function, we illustrate in the following example that asymptotic stability may not be obtained. We also illustrate the naive control (ua (t) = u∗ , ub (t) = u∗ ) may not provide asymptotic stability, and several key factors of the mechanisms involved for both controllers. Example 1: We first consider an example where the greedy controller fails to provide asymptotic stability, and where the brute force controller creates oscillations. Without loss of generality, we choose u∗ < m. Given 0 < ∆ < (a+b)/2 and k > 0, we consider the following initial datum on (a, b): u0 (x) = m ū −k ū + k ) if x ∈ (a, a+b 2 if x ∈ ( a+b + (2 p) ∆, 2 p ∈ {0, · · · , b−a − 1} 4∆ if x ∈ ( a+b + (2 p+ 2 p ∈ {0, · · · , b−a − 1} 4∆ a+b 2 1) ∆, + (2 p + 1) ∆) a+b 2 + (2 p + 2) ∆) where u∗ < m < ū and f (ū) = f (u∗ ). This case corresponds to the first row of equation (21), hence the applied boundary controls 6 is (ua (t) = m, ub (t) = u(t, b)). Since the characteristic speed of m is zero, the right boundary value converges towards m but never reaches it, hence the system remains in the configuration characterized by the first row of equation (21), and converges to the steady state m over the interval (a, b), not reaching the target u∗ < m. Hence we have stability but not asymptotic stability. For the same example, one may note that the brute force control (ua (t) = u∗ , ub (t) = u∗ ) has no action on the system when u(t, b) = ū + k since the control values are outside of the control space. When u(t, b) = ū − k, the brute force control induces slow backward moving shock waves (ū − k, u∗ ) from the right boundary which interact with fast forward moving shock waves (ū + k, ū − k) coming from the initial datum, and create slow forward moving shock waves (ū + k, u∗ ), hence we observe large oscillations at the right boundary (irrespective of the size of the oscillations in the initial datum). Eventually all the waves generating by oscillating initial datum exit the domain and the naive control produces backward moving shock waves (m, u∗ ) which yield convergence. VII. LYAPUNOV ASYMPTOTIC STABILITY In this section we design an improved controller which guarantees asymptotic stability and show that the associated solution to the corresponding IBVP remains in the class PWS+ for initial data u0 in the same class. A. Controller design The asymptotic convergence issue highlighted in the Example 1 stems from the fact that if the system reaches a configuration corresponding to the first row of equation (21), given the prescribed control values, it may remain in that configuration, which grants stability but prevents asymptotic stability for u∗ < m. Without loss of generality, we focus on the case u∗ < m. We define û > m by f (û) = f (u∗ ), ǔ > m by g(ǔ) = g(u∗) and ū = û+ǔ . We 2 propose the following values for the controller (ua (t), ub (t)): ua (t) = min{u(t, a), m}, ub (t) = ψ(u(t, b)), (22) where ψ(α) = χ]−∞,ū] (α) u∗ + (1 − χ]−∞,ū] (α))α and χI is the indicator function of the interval I. It is convenient to describe the control according to different cases as for the greedy control, thus we get: {m} × {u(t, b)} ∗ {m} × {u } if (u(t, a), u(t, b)) ∈ [m, +∞) × (ū, +∞) {u(t, a)} × {u∗ } {u(t, a)} × {u(t, b)} if (u(t, a), u(t, b)) ∈ (−∞, m) × (−∞, ū] if (u(t, a), u(t, b)) ∈ [m, +∞) × (−∞, ū] if (u(t, a), u(t, b)) ∈ (−∞, m) × (ū, +∞) (23) B. Existence of solution for asymptotically stabilizing controller To apply Theorem 4, we need to define boundary controls ua , ub , which will guarantee the solution to the corresponding IBVP to remain in the class PWS+ for initial data u0 in the same class. We prove such result for the boundary control providing asymptotic stability of the Lyapunov function (12). For simplicity we assume u∗ < m being the other case similar. Theorem 5: Consider an initial datum u0 in PWS+ and the boundary controls given by formula (21). Then the corresponding IBVP admits a unique solution which is in the class PWS+ for all times. Proof: We will show that for u0 in PWS+ , the boundary controls are well defined and piecewise smooth in time. Moreover the solution remains in the class PWS+ . First we enumerate the following cases: • Case 1. u(t, a) ≥ m, u(t, b) > ū. • Case 2. u(t, a) ≥ m, u(t, b) ≤ ū. • Case 3. u(t, a) < m, u(t, b) ≤ ū. • Case 4a. u(t, a) < m, u(t, b) > ū. In Case 1. the boundary controls generate a rarefaction wave from the left boundary, no wave on the right boundary and the traces verify u(t+, a) = m and u(t+, b) = u(t, b) > m (where u(t+, ·) indicates the limit at time t from the above). In Case 2. the boundary controls generate a rarefaction wave from the left boundary, a rarefaction wave or a shock on the right boundary and the traces verify u(t+, a) = m and u(t+, b) = u∗ . In Case 3. the boundary controls generate no wave from the left boundary, a rarefaction wave or a shock on the right boundary and the traces verify u(t+, a) = u(t, a) < m and u(t+, b) = u∗ . In Case 4. the boundary controls generate no wave from the left and right boundary and the traces verify u(t+, a) = u(t, a) < m and u(t+, b) = u(t, b) > ū > m. From this analysis we verify that the control is well defined for constant traces and the corresponding solution remains in the class PWS+ . It remains to verify that the solution is well defined and remains in the same class for waves interacting with the boundary and we proceed again by cases. In case 1: if wave interact with the left boundary, then no wave is produced, while for the right boundary a shock may arise and we transition to case 2. In case 2: if wave interact with the left or right boundary, then no wave is produced and we transition to case 1 or 3. In case 3: if wave interact with the left or right boundary, then no wave is produced and we remain in case 3 or transition to 4. In case 4: if wave interact with the left boundary, then no wave is produced, while for the right boundary a shock may arise and we transition to case 3. C. Non-local controls As we noticed the greedy control may not stabilize the system to u∗ , while the brute force control ua ≡ ub ≡ u∗ may overshoot (and not be in the control space) and produce oscillations. Finally, control (23) stabilizes the system, but the stabilization time can be far from optimal (and the same is for the brute force ones). Therefore, in this Section, we show a nonlocal control unl a,b which fast stabilizes the system to u∗ . We use the term nonlocal to indicate that this control will depend not only on the values of the traces u(t, a) and u(t, b) but on the value of the whole solution u. We focus again, for simplicity, on the case u∗ < m. Let A = supx∈[a,b] u0 (x) and  < m be such that f (Â) = f (A). For every U <  we define: T1 (U ) = (b − a) (A − U ) , f (U ) − f (A) and set unl a as in (23), while: U unl b (t) = u∗ T0 (U ) = T1 − (b − a) |f 0 (U )| 0 ≤ t ≤ T0 . T0 < t < +∞ (24) (25) The meaning of this control is as follows. First we send a large shock (u0 (a), U ) with negative speed to move the system in the zone u < m and then apply the stabilizing control. Notice that T1 is computed as the maximal time taken by the big shock to cross the interval [a, b], while T0 is the time at which the characteristic corresponding to u∗ should start from b to reach a 7 at time T1 . These choices will guarantee the desired effect. Notice also that T0 is a safe choice, but smaller values may give a better performance. In the following section we present numerical results of the implementation of the boundary control proposed. VIII. N UMERICAL EXAMPLES Time = 0.5 2 u0 (x) = 1 + 0.5 sin(10 ∗ x). 1.5 1 0.5 Density 1 0.5 0 -0.5 -1 -1 -1.5 -1.5 -2 -2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 Space Time = 0.5 2 1 0.5 Density Density 1 0.5 0 0.9 1 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1 0 -1 -1 -1.5 -1.5 -2 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 Space Time = 0.5 2 0.5 Space Time = 2 2 1.5 1.5 1 1 0.5 0.5 Density Density 0.8 -0.5 0 0 -0.5 0 -0.5 -1 -1 -1.5 -1.5 -2 -2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 Space Time = 0.5 2 0.5 Space Time = 2 2 1.5 1.5 1 1 0.5 0.5 Density Density 0.7 Time = 2 2 -2 0 -0.5 0 -0.5 -1 -1 -1.5 -1.5 -2 -2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 Space Time = 0.5 2 0.5 Space Time = 2 2 1.5 1.5 1 1 0.5 0.5 Density Density 0.6 1.5 -0.5 IX. C ONCLUSION 0.5 Space 1.5 In this article, we introduced a new technique for Lyapunov boundary stabilization results of weak entropy solutions to scalar conservation laws. We proved that under suitable regularity of the initial and boundary data, the solution to the initial-boundary value problem could be considered to be piecewise regular with a finite number of discontinuities. This allows the use of functional analysis tools available for smooth functions. We computed the derivative of the Lyapunov function candidate, and showed that the solution to the scalar conservation law with strictly convex flux is asymptotically stabilizable in the sense of Lyapunov, using the definition of the Oleinik entropy condition. We then showed how a greedy control, maximizing instantaneous decrease of the Lyapunov function, may not asymptotically stabilize the system. Therefore we designed a stabilizing control and a nonlocal ones which improve the stabilization time. 0 -0.5 (26) In figure 2 we present the evolution of the system under the greedy boundary control, the brute force control, the stabilizing control and the nonlocal ones with U = −2 and two different choices of the switching time T0 , see (25). We notice that the greedy control, as shown in Example 1, allows all oscillations to exit from right boundary but the solution does not converge to the steady state u ≡ u∗ = −1 (see Figure 2, first row). On the other side the brute force control does guarantee convergence to the steady state, but allows strong oscillations on the right boundary (see Figure 2, second row, first figure) and at time t = 2 the solution presents a large shock connecting a positive state to the steady state u∗ but still at space position x = 0.8 (see Figure 2, second row, second figure). The stabilizing control also stabilizes and, moreover, avoids oscillations (see Figure 2, third row, first figure). However, at time t = 2 the solution presents a large shock connecting a positive state to the steady state u∗ but still at space position x = 0.9 (see Figure 2, third row, second figure). The nonlocal control guarantees convergence and avoids oscillations. For the choice of switching time T0 = 1.5 (following (24)), at time t = 2 the solution presents a rarefaction wave connecting u = −2 to u∗ and traveling backward (see Figure 2, fourth row, second figure). For the choice of switching time T0 = 1, at time t = 2 the solution is already almost at the steady state, with a little rarefaction still exiting the left boundary (see Figure 2, fifth row, second figure). The downside of the nonlocal control is the fact that a large shock is introduced in the solution, which may be unfeasible for some specific application. Time = 2 2 1.5 Density In this section, we present numerical results obtained for a benchmark scenario. The numerical scheme used is the standard Godunov scheme [19] with 100 cells in space and a time discretization satisfying the tight Courant-Friedrich-Levy (CFL) condition [27]. We consider the flux function u 7→ u2 /2, the equilibrium state u∗ = −1, and the space domain [0, 1] with the oscillating initial condition: 0 -0.5 0 -0.5 -1 -1 -1.5 -1.5 -2 -2 0 0.1 0.2 0.3 0.4 0.5 Space 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 Space Fig. 2. Performance of the various controllers for an oscillating initial data. Plots report solutions at time t=0.5 and at time t=2. First row: greedy control. Second row: brute force control. Third row: stabilizing control. Fourth row: nonlocal control with T0 = 1.5. Fifth row: nonlocal control with T0 = 1. 8 ACKNOWLEDGMENTS The authors would like to thank Professor Miroslav Krstic for suggesting to work on this problem during his visit at UC Berkeley and for fruitful conversation leading to the genesis of this work. At this time the second author was a visiting scholar at UC Berkeley with financial support from CEMAGREF and from the France-Berkeley fund project. The authors would also like to thank their colleague Professor Daniel Work from UIUC for numerous fruitful discussions on a related topic. R EFERENCES [1] D. A MADORI. Initial-boundary value problems for nonlinear systems of conservation laws. NoDEA: Nonlinear Differential Equations and Applications, 4(1):1–42, 1997. [2] D. A MADORI and R. C OLOMBO. Continuous dependence for 2 × 2 conservation laws with boundary. Journal of Differential Equations, 138(2):229–266, 1997. [3] D. A MADORI and R. C OLOMBO. Viscosity solutions and standard Riemann semigroup for conservation laws with boundary. Rendiconti del Seminario Matematico della Universita di Padova, 99:219–245, 1998. [4] F. A NCONA and A. M ARSON. Asymptotic stabilization of systems of conservation laws by controls acting at a single boundary point. In Control methods in PDE-dynamical systems, volume 426 of Contemp. Math., pages 1–43. Amer. Math. Soc., Providence, RI, 2007. [5] F. A NCONA and A. M ARSON. Existence theory by front tracking for general nonlinear hyperbolic systems. Archive for rational mechanics and analysis, 185(2):287–340, 2007. [6] C. BARDOS, A.-Y. L EROUX, and J.-C. N EDELEC. First order quasilinear equations with boundary conditions. Communications in partial differential equations, 4(9):1017–1034, 1979. [7] G. BASTIN and J.-M. C ORON. Stability and boundary stabilization of 1-d hyberbolic systems. to appear, 2016. [8] S. B IANCHINI and A. B RESSAN. Vanishing viscosity solutions of nonlinear hyperbolic systems. Annals of Mathematics, pages 223–342, 2005. [9] S. B IANCHINI and L. S PINOLO. The boundary riemann solver coming from the real vanishing viscosity approximation. Arch. Rational Mech. Anal, page 2009. [10] A. B RESSAN. The unique limit of the Glimm scheme. Archive for Rational Mechanics and Analysis, 130(3):205–230, 1995. [11] A. B RESSAN, G. C RASTA, and B. P ICCOLI. Well-posedness of the Cauchy problem for n × n systems of conservation laws. American Mathematical Society, Providence, RI, 2000. [12] A. B RESSAN, T.-P. L IU, and T. YANG. L1 stability estimates for n×n conservation laws. Arch. Ration. Mech. Anal., 149(1):1–22, 1999. [13] P. C HRISTOFIDES. Nonlinear and robust control of PDE systems: methods and applications to transport-reaction processes. Birkhäuser, Boston, 2001. [14] J.-M. C ORON, S. E RVEDOZA, S. G HOSHAL, O. G LASS, and V. P ER ROLLAZ . Dissipative boundary conditions for 2 × 2 hyperbolic systems of conservation laws for entropy solutions in bv. preprint, http://arxiv.org/pdf/1512.05151.pdf. [15] C. DAFERMOS. Hyperbolic conservation laws in continuum physics. Springer-Verlag Berlin, Germany, 2010. [16] C. D ONADELLO and A. M ARSON. Stability of front tracking solutions to the initial and boundary value problem for systems of conservation laws. Nonlinear Differential Equations and Applications NoDEA, 14(5-6):569–592, 2007. [17] L. E VANS. Partial differential equations. American Mathematical Society, Providence, RI, 1998. [18] J. G LIMM. Solutions in the large for nonlinear hyperbolic systems of equations. Communication on pure and applied mathematics, 18:697– 715, 1965. [19] S. G ODUNOV. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. Matematicheskii Sbornik, 89(3):271–306, 1959. [20] E. H OPF. The partial differential equation ut + u ux = µ uxx . Communications on Pure and Applied Mathematics, 3(3):201–230, 1950. [21] M. K RSTIC. On global stabilization of Burgers equation by boundary control. Systems and Control Letters, 37(3):123–141, 1999. [22] M. K RSTIC and H. D ENG. Stabilization of nonlinear uncertain systems. Springer New York, 1998. [23] M. K RSTIC and A. S MYSHLYAEV. Boundary control of PDEs: A Course on Backstepping Designs. Society for Industrial and Applied Mathematics Philadelphia, PA, 2008. [24] S. K RUZHKOV. First order quasilinear equations in several independent variables. Mathematics of the USSR-Sbornik, 10(2):217–243, 1970. [25] P. L AX. Hyperbolic systems of conservation laws II. Communications on Pure and Applied Mathematics, 10:537–566, 1957. [26] P. L E F LOCH. Hyperbolic Systems of Conservation Laws: The theory of classical and nonclassical shock waves. Lectures in Mathematics, ETH Zürich, Birkhäuser Basel, 2002. [27] R. L EVEQUE. Finite volume methods for hyperbolic problems. Cambridge University Press, Cambridge, UK, 2002. [28] M. L IGHTHILL and G. W HITHAM. On kinematic waves, II: a theory of traffic flow on long crowded roads. Proceedings of the Royal Society of London, 229(1178):317–345, 1956. [29] O. O LEINIK. Discontinuous solutions of non-linear differential equations. Uspekhi Matematicheskikh Nauk, 12(3):3–73, 1957. [30] F. OTTO. Initial-boundary value problem for a scalar conservation law. Comptes rendus de l’Académie des sciences. Série 1, Mathématique, 322(8):729–734, 1996. [31] V. P ERROLLAZ. Asymptotic stabilization of entropy solutions to scalar conservation laws through a stationary feedback law. In Annales de l’Institut Henri Poincare (C) Non Linear Analysis, volume 30, pages 879–915. Elsevier, 2013. [32] P. R ICHARDS. Shock waves on the highway. Operations Research, 4(1):42–51, 1956. [33] N. S EGUIN and J. VOVELLE. Analysis and approximation of a scalar conservation law with a flux function with discontinuous coefficients. Mathematical Models and Methods in Applied Sciences, 13(02):221– 257, 2003. [34] D. S ERRE. Systems of conservation laws. Diderot, Paris, France, 1996. [35] L. S PINOLO. Vanishing viscosity solutions of a 2 x 2 triangular hyperbolic system with dirichlet conditions on two boundaries. 2006. [36] I. S TRUB and A. BAYEN. Weak formulation of boundary conditions for scalar conservation laws: an application to highway traffic modelling. International Journal of Robust Nonlinear Control, 16(16):733–748, 2006, doi: 10.1002/rnc.1099. [37] B. T EMPLE. Systems of conversation laws with invariant submanifolds. Transactions of the American Mathematical Society, 280(2):781–795, 1983.