# Numerical Methods for the Optimal Control of Scalar Conservation Laws

Sonja Steffensen, Michael Herty, Lorenzo Pareschi

► To cite this version:

Sonja Steffensen, Michael Herty, Lorenzo Pareschi. Numerical Methods for the Optimal Control of Scalar Conservation Laws. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.136-144, 10.1007/978-3-642-36062-6_14. hal-01347531

HAL Id: hal-01347531

https://inria.hal.science/hal-01347531

Submitted on 21 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Methods for the Optimal Control of Scalar Conservation Laws

Sonja Steffensen, Michael Herty*, and Lorenzo Pareschi **

RWTH Aachen University, Templergraben 55,
D-52065 Aachen, GERMANY

`{herty,steffensen}@mathc.rwth-aachen.de`

University of Ferrara, Department of Mathematics,
Via Machiavelli 35, I-44121 Ferrara, ITALY

lorenzo.pareschi@unife.it

**Abstract.** We are interested in a class of numerical schemes for the optimization of nonlinear hyperbolic partial differential equations. We present continuous and discretized relaxation schemes for scalar, one-conservation laws. We present numerical results on tracking type problems with nonsmooth desired states and convergence results for higher-order spatial and temporal discretization schemes.

**Keywords:** IMEX schemes, optimal control, conservation laws, Runge-Kutta methods

## 1 Introduction

We consider an optimal control problem for scalar conservation laws of the type

$$ \begin{array}{ll} \underset{u_0}{\text{minimize}} & J(u(T), u_0) \\ \text{subject to} & u_t + f(u)_x = 0, \quad u(0,x) = u_0(x), \end{array} \qquad (1) $$

Here, $J$ and $f$ are assumed to be smooth and possibly nonlinear functions. The initial value $u_0$ acts as control to the problem. It can be observed that the wave interactions that occur in the solution $u$ in the case of a nonlinear flux function $f$ pose the serious analytical challenges. Recently, the differentiability of $J$ with respect to $u_0$ could be proven in the sense of shift-differentiability. We refer to [6, 9–11, 28–30, 4, 32] for more details.

Here, a class of numerical methods applied to the optimal control problem (1) is studied. We only consider the case of smooth initial data and smooth solutions $u$ and refer to [4] for more details. For a numerical analysis including shock waves and in the case of the Lax-Friedrichs scheme we refer to [21, 32] and the references therein.

<footnote>* RWTH Aachen University, Templergraben 55, D-52065 Aachen, GERMANY. `{herty,steffensen}@mathc.rwth-aachen.de`</footnote>

<footnote>** University of Ferrara, Department of Mathematics, Via Machiavelli 35, I-44121 Ferrara, ITALY. `lorenzo.pareschi@unife.it`</footnote>
## 1.1 Relaxation Method

As motivation for a numerical scheme we follow the ideas of Jin and Xin [22]. Therein, a linear approximation (2) of the nonlinear hyperbolic equation

$$ \partial_t u + \partial_x f(u) = 0 $$

has been discussed. For initial conditions $u(x, 0) = u_0$ the approximation is

$$ \begin{aligned} \partial_t u + \partial_x v = 0, \qquad & u(x, 0) = u_0, \\ \partial_t v + a^2 \partial_x u = \frac{1}{\epsilon} (f(u) - v), \qquad & v(x, 0) = f(u_0) \end{aligned} \qquad (2) $$

where $\epsilon > 0$ is the relaxation rate and $a$ is a given constant satisfying the subcharacteristic condition $\max_u |f'(u)| \le a$. For $\epsilon$ being small, the solution $u$ of (2) satisfies $\partial_t u + \partial_x f(u) = \epsilon \partial_x ((a^2 - f(u)^2)\partial_x u)$ (cf. [22]). Applying the relaxation to the optimal control problem (1), we obtain

