Optimal control of the temperature in a catalytic converter

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Optimal control of the temperature in a catalytic converter Ingenuin Gasser ∗ , Martin Rybicki, Winnifried Wollner Universität Hamburg, Fachbereich Mathematik, Bundesstraße 55, 20146 Hamburg, Germany

article

info

Article history: Received 26 August 2013 Received in revised form 6 February 2014 Accepted 8 February 2014 Keywords: Catalytic converter Combustion Exhaust pipe Gas dynamics network Optimization Small Mach number

abstract This paper is concerned with the optimization of the dynamics of the gas flow in an exhaust pipe. As a prototypical question the start up heating of the catalytic converter is considered. Reaching the optimal temperature in the converter rapidly competes against the costs, i.e., releasing too much unburnt fuel in the exhaust gases. The underlying model is a one dimensional small Mach number asymptotic gas dynamic model for a mixture of burnt and unburnt gases, formulated on a network consisting of the various pieces of the exhaust tube. For the optimization, adjoint calculus is applied on the continuous level. The resulting system is then discretized. Finally, numerical examples show the validity of this approach. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction The control and the reduction of the emissions caused by vehicles is, and has been, an import issue over the last 5 decades. The first restrictions were introduced by the government of California (USA) in the early 1960s. In 1970 the European Community passed first laws regarding exhaust gas pollution. Today, there is the Euro 5 standard and the next Euro 6 standard will be compulsory in 2014 in Europe. Many other countries nowadays use or introduce similar restrictions. For the reduction of the concentration of CO, NOx and Cx Hy in the exhaust gas, there is a classical technical solution, namely in most cases two catalytic converters are installed in the exhaust pipe system. The function of the catalytic converters depends strongly on the temperature in the converters. There is a lower limit (about 600 °C) for a good function and an upper limit to avoid damages. In particular, right after engine start there is a critical time interval where the temperature in the converters is not high enough. A method of heating after the engine start is the combustion of unburnt gas in the catalytic converters. Modern exhaust systems can control the ratio of oxygen and fuel in the combustion chamber of the engine. By choosing a ratio with more fuel and less oxygen some unburnt fuel flows to the catalytic converters where it can be used for an exothermic reaction. Clearly, there is a competition between reaching fast the optimal converter temperature and using very little unburnt fuel in the exhaust gas. This is the main issue in this paper. We show, how to compute an optimal inflow distribution of unburnt gas (into the exhaust tube) with respect to a cost function that will be presented below. To face this problem, one has to model, appropriately, the gas dynamics in the exhaust pipe and the heat dynamics in the converters. In addition, we do not only need an appropriate model but also a model which allows fast direct simulations in order to be able to apply optimization tools, which typically need several simulations during the optimization. In this sense, this paper can be seen as a prototypical example for optimization in a compressible gas dynamic setting with the additional requirement of very fast simulation.

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (I. Gasser), [email protected] (M. Rybicki), [email protected] (W. Wollner). http://dx.doi.org/10.1016/j.camwa.2014.02.006 0898-1221/© 2014 Elsevier Ltd. All rights reserved.

2

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

There is not much mathematical literature on the modeling of gas dynamics in an exhaust pipe in the literature. We mention [1–3] on studies regarding the temperature in the converter based on more complex chemical models, with assumed homogeneous gas dynamics. Furthermore, we mention fully compressible fluid dynamic models for the exhaust tube, where the aim is not primarily the determination of the converter temperature but the computation of sound waves (see, e.g., [4,5]) or mechanical properties of the exhaust tube (see, e.g., [6,7]). A model which promises to be appropriate was developed recently in [8,9]. We shortly summarize the main features of this model:

• The model is one dimensional in space. The spacial axis is along the exhaust tube and all the values are mean values over the (non constant) cross sections.

• The basis is a fully compressible multi component gas dynamic model including surface friction and heat transfer through the surface, i.e., the exhaust tube walls.

• In this application the flow velocities are small (compared to the speed of sound), i.e., the Mach number is small. This fact is used to build a small Mach number asymptotic model in order to rule out sound waves. This reduces the simulation times significantly. • We use two gas components to model the underlying complex chemistry, namely burnt and unburnt gas. However, the crucial quantity for the heating in the catalytic converters is the total heat release of the chemical reaction which can easily be included in such a simple two component reaction model. • The exhaust tube is modeled as a network of tubes with constant cross section. Neither for the pre-asymptotic hyperbolic model – a system of nonlinear hyperbolic balance laws – nor for the asymptotic (hyperbolic–elliptic-like) model a well-posedness theory is available. This is a common problem in many complex, multi-physics, systems. However, assuming existence of a unique solution and differentiability with respect to the inflow of unburnt gas, we will derive a system of differential equations yielding the continuous adjoint to the calculated flow. These will be used to find an optimal inflow profile using a projected gradient method. We remark that, even in the pre-asymptotic regime, i.e., for hyperbolic systems – and even for scalar nonlinear hyperbolic equations – questions of uniqueness, differentiability, and optimality conditions are a field of active research. To give an impression of the different questions under consideration see, e.g., [10] for well-posedness of hyperbolic problems on networks, [11] for sensitivity analysis of hyperbolic conservation laws, [12,13] for optimization problems with hyperbolic equations on networks, [14] for control of Burgers equation with vanishing viscosity, or [15] for the control of unsteady compressible fluids. The asymptotic model, under consideration, allows direct simulations almost in real time on a standard laptop computer. The model derived in [8] was built on the basis of models proposed in [16,17]. In order to validate the accuracy of the asymptotic model, we compared its numerical solutions to the ones of the fully hyperbolic problem and found a good agreement (see [8]). Although it seems that we have drastically reduced the complexity of the model both from the application and mathematical point of view, we keep an appropriate, valid, and highly complex problem. From the application point of view, we deal with compressible low Mach number gas dynamics. Such problems arise in many application; ranging from astrophysics, meteorology to many engineering applications, e.g., gas pipeline networks [18]. From the mathematical point of view, we still have a multi component nonlinear system of PDE’s on a network. To do optimization on such a model is still a highly challenging topic. We mention only a few papers where similarly complex issues where studies. In [19] optimization of non reactive gas flow in a pipeline is considered, whereas the optimal control for traffic flow is discussed in [12]. In [20] the optimal control of glass cooling is studied. The paper is organized as follows. In Section 2, we introduce the model, whose main parts were derived in [8]. In Section 3, we define the optimization problem. In Section 4, we derive the optimality system. Section 5 is devoted to the discretization. Finally, in Section 6, we present numerical examples. 2. Model The following dimensionless model describes the gas dynamics in an exhaust pipe.

(Aρ)t + (Aρ u)x = 0, 

Aρ T + ε

R cv

Aρ

u

u|u| 1 − Cc Aχ ρ u, (Aρ u)t + (Aρ u2 )x + Apx = −Cf ρ ε 2  3 R u R + Aρ uT + ε Aρ + Aup = −h(T − TWall ) + q0 Aχ ρ zK (T ),

 2

2



t

cv

2

cv

(1)

x

(Aρ z )t + (Aρ uz )x = −Aχ ρ zK (T ), p = ρT , with the unknowns ρ, u, p, T , z as density, velocity, pressure, temperature, and ratio of unburnt gas, respectively. The terms on the left hand side originate from the one-dimensional Euler equations of gas dynamic in a pipe with variable cross section

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

3

A = A(x). By the terms on the right hand side, we describe the main physical effects which influence the dynamics of the gas. The first equation is the conservation of mass. On the right hand side of the second equation (momentum balance), we first have a wall friction (wall friction coefficient Cf ) and, secondly, a local friction (honeycomb structure friction coefficient Cc ) due to the honeycomb structure of the catalyst (locality is denoted by the indicator function χ ). In the third equation (energy balance), we have, on the one hand, energy loss due to heat transfer with the wall (heat transfer coefficient h, wall temperature TWall ), and, on the other hand, an energy gain due to the local exothermic reaction in the catalytic converter (heat release coefficient q0 ). The parameters cv and R denote the specific heat at constant volume and the ideal gas constant, respectively. The fourth equation (reaction balance) describes the dynamics of the ratio of unburnt gas z in the pipe. The temperature dependent reaction rate K (T ) is modeled by Arrhenius’ law. Lastly, we have the ideal gas law as a closing relation, in the fifth equation. The parameter ε := γ M 2 equals the product of the adiabatic exponent and the square of the Mach number and is known to be small (see [8]). The numerical simulation of the fully compressible model (1) is (although it is ‘‘only’’ 1D) costly. In [8] it is shown that for many purposes a simplified model on the basis of (1) can be used to simulate gas dynamics more efficiently. There are two major steps in the derivation of this simplified model: 1. Rather than introducing artificial intervals for a smooth change of the cross section A, one divides the exhaust pipe into pipes with constant cross sections and treat the whole system as a network. This has two major advantages: (a) Strong changes in A on small x-intervals are avoided. Thus, larger step sizes in space (and therefore, also in time) are possible. (b) The coupling conditions at the junctions allow to include models for pressure losses caused by the pipe’s geometry. 2. The flow in the exhaust pipe is of low Mach number (M) and (in this application) one is not interested in the propagation of sound waves, but mainly in the temperature in the catalyst. Thus, one can simplify the model by a small Mach number asymptotics (ε → 0). This leads to the simplified model (2), which is computationally much cheaper than the fully hyperbolic model (1) due to much larger step sizes in time given by the CFL condition. See [8] for the detailed derivation and comparison with the original model (1). The final version of this asymptotic model is the following:

ρti + (v i + Q i )ρxi = −qi ρ i , zti + (v i + Q i )zxi = −χ i z i K (T i ),

(2)

vti = Φ i /Ri (0), p0 = ρ i T i ,

with unknown density ρ i = ρ i (x, t ), time dependent velocity v i = v i (t ), ratio of unburnt gas z i = z i (x, t ), and temperature for the ith pipe T i = T i (x, t ). Before we give the definitions for the other quantities, let us first discuss the structure of the model. We see that instead of four coupled PDEs in (1) we have now only a system of two PDEs and one Integro-ODE (which is made precise below, see (9)). The super-index i ∈ {1, . . . , nP } denotes the pipe number, where nP is the amount of pipes in total. Since we have to distinguish between pipes that do and do not have a catalytic converter, we multiply terms, occurring only in pipes that have a catalyst, with the indicator function χ i , where

χ = i



1 0

pipe i has a catalyst, otherwise.

(3)

Furthermore, we have to define some terms in the Eqs. (2). All values for the parameters, can be found in Appendix A. 1. The velocity u is now split into the space independent velocity component v and the aggregated energy balance Q . u (x, t ) = v (t ) + Q [ρ , z ](x, t ) = v (t ) + i

i

i

i

i

i

x



qi [ρ i , z i ](y, t )dy,

(4)

0

qi [ρ i , z i ](x, t ) :=

1 

