© Springer Nature Switzerland AG 2020
E. Lindner et al. (eds.)Mathematical Modelling in Real Life ProblemsMathematics in Industry33https://doi.org/10.1007/978-3-030-50388-8_7

7. Optimal Heating of an Indoor Swimming Pool

Monika Wolfmayr1  
(1)
Faculty of Information Technology, University of Jyväskylä, Jyväskylä, Finland
 
 
Monika Wolfmayr
Email: monika.k.wolfmayr@jyu.fi
Keywords
Heat equationPDE-constrained optimizationControl constraintsProjected gradient methodFinite element methodImplicit Euler method

7.1 Introduction

Modeling the heating of an object is an important task in many applicational problems. Moreover, a matter of particular interest is to find the optimal heating of an object such that it has a desired temperature distribution after some given time. In order to formulate such optimal control problems and to solve them, a cost functional subject to a time-dependent partial differential equation (PDE) is derived. One of the profound works paving the way for PDE-constrained optimization’s relevance in research and application during the last couple of decades is Lion’s work [9] from 1971. Some recent published monographs discussing PDE-constrained optimization as well as various efficient computational methods for solving them are, e.g., [1, 6], and [11], where the latter one is used as basis for the discussion on solving the optimal heating problem of this work.

The goal of this work is to derive a simple mathematical model for finding the optimal heating of the air in a glass dome represented by a half sphere, where a swimming pool is located in the bottom of the dome and the heat sources (or heaters) are situated on a part of the boundary of the glass dome. The process from the model to the final numerical simulations usually involves several steps. The main steps in this work are the setting up of the mathematical model for the physical problem, obtaining some analytical results of the problem, presenting a proper discretization for the continuous problem and finally computing the numerical solution of the problem. The parabolic optimal control problem is discretized by the finite element method in space, and in time, we use the implicit Euler method for performing the time stepping. The used solution algorithm for the discretized problem is the projected gradient method, which is for instance applied in [5] as well as in more detail discussed in [4, 6, 8, 10].

We want to emphasize that the model and the presented optimization methods for the heating process of this work are only one example for a possible modeling and solution. In fact, the stated model problem has potential for many modeling tasks for students and researchers. For instance, different material parameters for the dome as well as for air and water could be studied more carefully. The optimal modeling of the heat sources could be stated as a shape optimization problem or instead of optimizing the temperature of the air in the glass dome, one could optimize the water temperature. This would correspond to a final desired temperature distribution corresponding to the boundary of the glass dome, where the swimming pool is located, for the optimal control problem. Another task for the students could be to compute many simulations with, e.g., Matlab’s pdeModeler to derive a better understanding of the problem in the pre-phase of studying the problem of this work. However, we only want to mention here a few other possibilities for modeling, studying and solving an optimal heating problem amongst many other tasks, and we are not focusing on them in the work presented here.

This article is organized as follows: First, the model of the heating process is formulated in Sect. 7.2. Next, Sect. 7.3 introduces the optimal control problem, which describes the optimal heating of the glass dome such that the desired temperature distribution is attained after a given time. In Sect. 7.4, proper function spaces are presented in order to discuss existence and uniqueness of the optimal control problem. We derive the reduced optimization problem in Sect. 7.5 before discussing its discretization and the numerical method for solving it, the projected gradient method, in more detail in Sect. 7.6. Numerical results are presented as well as conclusions are drawn in final Sect. 7.7.

7.2 Modeling

This section presents the modeling process. The physical problem is described in terms of mathematical language, which includes formulating an initial version of the problem, but then simplifying it in order to derive a version of the problem which is easier to solve. However, at the same time, the problem has to be kept accurate enough in order to compute an approximate solution being close enough to the original solution. That is exactly one of the major goals of mathematical modeling. In the following, we introduce the domain describing the glass dome, where an indoor swimming pool is located in the bottom of the dome, and the position of the heaters. The concrete equations describing the process of heating and the cost functional subject to them modeling the optimization task are discussed in next Sect. 7.3.