$$ \min_{u_0} J(u(\cdot, T), u_0) \quad \text{subject to} \quad \begin{cases} u_t + v_x = 0, \\ v_t + a^2 u_x = \frac{1}{\epsilon} (f(u) - v), \\ u(0, x) = u_0, \ v(0, x) = f(u_0) \end{cases} \qquad (3) $$

The corresponding adjoint equations for (3) are given by (cf. [?])

$$ \begin{aligned} -p_t - a^2 q_x &= \frac{q}{\epsilon} f'(u), \quad p(T, x) = p_T(x), \\ -q_t - p_x &= -\frac{q}{\epsilon}, \quad q(T, x) = q_T(x). \end{aligned} $$

For more information on the relaxation system, its limiting scheme for $\epsilon = 0$, further numerical analysis and extensions we refer to [1, 2, 5, 14, 3, 8, 22, 25, 27] and the references therein. Also, the computations are valid provided that all appearing functions are at least once differentiable. This is in general not the case for conservation laws.

## 2 IMEX-Runge-Kutta Discretization

Numerical discretization of the relaxation system using higher order temporal discretizations combined with higher order spatial discretization has been investigated in several recent publications as for example [22, 27]. We apply so called implicit-explicit Runge-Kutta methods [26, 27, 3] as temporal discretization (IMEX RK). Here, the explicit integration is used for the linear hyperbolic transport part and an implicit method is applied to the stiff source term. Implicit-explicit Runge-Kutta method have been studied in the context of control problems for example in [4, 19]. Define

$$ \mathbf{y} = (u, v)^T, \quad g(\mathbf{y}) = (v, a^2 u)^T \quad \text{and} \quad r(\mathbf{y}) := (0, -(v - f(u)))^T $$
then (2) becomes

$$ \mathbf{y}_t + g(\mathbf{y})_x = \frac{1}{\varepsilon} r(\mathbf{y}), \quad \text{and} \quad \mathbf{y}(0, x) = (u^0, f(u^0))^T(x) $$

Applying a suitable discretization $D_x$ of the spatial derivative yields the semi-discrete state equations

$$ \mathbf{y}' = -D_x g(\mathbf{y}) + \frac{1}{\epsilon} r(\mathbf{y}), \quad \mathbf{y}(0) = \mathbf{y}^0. \qquad (4) $$

*Remark 1.* Spatial discretizations for the linear transport part are well-known. The simplest possible is a first-order Upwind method:

$$ \frac{\partial}{\partial t} \mathbf{y}_j = -\frac{1}{\Delta x} \begin{pmatrix} 0 & 1 \\ a^2 & 0 \end{pmatrix} (\mathbf{y}_{j+1/2} - \mathbf{y}_{j-1/2}) + \frac{1}{\epsilon} r(\mathbf{y}_j), $$

where $\mathbf{y}_{j+1/2}$ is obtained by applying the first-order upwind method to characteristic variables $v \pm au$. Higher order MUSCL schemes, WENO schemes or central schemes have also been studied in this context.

The resulting semi-discrete optimal control problem is then given by:

$$ \begin{array}{ll} \text{minimize} & j(\mathbf{y}(T), \mathbf{y}^0) \\ \text{subject to} & \mathbf{y}' = -D_x g(\mathbf{y}) + \frac{1}{\epsilon} r(\mathbf{y}), \quad \mathbf{y}(0) = \mathbf{y}^0. \quad t \in [0, T] \end{array} \qquad (5) $$

In the context of relaxation schemes the semi-discrete problem is seen as a time-integration problem with stiff source which is discretized by an IMEX RK methods. For the numerical discretization we therefore consider the previous problem as an optimal control problem involving ordinary differential equations. Literature concerning the numerical analysis of Runge-Kutta methods for the optimality system of (5) have been studied in [17, 7, 24]. In [7, 17] partitioned Runge-Kutta methods for the optimality system are obtained using the *discretize-then-optimize* approach. The derived partitioned Runge-Kutta methods have been analysed with regard to symplecticity and order of convergence. In [19], Herty and Schleper, moreover, analysed the associated adjoint imex Runge-Kutta method that one obtains if an explicit method is applied to $D_x g(y)$ and a (diagonally) implicit method to $\frac{1}{\epsilon}r(y)$. In the following, we will analyse general partitioned Runge-Kutta methods using IMEX RK methods. More details can be found in [20]. Therein, the following IMEX Runge-Kutta discretization of (4) is studied.