γ p0

 i −h(T i (x, t ) − TWall (x, t )) + χ i q0 ρ i (x, t )z i (x, t )K (T i (x, t )) ,

(5)

with p0 given as zero order pressure in the small Mach number asymptotics and

 1 i T (x, t ) + Tout , (6) 2 where Tout is the scaled outside temperature. 2. The reaction rate K , which depends on the temperature/density of the gas mixture, is defined in the following way: i TWall (x, t ) :=

K (T i ) :=

K˜ 0 xr ur

 exp −

T˜ + Tr T i

 =

K˜ 0 xr ur

 exp −

T˜ + Tr p0

 ρi

=: Kρ (ρ i )

(7)

4

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

with the activation temperature T˜ + and reference values xr , ur , Tr for the spatial variable, velocity and temperature, respectively (see Table A.6). 3. The right hand side of the third equation is defined (with qi (x, t ) = qi [ρ i , z i ](x, t ) and Q i (x, t ) = Q i [ρ i , z i ](x, t )) as follows:

(Ri (x))(t ) :=

Li



ρ i (y, t )dy,

(8)

x

Φ [ρ , z , v ](t ) := i

i

i

i

pi1l

(t ) −

pi1r

(t ) −

Li



ρ (x, t ) (x, t )dx − i

Qti

Li



0

ρ i (x, t )(v i (t ) 0

+ Q i (x, t ))qi (x, t )dx − Cf



Li

ρ i ( x, t )

(v(t )i + Q i (x, t ))|v(t )i + Q i (x, t )|

0

− χ i Cc

2

dx

Li



ρ i (x, t )(v i (t ) + Q i (x, t )) dx.

(9)

0

The index set Icc ⊂ {1, . . . , nP } contains all pipes that have a catalytic converter. Since in a catalytic converter, the temperature is not considered to be known, we will have to augment this model. In order to describe the heating of the catalytic converter, in addition to the model in [8], we need to introduce an equation which models the evolution of the temperature of the catalyst. We model the temperature of the catalytic converter Tci as a i result of the heat exchange among catalyst temperature Tci and the average gas temperature TGas in the ith pipe. i TGas (t ) :=

1

Li



Li

T i (x, t )dx.

(10)

0

We choose the following ODE, for all pipes i ∈ Icc , to model this relation. i (Tci (t ))t = −hc (Tci (t ) − TGas (t )),

(11)

i

and we substitute q from (5) by qi [ρ i , z i , Tci ] :=

1 

γ p0

  i −h(T i − TWall ) + χ i −hc (T i (t ) − Tci (t )) + q0 ρ i z i K (T i ) ,

(12)

including the heat exchange with the catalytic converter in the energy balance qi . This has an impact on Q i and Φ i , i.e., Q i = Q i [ρ i , z i , Tci ] and Φ i = Φ i [ρ i , z i , v i , Tci ]. 2.1. The summarized model To summarize, we consider the following model to describe the physics in an exhaust pipe:

ρti + (v i + Q i )ρxi = −qi ρ i , zti + (v i + Q i )zxi = −χ i z i K (T i ),

(13)

vti = Φ i /Ri (0), j (Tcj )t = −hc (Tcj − TGas ),

for all pipes i = 1, . . . , nP , j ∈ Icc and (x, t ) ∈ (0, Li ) × (0, tend ), with closing relation

ρ i (x, t )T i (x, t ) = p0 ,

(14)

for all pipes i = 1, . . . , nP and (x, t ) ∈ [0, Li ] × [0, tend ], initial conditions

ρ i (x, 0) = ρici (x),

z i (x, 0) = zici (x),

v i (0) = vici ,

Tcj (0) = Tcj ic ,

(15)

for all pipes i = 1, . . . , nP , j ∈ Icc and x ∈ [0, Li ], boundary conditions for all times t ∈ [0, tend ]

ρ 1 (0, t ) = ρl (t ),

z 1 (0, t ) = zl (t ),

(16)

and coupling conditions for all pipes i = 2, . . . , nP

ρ i−1 (Li−1 , t ) = ρ i (0, t ),

z i−1 (Li−1 , t ) = z i (0, t ).

(17)

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

5

Remark 1. There are additional physical boundary conditions we need to prescribe, i.e., the pressure value at the exhaust pipe’s very left end very right end. p1 (0, t ) = pl (t ) = p0 + ε p1 l (t ),

pnP (LnP , t ) = pr (t ) = p0 + ε p1 r (t ).

(18)

For the mathematical model (13) these are just parameters (which appear in Φ , see (9)). Therefore it is not mentioned in the model description. Since we also need pressure values at the junctions, we need to prescribe a coupling condition for the pressure value. i

pi−1 (Li−1 , t ) = pi (0, t ) + fext .

(19)

See [8] or [21,22] for the definition of the pressure loss term fext . The numerical computation of all pressure values at each pipe end for each time step is not trivial and explained in detail in [8]. 3. The optimization problem 3.1. What do we want to optimize? As already mentioned in the introduction, the exothermic reaction of the unburnt gas has a major impact on the gas temperature in the catalyst. The desired reactions in the catalyst take place best if a certain temperature is achieved. More precisely, the catalyst is not functioning if the temperature is below a certain threshold (‘‘light of temperature’’), and may be damaged if a certain critical temperature is exceeded. It is thus natural to consider an optimization problem to reach a desired temperature Topt by controlling the inflow of unburnt gas zl . We will neglect the natural constraints induced by avoiding damage in this article. In order to formulate this optimization task, we first need a cost functional. As described above, on the one hand we want to reach the light of temperature in the catalyst, but on the other hand do not want to use too much fuel. We express this by the following cost functional.

J (Tci , zl ) :=

 1 2 i∈I cc

tend

(Tci (t ) − Topt )2 dt + σ

0

tend



zl (t )dt ,

(20)

0

with the additional upper and lower bound for the inflow of unburnt gas zl 0 ≤ zl (t ) ≤ zlmax

∀t ∈ [0, tend ].

(21)

The choice of an L1 -type term for the unburnt gas is physically well motivated by the linear costs associated with the amount of fuel used. However, there is no physical or chemical motivation for our choice of the tracking type term for the temperature. Measuring the costs by a L2 -type term is rather arbitrary and could be replaced by other functionals on the temperature. For the algorithmic realization, we assume that for any non-negative function zl ∈ L1 (0, tend ) there is a unique temperature Tci = T (zl ) ∈ L2 (0, tend ) given by the model (13)–(17). This gives rise to the reduced problem of finding an inflow distribution zl solving min j(zl ) = J (T (zl ), zl ) s.t. 0 ≤ zl ≤ zlmax .

(22)

To solve this numerically, we will apply a projected gradient method in this reduced setting, see Section 5. Before we can do so, we need to calculate the derivatives of the reduced cost functional which will be done in the following section. In the calculations, we follow the optimize-then-discretize approach, meaning that we will derive the necessary optimality conditions for the continuous problem and then discretize the optimality conditions to obtain a discrete problem, see our discussion in Section 5. 4. The optimality system and its derivation Remark 2. Up to now, there is no existence theory for the forward model (13)–(17). For a related model describing fire in tunnels existence of a unique solution is shown in [23]. However, our considered state Eqs. (13) are more complex. In contrast to the tunnel fire model, the heat term q is not a given function. Furthermore, we have two additional equations in the exhaust pipe model, namely the equation for ratio of unburnt gas (z) and the temperature of the catalytic converter (Tc ). Since the extension of the results in [23] to the model under consideration is beyond the scope of this paper, we assume the existence of a unique solution. Differentiability of this solution with respect to the inflow of unburnt gas could then be obtained by an application of the implicit function theorem, provided the linearized equations are solvable in the appropriate spaces. Since the precise setting for the function spaces would be provided by the nonlinear forward problem, we do not consider differentiability for the problem as well, but rather proceed formally in calculating the derivatives.

6

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

For the following computation, we need to assume that the velocity ui = v i + Q i is non-negative for all (x, t ) ∈ [0, Li ] × [0, tend ] in all pipes i = 1, . . . , nP . Then |v i + Q i | = v i + Q i in the term Φ . This assumption is fully in correspondence with the application. Since we want to determine the gradient of our cost functional J by an adjoint-based method, we need to calculate the adjoint. To do so, we formulate a Lagrangian functional. Let W := (w 1 , . . . , w nP ) be the vector of all state variables, i.e., w i := (ρ i , z i , v i , Tci ) and Λ := (λ1 , . . . , λnP ) the vector of all adjoint variables, i.e., λi := (ξρi , ξzi , ξvi , ξTi c , ηρ , ηz , νρi , νzi , νvi ,

νTi c , ζρi , ζzi ). Then the Lagrangian functional is    tend L(W , zl , Λ) = (Tci (t ) − Topt )2 dt + σ

−

nP  

nP  

0 tend

0

i∈Icc tend

ξρi (ρti + (v i + Q i )ρxi + qi ρ i )dxdt

(24)

ξzi (zti + (v i + Q i )zxi + χ i z i K (T i ))dxdt

(25)

Φi v − i R (0)

(26)

Li

 0

tend

 

0 tend

nP   i =1

−

ξv i



i t



i ξTi c ((Tci )t + hc (Tci − TGas ))dt

0

−

nP 

0

Li

νρ (ρ (x, 0) − ρ (x))dx − i

i

νvi (v i (0) − vici ) −

−

tend

0

nP   i=2

(28)

i ic



Li

nP  

νzi (z i (x, 0) − zici (x))dx

(29)

0

νTi c (Tci (0) − Tci ic )

(30)

i∈Icc

nP   i=2

ηz (z 1 (0, t ) − zl (t ))dt

i =1

i=1

−

tend



(27)

0

nP   i=1

dt

ηρ (ρ 1 (0, t ) − ρl (t ))dt −

−

−

(23)

Li



0

i =1

−

tend 0

i =1

−

zl (t )dt

0

0

i∈Icc

tend

tend

  ζρi ρ i−1 (Li−1 , t ) − ρ i (0, t ) dt

(31)

  ζzi z i−1 (Li−1 , t ) − z i (0, t ) dt .

(32)

0

Before starting with the computation, let us have a quick look at the lines of this large term. Line (23) is representing the cost functional. Lines (24)–(27) are the four state equations multiplied with Lagrangian multipliers ξ∗i , integrated over space and time and summed up over all pipes of the network. The last four lines are the boundary conditions multiplied with its Lagrangian multipliers η∗i and integrated over time (28), the initial conditions multiplied with its Lagrangian multipliers ν∗i integrated over space (29)–(30) and the coupling conditions multiplied with its Lagrangian multipliers ζ∗i and integrated over time (31)–(32). We derive the optimality system by computing the first variation with respect to Lagrangian multipliers, state variables, and the control quantity. Then, the stationarity requirement yields the following set of first order optimality conditions

δL = 0 ⇒ constraints or state equations, ∂λij ∂L = 0 ⇒ adjoint or co-state equations, ∂wji

(33)

