A third-order weighted essentially non-oscillatory scheme in optimal control problems governed by nonlinear hyperbolic conservation laws

The weighted essentially non-oscillatory (WENO) methods are popular and effective spatial discretization methods for nonlinear hyperbolic partial differential equations. Although these methods are formally first-order accurate when a shock is present, they still have uniform high-order accuracy right up to the shock location. In this paper, we propose a novel third-order numerical method for solving optimal control problems subject to scalar nonlinear hyperbolic conservation laws. It is based on the first-disretize-then-optimize approach and combines a discrete adjoint WENO scheme of third order with the classical strong stability preserving three-stage third-order Runge–Kutta method SSPRK3. We analyze its approximation properties and apply it to optimal control problems of tracking-type with non-smooth target states. Comparisons to common first-order methods such as the Lax–Friedrichs and Engquist–Osher method show its great potential to achieve a higher accuracy along with good resolution around discontinuities.


Introduction
We consider the optimal control problem u min 0 = arg min u0∈U ad J(y(T, •; u 0 ), y d ) (1) with the tracking-type functional where G(y(T, x; u 0 ), and y = y(T, x; u 0 ) is the scalar entropy solution at the final time T > 0 of the nonlinear hyperbolic conservation law (later referred to as state equation) Here, u 0 ∈ U ad ⊆ L ∞ (R) is the control and y d ∈ L 2 (R) denotes a given target towards which we strive to optimize.We assume that the flux function satisfies f ∈ C m (R) with sufficiently large m ∈ N and is convex, the admissible set U ad is non-empty, convex and closed, and the region of integration I in ( 2) is a bounded interval.Weak solutions to (4) are in general not unique, which implies that the physically relevant solution has to be chosen.As a fact we cite the well-known result from [21], which states that for u 0 ∈ L ∞ (R)∩BV (R) there exists a unique entropy solution in the sense of Krǔzkov in the class C([0, T ], L 1 loc (R))∩L ∞ (R×[0, T ]).In this work, we focus on the numerical treatment of optimal control problems (1) governed by hyperbolic conservation laws, which has been studied amongst others in [1,2,5,6,10,11,12,16,17,18,22,23,24,27,28,29].We will follow the first-disretize-then-optimize approach, i.e., equation ( 4) is first discretized in space and time by applying a weighted essentially non-oscillatory (WENO) scheme and a strong stability preserving Runge-Kutta (SSPRK) method.This leads to a finite dimensional optimal control problem, for which the first-order discrete optimality system can be derived and solved by existing optimization solvers such as nonlinear Newton-type algorithms.In spite of the large size of the resulting problems, the flexibility of this approach naturally allows the incorporation of additional constraints and bounds.Further advantages are the direct use of automatic differentiation techniques and the computation of discrete adjoints, which are consistent with the discrete optimal control problem.Symmetric approximations of Hessian matrices can be easily derived and result in a computational speedup.
The application of common methods from nonlinear optimization requires the computation of directional derivatives of the target functional J with respect to the control.An efficient computation of the gradient can be effectuated by using the so-called adjoint approach, in which the derivative is represented via the adjoint state.The crucial issue of hyperbolic conservation laws is the possible formation of shocks even for smooth initial data, for which reason the classical adjoint calculus does not apply.To overcome these difficulties, nonstandard variational concepts have been developed in [4,27,29], which incorporate the shock sensitivity in order to derive rigorous optimality conditions.The resulting nonconservative equation has been studied in [3,7,27].Their numerical resolution is intricate, since the interior boundary condition defined on a set of Lebesgue measure zero -existing for the continuous setting -is not present for the discrete counterpart.This inherent problem has been addressed in [1,10,11,12].The theory is, however, restricted to differentiable monotone schemes which have sufficiently large numerical diffusion and are of first order only.
To avoid unwanted smearing of the solution by large numerical diffusion and to overcome the lower order restriction of monotone schemes, we propose a novel approach based on WENO schemes introduced in [20,25,26].These schemes have been proven to approximate hyperbolic equations comprising both shocks and complex smooth solution structure with higher accuracy and adequate stability along with good resolution around discontinuities.Although these methods are formally first-order accurate when a shock is present, they still have uniform high-order accuracy right up to the shock location.By employing a global fluxsplitting, the numerical flux function becomes classically differentiable and therefore allows to develop discrete adjoint WENO methods of higher order.Since the third-order WENO method is often applied in applications, we consider this method in the context of optimal control in more detail.We prove that the discrete adjoint WENO3 method is third-order consistent in space for smooth solutions.A fully discrete method is derived by applying a third-order SSPRK method.We present numerical results and study the approximation behaviour of the adjoint WENO3 scheme.Finally, we solve an optimal control problem with discontinuous target and compare the performance of our novel scheme to common first-order schemes such as the modified Lax-Friedrichs and the Engquist-Osher scheme.Further examples can be found in [9].

Adjoint Equation and Reversible Solutions
In this section, we briefly recall some theoretical basics in order to set up appropriate adjoint equations for hyperbolic conservation laws.As pointed out in [4, Example 1], the solution operator S t : u 0 → y(t, •; u 0 ) is generically not differentiable in L 1 loc (R), for which reason the classical adjoint calculus does not apply.However, in [27] it has been shown that entropy solutions to hyperbolic conservation laws admit a generalized differentiable structure called shift-differentiability.Under suitable assumptions, a generalized Taylor expansion in L 1 loc of the form exists for all δu 0 ∈ L ∞ (R), where , is a bounded linear operator and S (xi) y is the shift variation defined by where Ω i = [min(x i , x i + δx i ), max(x i , x i + δx i )] and x 1 , ..., x N denote the locations of the down-jumps of the entropy solution.The important advantage of shift-variations is that this framework allows to develop an adjoint calculus for hyperbolic conservation laws by using an averaged sensitivity equation which avoids the linearization of (4) in the usual way, see [29] for further details.The directional derivative of J in (2) in the direction of δu 0 can then be represented by where p is the solution of the adjoint equation Here, p T (x) is given by where X s is the set of locations where y(T, •) possesses a shock and [w(x)] := w(x−)−w(x+), which naturally incorporates the shock sensitivity.Equation ( 8) is a linear transport equation with, in general, discontinuous coefficients.It admits multiple solutions, which requires the selection of the correct adjoint state.This is achieved by so-called reversible solutions that are defined along generalized characteristics [8].An illustrative demonstration is given in Example 2.1.Under suitable technical assumptions and for appropriate end data p T it can be shown that there exists a unique reversible solution to (8) that is bounded, L ∞ -stable, and T V -stable [27, Theorem 4.2.10 and Corollary 4.2 .11].In what follows, we will work with formulation (8) to derive a discrete adjoint WENO3 method.
Example 2.1.Let f (y) = 1 2 y 2 , u 0 (x) = −sign(x), T = 0.5, and y d (x) = 0 with x ∈ R. It is well-known that the unique entropy solution is given by y(t, x) = −sign(x), t ∈ [0, T ], and hence we have The area that is not occupied by the classical characteristics is called shock funnel.It is represented by the grey-coloured triangle in Fig. 1.In this region, p takes the constant value zero.The adjoint remains constant along the classical backwards characteristics outside this region.Hence, the reversible solution p is given by 3 Discrete Adjoint WENO3 Method In order to discretize (4) in space, we now consider solutions with compact support [a, b] for the entire time interval [0, T ].Thus, using as many points outside as needed and the Figure 1: Construction of the reversible solution p(t, x) from the end data p T (x) at T = 0.5.The shock funnel region is accentuated as grey-colored triangle.
compactness of the solution, the implementation of the zero boundary condition does not bear any difficulty.The interval [a, b] is partitioned into subintervals [x j−1/2 , x j+1/2 ] of the same size ∆x and with midpoints x j for j = 1, . . ., N .Setting u 0 := (u 0 (x 1 ), . . ., u 0 (x N )) T and defining spatial approximations y(t) := (y 1 (t), ...., y N (t)) T with y j (t) ≈ y(t, x j ), a spatial semi-discretization of (4) reads where the nonlinear operator F ∆x : R N → R N represents the discretization of ∂ x f (y).We choose a conservative finite difference where fj+1/2 : R m → R denotes the numerical flux at x j+1/2 , which is (at least) a Lipschitz continuous function of m neighboring values y i (t).In order to avoid the convergence of the scheme towards entropy violating solutions, we apply a global flux splitting Using the simple Lax-Friedrichs splitting f ± (y) = (f (y) ± αy)/2 with α := max u |f ′ (u)| yields the desired properties (f + ) ′ (y) ≥ 0 and (f − ) ′ (y) ≤ 0.Then, the numerical flux functions of the WENO3 method [26] are defined by where the weights are The smoothness indicators are given by and the linear weights are set to γ . Note that 0 < ε ≪ 1 is chosen in order to avoid that the denominator becomes zero.It is set to ε = 10 −6 in our numerical calculations.We would like to emphasize the observation that by construction the numerical fluxes f ± have the same smoothness dependency on its arguments as that of the physical flux function f (y).
Next we will derive the associated adjoint WENO3 scheme.Let f ∈ C 2 (R), i.e., there exists the Fréchet derivative of F ∆x defined in (13).The continuous optimal control problem is approximated by where U ad = {u ∈ R N : TV(u) ≤ C} is the discrete admissible set.Then, applying the common Lagrangian approach in R N with multipliers p(t) = (p 1 (t), . . ., p N (t)) T , the adjoint equation to (12) reads where ∇ y F ∆x is the Fréchet derivative of F ∆x and gradients are treated as row vectors.The initial condition (the adjoint equation works backwards in time) is the discrete counterpart to (9).Observe that the interior boundary condition does not appear here.A short calculation yields the componentwise description with the coefficients The indices of the numerical flux functions are directly related to their arguments, e.g.f + j+3/2 (y j , y j+1 , y j+2 ) due to (15).For later use, we note that i=−2,...,2 L i,j (y) = 0.
We will now study the consistency order of the adjoint WENO3 scheme, i.e., how accurate does the semi-discretization (20) approximate the continuous adjoint equation (8) in the case of smooth solutions.Inserting exact solution values p(t, x j ) and y(t, x j ) (still denoted by y j to simplify notation) in the semi-discrete scheme (21) gives the residual-type local spatial errors Taylor expansion around x j yields where we have already used that the sum of the L i,j disappears.The method is said to have adjoint consistency order q if r j (t) = O(∆x q ).In what follows, we will show that the adjoint WENO3 scheme satisfies all conditions for order q = 3.First, we have to calculate ∂ yj L i,j , i.e., particularly the derivatives of the numerical flux functions defined in (15), (16).Since ω ± 1 + ω ± 2 = 1 for all y(t), we deduce we find and We have the following three lemmata.
Proof.Taylor expansion gives f ± j = O(∆x 2 ).It remains to show that ∂ y k ω ± 1 = O(∆x).Indeed, we have The two quotients are bounded by (ω show O(∆x) for these terms and therefore also for ∂ y k ω ± 1 .
Due to the strict positivity of the weights, it remains to study the asymptotic behaviour of the nominator.Using the definitions ( 17) and ( 18), we have with the smoothness indicators Taylor expansion at x j yields in (32) which shows the assertion.
Proof.We first consider ω + 1 − γ + 1 .The difference can be expressed by The denominator is bounded from below by (γ Further, we deduce for the nominator Let h(x) := f + (y(x)) and define h j := f + (y(x j )).Inserting the smoothness indicators , Taylor expansion at x j yields Putting this together with the bound for the denominator stated above gives ω + 1 − γ + 1 = O(∆x 3 /ε), from which we can conclude the proof.The same arguments apply to the second difference ω − 1 − γ − 1 .
We are now ready to state the main result of this section.
In a last step, we use the asymptotic expressions for d i , i = 0, 1, 2, to calculate the residual-type local spatial error This concludes the proof.

Numerical Experiments
In this section, we will present some numerical examples for Burgers equation, i.e., we study problems with the nonlinear flux function f (y) = 1 2 y 2 in (4).The first example with smooth initial data and solution is chosen in order to check to third-order convergence of the discrete adjoint WENO3 method as stated in Theorem 3.1.In the second example, the approximation property of the discrete adjoint in the case of a shock in the initial solution is investigated and compared to approximations computed by means of the first-order modified Lax-Friedrichs (LF) and Engquist-Osher (EO) schemes.These schemes read with y n j ≈ y(n∆t, x j ), n T ∆t = T , and numerical fluxes given by Applying a standard Lagrangian approach and discrete adjoint calculus, the discrete adjoint schemes can be derived from [16,Prep. 3.1] as with the coefficients for the LF scheme and for the EO scheme.Convergence of these schemes has been intensively studied in [1,10,11,12,27].The choice γ = 1 leads to the classical LF method.Stability requirements for the adjoint LF and EO schemes yield the optimal value γ ⋆ = 0.5 together with the CFL-condition ∆t ≤ γ ⋆ ∆x/ sup |f ′ (y)|, see e.g.[16].Then, both schemes converge for Lipschitz continuous end data p T (x) in (8).The stronger condition ∆t ≤ γ ⋆ (∆x) 2−q / sup |f ′ (y)|, 0 < q < 1, ensures the convergence of the modified LF scheme for discontinuous end data, too [11,12].
Convergence for slightly modified end data and less numerical viscosity has been recently studied in [1].
In order to get a fully discrete scheme for WENO3, the differential equation ( 12) is numerically solved by the three-stage third-order strong-stability-preserving Runge-Kutta initial condition.The advantage of this approach is that it delivers a control whose entropy solution is close to the target and the location of the discontinuities almost coincide.Hence, the production of additional discontinuities within each iteration step is avoided, which improves the performance of the algorithm drastically.We will exemplify the influence of the choice of the initial guess in our optimal control problem.

Order Test for the Discrete Adjoint for Smooth Data
This section is devoted to numerically verify the third-order convergence of the adjoint WENO3 scheme.For this purpose, we choose the computational domain Ω T = (0, 0.5] × [−1.5, 1.5] and the objective functional with the smooth initial data The exact solution y(t, x) can be directly computed from the method of characteristics, i.e., y(t, x) = u 0 (x 0 (x(t), t)) with x 0 (x(t), t) being the solution of the nonlinear equation x(t) = x 0 + u 0 (x 0 ) t.A reference solution y T ≈ y(0.5, x) at the final time is computed by Newton's method with a high tolerance 10 −14 .Since shocks are not present, we find p T (x) = y(0.5,x) in (8).We also note that the characteristics curves of the adjoint problem coincide with the characteristic curves of the forward problem.Thus, the corresponding reversible solution p(0, x) at time t = 0 is given by u 0 (x), which serves as reference solution for the adjoint.
We use a sequence of spatial meshes with a number of grid points N = 150 • 2 i , i = 0, . . ., 6, and set ∆t = 0.5 ∆x.In order to keep the temporal error below O((∆x) 3 ), we apply the classical fourth-order four-stage explicit Runge-Kutta method (ERK4).Its adjoint time discretization has also order four [15] for smooth solutions and therefore the overall scheme is suitable to check the order three of the adjoint WENO3 method.We also present results for the forward WENO3 method to document the error of the approximated starting value p nT = y nT .The L ∞ -errors collected in Tab. 1 clearly show asymptotic order three of the spatial WENO3 discretization for both forward and adjoint numerical solution.

Approximation of the Discrete Adjoint in the Case of Shocks
We now consider discontinuous solutions with shocks.Our test case is taken from Example 2.1 with computational domain Ω T = (0, 0.5] × [−1, 1].The reversible solution p(0, x) at t = 0 is given by (11).We apply the above described forward and adjoint LF, EO, and WENO3 schemes with ∆x = 0.01, 0.002, and ∆t = 0.25 ∆x.The corresponding numerical approximations p 0 are shown in Fig. 2.
The first-order LF and EO schemes smear out the discontinuities, but deliver L ∞ -stable approximations and thus respect the analytical property of the adjoint.In the spirit of WENO schemes, the adjoint WENO3 delivers a quite sharp resolution of the shocks at the price of bounded over-and undershoots of around 5%.In Tab.∞ at time t = 0 for a sequence of spatial meshes with N = 150, 300, . . ., 9600 grid points.The convergence rates are computed from ln(E N /E 2N )/ ln (2), where E N stands for the corresponding error. in the shock funnel for x ∈ [−0.3, 0.3].All schemes converge quite rapidly.Note that convergence in the shock funnel is not always achieved since the interior boundary condition at shock positions as given in (9) does not appear on the discrete level, see the discussions in [1,11,12].

Optimal Control Problem with Discontinuous Target
We consider the optimal control problem (1) with the objective functional [16] J(y(0.5, and the discontinuous target y d defined by The optimal control u ⋆ 0 , which serves as a reference solution, is given in (11) for the adjoint Lax-Friedrichs (LF), Engquist-Osher (EO) and WENO3 scheme applied with ∆x = 0.01 (left), ∆x = 0.002 (right), and ∆t = 0.25 ∆x.
We will present results for two mesh sizes ∆x = 0.005, 0.002, and time steps ∆t = 0.25 ∆x.The initial guess for the control is computed from (52) with the individual method under consideration.For WENO3 and the coarser mesh size, it is shown in Fig. 3 together with the corresponding state solution.
In Fig. 4, the results of the gradient based optimization procedure described above for tolerances tol 1 = 10 −5 , tol 2 = 10 −7 , and mesh size ∆x = 0.005 for the adjoint WENO3 method are plotted.We can conclude that the adjoint WENO3 method allows to recover the initial data together with the final state solution adequately.The shock of the target is sharply resolved and the rarefaction of the initial data is also recovered.In order to compare these results with those obtained from the LF and EO schemes, we perform 50 iterations of the optimization algorithm for both mesh sizes.The calculated optimal controls and their corresponding final state solutions are collected in Fig. 5.The adjoint WENO3 method resolves the shock sharply.In contrast, the LF method is too diffusive and only provides an unsatisfactory shock resolution.The numerical artifacts around the shocks are huge.The optimized final state solution obtained by the EO scheme possesses very small numerical artifacts, but the shock is less sharply resolved and the spike of the target is slightly smeared out.In Tab. 3, we depict the iteration history for all runs of the optimization.In every case, the LF method performs poorer than the others.In terms of a low cost functional, the adjoint WENO3 method performs best.We also see the influence of the initial guess on the performance of the algorithm.This is due to the fact that the use of u 0 = 0 as starting control value produces artificial discontinuities within each iteration step.Table 3: Optimal control problem.Logarithmic values of the objective functional (55) at the beginning and after 50 iterations of the optimization algorithm, J 0 and J 50 , respectively.For comparison, values for an initial control u 0 = 0 for WENO3 are shown, too.

Summary
We have developed a novel adjoint WENO3 scheme to provide approximations of the gradient for optimal control problems governed by hyperbolic conservation laws and proved thirdorder consistency in space for sufficiently smooth solutions.The adjoint WENO3 method is able to sharply resolve discontinuities of reversible solutions.For an exemplary optimal control problem with discontinuous target, the method works very well and outperforms common first-order methods as the Lax-Friedrichs and Engquist-Osher schemes.

Acknowledgement
This work was supported by the Graduate School CE within the Centre for Computational Engineering at Technische Universität Darmstadt and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the collaborative research centre TRR154 "Mathematical modeling, simulation and optimisation using the example of gas networks" (Project-ID 239904186, TRR154/2-2018, TP B01).

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Optimal control problem.Initial control u 0 and optimal control u ⋆ 0 (left), initial state solution y nT at T = 0.5 and target y d (right), computed with the WENO3 scheme and mesh size ∆x = 0.005.

Table 1 :
2, we plot the L ∞ -error Burgers problem with smooth initial data and smooth solution: L ∞ -error of the forward solution y T − y nT ∞ at the final time T = 0.5 and adjoint solution u 0 −p 0