$$ \begin{array}{ll} Y_n^{(i)} = y_n + h \sum_{j=1}^{i-1} \tilde{a}_{ij} D_x g(Y_n^{(j)}) + h \sum_{j=1}^{i} a_{ij} \frac{1}{\epsilon} r(Y_n^{(j)}) & i = 1, \dots, s \\ y_{n+1} = y_n + h \sum_{i=1}^{s} \tilde{\omega}_i D_x g(Y_n^{(i)}) + h \sum_{i=1}^{s} \omega_i \frac{1}{\epsilon} r(Y_n^{(i)}), & n = 0, 1, 2, \dots \end{array} \qquad (6) $$
A nonlinear variable transformation and two intermediate states $\tilde{K}_n^{(i)}$ and $K_n^{(i)}$ give the equivalent system

$$ 
\begin{aligned}
\tilde{K}_n^{(i)} &= D_x g \left( y_n + h \sum_{j=1}^s \tilde{a}_{ij} \tilde{K}_n^{(j)} + h \sum_{j=1}^s a_{ij} K_n^{(j)} \right) & i = 1, \dots, s \\
K_n^{(i)} &= \frac{1}{\epsilon} r \left( y_n + h \sum_{j=1}^s \tilde{a}_{ij} \tilde{K}_n^{(j)} + h \sum_{j=1}^s a_{ij} K_n^{(j)} \right) & i = 1, \dots, s \\
y_{n+1} &= y_n + h \sum_{i=1}^s \tilde{\omega}_i \tilde{K}_n^{(i)} + h \sum_{i=1}^s \omega_i K_n^{(i)}, & n = 0, 1, 2, \dots
\end{aligned}
\quad (7) $$

The associated optimality systems for the two previous optimization problems then coincide and we refer to [20] for more details. It is proven that the adjoint schemes are equivalent to

$$ \tilde{P}^{(i)} = p_n - h \sum_{j=1}^{s} \tilde{\alpha}_{ij} g'(Y_n^{(j)})^T \bar{D}_x \tilde{P}^{(j)} - h \sum_{j=1}^{s} \alpha_{ij} \frac{1}{\epsilon} r'(Y_n^{(j)})^T P^{(j)} \quad i = 1, \dots, s $$

$$ P^{(i)} = p_n - h \sum_{j=1}^{s} \tilde{\beta}_{ij} g'(Y_n^{(j)})^T \bar{D}_x \tilde{P}^{(j)} - h \sum_{j=1}^{s} \beta_{ij} \frac{1}{\epsilon} r'(Y_n^{(j)})^T P^{(j)} \quad i = 1, \dots, s \quad (8) $$

$$ p_{n+1} = p_n - h \sum_{i=1}^{s} \tilde{\omega}_i g'(Y_n^{(i)})^T \bar{D}_x \tilde{P}^{(i)} - h \sum_{i=1}^{s} \omega_i \frac{1}{\epsilon} r'(Y_n^{(i)})^T P^{(i)} \quad n = 0, 1, \dots, N-1 $$

Here, the coefficients of the Runge-Kutta method $\tilde{\alpha}_{ij}$, $\alpha_{ij}$, $\tilde{\beta}_{ij}$ and $\beta_{ij}$ are given by

$$ \tilde{\alpha}_{ij} := \tilde{\omega}_j - \frac{\tilde{\omega}_j}{\tilde{\omega}_i} \tilde{a}_{ji}, \quad \alpha_{ij} := \omega_j - \frac{\omega_j}{\tilde{\omega}_i} \tilde{a}_{ji}, \quad \tilde{\beta}_{ij} := \tilde{\omega}_j - \frac{\tilde{\omega}_j}{\omega_i} a_{ji}, \quad \beta_{ij} := \omega_j - \frac{\omega_j}{\omega_i} a_{ji}. $$