∂L ∗ (z − zl ) ≥ 0 ∀ zl∗ ∈ [0, zlmax ] ⇒ reduced optimality condition. ∂ zl l 4.1. Derivation of the optimality condition

Let us start with the derivation of the optimality condition. This computation is the easiest, since zl only appears in the terms (23) and (28). Let δ zl be an arbitrarily chosen L1 -function, such that 0 ≤ zl + ϵδ zl ≤ zlmax , for sufficient small ϵ > 0.

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

7

Then the first variation of the Lagrangian L with respect to the control zl in direction δ zl is

∂ L(W , zl , Λ) (δ zl ) = ∂ zl

tend



δ zl (σ + ηz )dt .

0

Given that (33) holds for all feasible test-functions we assert that the inequality in the reduced optimality condition holds pointwise almost everywhere in [0, tend ] due to the fundamental lemma of calculus of variation (see, e.g., [24, Lemma 2.26]). This yields

(σ + ηz )(zl∗ − zl ) ≥ 0 ∀ 0 ≤ zl∗ ≤ zlmax ,

(34)

which in turn is equivalent to

 ≥0 σ + ηz = 0 ≤0

zl = 0, 0 < zl < zlmax , zl = zlmax .

(35)

4.2. Derivation of the constraints or state equations In order to obtain the state equations, we have to compute the first variation of the Lagrangian functional L with respect to the Lagrangian multipliers. The computation works just like in the above case. We then obtain the state system consisting of our differential equations (13) as well as our boundary (16), initial (15) and coupling conditions (17). 4.3. Derivation of the adjoint or co-state equations In this subsection, we will refuse to write super-indices, but compute the first variations for a pipe of length L = 1 with a catalytic converter (χ = 1). The computation of the adjoint equations for a pipe without a catalyst is just a special case (no equation for Tc and χ = 0). Note that, in the case of a single pipe there are no coupling conditions required. We will deal with this problem in Section 4.4. The following part is the most complicated one when it comes to computation of the first variation. This is the case, because all state variables except v appear in the functional Q = Q [ρ, z , Tc ]. We will show the computation of the first variation with respect to the space-independent velocity component v and the density ρ in this section. The remaining terms can be calculated in a similar manner, for details we refer to Appendix B. The first variation of L with respect to the density involves most of the terms in the Lagrangian and includes all techniques that are required for the computation of the variations with respect to the other two variables z and Tc . We start with ∂ L/∂v , since this is the easier one of the two cases. 4.3.1. First variation of L with respect to v The space independent velocity component v appears in the Lagrangian L in the density PDE (24), the ratio of unburnt gas PDE (25) and its own ODE (26) such as its initial condition (30). One has to keep in mind, that the functional Φ (9) has also a v -dependence. So before computing the first variation of L with respect to v , let us have a closer look at the functional Φ.

Φ [ρ, z , v + ϵδv, Tc ] − Φ [ρ, z , v, Tc ] = −ϵδv





1

  ϵ 2 δv 2 ρ q + Cf ρ(v + Q ) + χ Cc ρ dx −

2

0



   ϵδv ρ q + Cf ρ (v + Q ) + + χ Cc ρ dx

1

= −ϵδv 0

2

1



Cf ρ dx. 0

We already isolated the term of order ϵ 2 , since this term is of higher order and will vanish in the limit ϵ → 0.

   tend  1  ∂ L(W , zl , Λ) 1 (δv) = lim − ϵδv(ξρ ρx + ξz zx )dxdt − νv ϵδv  ϵ→0 ϵ ∂v 0 0 t =0  1    tend    ϵδv − ξv ϵδvt + ρ q + Cf ρ(v + Q ) + χ Cc ρ dx + o(ϵ) . R(0) 0 0 Using integration by parts, we shift the time derivative of δv to the co-state. This also leads to an evaluation at the boundaries. tend

 − 0

ξv δvt dt =

tend

 0

 t =t (ξv )t δv dt − ξv δv t =0end .

8

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

The limit ϵ → 0 gives

∂ L(W , zl , Λ) (δv) = ∂v

tend



1

  δv (ξv )t −

ξρ ρx + ξz zx dx

0

0

− ξv

1



1 R(0)

ρ q + Cf ρ(v + Q ) + χ Cc ρ dx



0

t =tend

dt − ξv δv t =0



− νv δv(0).

Now, we assume the derivative of the Lagrangian with respect to v , in the direction δv , to vanish. Since the variation δv is arbitrary, we choose a variation that vanishes at the boundaries, i.e., at t = 0 and t = tend . Then we deduce, by the fundamental theorem of variational calculus, that

−(ξv )t = −ξv S (0) −

1



ξρ ρx + ξz zx dx ,

(36)

0

with S (x) :=

1 R(0)

1



  ρ q + Cf ρ(v + Q ) + Cc ρ dy.

(37)

x

By choosing a variation that vanishes only at one of the end points t = 0, t = tend , we obtain, with the knowledge that (36) holds,

ξv = 0 for t = tend .

(38)

4.3.2. First variation of L with respect to ρ We have to compute the first variation of L with respect to ρ , i.e., ∂ L/∂ρ . Since the density appears in many terms in the Lagrangian, this computation is very long and technical. Therefore, we will divide this task in several subtasks, i.e., computation of the first variation of (a) (b) (c) (d)

the integral (24) with test function ξρ , the integral (25) with test function ξz , the integral (26) with test function ξv , the integral (25) with test function ξTc .

The remaining integrals depend linearly on ρ , hence the first variation is trivial to calculate. In the end, we will put all results of the subtasks back together and finalize with ∂ L/∂ρ . Before we start the computation of the derivative, let us focus on the dependence of q and Q upon the density ρ . We assume ρ ≥ c > 0, which is in correspondence with the physics. q[ρ + ϵδρ, z , Tc ] =

1

γ p0



 −h



p0

ρ + ϵδρ



− TWall − χ hc

p0

ρ + ϵδρ





− Tc + χ q0 (ρ + ϵδρ)zKρ (ρ + ϵδρ) .

Note that, we will use the notation, given in (7), for the reaction rate function. We observe that the density variable and its variation do not appear in a linear way, but as the denominator of a friction and the argument of the reaction rate. In order to subtract the evaluation of q at ρ + ϵδρ , we have to get rid of these terms. We do this by Taylor expansion: 1

1

1

= − ϵδρ 2 + o(ϵ), ρ + ϵδρ ρ ρ Kρ (ρ + ϵδρ) = Kρ (ρ) + ϵδρ Kρ′ (ρ) + o(ϵ). By replacing these terms by their Taylor expansions, we obtain q[ρ + ϵδρ, z , Tc ] = q[ρ, z , Tc ] + ϵδρ



1

γ p0

p0 (h + χ hc )

ρ2

 + χ q0 z (ρ Kρ′ (ρ) + Kρ (ρ)) + o(ϵ).

By definition of Q , this yields Q [ρ + ϵδρ, z , Tc ] − Q [ρ, z , Tc ] =

x



ϵδρ qρ dy + o(ϵ),

(39)

  ϵδρ qρ t dy + o(ϵ)

(40)

0

Qt [ρ + ϵδρ, z , Tc ] − Qt [ρ, z , Tc ] =

x

 0

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

9

where qρ :=

p0 (h + χ hc )



1

γ p0

ρ2

+



χ q0 z (ρ Kρ′ (ρ)

+ Kρ (ρ)) .

(41)

In order to shorten the expressions, we define the following abbreviations for the computation of ∂ L/∂ρ : q := q[ρ, z , Tc ],

qϵ := q[ρ + ϵδρ, z , Tc ],

Q := Q [ρ, z , Tc ],

Qϵ := Q [ρ + ϵδρ, z , Tc ].

(a) The ξρ -integral (24) tend



1



  ξρ (ρ + ϵδρ)t + (v + Qϵ )(ρ + ϵδρ)x + qϵ (ρ + ϵδρ) dxdt

− 0



0 tend



tend



1

  ξρ ρt + (v + Q )ρx + qρ dxdt

+ 0

0



1

=− 0

   x ϵδρ qρ dy dxdt + o(ϵ) =: Iξρ . ξρ ϵδρt + (v + Q )ϵδρx + qϵδρ + ρ qρ ϵδρ + ρx

(42)

0

0

Our aim is to isolate the variation δρ , such that we can apply the fundamental lemma of variational calculus again. Therefore we need to integrate by parts two times (one with respect to the spatial and once with respect to the time component). Furthermore, we need to use the following identity 1



f (x)

x



0



1



g (y)dy dx = 0

g (x)

1



0



f (y)dy dx x

which is a consequence of Fubini’s theorem. With this we manipulate (42), integrate by parts and obtain using (v + Q )x = q

 I ξρ =

tend





tend

0

1

1

   ϵδρ (ξρ )t + (v + Q )(ξρ )x − qρ ρξρ +

0

 ξρ ρx dy dxdt

x

−

ϵδρ(v + Q )ξρ dt

0

x=1

t =tend

1



ϵδρξρ dx

− 0

x=0

+ o(ϵ).

(43)

t =0

(b) The ξz -integral (25). Analog calculations give tend

 tend  1    − ξz zt + (v + Qϵ )zx + χ zKρ (ρ + ϵδρ) dxdt + ξz zt + (v + Q )zx + χ zKρ (ρ) dxdt 0 0 0 0   tend  1   x ′ =− ξz zx ϵδρ qρ dy + χ ϵδρ zKρ (ρ) dxdt + o(ϵ) 



1



0

0

tend





ϵδρ qρ

=− 0

0

1



1



0

x

 ξz zx dy + χ ξz zKρ′ (ρ) dxdt + o(ϵ).

(44)

(c) The ξv -integral (26). This is the most technical part, since ρ appears in this term so often. Let us develop this step by step. First, we have a look at the density integral in the denominator. By Taylor expansion we obtain: 1

1 0

ρ + ϵδρ dx

= 1 0

1



1

ρ dx

− 0

ϵδρ dx 

1

1 0

2 + o(ϵ) =

ρ dx

1 R(0)

−

1 R(0)2

1



ϵδρ dx + o(ϵ). 0

The term Φ evaluated at ρ + ϵδρ , will have the following form:

Φ [ρ + ϵδρ, z , v, Tc ] = Φ [ρ, z , v, Tc ] + ϵ Φρ [δρ; ρ, z , v, Tc ] + o(ϵ). To shorten the integrals, we introduce the following abbreviations:

Φ := Φ [ρ, z , v, Tc ],

Φϵ := Φ [ρ + ϵδρ, z , v, Tc ].

Please note the different meaning of Φϵ and the derivative Φρ . We will compute Φρ in detail later. Let us first have a look on the whole integral containing ξv . tend

 − 0



ξv vt −  1 0

tend

 = 0

ξv



1 R(0)



1

ρ + ϵδρ dx

ϵ Φρ −

Φϵ dt +

1 R(0)2

tend

 0

1

 0

 ξv vt −  1 0

1

