key: cord-0045679-pjzelxs9
authors: Herty, Michael; Salhi, Loubna; Seaid, Mohammed
title: A Relaxation Algorithm for Optimal Control Problems Governed by Two-Dimensional Conservation Laws
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50426-7_10
sha: a43e966ae116f50ffc77d704ef3eb803d7165a35
doc_id: 45679
cord_uid: pjzelxs9

We develop a class of numerical methods for solving optimal control problems governed by nonlinear conservation laws in two space dimensions. The relaxation approximation is used to transform the nonlinear problem to a semi-linear diagonalizable system with source terms. The relaxing system is hyperbolic and it can be numerically solved without need to either Riemann solvers for space discretization or a non-linear system of algebraic equations solvers for time discretization. In the current study, the optimal control problem is formulated for the relaxation system and at the relaxed limit its solution converges to the relaxed equation of conservation laws. An upwind method is used for reconstruction of numerical fluxes and an implicit-explicit scheme is used for time stepping. Computational results are presented for a two-dimensional inviscid Burgers problem.

In many applications, optimal control problems consist of a class of differential equations whose evolution and the behavior of their solutions can be controlled by involving external control laws. In the current study, we are interested in optimal control problems subject to the following two-dimensional nonlinear conservation law ∂ t u + ∇ · F(u) = 0, (x, y) ∈ Ω, t > 0, (1a) u(0, x, y) = u 0 (x, y),

where Ω is an open bounded domain in R 2 , (x, y) the space coordinates, t the time, u(t, x, y) is the control function, u 0 (x, y) the initial state and the flux F(u) = (f (u), g(u)) T , with f (u) and g(u) are nonlinear functions. In practice, optimal control problems require minimizing a cost functional J (u(T, x, y); u d (x, y)) based on the least-square method that associates a cost value to each possible behavior. Thus, the problem statement is min u0 J (u(T, x, y); u d (x, y)) := min

subject to the conservation law (1a). In (1b), u d (x, y) is the desired state at the final time T . Optimal control problems of type (1b) have received growing attention in both theoretical and numerical studies over recent decades. In most of these studies, control problems governed by hyperbolic equations have been less extensively treated compared to elliptic and parabolic control problems. This is mainly due to the fact that the semi-group generated by the hyperbolic conservation law is non-differentiable in L 1 whereas its domain of definition is an L 1 closed subset of BV . In the case of nonlinear conservation laws in one space dimension, a differential structure on general BV solutions has been presented and discussed in [4, 19] among others. The first-order optimality conditions for hyperbolic systems have been introduced in [5] based on the derived calculus. It turned out that the resulting adjoint equations are non-conservative which fail to recover stable solutions for problems with shocks. In [13, 14, 19] , numerical results for one-dimensional scalar problems with distributed control have been presented. More results for the case of a one-dimensional linear hyperbolic systems can be found in [8, 10, 15] . In [11] , a TVD Runge-Kutta method for the time discretization of such problems has been employed. It was shown that requiring high stability for both the discrete and adjoint states is too strong, limiting the method to first-order, regardless of the number of stages used in the method. Using the same discretization, authors in [11] have studied other conditions for the discrete adjoint such that the numerical approximation is of the best possible order. In [1] , the emphasis was placed on high-order linear multistep schemes for the time discretization of adjoint equations arising within optimal control problems. The authors reported that the so-called Adams methods may reduce to the first-order accuracy and that only BDF schemes may be used as higher order discretization for the hyperbolic relaxation systems in combination with a Lagrangian scheme. Theoretical and numerical methods using finite difference schemes combined with an immersed boundary method have been developed in [9] for a special class of optimal control problems namely, problems involving the shallow water equations and a geometric parameter to be optimized in the terminal cost. More recently, theoretical studies including a posteriori error estimates have been carried out for numerical schemes to solve multi-dimensional problems, based on adjoint equations, see for instance [16, 18] .

In the present work, we are interested in developing numerical algorithms for control problems of two-dimensional nonlinear conservation laws to achieve numerical stability without need to inclusion of extra artificial diffusion in the problem under study. For this purpose we consider the relaxation approximation of nonlinear conservation laws in the same manner as introduced in [12] . This approach approximates the nonlinear problem to semi-linear system with linear characteristic speeds, while preserving the hyperbolic structure on the expense of an additional equation and stiff source terms. Thus, the resulting relaxation system is semi-linear which allows for a Riemann-solver free treatment. The relaxation methods have been investigated by many authors, see [3] among others. First studies of relaxation systems with respect to control problems have been reported in [2] in case of one-dimensional scalar conservation laws. Numerical results are still very limited in the multi-dimensional cases and we therefore restrict ourselves to a numerical study including a first-order relaxation approximation. For the space discretization, we consider an upwind reconstruction of the numerical fluxes and an implicit-explicit method is used for the time integration.