## 2.1 Properties of Discrete IMEX-RK Optimality System

For the resulting scheme (6), (8) order conditions can be stated [20]. To this end we add a suitable equation for $\tilde{p}$ to the previous system.

$$ \tilde{p}_{n+1} = \tilde{p}_n - h \sum_{i=1}^{s} \tilde{\omega}_i f_y (Y_n^{(i)})^T \tilde{P}^{(i)} - h \sum_{i=1}^{s} \omega_i g_y (Y_n^{(i)})^T P^{(i)}. \quad (9) $$

The full method therefore is a standard additive Runge-Kutta scheme for

$$ 
\begin{aligned}
\mathbf{y}' &= -D_x g(\mathbf{y}) + \frac{1}{\epsilon} r(\mathbf{y}) \\
\tilde{\mathbf{p}}' &= g'(\mathbf{y})^T D_x \tilde{\mathbf{p}} + \frac{1}{\epsilon} r'(\mathbf{y})^T \mathbf{p} \\
\mathbf{p}' &= g'(\mathbf{y})^T D_x \tilde{\mathbf{p}} + \frac{1}{\epsilon} r'(\mathbf{y})^T \mathbf{p}
\end{aligned}
$$

If we define

$$ c_i := \sum_{j=1}^{s} a_{ij}, \quad \text{and} \quad \tilde{c}_i := \sum_{j=1}^{s} \tilde{a}_{ij}, $$
$$ \gamma_i := \sum_{j=1}^{s} \alpha_{ij}, \quad \text{and} \quad \tilde{\gamma}_i := \sum_{j=1}^{s} \tilde{\alpha}_{ij}, $$

$$ \delta_i := \sum_{j=1}^{s} \beta_{ij}, \quad \text{and} \quad \tilde{\delta}_i := \sum_{j=1}^{s} \tilde{\beta}_{ij} $$

then (1) holds true.

**Theorem 1.** Consider the Runge-Kutta scheme (6), (8), (9). This scheme is of

- **First-Order**: if (SRK1) is of first order

- **Second-Order**: if (SRK1) is of second order

- **Third-Order**: if (SRK1) is of third order and either

$$ \sum_{i=1}^{s} \omega_i \gamma_i^2 = \frac{1}{3}, \quad \sum_{i=1}^{s} \omega_i \tilde{\gamma}_i^2 = \frac{1}{3}, \quad \sum_{i=1}^{s} \omega_i \gamma_i \tilde{\gamma}_i = \frac{1}{3}, $$

are satisfied or if

$$ \sum_{i=1}^{s} \omega_i a_{ij} \gamma_i = \frac{1}{6}, \quad \sum_{i=1}^{s} \omega_i \tilde{a}_{ij} \tilde{\gamma}_i = \frac{1}{6} $$

and if

$$ \sum_{i=1}^{s} \omega_i a_{ij} \tilde{\gamma}_i = \frac{1}{6} \quad \text{or} \quad \sum_{i=1}^{s} \omega_i \tilde{a}_{ij} \gamma_i = \frac{1}{6} $$

are satisfied.

Note that the system (6) and (8) is not completely coupled, since the forward scheme (6) is solved independently of the adjoint scheme (8). General order conditions can be found e.g. in [23]. The proof of Theorem 1 and together with more details are discussed in [20].

# 3 Numerical Results

## 3.1 Scalar Example

As a simple example, we use a tracking type functional $J(u)$ together with Burgers' equation

$$ u_t + \left( \frac{u^2}{2} \right)_x = 0, $$

and the desired state $u_d$ at final time $T = 2.0$, that belongs to the initial condition $u_d(0, x) = \frac{1}{2} + \sin(x)$ and we start the optimization with the initial data
$u^s(0, x) \equiv 0.5$. Moreover, the spatial interval is given by $x \in [0, 2\pi]$, As discretization of the objective functional, we use