ρ dx

 Φ dt

 ϵδρ dx Φ dt + o(ϵ) =: Iξv .

(45)

10

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

The term ϵ Φρ consists of all the terms in which ϵ appears linearly, when evaluating Φ at ρ + ϵδρ . We will now compute those terms. Therefore, we split Φ in several parts.

Φ [ρ, z , v, Tc ] = p1 l − p1 r −

1

 0

 1  1  1 (v + Q )2 ρ Qt dx − ρ(v + Q )qdx − Cf ρ dx − χ Cc ρ(v + Q )dx . 2 0    0    0     

=:Φ1 [ρ,z ,Tc ]

=:Φ2 [ρ,z ,v,Tc ]

=:Φ3 [ρ,z ,v,Tc ]

=:Φ4 [ρ,z ,v,Tc ]

We compute the difference −Φi [ρ + ϵδρ, z , v, Tc ] + Φi [ρ, z , v, Tc ] =: −Φi;ϵ + Φi

− Φ1;ϵ + Φ1 = −

1



(ρ + ϵδρ)(Qϵ )t dx +

1



ρ Qt dx = −

ϵδρ Qt + [ϵδρ qρ ]t R(x)dx + o(ϵ).

(46)

0

0

0

1



For the remaining terms, we obtain in a similar fashion

−Φ2;ϵ + Φ2 = −

1



1



ϵδρ



(v + Q )2 2

0

−Φ4;ϵ + Φ4 = −χ Cc

 ρ qdy dx + o(ϵ),

x

0

−Φ3;ϵ + Φ3 = −Cf

1

  ϵδρ (v + Q )q + ρ(v + Q )qρ + qρ

1





1

 ρ(v + Q )dy dx + o(ϵ),

1

 ρ dy dx + o(ϵ).

+ qρ x

  ϵδρ (v + Q ) + qρ

0

x

Then, the directional derivative Φρ is given as

Φρ [δρ; ρ, z , v, Tc ] = −

1

 0

   1 φ δρ ρ(v + Q )qρ + qρ R(0)S (x) − dx − R(x)[δρ qρ ]t dx ρ 0

with



φ[ρ, z , v, Tc ] := − ρ Qt + ρ(v + Q )q + Cf ρ

(v + Q )2 2

 + χ Cc ρ(v + Q ) .

(47)

Before finalizing the ξv -integral, we have to integrate by parts in order to get rid of the time derivative of the perturbation δρ in the term (46). This yields  t =tend  tend  1  1 R(x) 1 − ξv [ϵδρ qρ ]t dxdt + ξv ϵδρ qρ R(x)dx R(0) R(0) 0 0 0 t =0    tend  1 R(x) = ϵδρ qρ ξv dxdt R(0) t 0 0    tend  1 R(x) Rt (x) Rt (0)R(x) = ϵδρ qρ (ξv )t + ξv − ξv dxdt R(0) R(0) R(0)2 0 0     tend  1  1 Rt (x) Rt (0)R(x) R(x) + ξv − ξv dxdt = ϵδρ qρ ξv S (0) + ξρ ρx + ξz zx dx R(0) R(0) R(0)2 0 0 0     tend  1  1 1 Rt (0)R(x) R(x) = R(x)S (0) + Rt (x) − dx. ϵδρ qρ ξρ ρx + ξz zx dx + ξv R(0) 0 R(0) R(0) 0 0 Thus, our ξv -integral (45) equals

   1 Rt (0)R(x) I ξv = −ρ(v + Q ) − R(0)S (x) + R(x)S (0) + Rt (x) − ϵδρ qρ ξv R(0) R(0) 0 0     1  tend  1 R(x) φ 1 Φ + ξρ ρx + ξz zx dx dxdt − − dxdt ϵδρ ξv R(0) 0 R ( 0 ) R ( 0 ) ρ 0 0  t =tend  1 1 − ξv ϵδρ qρ R(x)dx + o(ϵ). R(0) 0 t =0 (d) The ξTc -integral (27). 

tend



1

(48)

i Using the definition of TGas in (10) together with the relation (14), we obtain



tend

−

1

   ξTc (Tc )t + hc Tc −

0

0 tend



1



ϵδρ

=− 0

0

p 0 hc

ρ2

p0

ρ + ϵδρ

ξTc dxdt + o(ϵ).

 dx

tend

 dt + 0

1

   ξTc (Tc )t + hc Tc − 0

p0

ρ

 dx

dt

(49)

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

11

(e) The summary of the first variation of L with respect to ρ . Everything is prepared and we can conclude the calculations. Using the terms (43), (44), (48), (49), and the derivatives for the Eqs. (28) and (29), we obtain:

 tend  1   1 ∂ L(W , zl , Λ) ϵδρ (ξρ )t + (v + Q )(ξρ )x − qρ F dxdt (δρ) = lim ϵ→ 0 ∂ρ ϵ 0 0      tend  1 1 Φ φ p0 hc ϵδρ −ξv + − − χ zKρ′ (ρ)ξz − χ 2 ξTc dxdt R(0) R(0) ρ ρ 0 0 x=1  t =tend  tend  1  1 1 − ξv ϵδρ(v + Q )ξρ dt ϵδρξρ dx − ϵδρ qρ R(x)dx + R(0) 0 0 0 x =0 t =0   1  tend ϵδρ(x, 0)νρ dx + o(ϵ) , ϵδρ(0, t )ηρ dt − −

(50)

0

0

where F is defined as F [ρ, z , v, Tc , ξρ , ξz , ξv ] := ξρ ρ +

ξρ ρx + ξz zx dy − x

+ ξv

R(x)

1



R(0)

1



ξρ ρx + ξz zx dx 0

  Rt (0)R(x) ρ(v + Q ) + R(0)S (x) − R(x)S (0) − Rt (x) + . R(0) R(0) 1

(51)

Using that δρ is arbitrary, we end up with

−(ξρ )t − (v + Q )(ξρ )x = −qρ F − ξv

1 R(0)



Φ φ − R(0) ρ



− χ zKρ′ (ρ)ξz − χ

p0 hc

ρ2

ξTc ,

(52)

the terminal condition

ξρ = 0 for t = tend ,

(53)

and the boundary condition

ξρ (v + Q ) = 0 for x = 1.

(54)

Remark 3. Note that we will have a problem evaluating the boundary condition for x = 1, if the velocity u = v + Q at x = 1 would vanish at some time. We therefore have to assume positive velocities at x = 1 for all times t ∈ (0, tend ]. The only exception will be at t = 0, where we will have u(x, 0) = 0 in the case of an engine start. The boundary condition for ξρ at time t = 0 is then given by continuation. 4.4. Coupling conditions for the adjoint equations Now, let us discuss the derivation of coupling conditions for the adjoint equations. The above derivation was valid for a single pipe, not connected to any other pipe, i.e., we neglected the coupling conditions in the Lagrangian (31). Let us consider the first variation of the Lagrangian for the whole network (23)–(32) with respect to the density in an inner pipe i. Similarly, as outlined in Section 4.3.2, we would deduce (52) and (53). Since i ∈ {2, . . . , nP − 1} we have to neglect the boundary condition for the density (28). Finally, we would be left only with the following integrals: !

0 =

 tend   ∂ L(W , zl , Λ) (δρ) = − δρ i (Li , t ) (v i (t ) + Q i (Li , t ))ξρi (Li , t ) + ζρi+1 (t ) dt i ∂ρ 0  tend   + δρ i (0, t ) (v i (t ) + Q i (0, t ))ξρi (0, t ) + ζρi (t ) dt . 0

Using the fact that δρ i is arbitrary, we obtain

(v i (t ) + Q i (Li , t ))ξρi (Li , t ) = −ζρi+1 (t )

(55)

(v i (t ) + Q i (0, t ))ξρi (0, t ) = −ζρi (t ).

(56)

and

So increasing the indices of all terms in (56), we can write together with (55):

(v i (t ) + Q i (Li , t ))ξρi (Li , t ) = −ζρi+1 (t ) = (v i+1 (t ) + Q i+1 (0, t ))ξρi+1 (0, t ).

(57)

12

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

4.5. Summary—the optimality system We want to summaries the results of all computation in this section. Our optimality system now consists of (a) constraints or state equations, (b) adjoint or co-state equations and (c) the optimality condition. (a) Constraints or state equations. The constraints consist of the governing system of Eqs. (13), the closing relation (ideal gas law) (14) and the initial, boundary and coupling conditions (15)–(17). (b) Adjoint or co-state equations.1

−(ξρ )t − (v + Q )(ξρ )x = − qρ F − χ i

i

i

i

i

i

i



i

Kρ′ (ρ i )ξzi

p0 hc

ξ

i i 2 Li T c

z + (ρ )   i i i i i i i i i i −(ξz )t − (v + Q )(ξz )x = − qz F − ξz χ K (T ) − q ,  Li ξρi ρxi + ξzi zxi dx, −(ξvi )t = − ξvi S i (0) −



− ξv i

1 Ri (0)



Φi φi − Ri (0) ρi



,

(58)

0

−(ξTjc )t = −

Lj



j

0

j

qTc F j dx − hc ξTc + (Tcj − Topt ),

for all pipes i = 1, . . . , nP , j ∈ Icc and (x, t ) ∈ (0, Li ) × (0, tend ), the terminal conditions

ξρi (x, tend ) = 0,

ξzi (x, tend ) = 0,

ξvi (tend ) = 0,

ξTjc (tend ) = 0,

(59)

for all pipes i = 1, . . . , nP , j ∈ Icc and x ∈ [0, Li ], boundary conditions2

ξρnP (LnP , t )(v nP (t ) + Q nP (LnP , t )) = 0, ξznP (LnP , t )(v nP (t ) + Q nP (LnP , t )) = 0,

(60)

for all times t ∈ [0, tend ], and coupling conditions3

(v i (t ) + Q i (Li , t ))ξρi (Li , t ) = (v i+1 (t ) + Q i+1 (0, t ))ξρi+1 (0, t ), (v i (t ) + Q i (Li , t ))ξzi (Li , t ) = (v i+1 (t ) + Q i+1 (0, t ))ξzi+1 (0, t ),

(61)

for all pipes i = 1, . . . , nP and times t ∈ [0, tend ]. The variables qiρ , qiz , qiTc , φ i , S i and F i were defined in (41), (B.14), (B.1), (47), (37) and (51), respectively. Furthermore, since we will be interested in the quantity ηz , the relation4

ξz1 (0, t )(v 1 (t ) + Q 1 (0, t )) = ηz (t )

(62)

for all t ∈ [0, tend ] is needed. (c) Optimality condition

(σ + ηz )(zl∗ − zl ) ≥ 0 ∀ 0 ≤ zl∗ ≤ zlmax

(63)