The remainder of this paper is structured as follows. In Sect. 2, the relaxation approximation for the coupled optimal control problem and the nonlinear conservation laws is formulated. The space and time discretizations along with the approximation procedure of the solution is presented in Sect. 3. In Sect. 4, numerical results are presented for a test example of inviscid Burgers equation. Section 5 contains concluding remarks.

Following [12] , the relaxation approximation for (1a) allows to construct a corresponding linear hyperbolic system with a stiff source term that approximates the original problem with a small dissipative correction. Thus, the relaxation associated with (1a) reads

where τ is a small positive parameter that measures the relaxation rate, v and w are the relaxation variables, a 2 and b 2 are the characteristic speeds satisfying the sub-characteristic condition [12] f (u) 2

The initial conditions for the relaxation system (2a) are selected as

It is clear that, when τ tends to 0, the relaxation system (2) converges to the system of conservation law (1a). Note that the main advantage of numerically solving the relaxation system (2) over the original conservation law (1a) lies in the special structure of the linear characteristic fields and localized lower order terms. Indeed, the linear hyperbolic nature of (2) allows to approximate its solution easily by underresolved stable numerical discretization that uses neither Riemann solvers spatially nor nonlinear system of algebraic equations solvers temporally. Hence, using the relaxation approximation, the optimal control problem (1b) becomes

subject to the relaxation system (2a). Notice that a formal adjoint calculus leads to a first-order optimality conditions for the function u 0 . The calculations are rigorous provided that the solutions have sufficient regularity which however in general is not the case. Hence, we formulate the adjoint equations for the system (2) as

with terminal conditions given by

It should be stressed that the adjoint equations (4) have to be solved backwards in time and the gradient of the reduced cost functional is defined as

Again, when τ tends to 0, the system (4) converges to the adjoint problem associated with the conservation law (1a). Then, from the second and third equations in (4a), an expansion in terms of τ gives

Inserting these terms in the first equation of (4a) leads to

which is a viscous approximation to the formal adjoint of (4). Note that the gradient eventually vanishes at the minimum of the cost functional. Since u might develop discontinuities we will have to scope with discontinuous derivatives of the flux functions f (u) and g (u). However, since we use the relaxation approximation, the derivative functions f (u) and g (u) appear as source terms and not as a discontinuous transport coefficient as in (6). This problem has been investigated for one-dimensional problems in [19] . However, as pointed out in [2, 19] , the problem reappears in the small τ limit. In the one-dimensional case it can be shown that the first-order relaxation discretization converges to the reversible solution of a transport equation with discontinuous coefficient. Here, we focus on a numerical study of the optimality system (1a). Using the characteristic

an equivalent system associated with (2) can be reformulated as

The adjoint equations in characteristic form are therefore given by

This system is equivalent to a spatial splitting approximation of the adjoint equations (4). Introducing

we obtain from the equations in s that the solutions (p, q) satisfy

Similarly, for

Hence, the formulation (9) is precisely the spatial splitting applied to (4). Therefore, the adjoints in characteristic form are the same as the adjoint of the characteristic form when applying a dimensional splitting in the spatial variable. For the optimize-then-discretize approach discussed below it is therefore sufficient to state the discretization of the forward equations in characteristic form. A rigorous discussion of the relation between discrete adjoints, the characteristic variables and higher-order schemes can be found in [2] in the case of one-dimensional scalar advection equations.

Relaxation schemes are in fact a combination of non-oscillatory upwind space discretization and an implicit-explicit time integration of the resulting semidiscrete system, see for instance [3, 12] . The fully discrete system of the equations (2a) is referred to as a relaxing system, while that of the limiting system as the relaxation rate τ tends to zero is called a relaxed system. In this section, we formulate the space and time discretizations used for the numerical solution of optimal control problems and also formulate the algorithm used for the discrete gradient.

For the space discretization of the equations (2a), we cover the spatial domain with rectangular cells 1 2 ] of uniform sizes Δx and Δy for simplicity. The cells, C i,j , are centered at (x i = iΔx, y j = jΔy). We use the notations ω i± 1 2 

and