$$J(u(\cdot, T), u_0, u_d) = \frac{\Delta x}{2} \sum_{i=1}^{K} \|u_i - u_{di}\|^2.$$

Moreover, the discrete gradient of the reduced cost functional is given by

$$\nabla_{u_{0,i}} \tilde{J} = p_{0,i} + (Df(u_0)^T \mathbf{q}_0)_i.$$

In order to solve the optimal control problem, we apply a steepest descent method (with respect to the reduced cost functional) with fixed stepsize $0 < \alpha < 1$, i.e. we set $u_0^{k+1} = u_0^k + \alpha \nabla_{u_{0,i}} \tilde{J}$. As stopping criterion for the optimization process we test $|\tilde{J}(u_0, u_d)| < tol$ where $tol = 1E-2$ denotes a predefined stopping tolerance. We observe grid independence in the case where $u$ and $u_0$ are differentiable in space and time.

As first-order scheme, we test the Implicit-Explicit Euler scheme

$$\begin{aligned} u_i^* &= u_i^n \\ v_i^* &= v_i^n - \frac{\Delta t}{\epsilon} (v_i^* - f(u_i^*)) \\ u_i^{n+1} &= u_i^* - \Delta t D_x v_i^* \\ v_i^{n+1} &= v_i^* - \Delta t a^2 D_x v_i^* \end{aligned}$$

for the forward, as well as for the backward

$$\begin{aligned} q_i^* &= q_i^{n+1} - \Delta t D_x^* p_i^{n+1} \\ p_i^* &= p_i^{n+1} - \Delta t a^2 D_x^* q_i^{n+1} \\ q_i^n &= q_i^* - \frac{\Delta t}{\epsilon} q_i^n \\ p_i^n &= p_i^* + \frac{\Delta t}{\epsilon} q_i^n f'(u_i^n) \end{aligned}$$

The spatial gridsize is chosen to be $N_x = 300$, whereas the time discretization is done according to the CFL condition with constant $c_{CFL} = 0.5$.

<figure><img src="image_2.png" /></figure>
<table><thead><tr><td>N</td><td>Nr. of It.</td><td>CPU time (in sec.)</td></tr></thead><tbody><tr><td>100</td><td>44</td><td>1.795572e+01</td></tr><tr><td>150</td><td>43</td><td>3.768984e+01</td></tr><tr><td>200</td><td>42</td><td>6.380441e+01</td></tr><tr><td>300</td><td>41</td><td>1.491838e+02</td></tr></tbody></table>

# 4 Summary

We briefly discussed a class of numerical methods applied to an optimal control problem for scalar, hyperbolic partial differential equations. Order conditions for the temporal numerical discretization in the case of differentiable functions have been stated. Future work includes the analysis of additional properties of the derived numerical discretizations as for example strong stability and asymptotic preservation properties.

## Acknowledgments.

This work has been supported by DFG HE5386/7-1, HE5386/8-1 and by DAAD 50727872, 50756459 and 54365630. We also acknowledge the support of Ate- neo Italo-Tedesco (AIT) under the Vigoni project 2010-2012 'Adjoint implicit- -explicit methods for the numerical solution to optimization problems'.

## Appendix.

The discrete adjoint equations that correspond to the discrete optimization problem associated with (7) are

$$ \tilde{\xi}_n^{(i)} = h \tilde{\omega}_i p_{n+1} + h \sum_{j=1}^{s} \tilde{a}_{ji} g'(Y_n^{(j)})^T \bar{D}_x \tilde{\xi}_n^{(j)} + h \sum_{j=1}^{s} \tilde{a}_{ji} \frac{1}{\epsilon} r'(Y_n^{(j)})^T \xi_n^{(j)} $$