for all t ∈ (0, tend ). 5. Discretization In order to discretize the optimization problem given in the section above, there are two possibilities. First, we could consider any discretization of the state equation, the control, and the cost functional to obtain a finite dimensional optimization problem. For this problem one can derive optimality conditions, using adjoint calculus, similar to the continuous case above, but doing all calculations for the discrete equation. When following this approach, we have, by construction, the discrete derivative information on any given fixed mesh. Thus, the discrete optimization problem can be solved by any discrete derivative based optimization algorithm. However, depending on the choice of discretization the

1 See Appendix B for the derivation of equations for ξ i and ξ i . z Tc 2 The boundary condition for ξ nP is derived in Appendix B.2. z

3 See Appendix B.2.1 for the derivation of the coupling condition for ξ i . z 4 See Appendix B.2 and Eq. (B.22) for the derivation.

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

13

derivative information need not converge to a limit, and thus the sequence of discrete minimizers need not converge to a stationary point for the continuous limit problem, see, e.g., [25]. The second approach is to consider a discretization of the optimality system consisting of state and adjoint equations as well as the representation of the gradient (63). The advantage of this approach is that the discrete adjoint solutions converge to the solution of the adjoint (58) equation, provided the discretization is chosen in a reasonable way, due to construction. The drawback is also evident, the discretized adjoint need not coincide with the discrete adjoint obtained in the first approach. Thus the discretized adjoint provides only approximate derivative information on any discretization level. To deal with the discrepancy of the two approaches several possibilities exist. One is the construction of adjoint consistent schemes, i.e., discretization methods that ensure that both approaches yield the same adjoints. Such a reformulation is analyzed, for instance in [26] for Runge–Kutta schemes applied to optimal control with ODEs. In the context of PDEs, such reformulations may become more involved, since spatial irregularities may influence the choice of the temporal discretization, see, e.g., [27]. The second alternative, which we are going to pursue in the following is to consider the second approach while taking care, using appropriate mesh refinement, that the error in the discrete derivatives is sufficiently small. For a similar approach, where the error in the derivatives is assessed by a posteriori error estimation, see [28]. Thus, we follow the second approach to numerically calculate a solution to the necessary optimality conditions. In particular, the calculated discretized adjoint converges to the continuous adjoint, thus smallness of the calculated gradient is related to smallness of the real gradient. From the discussion above, we already know, that we only have an approximation to the discrete derivatives. If this approximation is not good enough, the calculated approximate negative gradient direction need not be a ‘‘descent directions’’ and thus fail to give descent for the discretized functional. However, if this happens any further iteration on the given discretization is misleading anyway, since discretization errors become dominant so that a refinement of the discretization is warranted. Hence failure of convergence without nearly satisfied optimality conditions serves us as a cheap estimate for the accuracy of the applied discretization; for more details see Table 2 and the discussion in Section 5.3. For the discretization of the state equation (13) an explicit upwind scheme for the spatial differential operator and explicit Euler for the time derivative is used (see [8] for details). Since the adjoint system (58) is posed backwards in space and time, we first substitute the time and space variables by tˆ := tend − t and xˆ = Li − x. Then the same numerical schemes that we use for solving the state system is used for the solution of (58). By the operation ‘‘·’’ between two elements of RM , we denote the weighted scalar product, which is an approximation for the L2 scalar product, i.e., for two mappings ϕ, ψ : [0, 1] → R and their discretizations ϕh , ψh ∈ RM we have

ϕh · ψh =

M 1 

M

(ϕh )m (ψh )m ≈ (ϕ, ψ)L2 =

1



m=1

ϕψ dy,

(64)

0

noting that we have uniform step sizes. In order to avoid too many sub-indices, we will not discriminate between the continuous and discrete state, adjoint and control variables, as we will always work with the discrete quantities in this section. Discretization of functionals such as j : L1 → R, will be denoted by the sub-index h, i.e., j ≈ jh . 5.1. Algorithm We will have two stopping criteria. 1.

STOP if the optimality condition ∥P (j′ )∥ ≈ projected gradient, i.e.,

1 M

∥P (j′h )∥2 < TOLopt with a tolerance TOLopt > 0 and P (j′h ) the discretized

σ + (ηz )m (P (j′h ))m = (P (σ (1, . . . , 1)T + ηz ))m = min(0, σ + (ηz )m ) max(0, σ + (ηz )m ) 

2.

0 < (zl )m < zlmax , (zl )m = 0, (zl )m = zlmax .

(65)

STOP if the value of the cost functional does not change any more, i.e., |jh (zlk ) − jh (zlk+1 )| < TOLdiff , with a tolerance TOLdiff > 0.

The first criterion is related to almost satisfied KKT conditions. The second criterion, however, can occur whenever the step (k) (k+1) size zl −zl tends to zero. This is the case, in particular, when the computed continuous gradient and the discrete gradient are too far apart. Thus if the algorithm stops due to the second criteria a refinement of the discretization is reasonable to assert convergence of the gradient used to determine the search direction. (0)

Algorithm 4. Pick an initial control zl ∈ RM , where M is the number of grid points in time. For k = 0, 1, 2, . . . repeat the following steps until one of the above stopping criteria is fulfilled:

14

I. Gasser et al. / Computers and Mathematics with Applications (

(k)

1. solve the constraints with control zl (k)

(k)

v(zl ), and Tc

(k)

)

–

(k)

(k)

to obtain the corresponding state variables ρ (k) = ρ(zl ), z (k) = z (zl ), v (k) =

= Tc (zl );

(k)

2. solve the adjoint system with state variables ρ (k) , z (k) , v (k) , Tc (k)

(k)

(k)

(k)

(k)

(k)

to obtain the adjoint variables ξρ(k) , ξz , ξv(k) , ξTc , ηz ;

3. use ηz to compute the reduced gradient j′h (zl ) = σ + ηz ∈ RM ; 4. compute step length α via projected line search (Armijo rule applied to jh (zl )); (k+1) 5. set zl = min(zlmax , max(0, zl(k) − α j′h (zl(k) ))) pointwise.

Remark 5. If we would have a consistent approximation of the adjoint equations the above algorithm would simply be a standard projected gradient algorithm, as it has been proposed by [29,30] and we would be content with stopping criterion 1. Its convergence, including inexact line search steps [31], has been analyzed, e.g., in [32]. The influence of discretization on the convergence has been considered, e.g., in [33] In our case, if the algorithm terminates due to stopping criterion 1, and h is sufficiently small (or M sufficiently large) so that the adjoint state is well enough approximated, the norm of the projected gradient is small. It is easy to see that the algorithm is always terminating on any given fixed mesh, since jh is bounded from below, and by construction jh is non increasing. Thus after finitely many iterations stopping criterion 2 must be satisfied. Once the mesh is refined, we can restart the algorithm on the new mesh. To avoid stopping of the algorithm due to too small slope of the cost functional it is advisable to pick TOLdiff = o(1) as M → ∞. 5.2. Numerical test: continuous VS discrete gradient In a first step, we test whether our implementation is correct. In particular, we test the implementation of the derivatives of jh , as we use them as stopping criteria in our algorithm. To do so, we compare directional derivatives with difference approximations, i.e., we check

  jh (zl + ϵδ zl ) − jh (zl )

Eϵ = DQh − j′h (zl ) · δ zl  = 





ϵ

  − j′h (zl ) · δ zl  → 0

for various values of ϵ . Before this, we first need to check the influence of the chosen discretization on Eϵ . In Table 1, we calculated DQh , j′h (zl ) · δ zl and Eϵ for the values (δ zl )m = 1 for all m = 1, . . . , M (unscaled (δ˜zl )m = 0.1) for various values of spatial grid points N. Due to the CFL-condition, this leads also to a refinement of the time mesh, i.e., M = O (N ), where M is the number of time grid points. As we can see, in Table 1, the difference quotient for ϵ = 1 is relatively stable with respect to the mesh size. However, the calculated derivatives are still sensitive to mesh refinement. This implies that even at N = 1600, we will have to expect effects of unresolved derivatives in our optimization algorithms. On the other hand, by comparing the subtables for σ = 0.01 and σ = 10, it is clear, that any numerical test for the correct implementation of the derivative will require a much more refined mesh in space and time. Table 1 Section 5.2: Difference quotient (DQh ) and discretized analytic gradient (j′h (zl )· δ zl ) for ϵ = 1 and zl = 0.1. N

50 100 200 400 800 1600

σ = 0.01

σ = 10

DQh

j′h (zl ) · δ zl

Eϵ

DQh

j′h (zl ) · δ zl

Eϵ

−85.3 −91.8 −93.9 −94.9 −95.4 −95.7

−2052.5 −1964.9 −1876.2 −1835.9 −1810.9 −1798.5

1967.1 1873.1 1782.3 1741.1 1715.5 1702.7

1579.7 1573.2 1571.1 1570.1 1569.6 1569.3

−387.4 −299.9 −211.2 −170.9 −145.9 −133.5

1967.1 1873.1 1782.3 1741.1 1715.5 1702.7

Remark 6. The reason why the error Eϵ in Table 1 is constant for different values of σ , is that for all ϵ > 0 it holds

  ϵ Eϵ = jh (zl + ϵδ zl ) − jh (zl ) − ϵ j′h (zl ) · δ zl    M M   ϵ  ϵ    Tc Tc (δ zl )m − (σ + (ηz )m )(δ zl )m  = jh (zl + ϵδ zl ) − jh (zl ) + σ   M m=1 M m=1   M   ϵ    = jThc (zl + δ zl ) − jThc (zl ) − (ηz )m (δ zl )m  ,   M m=1 and thus, Eϵ is independent of σ . Above, we denoted

 T jhc

(zl ) :=

 i∈Icc

M 1 

2M

 (( ) − Topt )

m=1

Tci m

2

≈

tend

1 i∈Icc

2

0

(Tci − Topt )2 dt .

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

15

Fig. 1. Section 5.2: Error between gradient and difference quotient.

To avoid large influence of the discretization onto Eϵ we note that the discretization error gets smaller if the end time tend is chosen smaller. Now we have to look at the behavior of Eϵ where t˜end is chosen between 1 and 2 s. In Fig. 1, we see the behavior of the error between directional derivatives and difference quotients for various choices of simulation times tend . As it is to be expected, for all values of tend as ϵ decreases so does the error as Eϵ = O (ϵ). As standard numerical analysis reveals, at some point round-off errors become dominant, leading to a behavior Eϵ = O (ϵ −1 ) as it can be seen in the graphic. We can see clearly that the point where round-off errors become dominate travels to larger ϵ as tend growth. However, at small final times, we can see that the error is small. Since the only change in the program is switching the value for the final time we conclude that our implementation yields correct values for the derivatives. 5.3. Numerical test: convergence failure and refinement of discretization As a next test, we come back to our statement at the beginning of Section 5. Namely, we investigate the effect of the inconsistent discretization onto the behavior of the gradient projection algorithm. To this end, we consider the behavior of (0) Algorithm 4 with the same initial value zl = 0 and σ = 0.01 for two different spatial (and thus also temporal) refinements. As we can see from Table 2, already in the first iteration slight differences in the value of the cost functional are visible, this has to be expected from what we have seen in the previous test case, as the discretization error on the interval (0, tend ) is again significant. More importantly for the stopping criterion, the norm of the projected gradient differs significantly between the two meshes. Table 2 Section 5.3: Results of the optimization algorithm for different numbers of spatial grid points N (rounded to five digits). Iteration