We have an indoor swimming pool which is located under a glass dome. For simplicity, it is assumed that we have an isolated system in the glass dome, so no heat can leak from the domain. The swimming pool covers the floor of the glass dome and we assume that the heaters are placed next to the floor up on the glass all around the dome. The target of the minimization functional is to reach a desired temperature distribution at the end of a given time interval (0, T), where T > 0 denotes the final time, with the least possible cost.

Due to the symmetry properties of the geometry as well as the uniform distribution of the water temperature, we reduce the three dimensional (3d) problem to a two dimensional (2d) one. The dimension reduction makes the numerical computations more simple. In the following, the 2d domain is denoted by Ω and its boundary by Γ = ∂Ω. We assume that $$\varOmega \subset \mathbb {R}^2$$ is a bounded Lipschitz domain. We subdivide its boundary Γ into four parts: the glass part Γ 1, the floor which is the swimming pool Γ 2 and the heaters Γ 3 and Γ 4. The domain Ω and its boundaries are illustrated in Fig. 7.1.
../images/475466_1_En_7_Chapter/475466_1_En_7_Fig1_HTML.png
Fig. 7.1

The domain Ω reduced from 3d to 2d due to symmetry properties describing the glass dome and its heaters placed at the ground of the glass dome next to the floor, hence subdividing the boundary Γ = ∂Ω into four parts Γ 1, Γ 2, Γ 3 and Γ 4

7.3 Optimal Control Problem

In this section, the optimal control problem is formulated, where an optimal control function u has to be obtained corresponding to the heating of the heat sources on the boundary Γ R := Γ 3 ∪ Γ 4 such that the state y reaches a desired temperature distribution y d after a given time T. This problem can be formulated in terms of a PDE-constrained optimization problem, which means minimizing a cost functional subject to a PDE and with u being the control function.

Let Q T := Ω × (0, T) denote the space-time cylinder with the lateral surface Σ := Γ × (0, T), where T > 0 denotes the final time. The optimal control problem is given as follows:
$$\displaystyle \begin{aligned} \min_{(y,u)} J(y,u) = \frac{1}{2} \int_{\varOmega} (y(\boldsymbol{x},T) - y_d(\boldsymbol{x}))^2 \, d\boldsymbol{x} + \frac{\lambda}{2} \int_0^T \int_{\varGamma_R} u(x,t)^2 \, ds \, dt\end{aligned} $$
(7.1)
such that
../images/475466_1_En_7_Chapter/475466_1_En_7_Equ2_HTML.png
(7.2)
$$\displaystyle \begin{aligned} \frac{\partial y}{\partial n} &= 0 &&\quad \text{on } \varSigma_1 := \varGamma_1 \times (0,T), {} \end{aligned} $$
(7.3)
$$\displaystyle \begin{aligned} y &= g &&\quad \text{on } \varSigma_2 := \varGamma_2 \times (0,T), {} \end{aligned} $$
(7.4)
$$\displaystyle \begin{aligned} \frac{\partial y}{\partial n} + \alpha y &= \beta u &&\quad \text{on } \varSigma_R := \varGamma_R \times (0,T), {} \end{aligned} $$
(7.5)
$$\displaystyle \begin{aligned} y(0) &= y_0 &&\quad \text{in } \varOmega, {} \end{aligned} $$
(7.6)
where g is the given constant water temperature, λ ≥ 0 is the cost coefficient or control parameter, and α and β are constants describing the heat transfer, which are modeling parameters and have to be chosen carefully. The control u denotes the radiator heating, which has to be chosen within a certain temperature range. Hence, we choose the control functions from the following set of admissible controls:
$$\displaystyle \begin{aligned} u \in U_{\text{ad}} = \{v \in L^2(\varSigma_R) : u_a(\boldsymbol{x},t) \leq u(\boldsymbol{x},t) \leq u_b(\boldsymbol{x},t) \text{ a.e. on } \varSigma_R\}, \end{aligned} $$
(7.7)
which means that u has to fulfill so called box constraints. Equations (7.3)–(7.5) are called Neumann, Dirichlet and Robin boundary conditions, respectively. Equation (7.3) characterizes a no-flux condition in normal direction. Equation (7.4) describes a constant temperature distribution. Regarding Eq. (7.5), α = β would be a reasonable choice from the physical point of view because this would mean that the temperature increase at this part of the boundary is proportional to the difference between the temperature there and outside. However, a decoupling of the parameters makes sense too, see [11], and does not change anything for the actual discussion and computations, since α = β could be chosen at any point.