to denote the point-values and the approximate cell-average of a generic function ω at (x i± 1 2 , y j , t n ), (x i , y j± 1 2 , t n ), and (x i , y j , t n ), respectively. We define the following finite differences

Then, the semi-discrete approximation of (2a) is

Similarly, the semi-discrete approximation of the adjoint equations (4a) is

Most relaxation schemes can be described as fractional step methods, in which the relaxation step is just a projection of the system into the local equilibrium.

The fully-discrete formulation of systems (11) and (12) can be obtained by the well-established IMEX methods, see for instance [17] . Indeed, the special structure of the nonlinear terms in (11) and (12) makes it trivial to evolve the flux terms explicitly and the stiff source terms implicitly. The semi-discrete formulations (11) or (12) can be rewritten in common ordinary differential equations notation as

where the time-dependent vector functions are defined accordingly for the forward problem (11) or for the backward problem (12) . Due to the presence of stiff terms in (13), one can not use fully explicit schemes to integrate the equations (13), particularly when τ tends to 0. On the other hand, integrating the equations (13) by fully implicit scheme, either linear or nonlinear algebraic equations have to be solved at every time step of the computational process. To find solutions of such systems is computationally very demanding. In this paper we consider an alternative approach based on the implicit-explicit (IMEX) Euler method. The non stiff stage of the splitting for F is straightforwardly treated by an explicit scheme, while the stiff stage for G is approximated by a diagonally implicit scheme. Let Δt = t n+1 −t n be the time step and Y n denotes the approximate solution at t = nΔt. We formulate the first-order IMEX scheme for the forward system (13) as

For the backward system (13) , the IMEX scheme is implemented as

Note that, using the above relaxation scheme neither linear algebraic equation nor nonlinear source terms can arise. In addition the relaxation schemes are stable independently of τ , so that the choice of Δt is based only on the usual CFL condition

where δ denotes the maximum cell size, δ = max(Δx, Δy). For the space discretization, a first-order upwind scheme is applied to the characteristic variables in (11) to obtain the numerical fluxes as

Thus, a first-order reconstruction of the numerical fluxes in the forward problem (11) yields

The numerical fluxes in the backward problem (12) are obtained by applying first order upwind scheme to the characteristic variables

Thus, a first-order reconstruction of the numerical fluxes in the backward problem (12) yields

In this study, the characteristic speeds a and b in the relaxation systems (2) and (4) are calculated locally at every cell as

It is worth saying that, larger a and b values usually add more numerical dissipation.

The implementation of the iterative optimization along with the Eulerian-Lagrangian numerical approach used in the implementation are performed in the same way as detailed in [7] . Thus, starting from the basic optimal control problem formulated as follows: Given a terminal state u d (x, y), find an initial datum u 0 (x, y) which by time t = T will either evolve into u(T, x, y) = u d (x, y) or will be as close as possible to u d in the L 2 -norm. To solve the problem iteratively, we implement the Algorithm 1 and generate a sequence of solutions u 0 (x, y) and w(x, y, 0) = g u (0) 0 (x, y) forward in time from t = 0 to t = T by using the relaxation method to obtain u (0) (T, x, y) , v (0) (T, x, y) = f u (0) 0 (T, x, y) and w (0) (T, x, y) = g u (0) 0 (T, x, y) .

-Solve the linear system (4a) backward in time from t = T to t = 0 using the relaxation method to obtain p (m) (0, x, y), q (m) (0, x, y) and r (m) (0, x, y). -Update the control u 0 , v 0 and w 0 using either a gradient descent or quasi-Newton method as described in [7] . -Solve the problem (2) subject to u(0, x, y) = u does not have to be stored during the iterations by using the developed method. In addition, although Algorithm 1 is similar to the continuous approach used in [6] , the focus is on the proposed numerical method to solve the problem (2) and thus, we do not need an approximation to the generalized tangent vectors to improve the gradient descent method.

To examine the performance of the relaxation algorithm to solve optimal control we present numerical results for a two-dimensional inviscid Burgers problem. In all the computational results presented in this section, the characteristic speeds a and b are locally chosen as in (21), the CFL number is fixed to 0.5 and time steps Δt are calculated according to the condition (16) . Here, the flux functions are defined by

The optimal control problems are solved in the domain [0, 1] × [0, 1] subject to period boundary conditions and equipped with the following initial data u(0, x, y) = sin 2 (πx) sin 2 (πy).

