Schl¨ogl Equation

The first considered system was discussed by Schl¨ogl [83] as a model for an auto-catalytic reaction scheme

A1+ 2X ^k

+

−⇀1

↽−

k⁻₁

3X, X ^k

+

−⇀2

↽−

k⁻₂

A2.

with velocity constantsk⁺₁, k₁⁻, k₂⁺, k₂⁻. A1andA2are chemicals with known constant concentrations c1 and c2 respectively. The constancy of the concentrations can be assured by continuously feeding the reactor where the reaction takes place. Then the reaction kinetics of the chemical X is given by the cubic polynomial

f(y) =−k⁻₁y³+k₁⁺y²−k⁺₂y+k₂⁻c2

where y denotes the concentration of X. For suitable velocity constants, f(y) = 0 has three real-valued roots and we can rewrite f(y) as

f(y) =−k(y−y1)(y−y2)(y−y3)

with constant parameters k and y_i (i = 1,2,3). By combining the reaction term with the diffusion of X we obtain the PDE for the concentration field y(x, t)

yt(x, t) =D∆y(x, t)−k(y(x, t)−y1)(y(x, t)−y2)(y(x, t)−y3) (5.1)

5 Numerical Implementation

whereD denotes the diffusion coefficient. For unbounded domains or large domains with Neumann boundary conditions it is known, cf. [20], that the solution of (5.1) has a wave type behaviour. This can lead to numerical difficulties in the optimal control problem. In our examples with small domains the travelling waves do not play a role.

In the following we consider the case y1 = −1, y2 = 0 and y3 = 1. Then the nonlinearity reduces to f(y) = −µy(y+ 1)(y−1) = µ(y−y³) where k = µ is the reaction parameter. The resulting reaction diffusion equation is in the literature also known as Chaffee Infante equation. Since this model is an easy representative of an equation with a non monotone nonlinearity, it is well suited as a benchmark example for optimization algorithms. Our first test example is a one dimensional system with distributed control. We look at

yt(x, t) =yxx(x, t) +µ(y(x, t)−y(x, t)³) +u(x, t) in (0,1)×(0,∞) (5.2a)

y(0, t) =y(1, t) = 0 in (0,∞) (5.2b)

y(x,0) =y₀(x) in (0,1) (5.2c)

where we impose a homogeneous Dirichlet condition on the boundary. As already discussed in Chapter 3 the origin (y ≡ 0) is an unstable equilibrium for µ ≥ π². Moreover, as we can see in Figure 5.1 (for µ= 15), there are in addition two stable equilibria (solid lines). It can be proven, cf. [21], that for a positive initial function y0(x) the solution converges to the positive equilibrium for t → ∞ and vice versa for negative initial data. Furthermore, any solution is global and there is no blow up in finite time, see also Remark 3.2. The solution trajectory of the uncontrolled equation can be found in Figure 6.1 with y0(x) = 0.2 sin(πx).

Although there exist models in application where a distributed control in the whole spatial domain makes sense, it is often more reasonable that the control can only act on parts of the boundary. Therefore, our next test example is the one dimensional Schl¨ogl equation with Neumann boundary control on the right side

y_t(x, t) =y_xx(x, t) +µ(y(x, t)−y(x, t)³) in (0,1)×(0,∞) (5.3a)

y(0, t) = 0 in (0,∞) (5.3b)

yx(1, t) =u(t) in (0,∞) (5.3c)

y(x,0) =y0(x) in (0,1). (5.3d)

On the left side we impose an homogeneous Dirichlet condition. Since in many chemical engineering problems only the flux and not the value on the boundary can be controlled, the use of Neumann control seems to be reasonable. In Chapter 3 we have seen that the uncontrolled equation is unstable for µ≥ ^π₄. In this case we have again two stable and one unstable equilibrium. The three equilibria of the PDE (5.3) are displayed in Figure 5.2 for µ= 15.

In the last example for the Schl¨ogl model we consider a two dimensional equation

5.1 Examples

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−0.6

−0.4

−0.2 0 0.2 0.4 0.6

x

y

Figure 5.1: Stable equilibria (solid lines) and unstable equilibrium (dashed line) of the uncontrolled one dimensional Schl¨ogl equation (5.2) withµ= 15.

with Neumann control on some parts of the boundary

yt(x, t) = ∆y(x, t) +µ(y(x, t)−y(x, t)³) in Ω×(0,∞) (5.4a)

∂νy(x, t) =u(x, t) on ΓNc×(0,∞) (5.4b)

∂νy(x, t) = 0 on ΓN0 ×(0,∞) (5.4c)

y(x,0) =y0(x) in Ω (5.4d)

The domain is given by the unit square with circular hole, see Figure 5.3 (a) for the detailed description. The boundary ∂Ω is divided into two parts: The boundary of the square ΓN0 and the boundary of the inner circle ΓNc. On ΓN0 we impose a homogeneous Neumann condition while the Neumann control acts on ΓNc.

If we interpret the solutiony(x, t) as a temperature field, the homogeneous Neumann condition on the boundary of the square can be seen as a perfect isolation of the domain. The circle inside the domain can be considered as a cooling device to control the temperature in the reactor by the heat flux.

For µ 6= 0 it is obvious that the uncontrolled equation has the three equilibria y ≡ −1, y ≡ 0 and y ≡ 1. (In the special case µ = 0 we get that y ≡ c is an equilibrium for each c∈ R.) The equilibrium y ≡0 is again unstable for arbitrary µ >0 while y≡ −1 andy ≡1 are stable.

In order to solve the PDE (5.4) we discretize the spatial variable by the finite element method introduced in Section 2.3.1. We use the Matlab PDE Toolbox to generate a mesh of triangular piecewise linear finite elements. The assembling is also done by Matlab. Since we do not utilize adaptive algorithms in space these steps can be done offline before the MPC algorithm starts. For the simulations in Chapter 6 we consider the three different meshes displayed in Figure 5.3 (b)-(d).

5 Numerical Implementation

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

x

y

Figure 5.2: Stable equilibria (solid lines) and unstable equilibrium (dashed line) of the uncontrolled one dimensional Schl¨ogl equation (5.3) withµ= 15.

The corresponding information about the triangulation can be found in Table 5.1.

With hmax we denote the maximum diameter of the triangles. The corresponding number of nodes is given by Mx. Since we use linear finite elements, Mx coincides with the number of state variables. The number of control variables (in this case the number of nodes on the circle) is denoted by Mu. Obviously, the control dimension grows much slower than the dimension of the state. This is due to the fact that we consider boundary control.

h_max M_x M_u Grid I 0.10 183 16 Grid II 0.05 652 24 Grid III 0.02 4089 48

Table 5.1: Parameter values for the triangulations in Figure 5.3 (b)-(d) with max-imum triangle diameter h_max, number of state variables M_x and control variablesMu.