$$ \xi_n^{(i)} = h \omega_i p_{n+1} + h \sum_{j=1}^{s} a_{ji} g'(Y_n^{(j)})^T \bar{D}_x \tilde{\xi}_n^{(j)} + h \sum_{j=1}^{s} a_{ji} \frac{1}{\epsilon} r'(Y_n^{(j)})^T \xi_n^{(j)} $$

$$ p_n = p_{n+1} + \sum_{i=1}^{s} g'(Y_n^{(i)})^T \bar{D}_x \tilde{\xi}_n^{(i)} + \sum_{i=1}^{s} \frac{1}{\epsilon} r'(Y_n^{(i)})^T \xi_n^{(i)} $$

$$ p_N = j'(y_N, y^0). $$

Moreover, the variable transformation that is needed to obtain (8) is given by

$$ \tilde{P}_n^{(i)} := \frac{\tilde{\xi}_n^{(i)}}{h \tilde{\omega}_i} \quad \text{and} \quad P_n^{(i)} := \frac{\xi_n^{(i)}}{h \omega_i} \quad (i = 1, \dots, s; \ n = 0, \dots, N-1). $$
On the other hand, using (6) the associated discrete adjoint equations are

$$ \zeta_n^{(i)} = h \left( \tilde{\omega}_i f_y(Y_n^{(i)}) + \omega_i g_y(Y_n^{(i)}) \right)^T p_{n+1} + \sum_{j=1}^{s} \tilde{a}_{ji} f_y(Y_n^{(i)})^T \zeta^{(j)} \\ + h \sum_{j=1}^{s} a_{ji} g_y(Y_n^{(i)})^T \zeta^{(j)} \qquad i = 1, \dots, s $$

$$ p_n = p_{n+1} + \sum_{i=1}^{s} \zeta_n^{(i)}, \qquad i = 1, \dots, N-1, \qquad p_N = j'(y_N) $$

which can be transformed into the scheme (8) using the variable transformation

$$ \tilde{P}_n^{(i)} := p_{n+1} + \sum_{j=1}^{s} \frac{\tilde{a}_{ji}}{\tilde{\omega}_i} \zeta_n^{(j)} \quad \text{and} \quad P_n^{(i)} := p_{n+1} + \sum_{j=1}^{s} \frac{a_{ji}}{\omega_i} \zeta_n^{(j)}. $$

## References

1. D. Aregba-Driollet and R. Natalini, Discrete kinetic schemes for multidimensional systems of conservation laws, SIAM J. Numer. Anal., 37 (2000), pp. 1973–2004 (electronic).

2. Areba-Driollet, D., Natalini, R.: Convergence of relaxation schemes for conservation laws. Applicable Analysis 61(1), 163–193 (1996).

3. U. Ascher, S. Ruuth, and R. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics, 25 (1997), 151–167

4. M. K. BANDA AND M. HERTY, Adjoint imex-based schemes for control problems governed by hyperbolic conservation laws, Computational Optimization and Applications, (2010).

5. Banda, M.K., Seaïd, M.: Higher-order relaxation schemes for hyperbolic systems of conservation laws. J. Numer. Math. 13(3), 171–196 (2005).

6. S. BIANCHINI, On the shift differentiability of the flow generated by a hyperbolic system of conservation laws, Discrete Contin. Dynam. Systems, 6 (2000), pp. 329–350.

7. J. F. Bonnans and J. Laurent-Varin, Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control, Numerische Mathematik, 103 (2006), 1–10.

8. S. Boscarino, L. Pareschi and G. Russo, Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit, Preprint, (2011).

9. A. BRESSAN AND G. GUERRA, Shift-differentiability of the flow generated by a conservation law, Discrete Contin. Dynam. Systems, 3 (1997), pp. 35–58.

10. A. BRESSAN AND M. LEWICKA, Shift differentials of maps in BV spaces, in Non-linear theory of generalized functions (Vienna, 1997), vol. 401 of Chapman & Hall/CRC Res. Notes Math., Chapman & Hall/CRC, Boca Raton, FL, 1999, pp. 47–61.