The goal is to find the optimal set of state and control (y, u) such that the cost is minimal.

7.4 Existence and Uniqueness

In this section, we discuss some basic results on the existence and uniqueness of the parabolic initial-boundary value problem (7.2)–(7.6), whereas we exclude the details. They can be found in [11]. We first introduce proper function spaces leading to a setting, where existence and uniqueness of the solution can be proved.

Definition 7.1
The normed space $$W^{1,0}_2(Q_T)$$ is defined as follows
$$\displaystyle \begin{aligned} W^{1,0}_2(Q_T) = \{ y \in L^2(Q_T) : D_i y \in L^2(Q_T) \, \forall i = 1, \ldots , d \} \end{aligned} $$
(7.8)
with the norm
$$\displaystyle \begin{aligned} \|y\|{}_{W^{1,0}_2(Q_T)} = \left( \int_0^T \int_{\varOmega} (|y(\boldsymbol{x},t)|{}^2 + |\nabla y(\boldsymbol{x},t)|{}^2 ) \, d\boldsymbol{x} \, dt \right)^{1/2}, \end{aligned} $$
(7.9)
where D iy denotes the spatial derivative of y in i-direction and d is the spatial dimension.

For the model problem of this work, the dimension is d = 2. In the following, let {V, ∥⋅∥V} be a real Banach space. More precisely, we will consider V = H 1(Ω) in this work.

