Conservation laws in the presence of shocks

The problem

We analyze a model tracking problem for a 1D scalar conservation law. It consists in optimizing the initial datum so to minimize a weighted distance to a given target during a given finite time horizon. To be more precise, given a finite time , a target function , and a positive weight function with compact support in , we consider the functional cost to be minimized , over a suitable class of initial data , defined by

(1)

where is the unique entropy solution of the scalar conservation law

(2)

Thus, the problem under consideration reads: To find such that

(3)

Here the flux is assumed to be smooth: . The initial datum will be assumed to belong to a suitable admissible class to ensure the existence of a minimizer.

Sensitivity of the state in the presence of shocks

Inspired in several results on the sensitivity of solutions of conservation laws in the presence of shocks in one-dimension (see [7, 4, 5, 6, 17, 13]), we focus on the particular case of solutions having a single shock. But the analysis can be extended to consider more general one-dimensional systems of conservation laws with a finite number of noninteracting shocks. We introduce the following hypothesis:

Hypothesis 1.1:Assume that is a weak entropy solution of (2) with a discontinuity along a regular curve , which is Lipschitz continuous outside . In particular, it satisfies the Rankine-Hugoniot condition on

Here we have used the notation: , for the jump at of any piecewise continuous function with a discontinuity at .

Note that divides into two parts: and , the sub-domains of to the left and to the right of respectively.
As we will see, in the presence of shocks, to deal correctly with optimal control and design problems, the state of the system needs to be viewed as constituted by the pair combining the solution of (2) and the shock location . This is relevant in the analysis of sensitivity of functions below and when applying descent algorithms.
We adopt the functional framework based on the generalized tangent vectors (see [7] and Definition 4.1 in [8]).
Let be the initial datum, that we assume to be Lipschitz continuous to both sides of a single discontinuity located at , and consider a generalized tangent vector for all . Let be a path which generates . For sufficiently small, the solution of (2) is Lipschitz continuous with a single discontinuity at , for all . Therefore, generates a generalized tangent vector . Moreover, in [8] it is proved that it satisfies the following linearized system:

(4)

Sensitivity of the cost in the presence of shocks

In this section we study the sensitivity of the functional with respect to variations associated with the generalized tangent vectors defined in the previous section. We first define an appropriate generalization of the Gateaux derivative of .

Definition 1.2:Let be a functional and be Lipschitz continuous with a discontinuity in , an initial datum for which the solution of (2) satisfies hypothesis (1.1). is Gateaux differentiable at in a generalized sense if for any generalized tangent vector and any family associated to the following limit exists,

and it depends only on and , i.e. it does not depend on the particular family which generates . The limit is the generalized Gateux derivative of in the direction .

The following result easily provides a characterization of the generalized Gateaux derivative of in terms of the solution of the associated adjoint system (6).

Proposition 1.3:The Gateaux derivative of can be written as follows

(5)

where the adjoint state pair satisfies the system

(6)

Let us briefly comment the result of Proposition 1.3 before giving its proof.
System (6) has a unique solution. In fact, to solve the backward system (6) we first define the solution on the shock from the conditions for . This determines the value of along the shock. We then propagate this information to both sides of , by characteristics (see Figure 1 where we illustrate this construction).

Figure 1: Characteristic lines entering on a shock and how they may be used to build the solution of the adjoint system both away from the shock and on its region of influence.

Formula (5) provides an obvious way to compute a first descent direction of at . We just take

(7)

Here, the value of must be interpreted as the optimal infinitesimal displacement of the discontinuity of .

In [8], when considering the inverse design problem, it was observed that the solution of the corresponding adjoint system at was discontinuous, with two discontinuities, one in each side of the original location of the discontinuity at . This was a reason not to use this descent direction and for introducing the alternating descent method. In the present setting, however, the adjoint state obtained is typically continuous. This is due to the fact that at both side of the discontinuity is defined by the method of characteristics and that, on the region of influence of the characteristics emanating from the shock, the continuity is preserved by the fact that, on one hand, itself is continuous as the primitive of an integrable function and that the data for and at are continuous too. Despite of this, as we shall see, the implementation of the alternating descent direction method is worth since it significantly improves the results obtained by the purely discrete approach.

Numerical Experiments

In this section we present some numerical experiments which illustrate the results obtained in an optimization model problem with each one of the numerical methods described in the previous section.
We have chosen as computational domain the interval and we have taken as boundary conditions in the numerical scheme, at each time step , the value of the initial data at the boundary. This can be justified if we assume that the initial datum is constant in a sufficiently large inner neighborhood of the boundary (which depends on the size of the -norm of the data under consideration and the time horizon ), due to the finite speed of propagation. A similar procedure is employed for the adjoint equation.
We underline once more that the solutions obtained with each method may correspond to global minima or local ones since the gradient algorithm does not distinguish them.
In the experiments we consider the Burgers’ equation, i.e. ,

(8)

The weight function , under consideration in the experiments is given by

And the time horizon .
To compare the efficiency of the different methods we consider a fixed (which satisfies the CFL condition). We then analyze the number of iterations that each method needs to attain a prescribed value of the functional.

Experiment 1

We first consider a piecewise constant target profile given by the solution of 8 with the initial condition given by

(9)

Note that, in this case, 9 yields a particular solution of the optimization problem and the minimum value of vanishes.
We solve the optimization problem 3 with the above described different methods starting from the following initialization for :

(10)

which also has a discontinuity but located on a different point.

Figure 2: Experiment 1. versus the number of iterations in the descent algorithm for the discrete and the alternating descent methods.

In Figures 3 and 4, we present the minimizers obtained by the methods above, and the associated solutions, Figures 5 and 6.
The initial datum obtained by the alternating descent method (Figures 4) is a good approximation of 9. The solution given by the discrete approach (Figure 3) presents added spurious oscillations. Furthermore, the discrete method is much slower and does not achieve the same level of accuracy since the functional does not decrease so much.

Experiment 2

The previous experiment indicates that the alternating descent method performs significantly better. In order to show that this is a systematic fact, which arises independently of the initialization of the method, we consider the target given by the solution of 8 with the initial condition given by

(11)

but this time we compare the performance of both methods starting from different initializations.
The obtained numerical results are presented in Figure 7.
We see that, regardless the initialization considered, the alternating descent method performs significantly better.
We observe that in the five experiments the alternating descent method perfumes better ensuring the descent of the functional in much fewer iterations and yielding smoother, less oscillatory approximation of the minimizer.
Note also that the discrete method, rather than yielding discontinuous approximations of the minimizer as the alternating descent method does, it produces an initial datum with a Lipschitz front. Observe that these are two different configurations that can lead to the same evolution for the Burgers equation after some time, once the front develops the discontinuity. This is in agreement with the fact that the functional to be minimized is only active in the time-interval .