11. A. BRESSAN AND A. MARSON, A variational calculus for discontinuous solutions to conservation laws, Communications Partial Differential Equations, 20 (1995), pp. 1491–1552.
12. A. BRESSAN AND W. SHEN, *Optimality conditions for solutions to hyperbolic balance laws*, Control methods in PDE-dynamical systems, Contemp. Math., 426 (2007), pp. 129–152.

13. G. Dimarco and L. Pareschi, *Asymptotic-Preserving IMEX Runge-Kutta methods for nonlinear kinetic equations*, preprint, (2011)

14. Gottlieb, S., Shu, C.W., Tadmor, E.: Strong stability preserving high-order time discretization methods. SIAM rev. **43**, 89–112 (2001)

15. A. L. Dontchev and W. W. Hager, *The Euler approximation in state constrained optimal control* Math. Comp., **70** (2001), 173–203

16. A. L. Dontchev and W. W. Hager and V. M. Veliov, *Second-order Runge-Kutta approximations in control constrained optimal control* SIAM J. Numer. Anal., **38** (2000), 202–226

17. W. W. Hager, *Runge-Kutta methods in optimal control and the transformed adjoint system*, Numerische Mathematik, **87** (2000), 247–282.

18. E. Hairer, S. P. Nørsett and G. Wanner, *Solving Ordinary Differential Equations, Part I, Nonstiff Problems*, Springer Series in Computational Mathematics, second edition (1993)

19. M. Herty and V. Schleper, *Time discretizations for numerical optimization of hyperbolic problems*, Applied Mathematics and Computation (2011)

20. M. Herty and L. Pareschi and S. Steffensen, *Implicit-Explicit Runge-Kutta schemes for numerical discretization of optimal control problems*, preprint available at University Ferrara, (2012)

21. M. Giles and S. Ulbrich, *Convergence of linearized and adjoint approximations for discontinuous solutions of conservation laws. Part 2: Adjoint approximations and extensions* SIAM J. Numer. Anal., **48**, (2010), pp. 905–921.

22. S. Jin and Z. P. Xin, *The relaxation schemes for systems of conservation laws in arbitrary space dimensions* Comm. Pure Appl. Math., **48** (1995), pp. 235–276.

23. C. A. Kennedy and M. H. Carpenter, *Additive Runge-Kutta schemes for convection-diffusion-reaction equations*, Appl. Num. Math., **44** (2003), 139-181.

24. J. Lang and J. Verwer, *W-Methods in optimal control* Preprint 2011, TU Darmstadt

25. Natalini, R., Terracina, A.: Convergence of a relaxation approximation to a boundary value problem for conservation laws. Comm. Partial Differential Equations **26**(7-8), 1235-1252 (2001).

26. L. Pareschi and G. Russo, *Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation*, J. Sci. Comput., **25** (2005), 129-155.

27. L. Pareschi and G. Russo, *Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations* Recent Trends in Numerical Analysis, Edited by L.Brugnano and D.Trigiante, **3** (2000), pp. 269–289.

28. S. ULBRICH, *Optimal Control of Nonlinear Hyperbolic Conservation Laws with Source Terms*, Technische Universitaet Muenchen, 2001.

29. S. ULBRICH, *On the superlinear local convergence of a filer-sqp method*, Technical Report, (2002).

30. S. ULBRICH, *Adjoint-based derivative computations for the optimal control of discontinuous solutions of hyperbolic conservation laws.*, Syst. Control Lett., 48 (2003), pp. 313–328.

31. A. Walther *Automatic differentiation of explicit Runge-Kutta methods for optimal control* J. Comp. Opt. Appl., **36** (2007), 83–108

32. C. CASTRO, F. PALACIOS, AND E. ZUAZUA, *An alternating descent method for the optimal control of the inviscid Burgers equation in the presence of shocks*, Math. Models Methods Appl. Sci., 18 (2008), pp. 369–416.