0 1 2 3 4 5 6 7 8 9

N = 50

N = 1600

∥P (j′h )∥2

jh

1 M

198.2154 109.9411 39.3864 39.2464 39.2349 39.2349 – – – –

0.0613 0.0241 0.0063 0.0060 0.0057 0.0057 – – – –

∥P (j′h )∥2

jh

1 M

169.9120 98.0908 28.5489 26.8061 26.1496 25.8218 24.2810 23.8705 23.8574 23.8574

0.0093 0.0035 0.0005 0.0012 0.0003 0.0014 0.0009 0.0006 0.0004 0.0004

As predicted since we do not calculate the discrete derivatives, the algorithm becomes stagnant once the error in the calculated gradient approximation becomes too large. Already after the third iteration the value of the cost functional is almost unchanged during the application of Algorithm 4 for N = 50. However, the projected gradient is still large, i.e., 1 ∥P (j′h )∥2 ≥ 5·10−3 . On the other hand when N = 1600 we can continue until M1 ∥P (j′h )∥2 ≈ 5·10−4 with significantly lower M

16

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

Fig. 2. Geometry of the exhaust pipe.

value of j before the cost functional is again stagnant. This confirms our expectation on the convergence of the algorithm and the possible cure for a lack of convergence by means of refinement. 6. Numerical examples We want to show results of two optimization problems we have simulated. 1. Setting: High cost of control, high starting control. Expectation: A decrease of ratio of unburnt gas zl is more important than achieving an optimal temperature in the catalysts. 2. Setting: Low cost of control, low starting control. Expectation: The ratio of unburnt gas zl should be increased in order to reach optimal temperature in the catalysts. For all simulations, that will be presented in this section, the following holds:

• The geometry of the pipe which is used in these simulations, is presented in the Fig. 2 (for data of the pipe’s geometry see Table A.5 in Appendix A).

• We choose N = 1600 to be the number of spatial grid points. The number of grid points in time M depend upon the CFL-condition, and therefore upon the reached velocities in the simulation.

• All quantities presented in this section are unscaled, i.e., they have a dimension. We will denote the unscaled quantities by the tilde symbol ‘‘ ˜· ’’. The table of reference values, used for the scaling, can be found in the Table A.6 in the Appendix. • The value of the parameters used for the simulation are listed in the Table A.7 in the Appendix. • For the verification of the simulation results, one should keep in mind the optimal temperature in the catalytic converters T˜opt = 800 K. • We will observe the temperature in the catalytic converters during the first t˜end = 60 s after the engine start.

• The upper bound for the control variable z˜l will be z˜lmax = 0.5. • The initial condition below corresponds to an engine start. For all pipes i = 1, . . . , 9 and for all x˜ ∈ [0, L˜ i ] ρ˜ ici (˜x, 0) = 1.2 u˜ iic (˜x, 0) = 0

m

kg m3

,

,

z˜ici (˜x, 0) = 0, (66) T˜ci ic (0) = 290.28 K.

s Note that these are the physical and not mathematical initial condition for this problem. For our mathematical model (13), we need an initial condition for v˜ i . From the above condition we can derive it by integrating Eq. (4) we obtain

vici =

Li



uiic (x)dx − 0

Li

 0

x

 0

1

γ p0

  i −h(Tici (y) − TWall ) + χ i (−hc (Tici (y) − Tci ic ) + q0 ρici (y)zici (y)K (Tici (y))) dydx.

By applying the ideal gas law to obtain Tic , we can deduce, with (66),

v˜ ici (0) = 0.

(67)

• The boundary condition for the density for all pipes i = 1, . . . , 9 and for all t˜ ∈ [0, t˜end ] is given by kg . (68) m3 The pressure boundary conditions p˜ l and p˜ r are only physical boundary conditions. For the mathematical model (13) they are only parameters (see Appendix A). The initial guess for the boundary condition for the ratio of unburnt gas (the quantity we want to control optimally) for the two simulations will be shown in the following subsections.

ρ˜ l (t ) = 0.4

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

17

Fig. 3. Example 1: Boundary condition for ratio of unburnt gas z˜ and temperatures in the catalytic converters T˜ci in (K) for some iterations (compare Table 3 for corresponding evaluation of the cost functional).

6.1. Example 1: high cost of control, high starting control The cost of control σ and the value for the first guess of the control variable zl (the boundary condition for the ratio of unburnt gas) for this simulation were

σ = 20,

z˜l (t˜) = 0.15 ∀t˜ ∈ [0, t˜end ].

(69)

The results of the simulation are illustrated in the Table 3 and Fig. 3. The first figure shows the mapping t˜ → z˜l (t˜), where z˜l is the boundary condition for the ratio of unburnt gas. The two other mappings show the temperature development over time in the two catalytic converters (pipe 2 and 4). The table shows the evaluation of the cost functional for the iterations done by the algorithm. J denotes the total cost functional, whereas Jz and JT i denote the cost of the control variable zl and c the cost of missing the optimal temperature Topt in the catalysts, respectively. We split the cost functional as follows jh = σ jzh +



Tl

jhc ,

(70)

l∈Icc

jzh Ti

:=

jhc :=

M tend 

M

(zl )m ≈

2M

zl (t )dt ,

(71)

0

m=1

M tend 

tend



((Tci )m − Topt )2 ≈

m=1

1 2

tend



(Tci (t ) − Topt )2 dt .

(72)

0

The tables contain the scaled quantities. The key quantities like relation between initial and final cost and the optimality conditions, will retain their informative value, despite scaling. Table 3 Example 1: Evaluation of the scaled cost functional (rounded to three decimal places) and optimality condition (rounded to five decimal places) for different control variables z˜l , computed by the optimization algorithm. (Compare Fig. 3 for corresponding control.) T2

T4

∥P (j′h )∥

Iteration

jh

jzh

jhc

jhc

1 M

0 1

5042.791 169.912

250 0

17.842 36.764

24.948 133.149

0.01698 0

We observe that after one iteration the optimization algorithm stops, since the optimality condition is fulfilled. Due to the high cost of the control, the optimal solution is z˜l (t˜) = 0 for all t˜ ∈ [0, t˜end ]. This means for our application, that one should use at least the stoichiometric amount of air in the combustion chamber of the engine, i.e., enough air for a complete combustion of the fuel, such that no unburnt gas enters the pipe. In other words, fuel is so expensive, it should not be used for heating up the catalytic converters. In Fig. 4, we can see the steady state solutions of the state variables for the initial control (black dotted lines) and the optimal control (green dashed lines). The solid black lines show the geometry of the exhaust pipe and the blue rectangles denote the catalysts.

18

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

Fig. 4. Example 1: Results of numerical simulation of the state variables velocity u˜ in ( ms ), ratio of unburnt gas z˜ and gas temperature T˜ in (K) at time t˜∗ = t˜end = 60 s (compare colors given in Table 3).

In the first plot of Fig. 4, we see the decreasing velocity profile. The reason for this behavior is the increasing density/decreasing temperature. The changes in diameter of the exhaust pipe lead to changes of the velocity at the junctions. The second plot shows the ratio of unburnt gas in the exhaust pipe. In the first iteration (black dotted line), we have a concentration of 0.15 as the boundary condition for the unscaled ratio of unburnt gas, which decreases in the two catalytic converters during the exothermic reaction. The temperature (third plot) increases in the catalytic converters, in the case in which we have a positive concentration of unburnt gas, and decreases over the whole exhaust pipe due to the heat transfer with the (colder) wall. 6.2. Example 2: low cost of control, low starting control The cost of control σ and the value for the first guess of the control variable zl (the boundary condition for the ratio of unburnt gas) for this simulation are

σ = 0.01,

z˜l (t˜) = 0

∀t˜ ∈ [0, t˜end ].

(73)

The results of the simulation are illustrated in the Table 4 and Fig. 5. In this scenario (fuel is cheap), the optimization algorithm suggests to use more of fuel. After 9 iterations this yields to our ‘‘optimal’’ control, although the stopping criteria which led to the abortion of the algorithm, was the second criterion (no change of the cost functional due to 20 line search attempts). Nevertheless, the obtained control leads to a fast heating a temperature close to the optimal T˜opt = 800 K in both catalytic converters, i.e., from the application’s point of view: a satisfying result. Table 4 Example 2: Evaluation of the scaled cost functional (rounded to three decimal places) and optimality condition (rounded to five decimal places) for different control variables zl , computed by the optimization algorithm. (Compare Fig. 5 for corresponding control.) T2

T4

∥P (j′h )∥

Iteration

jh

jzh

jhc

jhc

1 M

0 1 2 3 4 5 6 7 8 9

169.912 98.091 28.549 26.806 26.150 25.822 24.281 23.870 23.857 23.857

0 482.276 317.238 282.716 333.782 266.863 283.835 295.016 302.555 302.555

36.765 51.327 15.390 10.502 14.747 8.843 9.858 10.702 11.346 11.346

133.148 41.941 9.986 13.477 8.065 14.310 11.584 10.219 9.486 9.486

0.00933 0.00347 0.00058 0.00130 0.00035 0.00139 0.00091 0.00061 0.00042 0.00042

7. Conclusion In this paper, we were able to answer the question, how to ensure reaching an optimal temperature in a catalytic converter of an exhaust pipe after the engine start by controlling ratio of unburnt gas in the gas mixture. This was achieved using the formal continuous adjoint for the calculation of derivatives in a projected gradient algorithm. Numerical examples demonstrate the feasibility of the approach, and show that stagnation of the projected gradient algorithm is due

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

19

Fig. 5. Example 2: Boundary condition for ratio of unburnt gas z˜ and temperatures in the catalytic converters T˜ci in (K) for some iterations (compare Table 4 for corresponding evaluation of the cost functional). Table A.5 Pipe geometry. Pipe number

Length (m)

Radius (m)

CC

1 2 3 4 5 6 7 8 9

0.415 0.12 0.93 0.1 0.45 0.47 0.17 0.43 0.515

0.021 0.04 0.021 0.06 0.021 0.06 0.021 0.095 0.021

0 1 0 1 0 0 0 0 0

Sum

3.6

–

2

Table A.6 Reference values. Quantity

Meaning

Unit

Reference quantity

Reference value

t˜ x˜ ρ˜ u˜ p˜ T˜ z˜

Time Space Density Velocity Pressure Temperature Ratio of unburnt gas

s m kg m−3 m s−1 kg m−1 s−2 K

tr = xr /ur xr = L˜

0.36 s 3.6 m 1.2 kg m−3 10 m s−1 105 kg m−1 s−2 290.28 K 0.1

ρr