We solve the optimization problem for terminal time T = 0.2 using a relaxation rate τ = 10 −6 on three different meshes with 100 × 100, 200 × 200 and 400 × 400 control volumes. For each of these runs, we display the initial data u 0 , reference solution and the optimized solution u t along with the gradient of the reduced cost functional defined in (5). In Fig. 1 we present numerical results obtained on a mesh with 100 × 100 control volumes. Those results obtained on meshes with 200 × 200 and 400 × 400 control volumes are displayed in Fig. 2 and Fig. 3 , respectively. It is clear that the proposed algorithm resolves the desired solution for this problem and it captures all small features appearing in computational domain. The reference solution and the initial condition appear to be similar confirming the convergence of the proposed numerical techniques. As can be seen in the presented results, a shock is formed in the solution u t propagating along the main diagonal in the domain. The effect of mesh refinement on the computed solutions is noticeable in these figures. It is also clear that our relaxation methods accurately capture the shock and its propagation along the diagonal. However, due to the numerical dissipation, the resolved shock has been smeared out in the results obtained on a mesh with 100 × 100 control volumes. As expected, the numerical results obtained on this mesh are more diffusive than those computed using meshes with 200 × 200 and 400 × 400 control volumes. To further visualize this effect we display in Fig. 4 the cross-sections along the main diagonal y = x for the results on the considered meshes. It is apparent that the gradient resolution is deteriorated with the excessive dissipation included by the coarse mesh with 100 × 100 control volumes. On the other hand, the solutions are completely free of spurious oscillations and the shocks are well resolved by the proposed method without nonlinear computational tools. It should be that the number of iterations in the optimal control problem does not overpass 23 iterations for all considered meshes. These features clearly demonstrate the efficiency achieved by the proposed method for solving optimal control problems for the inviscid Burgers equation. The performance of the method is very attractive since the computed solution remains stable and accurate even when coarse meshes are used without requiring Riemann solvers or complicated techniques to reconstruct the numerical fluxes.

A class of numerical methods for solving optimal control problems governed by nonlinear conservation laws in two space dimensions has been presented and assessed. As solvers for the forward and backward problems we implement a relaxation method combining the upwind reconstruction for space discretization and implicit-explicit scheme for time integration. These techniques solve the nonlinear conservation laws without relying on Riemann solvers or linear solvers of algebraic equations. The optimal control problem is formulated for the relaxation system and at the relaxed limit its solution converges to the relaxed equation of conservation laws. The proposed method has been tested on an optimal control problem for the two-dimensional inviscid Burgers. The obtained results indicate good shock resolution with reasonable accuracy in smooth regions and without any nonphysical oscillations near the shock areas. Although, we have studied only the case of first-order relaxation methods, the extension to high-order reconstructions would be an encouraging next step and requires an in-depth study on optimal control problems to deal with the nonlinear structure of hyperbolic systems of conservation laws. Finally, we should point out that d the algorithm presented in this paper can be highly optimized for vector computers, because it does not require nonlinear solvers and contain no recursive elements. Some difficulties arise from the fact that for efficient vectorization the data should be stored contiguously within long vectors rather than two-dimensional arrays.

Linear multistep methods for optimal control problems and applications to hyperbolic relaxation systems

Adjoint IMEX-based schemes for control problems governed by hyperbolic conservation laws

Higher-order relaxation schemes for hyperbolic systems of conservation laws

On the shift differentiability of the flow generated by a hyperbolic system of conservation laws

Optimality conditions for solutions to hyperbolic balance

An alternating descent method for the optimal control of the inviscid burgers equation in the presence of shocks

An Eulerian-Lagrangian method for optimization problems governed by multidimensional nonlinear hyperbolic PDEs. Comput

Optimal time for the controllability of linear hyperbolic systems in one-dimensional space

Optimal control problem for viscous systems of conservation laws, with geometric parameter, and application to the Shallow-Water equations

Solutions of lp-norm-minimal control problems for the wave equation

Total variation diminishing schemes in optimal control of scalar conservation laws

The relaxation schemes for systems of conservation laws in arbitrary space dimensions

Optimal, globally constraint-preserving, DG (TD) 2 schemes for computational electrodynamics based on two-derivative Runge-Kutta timestepping and multidimensional generalized riemann problem solvers-a von Neumann stability analysis

On the properties of discrete adjoints of numerical methods for the advection equation

Optimal boundary control of hyperbolic equations with pointwise state constraints

The relation between primal and dual boundary conditions for hyperbolic systems of equations

Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations

The group-theoretical analysis of nonlinear optimal control problems with hamiltonian formalism

A sensitivity and adjoint calculus for discontinuous solutions of hyperbolic conservation laws with source terms