Definition 7.2
The space L p(0, T;V ), 1 ≤ p < ∞, denotes the linear space of all equivalence classes of measurable vector valued functions y : [0, T] → V  such that
$$\displaystyle \begin{aligned} \int_0^T \|y(t)\|{}_V^p \, dt &lt; \infty. \end{aligned} $$
(7.10)
The space L p(0, T;V ) is a Banach space with respect to the norm
$$\displaystyle \begin{aligned} \|y\|{}_{L^p(0,T;V)} := \left(\int_0^T \|y(t)\|{}_V^p \, dt \right)^{1/p}. \end{aligned} $$
(7.11)
Definition 7.3
The space W(0, T) = {y ∈ L 2(0, T;V ) : y′∈ L 2(0, T;V ∗)} is equipped with the norm
$$\displaystyle \begin{aligned} \|y\|{}_{W(0,T)} = \left( \int_0^T (|y(t)|{}_V^2 + |y'(t)|{}_{V^*}^2 ) \, dt \right)^{1/2}. \end{aligned} $$
(7.12)
It is a Hilbert space with the scalar product
$$\displaystyle \begin{aligned} (u,w)_{W(0,T)} = \int_0^T (u(t),w(t))_V \, dt + \int_0^T (u'(t),w'(t))_{V^*} \, dt. \end{aligned} $$
(7.13)

The relation V ⊂ H = H ∗⊂ V ∗ is called a Gelfand or evolution triple and describes a chain of dense and continuous embeddings.

The problem (7.2)–(7.6) has a unique weak solution $$y \in W^{1,0}_2(Q_T)$$ for a given u ∈ U ad. Moreover, the solution depends continuously on the data, which means that there exists a constant c > 0 being independent of u, g and y 0 such that
$$\displaystyle \begin{aligned} \max_{t \in [0,T]} \|y(\cdot,t)\|{}_{L^2(\varOmega)} + \|y\|{}_{W^{1,0}_2(Q_T)} \leq c (\|u\|{}_{L^2(\varSigma_R)} + \|g\|{}_{L^2(\varSigma_2)} + \|y_0\|{}_{L^2(\varOmega)}) \end{aligned} $$
(7.14)
for all u ∈ L 2(Σ R), g ∈ L 2(Σ 2) and y 0 ∈ L 2(Ω). Hence, problem (7.2)–(7.6) is well-posed in $$W^{1,0}_2(Q_T)$$ . Furthermore, since $$y \in W^{1,0}_2(Q_T)$$ and it is a weak solution of problem (7.2)–(7.6), y also belongs to W(0, T). The following estimate holds:
$$\displaystyle \begin{aligned} \|y\|{}_{W(0,T)} \leq \tilde c (\|u\|{}_{L^2(\varSigma_R)} + \|g\|{}_{L^2(\varSigma_2)} + \|y_0\|{}_{L^2(\varOmega)}) \end{aligned} $$
(7.15)
for some constant $$\tilde c &gt; 0$$ being independent of u, g and y 0. Hence, problem (7.2)–(7.6) is also well-posed in the space W(0, T).

Note that the Neumann boundary conditions on Γ 1 are included in both estimates (7.14) and (7.15) related to the well-posedness of the problem (as discussed in [11]). However, the Neumann boundary conditions (7.3) are equal to zero.

Under the assumptions that $$\varOmega \subset \mathbb {R}^2$$ is a bounded Lipschitz domain with boundary Γ, λ ≥ 0 is a fixed constant, y d ∈ L 2(Q T), α, β ∈ L ∞(Σ R), and u a, u b ∈ L 2(Σ R) with u a ≤ u b a.e. on Σ R, together with the existence and uniqueness result on the parabolic initial-boundary value problem (7.2)–(7.6) in W(0, T), the optimal control problem (7.1)–(7.7) has at least one optimal control $$\bar u \in U_{\text{ad}}$$ . In case of λ > 0 the optimal control $$\bar u$$ is uniquely determined.

7.5 Reduced Optimization Problem

In order to solve the optimal control problem (7.1)–(7.7), we derive the so called reduced optimization problem first.

Since the problem is well-posed as discussed in the previous section, we can formally eliminate the state equation (7.2)–(7.6) and the minimization problem reads as follows
$$\displaystyle \begin{aligned} \min_{u} \bar{J}(u) = \frac{1}{2} \int_{\varOmega} (y_u(\boldsymbol{x},T) - y_d(\boldsymbol{x}))^2 \, d\boldsymbol{x} + \frac{\lambda}{2} \int_0^T \int_{\varGamma_R} u(x,t)^2 \, ds \, dt \end{aligned} $$
(7.16)
such that (7.7) is satisfied. The problem (7.16) is called reduced optimization problem. Formally, the function y u denotes that the state function is depending on u. However, for simplicity we can set again y u = y. In order to solve the problem (7.16), we apply the projected gradient method. The gradient of $$\bar {J}$$ has to be calculated by deriving the adjoint problem which is given by
../images/475466_1_En_7_Chapter/475466_1_En_7_Equ17_HTML.png
(7.17)
$$\displaystyle \begin{aligned} \frac{\partial p}{\partial n} &amp;= 0 &amp;&amp;\quad \text{on } \varSigma_1, \end{aligned} $$
(7.18)
$$\displaystyle \begin{aligned} p &amp;= 0 &amp;&amp;\quad \text{on } \varSigma_2, \end{aligned} $$
(7.19)
$$\displaystyle \begin{aligned} \frac{\partial p}{\partial n} + \alpha p &amp;= 0 &amp;&amp;\quad \text{on } \varSigma_R, \end{aligned} $$
(7.20)
$$\displaystyle \begin{aligned} p(T) &amp;= y(T) - y_d &amp;&amp;\quad \text{in } \varOmega. \end{aligned} $$
(7.21)
The gradient of $$\bar {J}$$ is given by
$$\displaystyle \begin{aligned} \nabla \bar{J}(u(\boldsymbol{x},t)) = \beta \chi_{\varGamma_R} p(\boldsymbol{x},T-t) + \lambda u(\boldsymbol{x},t) \end{aligned} $$
(7.22)
with $$\chi _{\varGamma _R}$$ denoting the characteristic function on Γ R. The projected gradient method can be now applied for computing the solution of the PDE-constrained optimization problem (7.1)–(7.7). We denote by
$$\displaystyle \begin{aligned} \mathcal{P}_{[u_a,u_b]}(u) = \max\{u_a, \min\{u_b,u\}\} \end{aligned} $$
(7.23)
the projection onto the set of admissible controls U ad.

Now putting everything together, the optimality system for (7.1)–(7.7) and a given λ > 0 reads as follows

../images/475466_1_En_7_Chapter/475466_1_En_7_Equ24_HTML.png
(7.24)
In case that λ = 0, the projection formula changes to
$$\displaystyle \begin{aligned} u(\boldsymbol{x},t) &amp;= u_a(\boldsymbol{x},t), \qquad \text{ if } \beta(\boldsymbol{x},t) p(\boldsymbol{x},t) &gt; 0, \\ u(\boldsymbol{x},t) &amp;= u_b(\boldsymbol{x},t), \qquad \text{ if } \beta(\boldsymbol{x},t) p(\boldsymbol{x},t) &lt; 0. \end{aligned} $$
(7.25)
Remark 7.1

In case that there are no control constraints imposed, the projection formula simplifies to u = −λ −1βp.

7.6 Discretization and Numerical Method

In order to numerically solve the optimal control problem (7.1)–(7.7), which is equivalent to solving (7.24), we discretize the heat equation in space by the finite element method and in time. We use the implicit Euler method for performing the time stepping.

We approximate the functions y, u and p by finite element functions y h, u h and p h from the conforming finite element space V h = span{φ 1, …, φ n} with the basis functions {φ i(x) : i = 1, 2, …, n h}, where h denotes the discretization parameter with n = n h = dim V h = O(h −2). We use standard, continuous, piecewise linear finite elements and a regular triangulation $$\mathcal {T}_h$$ to construct the finite element space V h. For more information, we refer the reader to [3] as well as to the newer publications [2, 7]. Discretizing problem (7.24) by computing its weak formulations and then inserting the finite element approximations for discretizing in space leads to the following discrete formulation:
$$\displaystyle \begin{aligned} M_h \underline{y}_{h,t} + K_h \underline{y}_h + \alpha M_h^{\varGamma_R} \underline{y}_h &amp;= \beta M_h^{\varGamma_R} \underline{u}_h, \qquad \underline{y}_h(0) = y_0, \end{aligned} $$
(7.26)
$$\displaystyle \begin{aligned} -M_h \underline{p}_{h,t} + K_h \underline{p}_h + \alpha M_h^{\varGamma_R} \underline{p}_h &amp;= 0, \qquad \underline{p}_h(T) = \underline{y}_h(T) - y_d, \end{aligned} $$
(7.27)
together with the projection formula
$$\displaystyle \begin{aligned} \underline{u}_h = \mathcal{P}_{[u_a,u_b]}(-\frac{1}{\lambda} \beta \underline{p}_h) \end{aligned} $$
(7.28)
for λ > 0. The problem (7.26)–(7.28) has to be solved with respect to the nodal parameter vectors
$$\displaystyle \begin{aligned} \underline{y}_h = (y_{h,i})_{i=1,\dots,n}, \quad \underline{u}_h = (u_{h,i})_{i=1,\dots,n}, \quad \underline{p}_h = (p_{h,i})_{i=1,\dots,n} \,\, \in \mathbb{R}^n \end{aligned} $$
of the finite element approximations $$y_{h}(\boldsymbol {x}) = \sum _{i=1}^n y_{h,i} \varphi _i(\boldsymbol {x})$$ , $$u _{h}(\boldsymbol {x}) = \sum _{i=1}^n u_{h,i} \varphi _i(\boldsymbol {x})$$ and $$p _{h}(\boldsymbol {x}) = \sum _{i=1}^n p_{h,i} \varphi _i(\boldsymbol {x})$$ . The values for $$ \underline {y}_h$$ are set to g in the nodal values on the boundary Γ 2 and the problems are solved only for the degrees of freedom. The matrices M h, $$M_h^{\varGamma _R}$$ and K h denote the mass matrix, the mass matrix corresponding only to the Robin boundary Γ R and the stiffness matrix, respectively. The entries of the mass and stiffness matrices are defined by the integrals
../images/475466_1_En_7_Chapter/475466_1_En_7_Equb_HTML.png
For the time stepping we use the implicit Euler method. After implementing the finite element discretization and the Euler scheme, we apply the following steps of the projected gradient method:
../images/475466_1_En_7_Chapter/475466_1_En_7_Figa_HTML.png
Remark 7.2

For a first implementation, the step length γ = γ k can be chosen constant for all k. However, a better performance is achieved by applying a line search strategy as for instance the Armijo or Wolfe conditions to obtain the best possible γ k in every iteration step k. We refer the reader to the methods discussed for instance in [5]. However, these strategies are not subject to the present work.

7.7 Numerical Results and Conclusions

In this section, we present numerical results for solving the type of model optimization problem discussed in this article and draw some conclusions in the end. The numerical experiments were computed in Matlab. The meshes were precomputed with Matlab’s pdeModeler. The finite element approximation and time stepping as well as the projected gradient algorithm were implemented according to the discussions in the previous two sections.

In Fig. 7.2, the nodes corresponding to the interior nodes (‘g.’), the Neumann boundary Γ 1 (‘kd’), the Dirichlet boundary Γ 2 (‘bs’) and the Robin boundary Γ R (‘r*’) are illustrated.
../images/475466_1_En_7_Chapter/475466_1_En_7_Fig2_HTML.png
Fig. 7.2

The nodes marked corresponding to different boundaries and the interior of the domain Ω on a mesh with 4073 nodes

In the numerical experiments, we choose the following given data: the water temperature g = 20, the parameters α = β = 102, the final time T = 1, the box constraints u a = 20 and u b = 60, the desired final temperature y d = 30 and the initial value y 0 = 0 satisfying the boundary conditions. For the step lengths γ k of the projected gradient algorithm, we choose the golden ratio γ k = γ = 1.618 constant for all iteration steps k. The stopping criteria include that the norm of the errors $$e_{k+1} := \|u_h^{k+1}-u_h^k\|/\|u_h^k\| &lt; \epsilon _1$$ or |e k+1 − e k| < 𝜖 2 with 𝜖 1 = 10−1 and 𝜖 2 = 10−2 have to be fulfilled as well as setting a maximum number of iteration steps k max = 20 with k < k max.

In the first numerical experiment, we choose a fixed value for the cost coefficient λ = 10−2 and compute the solution for different mesh sizes n ∈{76, 275, 1045, 4073, 16081}. Table 7.1 presents the number of iterations needed until the stopping criteria were satisfied, for different mesh sizes and time steps. The number of time steps was chosen corresponding to the mesh size in order to guarantee that the CFL (Courant-Friedrichs-Lewy) condition is fulfilled.
Table 7.1

Number of iterations needed to satisfy the stopping criteria for different mesh sizes and numbers of time steps for a fixed cost coefficient λ = 10−2

Mesh size

Time steps

Iteration steps

76

125

7

275

250

5

1045

1000

19

4073

4000

19

16,081

16,000

4

In the set of Figs. 7.3, 7.4, 7.5, 7.6, 7.7, and 7.8, the approximate solutions y h defined in Matlab as y are presented for the final time t = T = 1 computed on the different meshes including one figure, Fig. 7.6, presenting the adjoint state p for the mesh with 1045 nodes. We present only one figure for the adjoint state, since for other mesh sizes the plots looked similar.
../images/475466_1_En_7_Chapter/475466_1_En_7_Fig3_HTML.png
Fig. 7.3

The approximate solution y for final time t = T = 1 on a mesh with 76 nodes

../images/475466_1_En_7_Chapter/475466_1_En_7_Fig4_HTML.png
Fig. 7.4

The approximate solution y for final time t = T = 1 on a mesh with 275 nodes

../images/475466_1_En_7_Chapter/475466_1_En_7_Fig5_HTML.png
Fig. 7.5

The approximate solution y for final time t = T = 1 on a mesh with 1045 nodes

../images/475466_1_En_7_Chapter/475466_1_En_7_Fig6_HTML.png
Fig. 7.6

The approximate adjoint state p for time t = 0 on a mesh with 1045 nodes

../images/475466_1_En_7_Chapter/475466_1_En_7_Fig7_HTML.png
Fig. 7.7

The approximate solution y for final time t = T = 1 on a mesh with 4073 nodes

../images/475466_1_En_7_Chapter/475466_1_En_7_Fig8_HTML.png
Fig. 7.8

The approximate solution y for final time t = T = 1 on a mesh with 16,081 nodes

In the second set of numerical experiments, we compute solutions for different cost coefficients λ on two different grids: one with mesh size n = 275 and 250 time steps, and another one with mesh size n = 1045 and 1000 time steps. The numerical results including the number of iteration steps needed are presented in Table 7.2. It can be observed that the numbers of iteration steps for the first case (275/250) are all very similar for different values of λ, even for the lower ones. In the second case (1045/1000), the numbers of iteration steps are getting higher the lower the values of λ. However, the results are satisfactory for these cases too. As example the approximate solution y for the final time t = T = 1 computed for λ = 10−4 is presented in Fig. 7.9. The approximate solution for λ = 10−2 has already been presented in Fig. 7.5.
../images/475466_1_En_7_Chapter/475466_1_En_7_Fig9_HTML.png
Fig. 7.9

The approximate solution y for final time t = T = 1 on a mesh with 1045 nodes for the value λ = 10−4

Table 7.2

Number of iterations needed to satisfy the stopping criteria for different cost coefficients λ ∈{10−4, 10−2, 1, 102, 104} on grids with mesh sizes n = 275 and n = 1045 with 250 and 1000 time steps, respectively

λ

Iteration steps (275/250)

Iteration steps (1045/1000)

10−4

5

19

10−2

5

19

1

7

19

102

4

4

104

4

4

The results of Tables 7.1 and 7.2 were included as example how one can perform different tables for different parameter values or combinations. Students or researchers could compute exactly these different kinds of numerical experiments in order to study the practical performance of the optimization problem.

After presenting the numerical results, we have to mention again that the step length γ = γ k of the projected gradient method has been chosen constant for all k in all computations. However, better results should be achieved by applying a suitable line search strategy, see Remark 7.2. With this we want to conclude that the optimal control problem discussed in this work is one model formulation for solving the optimization of heating a domain such as the air of a swimming pool area surrounded by a glass dome. Modeling and solving the optimal heating of a swimming pool area has potential for many different formulations related to mathematical modeling, discussing different solution methods and performing many numerical tests including “playing around” with different values for the parameters, constants and given functions, and finally choosing proper ones for the model problem. All of these tasks can be performed by students depending also on their previous knowledge, interests and ideas.

Acknowledgements

I would like to thank my students T. Bazlyankov, T. Briffard, G. Krzyzanowski, P.-O. Maisonneuve and C. Neßler for their work at the 26th ECMI Modelling Week which provided part of the starting point for this work. I gratefully acknowledge the financial support by the Academy of Finland under the grant 295897.