ur pr Tr = pr /(Rρr ) zr

to insufficient resolution of the differential equations and can thus be healed by refinement of the discretization. To resume, this paper can be seen as a successful example of how optimization of a complex system under the constraint of a set of a nonlinear system of PDE’s of mixed type, i.e., a compressible low Mach number gas dynamic model, can be realized in reasonable simulation times. Appendix A. Data

h :=

˜ r 4hx d˜ ρr ur cv

,

hc :=

h˜ c xr

ρr ur cv

,

q0 :=

ρr zr q˜ 0 R pr cv

,

Cf :=

C˜f xr d˜

,

Cc :=

C˜c xr ur

.

(A.1)

The physical boundary conditions for the pressure used in all simulations are p˜ l = 1.01 bar,

p˜ r = 1 bar.

(A.2)

The following Tables A.5–A.7 contain the information of the pipe’s geometry, reference values and parameters used for the simulations, respectively.

20

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

Table A.7 Parameters, units, parameter values. Parameter

Meaning

C˜f

Unit

cv h˜ p0 R T˜out

γ

Wall friction Specific heat at constant volume Heat transfer rate with wall Scaled pressure in leading order Ideal gas constant External temperature Adiabatic exponent

C˜c T˜ + h˜ c K˜ 0 q˜ 0 T˜opt

Catalyst friction Activation temperature Heat transfer rate with catalyst Reaction rate Heat release rate in catalyst Optimal temperature in catalyst

m2 s−2 K−1 kg s−3 K−1 m2 s−2 K−1 K

s−1 K kg s−3 K−1 s−1 m2 s−2 K

Parameter value 0.0241 717.7 100 1 287.08 290.28 1.4 800 600 100 100 5 · 106 800

Appendix B. First variation of L with respect to Tc and z For the sake of completeness in the derivation of the adjoint system, we will compute the missing first variation of the Lagrangian L (23)–(32) with respect to Tc and z. We will derive the adjoint equations like in Section 4.3, i.e., for a pipe of length L = 1 with a catalytic converter (χ = 1). B.1. First variation of L with respect to Tc Before we start with the difference quotient, we first want to determine the differences of the q, Q and Qt , evaluated at Tc + ϵδ Tc and Tc : q[ρ, z , Tc + ϵδ Tc ](x, t ) − q[ρ, z , Tc ](x, t ) = ϵδ Tc (t )qTc , Q [ρ, z , Tc + ϵδ Tc ](x, t ) − Q [ρ, z , Tc ](x, t ) = ϵδ Tc (t )

x



qTc dy, 0

Qt [ρ, z , Tc + ϵδ Tc ](x, t ) − Qt [ρ, z , Tc ](x, t ) = [ϵδ Tc (t )]t

x



qTc dy, 0

with qTc :=

hc

γ p0

.

(B.1)

We will compute the first variation of the Lagrangian with respect to Tc step by step. Since Tc also appears in the cost functional, we want to compute the difference quotient of it as well. In order to shorten the expressions, we define the following abbreviations for this section: q := q[ρ, z , Tc ],

qϵ := q[ρ, z , Tc + ϵδ Tc ],

Q := Q [ρ, z , Tc ],

Qϵ := Q [ρ, z , Tc + ϵδ Tc ].

(a) The cost functional (23). It is straightforward to see

J [Tc + ϵδ Tc , zl ] − J [Tc , zl ] =

tend



ϵδ Tc (Tc − Topt )dt + o(ϵ).

(B.2)

0

(b) The ξρ -integral (24) tend

 tend  1     ξρ ρt + (v + Qϵ )ρx + qϵ ρ dxdt + ξρ ρt + (v + Q )ρx + qρ dxdt 0 0 0 0   tend  1   x =− ξρ ρx ϵδ Tc qTc dy + ϵδ Tc ρ qTc dxdt 



1

−

0

 =− 0

0 tend

ϵδ Tc

0 1



1

 qTc

0

x

 ξρ ρx dy + ξρ ρ dxdt .

(B.3)

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

21

(c) The ξz -integral (25) tend





1

ξz zt + (v + Qϵ )zx + χ zK (T ) dxdt + 

− 0 tend

0





ξz zx ϵδ Tc

tend



  ξz zt + (v + Q )zx + χ zK (T ) dxdt 0

x



qTc dydxdt 0

0

0

1



0 1



=− =−

tend



ϵδ Tc

1



1



ξz zx dydxdt .

qTc

(B.4)

x

0

0

(d) The ξv -integral (26). With

Φ := Φ [ρ, z , v, Tc ],

Φϵ := Φ [ρ, z , v, Tc + ϵδ Tc ],

we deduce tend

    tend  1 1 ξv vt − ξv vt − Φϵ dt + Φ dt R(0) R(0) 0 0  1  x  1  x  tend 1 =− ξv ρ[ϵδ Tc ]t qTc dydx + ϵδ Tc ρ(v + Q )qTc + ϵδ Tc ρ q qTc dydx R(0) 0 0 0 0 0   1  x  1  x + Cf ϵδ Tc ρ(v + Q ) qTc dydx + χ Cc ϵδ Tc ρ qTc dydx dt + o(ϵ) 

−

0

0

tend

 =−

ξv

0

1 R(0)

1



ϵδ Tc qTc

+ Cf



0

0

1



[ϵδ Tc ]t qTc R(x)dx +

ϵδ Tc ρ(v + Q )qTc dx +

0

0

1



ρ(v + Q )dydx + χ Cc

0

1



x

1

 0

1

ϵδ Tc qTc

1



0

1



ϵδ Tc qTc

ρ qdydx x



ρ dydx dt + o(ϵ) x

:= Iξv .

(B.5)

We have to deal with the time derivative in this equation. Since we finally want to isolate the variation δ Tc we have to integrate by parts. tend



1



ξv

− 0

0

R(x) R(0)

tend



1



ϵδ Tc qTc ξv

[ϵδ Tc ]t qTc dx = 0

R(x)





R(0)

0

1



ϵδ Tc qTc ξv

dx − 0

t

R(x) R(0)

t =tend dx

.

t =0

With (36) we can get rid of the time derivative5



ξv

R(x)



R(0)

= (ξv )t t

 =

R(x) R(0)

+ ξv

ξv S (0) +



Rt (x) R(0)

− ξv

Rt (0)R(x)

1

ξρ ρx + ξz zx dx



=

1



R(0)

ξρ ρx + ξz zx dx + ξv 0

R(x) R(0)

0

R(x)

(B.6)

R(0)2



1 R(0)

+ ξv

Rt (x) R(0)

− ξv

Rt (0)R(x)

(B.7)

R(0)2

R(x)S (0) + Rt (x) −

Rt (0)R(x) R(0)

1





ρ dy .

(B.8)

x

Finally, we obtain for the ξv —integral (B.5): tend

 I ξv =

ϵδ Tc

0

+

1



R(x) R(0)

0

1



1 R(0)

R(x)S (0) − R(0)S (x) + Rt (x) −

Rt (0)R(x)

 − ρ(v + Q )

R(0)  t =tend  1 1 ξρ ρx + ξz zx dx dxdt − ϵδ Tc ξv qTc R(x)dx + o(ϵ). R(0) 0 t =0 0





qTc ξv



(B.9)

(e) The ξTc -integral (27) tend

 −

  ξTc (Tc + ϵδ Tc )t − hc (Tgas − (Tc + ϵδ Tc )) dt +

0

tend



  ξTc (Tc )t − hc (Tgas − Tc ) dt

0 tend

 =−

  ξTc (ϵδ Tc )t + hc ϵδ Tc dt =

0

tend

 0

   t =t ϵδ Tc (ξTc )t − hc ξTc dt − ϵδ Tc ξTc t =0end .

(f) The summary of the computation for the first variation of L with respect to Tc . 5 The variable S is defined in (37).

(B.10)

22

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

Since everything is prepared, we can start with the derivative. Combining the terms (B.2), (B.3), (B.4), (B.9), and (B.10), we obtain:

   tend  1 ∂ L(W , zl , Λ) 1 qTc Fdx dt (δ Tc ) = lim ϵδ Tc (ξTc )t − hc ξTc + (Tc − Topt ) − ϵ→0 ϵ ∂ Tc 0 0    t =tend  1  t =tend  1  − ϵδ Tc ξv − ϵδ Tc ξTc t =0 − νTc ϵδ Tc  +o(ϵ) qTc R(x)dx R(0) 0 t =0 t =0

(B.11)

where F is defined (51). By using that δ Tc is arbitrary, we end up with

−(ξTc )t = −

1



qTc Fdx − hc ξTc + (Tc − Topt ),

(B.12)

0

and the terminal condition

ξTc = 0 for t = tend .

(B.13)

B.2. First variation of L with respect to z Before we start with the difference quotient, we first want to determine the differences of the q and Q , where z appears as a dependent variable. q[ρ, z + ϵδ z , Tc ] − q[ρ, z , Tc ] = ϵδ zqz Q [ρ, z + ϵδ z , Tc ] − Q [ρ, z , Tc ] = Qt [ρ, z + ϵδ z , Tc ] − Qt [ρ, z , Tc ] =

x



ϵδ z (y, t )qz (y, t )dy

0 x

[ϵδ z (y, t )qz (y, t )]t dy

0

with qz := χ

q0

γ p0

ρ K (T ).

(B.14)

We compute now the variation of the Lagrangian step by step. In order to shorten the expressions, we define only for the computation of ∂ L/∂ z the following terms: q := q[ρ, z , Tc ],

qϵ := q[ρ, z + ϵδ z , Tc ],

Q := Q [ρ, z , Tc ],

Qϵ := Q [ρ, z + ϵδ z , Tc ].

(a) The ξρ -integral (24) tend

 tend  1     ξρ ρt + (v + Qϵ )ρx + qϵ ρ dxdt + ξρ ρt + (v + Q )ρx + qρ dxdt 0 0 0 0   tend  1   x =− ξρ ρx ϵδ zqz dy + ϵδ z ρ qz dxdt 



1

−

0

0 tend



0 1



ϵδ z qz

=− 0

1

 

0



ξρ ρx dy + ξρ ρ qz dxdt .

(B.15)

x

(b) The ξz -integral (25) tend





1

  ξz (z + ϵδ z )t + (v + Qϵ )(z + ϵδ z )x + χ (z + ϵδ z )K (T ) dxdt

− 0

0 tend



1



  ξz zt + (v + Q )zx + χ zK (T ) dxdt 0 0   tend  1   x =− ξz ϵδ zt + (v + Q )ϵδ zx + zx ϵδ zqz dy + χ ϵδ zK (T ) dxdt + o(ϵ) +

0

0

tend



1



ϵδ z (ξz )t + (v + Q )(ξz )x + ξz q − qz

= 0

0



0 tend

 − 0



1

 ξz zx dy − χ ξz K (T ) dxdt

x

ϵδ z (v + Q )ξz dt

x=1

 −

x =0

0

1

ϵδ z ξz dx

t =tend t =0

+ o(ϵ).

(B.16)

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

23

(c) The ξv -integral (26). With

Φ := Φ [ρ, z , v, Tc ],

Φϵ := Φ [ρ, z + ϵδ z , v, Tc ],

we deduce tend





ξv vt −

− 0

tend

 =−

ξv

0



1 R(0)

1 R(0)

1

ρ

tend

=−

ξv

0

1 R(0)

1



ϵδ zqz

+ Cf 0

Φ dt

1

ϵδ z ρ(v + Q )qz + ρ q

x

ϵδ zqz dydx + χ Cc

1



ρ

1

[ϵδ zqz ]t R(x)dx +

x



ϵδ zqz dydx dt + o(ϵ)

1

ϵδ z ρ(v + Q )qz dx +

ρ(v + Q )dydx + χ Cc x

 0

0

1

1



ϵδ zqz dydx



0

0



x

 0

0

0



R(0)

[ϵδ zqz ]t dydx +

0





1



0



0



ξv vt −

x



0

ρ(v + Q )

+ Cf



0



1



Φϵ dt +

tend



ϵδ zqz



0

x

1

1

ϵδ zqz

1



ρ qdydx x

 ρ dydx dt + o(ϵ) =: Iξv .

(B.17)

Like in the case of Tc , we have to deal with the time derivative in this equation. Since we want to isolate the variation δ z we have to integrate by parts. tend



1



ξv

− 0

0

R(x) R(0)

tend





[ϵδ zqz ]t dxdt = 0

0

1

   1 t =tend R(x) R(x) ϵδ zqz ξv dxdt − ϵδ zqz ξv dx . R(0) t R(0) 0 t =0

We replace the time derivative like in (B.6)–(B.8). We obtain finally for the ξv —integral (B.17)

   Rt (0)R(x) 1 −ρ(v + Q ) − R(0)S (x) + R(x)S (0) + Rt (x) − ϵδ z ξv qz R(0) R(0) 0 0   1 t =tend  1 R(x) R(x) + qz ξρ ρx + ξz zx dx dxdt − ϵδ zqz ξv dx + o(ϵ). R(0) 0 R (0) 0 t =0

 I ξv =

tend



1

(B.18)

We can combine the terms (B.15), (B.16), and (B.18) to get:

 tend  1   ∂ L(W , zl , Λ) 1 (δ z ) = lim ϵδ z (ξz )t + (v + Q )(ξz )x − ξz (χ K (T ) − q) − qz F dxdt ϵ→0 ϵ ∂z 0 0  tend x=1  1 t =tend  1 t =tend R(x) − ϵδ z (v + Q )ξz dt − ϵδ z ξz dx − ϵδ zqz ξv dx R(0) 0 0 0 x =0 t =0 t =0   tend  1 − ϵδ z (0, t )η2 (t )dt − ϵδ z (x, 0)ν2 (x)dx + o(ϵ) . 0

0

Using that δ z is arbitrary, we end up with

−(ξz )t − (v + Q )(ξz )x = −qz F − ξz (χ K (T ) − q)

(B.19)

with terminal condition

ξz = 0 for t = tend ,

(B.20)

and boundary condition

ξz (v + Q ) = 0 for x = 1.

(B.21)

Furthermore, the spatial boundary condition for ξz at x = 0 plays an important role. By this conditions we compute the adjoint ηz , which also appears in the optimality condition, i.e., is used for the computation of the gradient.

ξz (v + Q ) = ηz for x = 0.

(B.22)

24

I. Gasser et al. / Computers and Mathematics with Applications (

)

–

B.2.1. Coupling conditions for the adjoint equation for ξz We follow the derivation from Section 4.4. Let us consider the first variation of the Lagrangian for the whole network (23)–(32) with respect to the ratio of unburnt gas in an inner pipe i. As outlined in Appendix B.2, we would deduce (B.19) and (B.20). Since i ∈ {2, . . . , nP − 1}, we have to neglect the boundary condition for the ratio of unburnt gas (28). Finally, we would be left only with the following integrals:

 tend   ∂ L(W , zl , Λ) δ z i (Li , t ) (v i (t ) + Q i (Li , t ))ξzi (Li , t ) + ζzi+1 (t ) dt 0 = (δ z ) = − ∂ zi 0  tend  i  i δ z (0, t ) (v (t ) + Q i (0, t ))ξzi (0, t ) + ζzi (t ) dt . + !

0

Using the fact that δ z i is arbitrary, we obtain

(v i (t ) + Q i (Li , t ))ξzi (Li , t ) = −ζzi+1 (t )

(B.23)

(v i (t ) + Q i (0, t ))ξzi (0, t ) = −ζzi (t ).

(B.24)

and

Increasing the indices of all terms in (B.24), we can write combine this with (B.23) to get:

(v i (t ) + Q i (Li , t ))ξzi (Li , t ) = −ζzi+1 (t ) = (v i+1 (t ) + Q i+1 (0, t ))ξzi+1 (0, t ).

(B.25)

References [1] S. Rjasanow, Heat transfer in an insulated exhaust pipe, J. Engrg. Math. 29 (1995) 33–49. [2] S.-M. Liang, S.-J. Tsai, S.-F. Wang, An effective approach for calculation of exhaust pipe flows, J. Mech. 25 (2009) 177–188. [3] G. D’Errico, G. Ferrari, A. Onorati, Numerical modeling of unsteady reacting flows in the exhaust system of a S.I. engine including the catalytic converter, in: Seoul 2000 FISITA World Automotive Congress, 2000. [4] M. Hornikx, W. De Roeck, W. Desmet, Flow and geometrical effects on radiated noise from exhausts computed by a hybrid extended Fourier PSTD method, in: 17th AIAA/CEAS Aeroacoustics Conference, 32nd AIAA Aeroacoustics Conference, 2011. [5] S.-M. Lianga, S.-M. Taib, Analysis and prediction of shock-induced near-field acoustics from an exhaust pipe, Comput. & Fluids 45 (2011) 222–232. [6] V. Jakubauskas, D. Weaver, Transverse vibrations of bellows expansion joints—part I and II, J. Fluids Struct. 12 (1998) 445–473. [7] S.-F. Ling, Vibration isolation of exhaust pipe under vehicle chassis, Int. J. Veh. Des. 15 (1994) 131–142. [8] I. Gasser, M. Rybicki, Modelling and simulation of gas dynamics in an exhaust pipe, Appl. Math. Model. 37 (2013) 2747–2764. [9] I. Gasser, M. Rybicki, W. Wollner, Modelling, simulation and optimization of gas dynamics in an exhaust pipe, in: Proceedings of HYP 2012. [10] M. Gugat, M. Herty, A. Klar, G. Leugering, V. Schleper, Well-posedness of networked hyperbolic systems of balance laws, in: Constrained Optimization and Optimal Control for Partial Differential Equations, Vol. 160, Springer, Basel, 2012, pp. 123–146. [11] S. Ulbrich, A sensitivity and adjoint calculus for discontinuous solutions of hyperbolic conservation laws with source terms, SIAM J. Control Optim. 41 (3) (2002) 740–797 (electronic). [12] M. Gugat, M. Herty, A. Klar, G. Leugering, Optimal control for traffic flow networks, J. Optim. Theory Appl. 126 (2005) 589–616. [13] C. Kirchner, M. Herty, S. Göttlich, A. Klar, Optimal control for continuous supply network models, Netw. Heterog. Media 1 (4) (2006) 675–688. [14] C. Castro, F. Palacios, E. Zuazua, Optimal control and vanishing viscosity for the Burgers equation, in: Integral Methods in Science and Engineering. Vol. 2, Birkhäuser Boston Inc., Boston, MA, 2010, pp. 65–90. [15] S. Collis, K. Ghayour, M. Heinkenschloss, M. Ulbrich, S. Ulbrich, Optimal control of unsteady compressible viscous flows, Internat. J. Numer. Methods Fluids 40 (11) (2002) 1401–1429. [16] R. Petrucci, Un modello matematico di combustione dei gas di scarico in marmitte catalitiche, Master’s Thesis, Roma Tre—Universit delgi Studi, 2007. [17] R. Natalini, L. Lacoste, Mathematical Modeling of Chemical Processes in Exhaust Pipe, Tech. Rep., IAC Rome, 2004. [18] J. Brouwer, I. Gasser, M. Herty, Gas pipeline models revisited: model hierarchies, nonisothermal models, and simulations of networks, Multiscale Model. Simul. 9 (2011) 601–623. [19] M. Herty, V. Sachers, Adjoint calculus for optimization of gas networks, Netw. Heterog. Media 2 (4) (2007) 733–750. [20] M. Frank, A. Klar, R. Pinnau, Optimal control of glass cooling using simplified PN theory, Transport Theory Statist. Phys. 39 (2–4) (2010) 282–311. [21] T.C. Company, Flow of Fluids Through Valves, Fittings, and Pipe—Technical Paper No. 410, Crane Valves, 1982. [22] R. Mulley, Flow of Industrial Fluids—Theory and Equations, CRC Press, 2004. [23] I. Gasser, H. Steinrück, An asymptotic one-dimensional model to describe fires in tunnels III, SIAM J. Appl. Math. 66 (6) (2006) 2187–2203. [24] F. Tröltzsch, Optimal Control of Partial Differential Equations: Theory, Methods, and Applications, in: Graduate Studies in Mathematics, American Mathematical Society, 2010. [25] E. Polak, Optimization: Algorithms and Consistent Approximations, in: Applied Mathematical Sciences, vol. 124, Springer-Verlag, 1997. [26] W.W. Hager, Runge–Kutta methods in optimal control and the transformed adjoint system, Numer. Math. 87 (2) (2000) 247–282. [27] C. Goll, R. Rannacher, W. Wollner, The damped Crank–Nicolson time-marching scheme for the adaptive solution of the Black–Scholes equation, Lehrstuhl Numerische Mathematik, 2011, Preprint. [28] J.C. Ziems, S. Ulbrich, Adaptive multilevel inexact SQP methods for PDE-constrained optimization, SIAM J. Optim. 21 (1) (2011) 1–40. [29] A.A. Goldstein, Convex programming in Hilbert space, Bull. Amer. Math. Soc. 70 (1964) 709–710. [30] E. Levitin, B. Polyak, Constrained minimization methods, USSR Comput. Math. Math. Phys. 6 (5) (1966) 1–50. [31] D.P. Bertsekas, On the Goldstein–Levitin–Polyak gradient projection methods, IEEE Trans. Automat. Control 21 (1976) 174–184. [32] J.C. Dunn, Global and asymptotic convergence rate estimates for a class of projected gradient processes, SIAM J. Control Optim. 19 (3) (1981) 368–400. [33] J.C. Dunn, On state constraint representations and mesh-dependent gradient projection convergence rates for optimal control problems, SIAM J. Control Optim. 39 (4) (2000) 1082–1111.