# On Nonlinear Inverse Problems of Heat Transfer with Radiation Boundary Conditions: Application to Dehydration of Gypsum Plasterboards Exposed to Fire

Authors: A. Belmiloudi; F. Mahé
Journal: Advances in Numerical Analysis
 (2014)

---

## Abstract

The paper investigates boundary optimal controls and parameter estimates to the well-posedness nonlinear model of dehydration of thermic problems. We summarize the general formulations for the boundary control for initial-boundary value problem for nonlinear partial differential equations modeling the heat transfer and derive necessary optimality conditions, including the adjoint equation, for the optimal set of parameters minimizing objective functionsJ. Numerical simulations illustrate several numerical optimization methods, examples, and realistic cases, in which several interesting phenomena are observed. A large amount of computational effort is required to solve the coupled state equation and the adjoint equation (which is backwards in time), and the algebraic gradient equation (which implements the coupling between the adjoint and control variables). The state and adjoint equations are solved using the finite element method.


---

## Body

## 1. Introduction and Mathematical Setting of the Problem

### 1.1. Motivation

Since several years, a considerable effort has been made to develop materials having a good fire resistance. Such materials must provide a sufficient mechanical resistance to avoid the premature collapse of a building structure undergoing a fire. Consequently this type of material must withstand significant heating without burning and keep its mechanical resistance sufficient. Criteria which permit appreciating the fire resistance of materials are given by several norms which define the minimum fire exposure duration that must support the building structure.One of the building materials presenting the best fire resistance is gypsum plasterboard, which in turn is due to the dehydration phenomenon. This material presents the particularity to undergoing two chemical reactions of dehydration during its heating. These two endothermic reactions considerably slow down the heating of the material (since the dehydration process consumes large amount of heat) and provide the plasterboard excellent fire resistance.The main particularity of gypsum plasterboard is that it contains21
% of chemically bound water by weight. When the temperature reaches 100°C in a point of the plasterboard, a reaction of dehydration occurs in the material. This chemical reaction dissociates a certain quantity of water which is combined to the crystal lattice. In terms of fire safety, the reaction of dehydration and the vaporization of free water absorb a certain amount of energy which significantly slows down the heating of the material and in particular the temperature rise on the unexposed side of plasterboard.Although necessary, experimental testing is not a convenient way to estimate the fire resistance of a material. Indeed, full-scale testing poses the problem of the high cost of the experimental setup and the difficulty to implement the experiment. In addition pilot-scale testing does not allow to accurately reproduce the real conditions of a fire exposure. Consequently, the development of a mathematical model and the numerical simulation of the heating of gypsum plasterboard exposed to fire appear as a suitable means to study the thermal behaviour of the material during a fire exposure.


### 1.2. Thermochemistry of Gypsum

Gypsum plasterboard is commonly used as construction material to improve fire resistance of building structures. The pure gypsum, existing at the natural state as a more or less compact rock, is composed of calcium sulphate with free water at equilibrium moisture (approximately3
%) and approximately 20
% per weight of chemically combined water of crystallization (see, e.g., [1–3]). Its chemical formula is CaSO4·2H2O (calcium sulphate di-hydrate). The industrials add various chemical elements (in small quantities) in order to increase their performance when exposed to elevated temperatures. The chemical reaction which consists in removing chemically combined water of crystallization is called calcination. During heating, gypsum plaster undergoes two endothermic decomposition reactions. The first dehydration reaction occurs at approximately 100–120°C when the calcium sulphate di-hydrate is converted to calcium sulphate hemihydrate (the reaction is always complete by 160°C) as shown by the following reaction:
(1)
CaSO
4
·
2
H
2
O
⟶
CaSO
4
·
1
2
H
2
O
+
3
2
H
2
O
The amount of energy required by this first dehydration is about 500 kJ per kg of gypsum (see [4]).The second dehydration reaction occurs when calcium sulphate hemihydrate is converted to calcium sulphate anhydrate as shown by the following reaction:(2)
CaSO
4
·
1
2
H
2
O
⟶
CaSO
4
+
1
2
H
2
O
The amount of energy corresponding to this second reaction is about 169 kJ per kg (see [5]).Remark 1.
Other amounts of energy can be found in the literature for these reactions; see [6] for a review.Both reactions are endothermic, produce water vapour, and absorb a large amount of energy. The effect of the endothermic reactions on the heating of the wall of plasterboard is taken into account by including the latent heats of reactions (1) and (2) in the specific heat evolution. The first dehydration reaction occurs at approximately 100–160°C; on the other hand, there is some controversy to when the second dehydration reaction occurs. Andersson et al. [7] (e.g.,) estimate that the second reaction occurs between 210°C and 300°C. That consists in introducing two peaks in the evolution of the specific heat according to the temperature, corresponding to the temperatures to which the reactions occur. The areas under the two peaks are equal to the latent heats of the two chemical reactions. Other experiences show that this second reaction occurs immediately after the first one [8]. In the numerical examples, we will choose this model with only one peak between 100 and 170°C in the evolution of the specific heat. The information on the thermophysical properties of gypsum plasterboard, at high temperatures, are difficult to measure and then are limited, because the derived results are always complicated by the dynamic nature of the (fire resistive) materials and vary considerably with the used method of measurement (a wide variety of experimental techniques exists for measuring these properties) and the rate of temperature change (for more details see, e.g., [6, 9–11]).Remark 2.
The model of dehydration can be completed by other reactions at high temperatures. For example, a third endothermic reaction, for decomposition of CaCO3, is introduced at 650°C in [8] or 800°C in [12]. In [13], the authors note an exothermic reaction around 400°C.


### 1.3. Outline of the Paper

The paper is organized as follows. In the next subsection, we give a sketch of the modeling leading to problem and we establish the governing equations. The model is presented for dehydration of gypsum plasterboards exposed to fire but the model is general and the method is valid for other applications. In Section2 we give a description of the parameter estimates (identification problems) as nonlinear optimal control problems with boundary control. This includes results concerning the existence of the optimal solutions, necessary optimality conditions (necessary to develop numerical optimization methods), the optimization problem, and adjoint model. Section 3 contains details of the computational algorithm and numerical simulation optimizations of the optimal control problems. Numerical results for several examples are presented and a realistic situation is analyzed. Section 4 contains a summary and a discussion of future work.


### 1.4. Modeling of the Wall of Plasterboard Heating and the Direct Forward Model

This section is devoted to an introduction of the derivation of dehydration of gypsum plasterboards (exposed to fire) model that we study. It is well known that the problem of heat and mass transfer in plasterboard exposed to fire is essentially one-dimensional, so the model is derived in one-dimensional formulation.Let us consider a wall of plasterboard exposed to fire, which is located vertically on a retaining structural frame. The left hand side is exposed to a heat source, as may occur in furnace in which fire tests are conducted. We suppose that the depth and the width of the wall are much bigger than the thickness (L). Therefore, heat fluxes in lateral (y) direction and vertical (z) direction can be neglected in front of the heat flux in direction x. On the other hand, the heat source is distributed uniformly on the heated side of the plasterboard. Applying these physical considerations, the heat transfer can be treated as 1-dimensional process in x-direction. The one-way heat transfer through a plane wall is described by the heat equation in its one-dimensional form
(3)
λ
(
u
)
∂
u
∂
t
(
x
,
t
)
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
(
x
,
t
)
=
0
,
llllllllllllllllllllllllllllllllll
in
𝒬
=
Ω
×
(
0
,
T
)
,
u
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
Ω
,
where Ω is the boundary subset (
0
,
L
), u is the temperature, λ
=
ρ
c, ρ is the density of body material, c is the specific heat, and κ is the thermal conductivity. The functions λ and κ are variables, positive, and bounded.The external surface of the wall of plasterboard exposed to fire receives a heat flux which consists of convective and radiative components. Consequently, the boundary condition on this side is written, fort in (
0
,
T
), as
(4)
-
κ
(
u
(
0
,
t
)
)
∂
u
∂
x
(
0
,
t
)
=
h
0
(
u
f
(
t
)
-
u
(
0
,
t
)
)
+
φ
0
(
|
u
f
(
t
)
|
3
u
f
(
t
)
-
|
u
(
0
,
t
)
|
3
u
(
0
,
t
)
)
,
where, for all t
∈
[
0
,
T
], h
0
(
t
)
>
0 is the convective heat transfer coefficient between the furnace and the plasterboard surface, u
f
(
t
)
>
0 is the furnace temperature, φ
0
(
t
)
=
σ
ϵ
0
(
t
)
>
0, σ
=
5.67
·
1
0
-
8
W
/
m
2
K
4 is Stefan-Boltzmann's constant, and ϵ
0 is effective emissivity of the surface.The external surface of the plasterboard, which is not exposed to fire, transfers heat to the surroundings by means of convection and radiation. As the previous case, the boundary condition on this side is written, fort in (
0
,
T
), as
(5)
κ
(
u
(
L
,
t
)
)
∂
u
∂
x
(
L
,
t
)
=
h
1
(
u
a
(
t
)
-
u
(
L
,
t
)
)
+
φ
1
(
|
u
a
(
t
)
|
3
u
a
(
t
)
-
|
u
(
L
,
t
)
|
3
u
(
L
,
t
)
)
,
where, for all t
∈
[
0
,
T
], h
1
(
t
)
>
0 is convective heat transfer coefficient between the surroundings and the plasterboard surface, u
a
(
t
)
>
0 is the surroundings temperature, φ
1
(
t
)
=
σ
ϵ
1
(
t
)
>
0, and ϵ
1 is effective emissivity of the surface.Remark 3.
In the sequel, to simplify the notations, we denoteu for u
(
x
,
t
), u
(
0
) for u
(
0
,
t
), u
(
L
) for u
(
L
,
t
), u
(
T
) for u
(
x
,
T
), u
f for u
f
(
t
), and u
a for u
a
(
t
).We assume that(H1)
functionsf
=
(
h
0
,
h
1
)
∈
(
L
2
(
0
,
T
)
)
2 and ϕ
=
(
φ
0
,
φ
1
)
∈
(
L
2
(
0
,
T
)
)
2 and satisfy the following pointwise constraint:
(6)
0
≤
a
i
≤
φ
i
≤
b
i
a
.
e
.
t
∈
(
0
,
T
)
,
for
i
=
0,1
,
0
≤
c
i
≤
h
i
≤
d
i
a
.
e
.
t
∈
(
0
,
T
)
,
for
i
=
0,1
,for some positive constantsa
i
,
b
i
,
c
i
,
d
i
,
i
=
0,1,(H2)
operatorsκ
(
v
), λ
(
v
) are sufficiently regular with 0
<
K
0
≤
κ
(
v
)
≤
K
1 and with 0
<
C
0
≤
λ
(
v
)
≤
C
1 for any v
>
0.Remark 4.
(
1
) Emissivity of a material is defined as the ratio of energy radiated by a particular material to energy radiated by a black body at the same temperature. It is a dimensionless quantity (i.e., a quantity without a physical unit).
(
2
) In the physical case there are not absolute values under the boundary conditions (since the temperature is nonnegative). For real physical data and operators λ and κ, we can prove by using the maximum principle that the temperature is positive and then we can remove the absolute values. It is useful for the mathematical study of the problem.
(
3
) It is clear that we can obtain in the same way the model in N-dimensional for N
≤
3 as follows:
(7)
λ
(
u
)
∂
u
∂
t
-
div
⁡
(
κ
(
u
)
∇
u
)
=
0
,
in
𝒬
,
u
(
·
,
0
)
=
u
0
,
in
Ω
,
κ
(
u
)
∇
u
·
n
=
h
0
(
u
f
-
u
)
+
φ
0
(
|
u
f
|
3
u
f
-
|
u
|
3
u
)
,
llllllllllllllllllllllllllllll
in
Σ
f
=
Γ
f
×
(
0
,
T
)
,
κ
(
u
)
∇
u
·
n
=
h
1
(
u
a
-
u
)
+
φ
1
(
|
u
a
|
3
u
a
-
|
u
|
3
u
)
,
lllllllllllllllllllllllllllll
in
Σ
c
=
Γ
c
×
(
0
,
T
)
,
u
=
u
a
,
in
Σ
s
=
Γ
s
×
(
0
,
T
)
,
where Ω
⊂
ℝ
N is an open bounded domain with boundary Γ
=
∂
Ω, Γ
=
Γ
c
∪
Γ
f
∪
Γ
s such that Γ
c
∩
Γ
f
∩
Γ
s
=
∅, and n is the outward unit normal vector on Γ. Boundaries Γ
c and Γ
f denote the cold side and the fire side, respectively, and Γ
s denotes the other surface of the plasterboard. The well-posedness of problem (7) in 3-dimensional can be obtained in similar way as in [14].A priori, most researchers who have worked on the modeling of the behavior of gypsum board (literature in the public domain in this field is sparse; see, e.g., [3–5, 7, 11, 15]), have assumed the convective heat transfer coefficients h
0 and h
1 and the relative emissivity ϵ
0 as constants or/and neglected the relative emissivity ϵ
1 on cold surface. The choice of a constant for these coefficients is not a good physical representation of plasterboard exposed to fire. Moreover, the convective heat transfer coefficients and emissivities depend, among other things, on the state of the external surfaces. During the exposition to fire, mechanical resistance of the external surfaces decreases and that causes appearance and growing of crazes (degradation). These modifications of surface states of the external surfaces modify convective and radiative heat transfers. Consequently, the resultant emissivities and convective heat transfer coefficients depend on temperature and large uncertainties exist in regard to the quality of the data reported. Moreover, the work of Belmiloudi and Le Meur [16] shows very clearly that the radiative heat transfer between the unexposed surface and the surrounding cannot be neglected. Then, it is necessary to estimate the convective heat transfers and the emissivity coefficients.To satisfy this requirement, we estimate these parameters by using inverse problem techniques as optimal control methods. It is clear that the accuracy of the parameter estimate from furnace, fire test data, and target observations (or measurement results) depends significantly on the thermal characteristics of a furnace, on the geometry of the studied element, and on the input thermal properties of the material. So, it is important to have a consistent set of values for these data.


## 1.1. Motivation

Since several years, a considerable effort has been made to develop materials having a good fire resistance. Such materials must provide a sufficient mechanical resistance to avoid the premature collapse of a building structure undergoing a fire. Consequently this type of material must withstand significant heating without burning and keep its mechanical resistance sufficient. Criteria which permit appreciating the fire resistance of materials are given by several norms which define the minimum fire exposure duration that must support the building structure.One of the building materials presenting the best fire resistance is gypsum plasterboard, which in turn is due to the dehydration phenomenon. This material presents the particularity to undergoing two chemical reactions of dehydration during its heating. These two endothermic reactions considerably slow down the heating of the material (since the dehydration process consumes large amount of heat) and provide the plasterboard excellent fire resistance.The main particularity of gypsum plasterboard is that it contains21
% of chemically bound water by weight. When the temperature reaches 100°C in a point of the plasterboard, a reaction of dehydration occurs in the material. This chemical reaction dissociates a certain quantity of water which is combined to the crystal lattice. In terms of fire safety, the reaction of dehydration and the vaporization of free water absorb a certain amount of energy which significantly slows down the heating of the material and in particular the temperature rise on the unexposed side of plasterboard.Although necessary, experimental testing is not a convenient way to estimate the fire resistance of a material. Indeed, full-scale testing poses the problem of the high cost of the experimental setup and the difficulty to implement the experiment. In addition pilot-scale testing does not allow to accurately reproduce the real conditions of a fire exposure. Consequently, the development of a mathematical model and the numerical simulation of the heating of gypsum plasterboard exposed to fire appear as a suitable means to study the thermal behaviour of the material during a fire exposure.


## 1.2. Thermochemistry of Gypsum

Gypsum plasterboard is commonly used as construction material to improve fire resistance of building structures. The pure gypsum, existing at the natural state as a more or less compact rock, is composed of calcium sulphate with free water at equilibrium moisture (approximately3
%) and approximately 20
% per weight of chemically combined water of crystallization (see, e.g., [1–3]). Its chemical formula is CaSO4·2H2O (calcium sulphate di-hydrate). The industrials add various chemical elements (in small quantities) in order to increase their performance when exposed to elevated temperatures. The chemical reaction which consists in removing chemically combined water of crystallization is called calcination. During heating, gypsum plaster undergoes two endothermic decomposition reactions. The first dehydration reaction occurs at approximately 100–120°C when the calcium sulphate di-hydrate is converted to calcium sulphate hemihydrate (the reaction is always complete by 160°C) as shown by the following reaction:
(1)
CaSO
4
·
2
H
2
O
⟶
CaSO
4
·
1
2
H
2
O
+
3
2
H
2
O
The amount of energy required by this first dehydration is about 500 kJ per kg of gypsum (see [4]).The second dehydration reaction occurs when calcium sulphate hemihydrate is converted to calcium sulphate anhydrate as shown by the following reaction:(2)
CaSO
4
·
1
2
H
2
O
⟶
CaSO
4
+
1
2
H
2
O
The amount of energy corresponding to this second reaction is about 169 kJ per kg (see [5]).Remark 1.
Other amounts of energy can be found in the literature for these reactions; see [6] for a review.Both reactions are endothermic, produce water vapour, and absorb a large amount of energy. The effect of the endothermic reactions on the heating of the wall of plasterboard is taken into account by including the latent heats of reactions (1) and (2) in the specific heat evolution. The first dehydration reaction occurs at approximately 100–160°C; on the other hand, there is some controversy to when the second dehydration reaction occurs. Andersson et al. [7] (e.g.,) estimate that the second reaction occurs between 210°C and 300°C. That consists in introducing two peaks in the evolution of the specific heat according to the temperature, corresponding to the temperatures to which the reactions occur. The areas under the two peaks are equal to the latent heats of the two chemical reactions. Other experiences show that this second reaction occurs immediately after the first one [8]. In the numerical examples, we will choose this model with only one peak between 100 and 170°C in the evolution of the specific heat. The information on the thermophysical properties of gypsum plasterboard, at high temperatures, are difficult to measure and then are limited, because the derived results are always complicated by the dynamic nature of the (fire resistive) materials and vary considerably with the used method of measurement (a wide variety of experimental techniques exists for measuring these properties) and the rate of temperature change (for more details see, e.g., [6, 9–11]).Remark 2.
The model of dehydration can be completed by other reactions at high temperatures. For example, a third endothermic reaction, for decomposition of CaCO3, is introduced at 650°C in [8] or 800°C in [12]. In [13], the authors note an exothermic reaction around 400°C.


## 1.3. Outline of the Paper

The paper is organized as follows. In the next subsection, we give a sketch of the modeling leading to problem and we establish the governing equations. The model is presented for dehydration of gypsum plasterboards exposed to fire but the model is general and the method is valid for other applications. In Section2 we give a description of the parameter estimates (identification problems) as nonlinear optimal control problems with boundary control. This includes results concerning the existence of the optimal solutions, necessary optimality conditions (necessary to develop numerical optimization methods), the optimization problem, and adjoint model. Section 3 contains details of the computational algorithm and numerical simulation optimizations of the optimal control problems. Numerical results for several examples are presented and a realistic situation is analyzed. Section 4 contains a summary and a discussion of future work.


## 1.4. Modeling of the Wall of Plasterboard Heating and the Direct Forward Model

This section is devoted to an introduction of the derivation of dehydration of gypsum plasterboards (exposed to fire) model that we study. It is well known that the problem of heat and mass transfer in plasterboard exposed to fire is essentially one-dimensional, so the model is derived in one-dimensional formulation.Let us consider a wall of plasterboard exposed to fire, which is located vertically on a retaining structural frame. The left hand side is exposed to a heat source, as may occur in furnace in which fire tests are conducted. We suppose that the depth and the width of the wall are much bigger than the thickness (L). Therefore, heat fluxes in lateral (y) direction and vertical (z) direction can be neglected in front of the heat flux in direction x. On the other hand, the heat source is distributed uniformly on the heated side of the plasterboard. Applying these physical considerations, the heat transfer can be treated as 1-dimensional process in x-direction. The one-way heat transfer through a plane wall is described by the heat equation in its one-dimensional form
(3)
λ
(
u
)
∂
u
∂
t
(
x
,
t
)
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
(
x
,
t
)
=
0
,
llllllllllllllllllllllllllllllllll
in
𝒬
=
Ω
×
(
0
,
T
)
,
u
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
Ω
,
where Ω is the boundary subset (
0
,
L
), u is the temperature, λ
=
ρ
c, ρ is the density of body material, c is the specific heat, and κ is the thermal conductivity. The functions λ and κ are variables, positive, and bounded.The external surface of the wall of plasterboard exposed to fire receives a heat flux which consists of convective and radiative components. Consequently, the boundary condition on this side is written, fort in (
0
,
T
), as
(4)
-
κ
(
u
(
0
,
t
)
)
∂
u
∂
x
(
0
,
t
)
=
h
0
(
u
f
(
t
)
-
u
(
0
,
t
)
)
+
φ
0
(
|
u
f
(
t
)
|
3
u
f
(
t
)
-
|
u
(
0
,
t
)
|
3
u
(
0
,
t
)
)
,
where, for all t
∈
[
0
,
T
], h
0
(
t
)
>
0 is the convective heat transfer coefficient between the furnace and the plasterboard surface, u
f
(
t
)
>
0 is the furnace temperature, φ
0
(
t
)
=
σ
ϵ
0
(
t
)
>
0, σ
=
5.67
·
1
0
-
8
W
/
m
2
K
4 is Stefan-Boltzmann's constant, and ϵ
0 is effective emissivity of the surface.The external surface of the plasterboard, which is not exposed to fire, transfers heat to the surroundings by means of convection and radiation. As the previous case, the boundary condition on this side is written, fort in (
0
,
T
), as
(5)
κ
(
u
(
L
,
t
)
)
∂
u
∂
x
(
L
,
t
)
=
h
1
(
u
a
(
t
)
-
u
(
L
,
t
)
)
+
φ
1
(
|
u
a
(
t
)
|
3
u
a
(
t
)
-
|
u
(
L
,
t
)
|
3
u
(
L
,
t
)
)
,
where, for all t
∈
[
0
,
T
], h
1
(
t
)
>
0 is convective heat transfer coefficient between the surroundings and the plasterboard surface, u
a
(
t
)
>
0 is the surroundings temperature, φ
1
(
t
)
=
σ
ϵ
1
(
t
)
>
0, and ϵ
1 is effective emissivity of the surface.Remark 3.
In the sequel, to simplify the notations, we denoteu for u
(
x
,
t
), u
(
0
) for u
(
0
,
t
), u
(
L
) for u
(
L
,
t
), u
(
T
) for u
(
x
,
T
), u
f for u
f
(
t
), and u
a for u
a
(
t
).We assume that(H1)
functionsf
=
(
h
0
,
h
1
)
∈
(
L
2
(
0
,
T
)
)
2 and ϕ
=
(
φ
0
,
φ
1
)
∈
(
L
2
(
0
,
T
)
)
2 and satisfy the following pointwise constraint:
(6)
0
≤
a
i
≤
φ
i
≤
b
i
a
.
e
.
t
∈
(
0
,
T
)
,
for
i
=
0,1
,
0
≤
c
i
≤
h
i
≤
d
i
a
.
e
.
t
∈
(
0
,
T
)
,
for
i
=
0,1
,for some positive constantsa
i
,
b
i
,
c
i
,
d
i
,
i
=
0,1,(H2)
operatorsκ
(
v
), λ
(
v
) are sufficiently regular with 0
<
K
0
≤
κ
(
v
)
≤
K
1 and with 0
<
C
0
≤
λ
(
v
)
≤
C
1 for any v
>
0.Remark 4.
(
1
) Emissivity of a material is defined as the ratio of energy radiated by a particular material to energy radiated by a black body at the same temperature. It is a dimensionless quantity (i.e., a quantity without a physical unit).
(
2
) In the physical case there are not absolute values under the boundary conditions (since the temperature is nonnegative). For real physical data and operators λ and κ, we can prove by using the maximum principle that the temperature is positive and then we can remove the absolute values. It is useful for the mathematical study of the problem.
(
3
) It is clear that we can obtain in the same way the model in N-dimensional for N
≤
3 as follows:
(7)
λ
(
u
)
∂
u
∂
t
-
div
⁡
(
κ
(
u
)
∇
u
)
=
0
,
in
𝒬
,
u
(
·
,
0
)
=
u
0
,
in
Ω
,
κ
(
u
)
∇
u
·
n
=
h
0
(
u
f
-
u
)
+
φ
0
(
|
u
f
|
3
u
f
-
|
u
|
3
u
)
,
llllllllllllllllllllllllllllll
in
Σ
f
=
Γ
f
×
(
0
,
T
)
,
κ
(
u
)
∇
u
·
n
=
h
1
(
u
a
-
u
)
+
φ
1
(
|
u
a
|
3
u
a
-
|
u
|
3
u
)
,
lllllllllllllllllllllllllllll
in
Σ
c
=
Γ
c
×
(
0
,
T
)
,
u
=
u
a
,
in
Σ
s
=
Γ
s
×
(
0
,
T
)
,
where Ω
⊂
ℝ
N is an open bounded domain with boundary Γ
=
∂
Ω, Γ
=
Γ
c
∪
Γ
f
∪
Γ
s such that Γ
c
∩
Γ
f
∩
Γ
s
=
∅, and n is the outward unit normal vector on Γ. Boundaries Γ
c and Γ
f denote the cold side and the fire side, respectively, and Γ
s denotes the other surface of the plasterboard. The well-posedness of problem (7) in 3-dimensional can be obtained in similar way as in [14].A priori, most researchers who have worked on the modeling of the behavior of gypsum board (literature in the public domain in this field is sparse; see, e.g., [3–5, 7, 11, 15]), have assumed the convective heat transfer coefficients h
0 and h
1 and the relative emissivity ϵ
0 as constants or/and neglected the relative emissivity ϵ
1 on cold surface. The choice of a constant for these coefficients is not a good physical representation of plasterboard exposed to fire. Moreover, the convective heat transfer coefficients and emissivities depend, among other things, on the state of the external surfaces. During the exposition to fire, mechanical resistance of the external surfaces decreases and that causes appearance and growing of crazes (degradation). These modifications of surface states of the external surfaces modify convective and radiative heat transfers. Consequently, the resultant emissivities and convective heat transfer coefficients depend on temperature and large uncertainties exist in regard to the quality of the data reported. Moreover, the work of Belmiloudi and Le Meur [16] shows very clearly that the radiative heat transfer between the unexposed surface and the surrounding cannot be neglected. Then, it is necessary to estimate the convective heat transfers and the emissivity coefficients.To satisfy this requirement, we estimate these parameters by using inverse problem techniques as optimal control methods. It is clear that the accuracy of the parameter estimate from furnace, fire test data, and target observations (or measurement results) depends significantly on the thermal characteristics of a furnace, on the geometry of the studied element, and on the input thermal properties of the material. So, it is important to have a consistent set of values for these data.


## 2. The Inverse Problem Formulation

### 2.1. Problem Formulation

We assume that there exists a unique solutionu of problem (3), (6) with boundary conditions (4), (5) under some hypotheses for the data and some regularity of the operators λ and κ, satisfying the following regularity (by using [16, 17]):
(8)
u
∈
L
2
(
0
,
T
;
H
1
(
Ω
)
)
∩
L
∞
(
0
,
T
;
L
2
(
Ω
)
)
,
∂
u
∂
t
∈
L
2
(
0
,
T
;
(
H
1
(
Ω
)
)
′
)
,
|
u
(
0
)
|
3
u
(
0
)
∈
L
5
/
4
(
0
,
T
)
,
|
u
(
L
)
|
3
u
(
L
)
∈
L
5
/
4
(
0
,
T
)
.Introduce now the mappingℱ which maps the source term h of (3)–(5) into the corresponding solution u
=
ℱ
(
h
), where h
=
(
h
0
,
h
1
). In this section we formulate the inverse problem as an optimal control problem. The control procedure consists of finding the optimal controls h
opt
=
(
h
0
opt
,
h
1
opt
) and the corresponding optimal temperature u
opt
=
ℱ
(
h
opt
) which minimize a cost criterion J
(
h
). The cost functional J measures the distance between a measured temperature (the observations u
obs, u
obs
0,    u
obs
L, and u
obs
T) and the corresponding predicted temperature obtained from the primal (or direct) model (3)–(5). Precisely we will study the following optimal control problem (
𝒫
c
).Find(
u
,
h
) such that the following objective function:
(9)
J
(
h
)
=
1
2
∫
‍
∫
𝒬
‍
|
𝒞
(
u
-
u
obs
)
|
2
d
x
d
t
+
1
2
∫
Ω
‍
|
𝒟
(
u
(
T
)
-
u
obs
T
)
|
2
d
x
+
1
2
∫
0
T
‍
[
γ
0
|
u
(
0
)
-
u
obs
0
|
2
+
γ
1
|
u
(
L
)
-
u
obs
L
|
2
]
d
t
+
1
2
∫
0
T
‍
|
ℳ
h
|
2
d
tis minimized with respect toh
∈
𝒴
ad subject to (3)–(5), where 𝒴
ad is the set of admissible controls, u
=
ℱ
(
h
), γ
0 and γ
1 are predefined nonnegative weights, and ℳ
=
diag
⁡
(
α
0
,
α
1
) are predefined nonnegative weights such that α
0
+
α
1
≠
0. The operators 𝒞 and 𝒟 are unbounded operators on L
2
(
Ω
) satisfying (∀
v
∈
H
1
(
Ω
)):
(10)
∥
𝒞
v
∥
L
2
(
Ω
)
2
≤
δ
1
∥
v
∥
L
2
(
Ω
)
2
+
δ
2
∥
v
∥
H
1
(
Ω
)
2
,
∥
𝒟
v
∥
L
2
(
Ω
)
2
≤
δ
1
∥
v
∥
L
2
(
Ω
)
2
+
δ
2
∥
v
∥
H
1
(
Ω
)
2
,
with δ
1
+
δ
2
>
0 and δ
i
≥
0, for i
=
1,2.Remark 5.
The general form (9) of the objective function allows considering different kinds of problems (we denote I as the identity and d as a positive scalar): (1)
distributed observations control problems with𝒞
=
d
I, γ
0
=
γ
1
=
0, and 𝒟
=
0,
(2)
boundary observations control problems with𝒞
=
0, γ
0
=
γ
1
=
d, and 𝒟
=
0,
(3)
final observations control problems with𝒞
=
0, γ
0
=
γ
1
=
0, and 𝒟
=
d
I,
(4)
single control function, for example,h
1, with α
1
≠
0 and α
0
=
0.
We consider here an optimal control problem with Tichonov regularization terms [18]. The nonnegative weights α
0 and α
1 play the role of regularization parameters and take small values in the numerical simulations.According to (6) the set of admissible controls describing the constraint is
(11)
𝒴
ad
=
{
h
=
(
h
0
,
h
1
)
∈
(
L
2
(
0
,
T
)
)
2
such
that
(
L
2
(
0
,
T
)
)
2
0
≤
c
i
≤
h
i
≤
d
i
a
.
e
.
in
(
0
,
T
)
,
for
i
=
0,1
}
.


### 2.2. First-Order Necessary Conditions

Assume that the nonlinear control problem (9) admits an optimal solution (for similar result see, e.g., [19, 20]); the necessary conditions for this optimum is given by the following theorems (see [19]).Theorem 6.
IfJ attains a (local) minimum at a point h
*
∈
𝒴
a
d, then the following first optimality conditions hold:
(12)
J
′
(
h
*
)
·
(
h
-
h
*
)
≥
0
,
∀
h
∈
𝒴
a
d
,
where J
′ is the directional derivative of J.In order to solve numerically the optimal control problem it is necessary to derive the gradient of the cost functionalJ with respect to the control h. For this we suppose that the operator solution ℱ is continuously differentiable on 𝒴
ad and its derivative w
=
ℱ
′
(
h
)
·
g
=
lim
⁡
ϵ
→
0
(
(
ℱ
(
h
+
ϵ
g
)
-
ℱ
(
h
)
)
/
ϵ
) is the unique solution of the following system (for h
=
(
h
0
,
h
1
) and g
=
(
g
0
,
g
1
)):
(13)
λ
(
u
)
∂
w
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
w
∂
x
)
+
F
(
u
,
w
)
=
0
,
in
𝒬
,
-
κ
(
u
(
0
)
)
∂
w
∂
x
(
0
)
=
-
[
h
0
-
G
(
u
)
(
0
)
+
4
φ
0
|
u
(
0
)
|
3
]
w
(
0
)
+
g
0
(
u
f
-
u
(
0
)
)
,
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
w
∂
x
(
L
)
=
-
[
h
1
+
G
(
u
)
(
L
)
+
4
φ
1
|
u
(
L
)
|
3
]
w
(
L
)
+
g
1
(
u
a
-
u
(
L
)
)
,
in
(
0
,
T
)
,
w
(
x
,
0
)
=
0
,
x
∈
Ω
,where operatorsF and G are defined as follows:(14)
F
(
u
,
w
)
=
∂
∂
t
(
λ
(
u
)
)
w
-
∂
∂
x
(
w
∂
∂
x
(
κ
(
u
)
)
)
,
G
(
u
)
=
∂
∂
x
(
κ
(
u
)
)
.Remark 7.
We verify easily that(15)
λ
(
u
)
∂
w
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
w
∂
x
)
+
F
(
u
,
w
)
=
∂
(
λ
(
u
)
w
)
∂
t
-
∂
2
∂
x
2
(
κ
(
u
)
w
)
,
κ
(
u
(
0
)
)
∂
w
∂
x
(
0
)
+
G
(
u
)
(
0
)
w
(
0
)
=
∂
(
κ
(
u
)
w
)
∂
x
(
0
)
,
κ
(
u
(
L
)
)
∂
w
∂
x
(
L
)
+
G
(
u
)
(
L
)
w
(
L
)
=
∂
(
κ
(
u
)
w
)
∂
x
(
L
)
.In order to derive (13), we write the two systems satisfied by u
ϵ
=
ℱ
(
h
+
ϵ
g
) and u
=
ℱ
(
h
) as
(16)
λ
(
u
ϵ
)
∂
u
ϵ
∂
t
-
∂
∂
x
(
κ
(
u
ϵ
)
∂
u
ϵ
∂
x
)
=
0
,
in
𝒬
,
-
κ
(
u
ϵ
(
0
)
)
∂
u
ϵ
∂
x
(
0
)
=
(
h
0
+
ϵ
g
0
)
(
u
f
-
u
ϵ
(
0
)
)
+
φ
0
(
|
u
f
|
3
u
f
-
|
u
ϵ
(
0
)
|
3
u
ϵ
(
0
)
)
,
+
φ
0
|
u
f
|
3
u
f
-
|
u
ϵ
(
0
)
|
3
in
(
0
,
T
)
,
κ
(
u
ϵ
(
L
)
)
∂
u
ϵ
∂
x
(
L
)
=
(
h
1
+
ϵ
g
1
)
(
u
a
-
u
ϵ
(
L
)
)
+
φ
1
(
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
u
ϵ
(
L
)
)
,
+
φ
1
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
in
(
0
,
T
)
,
u
ϵ
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
Ω
,
(17)
λ
(
u
)
∂
u
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
=
0
,
in
𝒬
,
-
κ
(
u
(
0
)
)
∂
u
∂
x
(
0
)
=
h
0
(
u
f
-
u
(
0
)
)
+
φ
0
(
|
u
f
|
3
u
f
-
|
u
(
0
)
|
3
u
(
0
)
)
,
+
φ
1
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
∂
x
(
L
)
=
h
1
(
u
a
-
u
(
L
)
)
+
φ
1
(
|
u
a
|
3
u
a
-
|
u
(
L
)
|
3
u
(
L
)
)
,
+
φ
1
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
in
(
0
,
T
)
,
u
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
Ω
.By subtracting (16) from (17), dividing the obtained system by ϵ, and letting ϵ tend to zero, the system satisfied by w, called the tangent linear model (TLM), is given by the system (13).We can now show the first-order necessary conditions (optimality conditions) and calculate the gradient ofJ by using TLM and by introducing an intermediate costate model. We introduce the following projection:
(18)
Π
[
a
,
b
]
(
Ψ
)
=
max
⁡
(
a
,
min
⁡
(
b
,
Ψ
)
)
,
where Ψ is an arbitrary function and a
,
b are given real constants.Theorem 8.
If J attains a (local) minimum at a pointh
*
∈
𝒴
a
d, then
(19)
h
0
*
=
Π
[
c
0
,
d
0
]
(
(
-
u
f
+
u
*
(
0
)
)
u
~
*
(
0
)
α
0
)
,
h
1
*
=
Π
[
c
1
,
d
1
]
(
(
-
u
a
+
u
*
(
L
)
)
u
~
*
(
L
)
α
1
)
,
or in the variational inequality formulation (for all h
∈
𝒴
a
d)
(20)
∫
0
T
‍
(
h
1
-
h
1
*
)
[
(
u
a
-
u
*
(
L
)
)
u
~
*
(
L
)
+
α
1
h
1
*
]
d
t
+
∫
0
T
‍
(
h
0
-
h
0
*
)
[
(
u
f
-
u
*
(
0
)
)
u
~
*
(
0
)
+
α
0
h
0
*
]
d
t
≥
0
,
where u
*
=
ℱ
(
h
*
) and u
~
*
=
ℱ
*
(
h
*
) with the function u
~
=
ℱ
*
(
h
) which is the unique solution of the adjoint (costate) problem (with initial value given at final time T) given by
(21)
-
λ
(
u
)
∂
u
~
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
~
∂
x
)
+
G
(
u
)
∂
u
~
∂
x
=
𝒞
*
𝒞
(
u
-
u
o
b
s
)
,
lllllllllllllllllllll
i
n
𝒬
,
-
κ
(
u
(
0
)
)
∂
u
~
∂
x
(
0
)
=
-
(
h
0
+
4
φ
0
|
u
(
0
)
|
3
)
u
~
+
γ
0
(
u
(
0
)
-
u
o
b
s
0
)
,
i
n
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
~
∂
x
(
L
)
=
-
(
h
1
+
4
φ
1
|
u
(
L
)
|
3
)
u
~
+
γ
1
(
u
(
L
)
-
u
o
b
s
L
)
,
i
n
(
0
,
T
)
,
λ
(
u
(
·
,
T
)
)
u
~
(
·
,
T
)
=
𝒟
*
𝒟
(
u
(
·
,
T
)
-
u
o
b
s
T
)
,
i
n
Ω
.
Moreover, the gradient ofJ at any element h of 𝒴
a
d can be given by
(22)
∂
J
∂
h
(
h
)
=
(
(
u
f
-
u
(
0
)
)
u
~
(
0
)
+
α
0
h
0
,
(
u
f
-
u
(
0
)
)
(
u
a
-
u
(
L
)
)
u
~
(
L
)
+
α
1
h
1
)
.Proof.
See the Appendix.We point out that the adjoint problem (21), which is backward in time, can be transformed into an initial-boundary value problem by the time transformation t
:
=
T
-
t, which allows to employ [16] for the existence of a unique solution of (21) for a sufficiently regular data.


### 2.3. Optimization Procedure

By solving successively the direct problem and the adjoint problem, we can therefore calculate the gradient of the objective function relative to the control parametersh
=
(
h
0
,
h
1
). Once the gradient of the objective function J, ∇
J, is known, we can seek a minimum of J. For a given observation (
u
obs
,
v
obs
), the optimization algorithm is summarized in Figure 1.Figure 1
Optimization algorithm:J is minimized until some convergence criteria are attained.


## 2.1. Problem Formulation

We assume that there exists a unique solutionu of problem (3), (6) with boundary conditions (4), (5) under some hypotheses for the data and some regularity of the operators λ and κ, satisfying the following regularity (by using [16, 17]):
(8)
u
∈
L
2
(
0
,
T
;
H
1
(
Ω
)
)
∩
L
∞
(
0
,
T
;
L
2
(
Ω
)
)
,
∂
u
∂
t
∈
L
2
(
0
,
T
;
(
H
1
(
Ω
)
)
′
)
,
|
u
(
0
)
|
3
u
(
0
)
∈
L
5
/
4
(
0
,
T
)
,
|
u
(
L
)
|
3
u
(
L
)
∈
L
5
/
4
(
0
,
T
)
.Introduce now the mappingℱ which maps the source term h of (3)–(5) into the corresponding solution u
=
ℱ
(
h
), where h
=
(
h
0
,
h
1
). In this section we formulate the inverse problem as an optimal control problem. The control procedure consists of finding the optimal controls h
opt
=
(
h
0
opt
,
h
1
opt
) and the corresponding optimal temperature u
opt
=
ℱ
(
h
opt
) which minimize a cost criterion J
(
h
). The cost functional J measures the distance between a measured temperature (the observations u
obs, u
obs
0,    u
obs
L, and u
obs
T) and the corresponding predicted temperature obtained from the primal (or direct) model (3)–(5). Precisely we will study the following optimal control problem (
𝒫
c
).Find(
u
,
h
) such that the following objective function:
(9)
J
(
h
)
=
1
2
∫
‍
∫
𝒬
‍
|
𝒞
(
u
-
u
obs
)
|
2
d
x
d
t
+
1
2
∫
Ω
‍
|
𝒟
(
u
(
T
)
-
u
obs
T
)
|
2
d
x
+
1
2
∫
0
T
‍
[
γ
0
|
u
(
0
)
-
u
obs
0
|
2
+
γ
1
|
u
(
L
)
-
u
obs
L
|
2
]
d
t
+
1
2
∫
0
T
‍
|
ℳ
h
|
2
d
tis minimized with respect toh
∈
𝒴
ad subject to (3)–(5), where 𝒴
ad is the set of admissible controls, u
=
ℱ
(
h
), γ
0 and γ
1 are predefined nonnegative weights, and ℳ
=
diag
⁡
(
α
0
,
α
1
) are predefined nonnegative weights such that α
0
+
α
1
≠
0. The operators 𝒞 and 𝒟 are unbounded operators on L
2
(
Ω
) satisfying (∀
v
∈
H
1
(
Ω
)):
(10)
∥
𝒞
v
∥
L
2
(
Ω
)
2
≤
δ
1
∥
v
∥
L
2
(
Ω
)
2
+
δ
2
∥
v
∥
H
1
(
Ω
)
2
,
∥
𝒟
v
∥
L
2
(
Ω
)
2
≤
δ
1
∥
v
∥
L
2
(
Ω
)
2
+
δ
2
∥
v
∥
H
1
(
Ω
)
2
,
with δ
1
+
δ
2
>
0 and δ
i
≥
0, for i
=
1,2.Remark 5.
The general form (9) of the objective function allows considering different kinds of problems (we denote I as the identity and d as a positive scalar): (1)
distributed observations control problems with𝒞
=
d
I, γ
0
=
γ
1
=
0, and 𝒟
=
0,
(2)
boundary observations control problems with𝒞
=
0, γ
0
=
γ
1
=
d, and 𝒟
=
0,
(3)
final observations control problems with𝒞
=
0, γ
0
=
γ
1
=
0, and 𝒟
=
d
I,
(4)
single control function, for example,h
1, with α
1
≠
0 and α
0
=
0.
We consider here an optimal control problem with Tichonov regularization terms [18]. The nonnegative weights α
0 and α
1 play the role of regularization parameters and take small values in the numerical simulations.According to (6) the set of admissible controls describing the constraint is
(11)
𝒴
ad
=
{
h
=
(
h
0
,
h
1
)
∈
(
L
2
(
0
,
T
)
)
2
such
that
(
L
2
(
0
,
T
)
)
2
0
≤
c
i
≤
h
i
≤
d
i
a
.
e
.
in
(
0
,
T
)
,
for
i
=
0,1
}
.


## 2.2. First-Order Necessary Conditions

Assume that the nonlinear control problem (9) admits an optimal solution (for similar result see, e.g., [19, 20]); the necessary conditions for this optimum is given by the following theorems (see [19]).Theorem 6.
IfJ attains a (local) minimum at a point h
*
∈
𝒴
a
d, then the following first optimality conditions hold:
(12)
J
′
(
h
*
)
·
(
h
-
h
*
)
≥
0
,
∀
h
∈
𝒴
a
d
,
where J
′ is the directional derivative of J.In order to solve numerically the optimal control problem it is necessary to derive the gradient of the cost functionalJ with respect to the control h. For this we suppose that the operator solution ℱ is continuously differentiable on 𝒴
ad and its derivative w
=
ℱ
′
(
h
)
·
g
=
lim
⁡
ϵ
→
0
(
(
ℱ
(
h
+
ϵ
g
)
-
ℱ
(
h
)
)
/
ϵ
) is the unique solution of the following system (for h
=
(
h
0
,
h
1
) and g
=
(
g
0
,
g
1
)):
(13)
λ
(
u
)
∂
w
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
w
∂
x
)
+
F
(
u
,
w
)
=
0
,
in
𝒬
,
-
κ
(
u
(
0
)
)
∂
w
∂
x
(
0
)
=
-
[
h
0
-
G
(
u
)
(
0
)
+
4
φ
0
|
u
(
0
)
|
3
]
w
(
0
)
+
g
0
(
u
f
-
u
(
0
)
)
,
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
w
∂
x
(
L
)
=
-
[
h
1
+
G
(
u
)
(
L
)
+
4
φ
1
|
u
(
L
)
|
3
]
w
(
L
)
+
g
1
(
u
a
-
u
(
L
)
)
,
in
(
0
,
T
)
,
w
(
x
,
0
)
=
0
,
x
∈
Ω
,where operatorsF and G are defined as follows:(14)
F
(
u
,
w
)
=
∂
∂
t
(
λ
(
u
)
)
w
-
∂
∂
x
(
w
∂
∂
x
(
κ
(
u
)
)
)
,
G
(
u
)
=
∂
∂
x
(
κ
(
u
)
)
.Remark 7.
We verify easily that(15)
λ
(
u
)
∂
w
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
w
∂
x
)
+
F
(
u
,
w
)
=
∂
(
λ
(
u
)
w
)
∂
t
-
∂
2
∂
x
2
(
κ
(
u
)
w
)
,
κ
(
u
(
0
)
)
∂
w
∂
x
(
0
)
+
G
(
u
)
(
0
)
w
(
0
)
=
∂
(
κ
(
u
)
w
)
∂
x
(
0
)
,
κ
(
u
(
L
)
)
∂
w
∂
x
(
L
)
+
G
(
u
)
(
L
)
w
(
L
)
=
∂
(
κ
(
u
)
w
)
∂
x
(
L
)
.In order to derive (13), we write the two systems satisfied by u
ϵ
=
ℱ
(
h
+
ϵ
g
) and u
=
ℱ
(
h
) as
(16)
λ
(
u
ϵ
)
∂
u
ϵ
∂
t
-
∂
∂
x
(
κ
(
u
ϵ
)
∂
u
ϵ
∂
x
)
=
0
,
in
𝒬
,
-
κ
(
u
ϵ
(
0
)
)
∂
u
ϵ
∂
x
(
0
)
=
(
h
0
+
ϵ
g
0
)
(
u
f
-
u
ϵ
(
0
)
)
+
φ
0
(
|
u
f
|
3
u
f
-
|
u
ϵ
(
0
)
|
3
u
ϵ
(
0
)
)
,
+
φ
0
|
u
f
|
3
u
f
-
|
u
ϵ
(
0
)
|
3
in
(
0
,
T
)
,
κ
(
u
ϵ
(
L
)
)
∂
u
ϵ
∂
x
(
L
)
=
(
h
1
+
ϵ
g
1
)
(
u
a
-
u
ϵ
(
L
)
)
+
φ
1
(
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
u
ϵ
(
L
)
)
,
+
φ
1
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
in
(
0
,
T
)
,
u
ϵ
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
Ω
,
(17)
λ
(
u
)
∂
u
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
=
0
,
in
𝒬
,
-
κ
(
u
(
0
)
)
∂
u
∂
x
(
0
)
=
h
0
(
u
f
-
u
(
0
)
)
+
φ
0
(
|
u
f
|
3
u
f
-
|
u
(
0
)
|
3
u
(
0
)
)
,
+
φ
1
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
∂
x
(
L
)
=
h
1
(
u
a
-
u
(
L
)
)
+
φ
1
(
|
u
a
|
3
u
a
-
|
u
(
L
)
|
3
u
(
L
)
)
,
+
φ
1
|
u
a
|
3
u
a
-
|
u
ϵ
(
L
)
|
3
in
(
0
,
T
)
,
u
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
Ω
.By subtracting (16) from (17), dividing the obtained system by ϵ, and letting ϵ tend to zero, the system satisfied by w, called the tangent linear model (TLM), is given by the system (13).We can now show the first-order necessary conditions (optimality conditions) and calculate the gradient ofJ by using TLM and by introducing an intermediate costate model. We introduce the following projection:
(18)
Π
[
a
,
b
]
(
Ψ
)
=
max
⁡
(
a
,
min
⁡
(
b
,
Ψ
)
)
,
where Ψ is an arbitrary function and a
,
b are given real constants.Theorem 8.
If J attains a (local) minimum at a pointh
*
∈
𝒴
a
d, then
(19)
h
0
*
=
Π
[
c
0
,
d
0
]
(
(
-
u
f
+
u
*
(
0
)
)
u
~
*
(
0
)
α
0
)
,
h
1
*
=
Π
[
c
1
,
d
1
]
(
(
-
u
a
+
u
*
(
L
)
)
u
~
*
(
L
)
α
1
)
,
or in the variational inequality formulation (for all h
∈
𝒴
a
d)
(20)
∫
0
T
‍
(
h
1
-
h
1
*
)
[
(
u
a
-
u
*
(
L
)
)
u
~
*
(
L
)
+
α
1
h
1
*
]
d
t
+
∫
0
T
‍
(
h
0
-
h
0
*
)
[
(
u
f
-
u
*
(
0
)
)
u
~
*
(
0
)
+
α
0
h
0
*
]
d
t
≥
0
,
where u
*
=
ℱ
(
h
*
) and u
~
*
=
ℱ
*
(
h
*
) with the function u
~
=
ℱ
*
(
h
) which is the unique solution of the adjoint (costate) problem (with initial value given at final time T) given by
(21)
-
λ
(
u
)
∂
u
~
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
~
∂
x
)
+
G
(
u
)
∂
u
~
∂
x
=
𝒞
*
𝒞
(
u
-
u
o
b
s
)
,
lllllllllllllllllllll
i
n
𝒬
,
-
κ
(
u
(
0
)
)
∂
u
~
∂
x
(
0
)
=
-
(
h
0
+
4
φ
0
|
u
(
0
)
|
3
)
u
~
+
γ
0
(
u
(
0
)
-
u
o
b
s
0
)
,
i
n
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
~
∂
x
(
L
)
=
-
(
h
1
+
4
φ
1
|
u
(
L
)
|
3
)
u
~
+
γ
1
(
u
(
L
)
-
u
o
b
s
L
)
,
i
n
(
0
,
T
)
,
λ
(
u
(
·
,
T
)
)
u
~
(
·
,
T
)
=
𝒟
*
𝒟
(
u
(
·
,
T
)
-
u
o
b
s
T
)
,
i
n
Ω
.
Moreover, the gradient ofJ at any element h of 𝒴
a
d can be given by
(22)
∂
J
∂
h
(
h
)
=
(
(
u
f
-
u
(
0
)
)
u
~
(
0
)
+
α
0
h
0
,
(
u
f
-
u
(
0
)
)
(
u
a
-
u
(
L
)
)
u
~
(
L
)
+
α
1
h
1
)
.Proof.
See the Appendix.We point out that the adjoint problem (21), which is backward in time, can be transformed into an initial-boundary value problem by the time transformation t
:
=
T
-
t, which allows to employ [16] for the existence of a unique solution of (21) for a sufficiently regular data.


## 2.3. Optimization Procedure

By solving successively the direct problem and the adjoint problem, we can therefore calculate the gradient of the objective function relative to the control parametersh
=
(
h
0
,
h
1
). Once the gradient of the objective function J, ∇
J, is known, we can seek a minimum of J. For a given observation (
u
obs
,
v
obs
), the optimization algorithm is summarized in Figure 1.Figure 1
Optimization algorithm:J is minimized until some convergence criteria are attained.


## 3. Numerical Analysis and Simulations

In this section, we seek to estimate separately the parametersh
0 and h
1. So, in order to facilitate the presentation, we denote by η the control function which plays the role of h
i, for i
=
0 or 1. Moreover, we consider two kinds of control problems:(i)
distributed observations control problems, where𝒞 is the identity function, γ
0
=
γ
1
=
0, 𝒟
=
0, and
(23)
J
(
η
)
=
1
2
∥
u
-
u
obs
∥
L
2
(
𝒬
)
2
,
(ii)
boundary observations control problems with𝒞
=
0, γ
0
=
γ
1
=
1, 𝒟
=
0, and
(24)
J
(
η
)
=
1
2
∥
u
(
0
)
-
u
obs
0
∥
L
2
(
0
,
T
)
2
+
1
2
∥
u
(
L
)
-
u
obs
L
∥
L
2
(
0
,
T
)
2
,where u is the solution of the direct problem corresponding to the control function η which will be denoted in the sequel by ℱ
(
η
). Then, the expression of the corresponding gradient can be given (for both kinds of cost functions) by
(25)
g
=
∇
J
(
η
)
=
[
u
f
-
u
(
0
,
·
)
]
u
~
(
0
,
·
)
if
η
=
h
0
,
g
=
∇
J
(
η
)
=
[
u
a
-
u
(
L
,
·
)
]
u
~
(
L
,
·
)
if
η
=
h
1
,
where u
~ is the solution of adjoint problem (21), corresponding to the direct model.


### 3.1. Numerical Implementation and Outline

As noted in the previous subsection, solving the nonlinear boundary control problem (9) by a gradient method requires, at each iteration of the optimization algorithm, solving the direct problem and its corresponding adjoint problem. In order to solve numerically these two problems, we use first the discretization in space with Lagrange finite elements of order 2 and the derived differential systems are integrated by using variable-order, variable-step-size backward differentiation formulas [21]. Finally, the obtained nonlinear algebraic systems are solved with Newton's method [22] and at each iteration a direct method is used to solve the considered linear system. In order to solve the optimization problems, we have used the gradient method (GM) and the conjugate gradient method (CGM). For the CGM method we have considered the following well-known descent direction methods: Fletcher-Reeves, Polak-Ribière, Hestenes-Stiefel, and the recent method of Dai and Yuan [23]. More precisely, for the iteration index k, we denote by g
(
k
) the numerical approximation of the gradient function g (given by (25)), d
(
k
) the descent direction, ρ
k the descent step, and η
k the numerical approximation of the control function η; at the kth iteration of the algorithm, the considered gradient schemes are as follows.(1)
Initialization:η
0 (given):
(a)
computeu
0
=
ℱ
(
η
0
) by solving the direct problem;
(b)
computeu
~
0 by solving the adjoint problem;
(c)
gradient ofJ at η
0, g
(
0
) is given by (25);
(d)
determine the direction:d
(
0
)
=
-
g
(
0
);
(e)
determineη
1
=
η
0
+
ρ
0
d
(
0
) and initialize k
=
1.
(2)
Computeu
k
=
ℱ
(
η
k
) by solving the direct problem.
(3)
Computeu
~
k by solving the adjoint problem.
(4)
Gradient ofJ at η
k, g
(
k
) is given by (25).
(5)
Determine the directiond
(
k
) by one of the expressions shown in Table 1 (where (
·
;
·
) and ∥
·
∥ are the scalar product and its associated norm).
(6)
Determine the descent stepρ
k by the methods shown in Table 2.
(7)
Determineη
k
+
1
=
η
k
+
ρ
k
d
(
k
).
(8)
If the gradient is sufficiently small then end; else setk
:
=
k
+
1 and go to (
2
). The approximation of the optimal solution η
* is η
k.Table 1
Method
Notation
Descent direction
Parameter
Gradient
Grad
d
(
k
)
=
-
g
(
k
)
—
Fletcher-Reeves CGM
CGFR
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
∥
g
(
k
)
∥
2
∥
g
(
k
-
1
)
∥
2
Polak-Ribieres CGM
CGPR
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
(
g
(
k
)
-
g
(
k
-
1
)
;
g
(
k
)
)
∥
g
(
k
-
1
)
∥
2
Hestenes-Stiefel CGM
CGHS
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
(
g
(
k
)
;
g
(
k
)
-
g
(
k
-
1
)
)
(
d
(
k
-
1
)
;
g
(
k
)
-
g
(
k
-
1
)
)
Dai-Yuan CGM
CGDY
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
∥
g
(
k
)
∥
2
(
d
(
k
-
1
)
;
g
(
k
)
-
g
(
k
-
1
)
)Table 2
Method
Notation
Description
Constant step
cst
ρ
ρ
k
=
ρ
Decreasing step
dim
⁡
ρ
ρ
0
=
ρ,  ρ
k
=
ρ
k
-
1
2 if J
(
η
k
)
>
J
(
η
k
-
1
) (not changed elsewhere)
Minimum step
min
⁡
ρ
k
=
min
⁡
⁡
(
1
,
1
∥
g
(
k
)
∥
∞
)
Optimal step
opt
ρ
k is an approximate solution of min
ρ
≥
0
⁡
J
(
η
k
+
ρ
d
(
k
)
)Remark 9.
An approximate value of the optimal step is obtained using the linear approximation ofℱ. For example, in the case of the distributed observations control problem, we have
(26)
ρ
k
≃
-
(
u
k
-
u
obs
,
w
k
)
L
2
(
𝒬
)
∥
w
k
∥
L
2
(
𝒬
)
,
where w
k is solution of (13) for the iteration k.Remark 10.
If the exact control functionη
ex is known, we can measure the efficiency of the method with the following relative error on η
k:
(27)
errel
(
η
k
)
=
∥
η
ex
-
η
k
∥
L
2
(
𝒬
)
∥
η
ex
∥
L
2
(
𝒬
)
.Then, to solve the optimal control problem (9), we have developed a specific software, based on Comsol and Matlab tools, taking into account the following: first the nature of the nonlinearity in the operators λ and κ and the nonlinear radiative term on the boundary of the domain and second the nature of the adjoint problem which is backward in time and coupled with the direct problem. In order to validate our approach, we have studied several examples in different situations and used different descent direction methods for the optimization algorithm. In this paper we present only two application examples and a realistic simulation model. For the first application, we consider the model without the radiative boundary conditions and we assume that the observation u
obs is an analytical given function. In the second one, the observation is computed by solving the model with the radiative boundary conditions, corresponding to given convective heat transfer parameters. For the realistic situation, we consider the real parameters and data given in [4, 16] and we construct the observation model as in the second application. With these examples, we validate the method and discuss two numerical questions: the improvement of the computed solution at final time and the stability of the algorithm for noisy data.The numerical simulations are performed on a computer with a processor 2.8 GHz Core 2 Duo and 4 Gb memory and taken between 1400 and 2200 s for 50 optimization iterations (depending on the descent direction method and on the treated example).


### 3.2. Numerical Examples and Validation

In this section, we denote byη
ex the exact value of the control function η that we want to approximate by η
k, at the kth iteration of the gradient algorithm. Moreover, we assume that T
=
L
=
1 and we fix the time step δ
t and space step δ
x to 1
0
-
3.To simplify the presentation, we detail only the results for the control on the convective heat transfer coefficienth
1 and give some comments for the control on the convective heat transfer coefficient h
0.


#### 3.2.1. Example without Radiation Conditions

In this example, the control functionη plays the role of h
1; the observation is the given functions u
obs
=
e
t
+
L
-
x and η
ex
(
=
h
1
)
=
1
+
cos
⁡
(
t
). The other operators and data of the model are given by h
0
=
1
,
φ
0
=
φ
1
=
0
,
λ
(
u
)
=
1
+
u, κ
(
u
)
=
1
+
u
2, u
0
(
x
)
=
e
L
-
x
,
u
f
(
t
)
=
[
2
+
e
2
(
t
+
L
)
]
e
t
+
L
,
u
a
(
t
)
=
[
1
-
(
1
+
e
2
t
)
/
(
1
+
cos
⁡
t
)
]
e
t and we have added the following function F
(
x
,
t
)
=
e
2
(
t
+
L
-
x
)
(
1
-
3
e
t
+
L
-
x
) such that (3) becomes
(28)
λ
(
u
)
∂
u
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
=
F
to ensure that u
obs is the analytical solution of the direct problem (corresponding to η
ex).Remark 11.
We want to emphasize on the fact that the addition of the artificial right-hand memberF does not change the formulation of the adjoint problem, the expression of the gradient of the objective function, and the optimization algorithm.To validate our approach, we have tested all the gradient methods described above. With all these gradient methods we get computed valuesη
k which converge towards the exact value η
ex: the relative error (27) decreases when k increases (see Figures 2, 3, and 4). But the convergence speed depends on the method and on the initial value. So, we are going to discuss these two points for the distributed control problem (23) and then compare with the results of the boundary control problems (24).Figure 2
Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k for different initializations with the CGDYOPT method.Example of Section3.2.1—η
k
(
t
) for some k with different initializations: (a) η
0
=
η
ex
(
0
)
=
2 and (b) η
0
=
η
ex
(
T
)
=
1
+
cos
⁡
(
T
).
(a)
(b)Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k with (a) CGPR method and different optimization steps and (b) an optimal step and different gradient methods.
(a)
(b)First, to use the optimization algorithm we have chosen an initial valueη
0 for the control function η. We present on Figure 2 the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for different initializations and on Figure 3 the function η
k
(
t
) for some k (computed with the CGDYOPT method: Dai-Yuan conjugate gradient with optimal step) and for these different initializations: (a) η
0
=
2.2, (b) η
0
=
η
ex
(
0
), (c) η
0
=
1.7, (d) η
0
=
η
ex
(
T
), (e) η
0
=
1, and (f) η
0
=
[
η
ex
(
T
)
-
η
ex
(
0
)
]
(
t
/
T
)
+
η
ex
(
0
). As expected, the method is convergent for each initialization but the convergence speed depends on the proximity of the initial value η
0 from the exact value η
ex but the best accuracy that can be obtained with the method depends on the gap between η
0
(
T
) and η
ex
(
T
).Indeed, for each numerical simulation, we note thatη
k
(
T
)
=
η
0 for each k (see Figure 3 for the function η
k
(
t
) with (a) η
0
=
η
ex
(
0
)
=
2 and (b) η
0
=
η
ex
(
T
)
=
1
+
cos
⁡
(
T
)): the last value of η
k is never modified. For this kind of optimization problem, there is no final observation and the gradient (like the descent direction) is always zero for t
=
T. Then, the parameter cannot be well approximated near the final time T. Two methods will be proposed in the next subsection to solve this problem. Of course, this has no effect when the initialization is a continuous function with η
0
(
T
)
=
η
ex
(
T
) (see Figure 3(b)). From a physical point of view, the good choice for the initial value can be η
0 such that η
0
(
0
)
=
η
ex
(
0
) because the parameters are assumed known at ambient temperature (corresponding to t
=
0). It is why we use this initialization for the next examples.To compare the different gradient methods, we fix the initialization toη
0
=
η
ex
(
T
) and we present on Figure 4(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for the CGPR method with different optimization steps and on Figure 4(b) the convergence curves for the different gradient methods with optimal step. We conclude that the choice of the optimal step gives a significant improvement compared to other possibilities. As expected the conjugate gradient methods have a better convergence speed than the gradient method. Among the conjugate gradient methods, the Dai-Yuan method seems to be the better choice because it gives the better result for almost all iterations. This is in agreement with the literature results in optimization. It is why we present mainly the results for this method for the next examples.Now, we want to compare the distributed observations control problem and the boundary observations control problem. We present on Figure5 the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for the CGDY method and η
0
=
η
ex
(
0
). In this case, the convergence is a little bit better for the boundary control problem but it depends on the initialization. The important point is that the results are close and it is sufficient to measure the values of u only in x
=
0 and x
=
L.Figure 5
Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k for the distributed and boundary control problems.Remark 12.
To explore the difference between the control on the cold side of the domain (x
=
L) and the control on the hot side of the domain (x
=
0), we have tested the control on the convective heat transfer coefficient h
0 (the function η plays here the role of h
0). We have considered the same data as in the previous example except for the values of h
0 and h
1 that are exchanged η
ex
(
=
h
0
)
=
1
+
cos
⁡
(
t
)
,
h
1
=
1, and the related data u
f
(
t
)
=
[
1
+
(
1
+
e
2
(
t
+
L
)
)
/
(
1
+
cos
⁡
t
)
]
e
t
+
L
,
u
a
(
t
)
=
-
e
3
t, and we obtained similar results to the previous ones.


#### 3.2.2. How to Improve the Result at Final Time?

We propose here two methods to overcome the difficulty to compute the solution of the inverse problem at the final time, the gradient of the functional being null for each step of the algorithm. The first possibility is to add a term with final observation to the functional and the second to deduce the value of the solution at final time by a linear approximation of previous values as presented in [24]. We present these methods for the distributed control problem.The first method consists of replacing the expression ofJ in (23) by
(29)
J
(
η
)
=
1
2
∥
u
-
u
obs
∥
L
2
(
𝒬
)
2
+
β
2
∥
u
(
T
)
-
u
ob
s
T
∥
L
2
(
Ω
)
2
,
where β is a positive constant. More precisely, 𝒟 is the identity function multiplied by β in (9). The expression of the gradient (25) is unchanged, but the definition of the adjoint problem (21) is modified because of the new value of 𝒟.The choice of the value ofβ is important. If it is too large the algorithm does not converge to the expected solution for the same initialization because, with only a final observation (problem of exact controllability [25, 26]), different values of η give solutions u of the direct problem with the same final value. If it is too small, there is no change with the new functional. As we can see in Figure 6, a good choice consists to give the final observation the same importance than each of the other observations with β being a value close to the time step δ
t
=
1
0
-
3. We use here the data of the example of Section 3.2.1 defined in the previous section.Example of Section3.2.1 with final observation: (a) log
⁡
10
(
errel
(
η
k
)
) versus k for different β and (b) η
k
(
t
) for some k with β
=
0.001.
(a)
(b)The second method consists of a linear approximation of the solution at final time by previous values of the solution. In [24], is proposed the following formula:
(30)
η
k
(
T
)
=
2
η
k
(
T
-
δ
t
)
-
η
k
(
T
-
2
δ
t
)
which is efficient if the time step is not too small. For δ
t
=
1
0
-
3, the value of η
k
(
T
) converges very slowly to the exact value; see Figure 7. This convergence can be accelerated with the following generalization of the formula:
(31)
η
k
(
T
)
=
2
η
k
(
T
-
n
l
c
δ
t
)
-
η
k
(
T
-
2
n
l
c
δ
t
)
with a small integer n
l
c which should not be too large to avoid instability. We present on Figure 8(a) the convergence curves for different values of n
l
c and on Figure 8(b)  η
k computed with n
l
c
=
2. We will use this value in the sequel.Figure 7
Example of Section3.2.1—linear approximation—η
k
(
t
) for some k with n
l
c
=
1.Example of Section3.2.1—linear approximation (a) log
⁡
10
(
errel
(
η
k
)
) versus k for different n
l
c and (b) η
k
(
t
) for some k with n
l
c
=
2.
(a)
(b)


#### 3.2.3. Example with Radiation Conditions

For this second kind of tests, there is no more analytical expression for the observationu
obs which is now computed as the solution of the direct problem (3)-(4)-(5) with the exact value η
ex of the control function η. Then, solving the control problem from an initialization η
0 of η, we compute approximations η
k of η
ex using the CGDYOPT method.As in the previous section we detail the results for the control on the convective heat transfer coefficienth
1 and give a remark for the control on the convective heat transfer coefficient h
0.With the datah
0
=
1
,
h
1
=
1
+
cos
⁡
(
t
)
,
φ
0
=
φ
1
=
1
0
-
4, λ
(
u
)
=
1
+
u
,
κ
(
u
)
=
1
+
u
2
,
u
a
=
10
,
u
f
(
t
)
=
20
+
ln
⁡
(
t
)
,
u
0
=
u
a, we compute the observation u
obs as numerical solution of (3)-(4)-(5). It is a problem with radiation boundary conditions.The observation being known on the boundary and because of the time dependence ofh
0 and h
1, we can consider an optimization problem without radiation condition on the boundary. Indeed, the boundary conditions of the direct problem become
(32)
-
κ
(
u
(
0
)
)
∂
u
∂
x
(
0
)
=
H
0
(
u
f
-
u
(
0
)
)
,
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
∂
x
(
L
)
=
H
1
(
u
a
-
u
(
L
)
)
,
in
(
0
,
T
)
,
with
(33)
H
0
=
h
0
+
φ
0
(
u
f
+
u
obs
)
(
u
f
2
+
u
obs
2
)
,
H
1
=
h
1
+
φ
1
(
u
a
+
u
obs
)
(
u
a
2
+
u
obs
2
)
.The functionsh
0 and H
0, h
1 and H
1 are compared in Figure 9.Example of Section3.2.3—comparing data (a) h
0 and H
0 and (b) h
1 and H
1.
(a)
(b)In this example, the control functionη plays the role of H
1 with η
ex
=
H
1 and we solve the same optimization problem as in Example of Section 3.2.1 with F
=
0 and (
h
0
,
h
1
) replaced by (
H
0
,
H
1
). The functions u
obs and u
a being known, we can deduce φ
1 from H
1 and h
1 or h
1 from H
1 and φ
1.For the initializationη
0
=
η
ex
(
0
), we present on Figure 10(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for both kinds of problems (with distributed observations and boundary observations) and (b) the function η
k
(
t
) for some k (corresponding to the boundary observation). The results are good but contrary to what is expected, we note again a better result with the boundary observations (it depends on the initialization, e.g., with η
0
=
η
ex
(
T
) the method with distributed observations is a little bit better).Example of Section3.2.3 with n
l
c
=
2: (a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
k
(
t
) for some k with boundary observation.
(a)
(b)Remark 13.
We observed that the result is worse if we initialize our algorithm withη
0
=
η
ex
(
T
) (instead of η
0
=
η
ex
(
0
)). It can be explained by the fact that for η
0
=
η
ex
(
T
), η
k is not changed for t small because without a source term in (3), the value of H
1, used in the boundary condition at the point x
=
L, has no significant effect immediately on u. This comes from the fact that the source term u
f is only on the boundary at the point x
=
0, and there is a delay depending on the parameters of the equation. We checked that there is no more delay when the source term and the control are at the same point x
=
0.


#### 3.2.4. Numerical Stability of the Algorithm

In real problems, the observationu
obs comes from experimental data with measurement errors and it is interesting to study the influence of these errors on the result given by the algorithm. As in the stability study presented in [24] for a linear problem, the noisy data are generated at each point number i of the mesh by the formula
(34)
u
~
obs
,
i
=
u
obs
,
i
(
1
+
δ
ζ
i
)
,
where δ gives the relative noise level and ζ
i is a random number obtained by a normal distribution of mean zero and unit standard deviation (Matlab function randn).On Figure11, we present the results for the example of Section 3.2.1 with boundary observation and final observation (β
=
10
-
2) (the results are similar with n
l
c
=
2 instead of final observation) and initialization η
0
(
t
)
=
η
ex
(
0
), δ
t
=
1
0
-
3
,
δ
x
=
1
0
-
2. For the 3 noise levels (1
%, 3
%, and 5
%), the computed solution remains close to the exact solution with a gap comparable to the noise level.Example of Section3.2.1 with β
=
0.01 and boundary observation—(a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
20
(
t
), for different values of noise level δ.
(a)
(b)The results are completely different for the example of Section3.2.3 (see Figure 12). For noise level δ
≥
0.0005, the computed solution does not remain close to the exact solution. The radiation condition on the boundary, with term (
|
u
a
|
3
u
a
-
|
u
|
3
u
), is very sensitive to the perturbations of the data.Figure 12
Example of Section3.2.3 with β
=
0.01 and boundary observation—η
7
(
t
) for different values of noise level δ.


### 3.3. Simulation on Real Data

We now apply our method to a realistic problem. To that purpose, we take similarly data as in [4, 16]. We consider a gypsum plasterboard, with thickness L
=
4cm, exposed to fire in x
=
0, during T
=
1200 s. The surroundings temperature is u
a
=
2
0
∘C and the furnace temperature is u
f
=
u
a
+
345
log
⁡
10
(
1
+
8
t
/
60
) (see Figure 13). The initial temperature is u
0
=
u
a. The thermal conductivity is given by
(35)
κ
(
u
)
=
{
0.25
if
u
<
112
.
5
°
C
,
0.1116
+
6.578
·
1
0
-
5
u
if
112
.
5
∘
C
≤
u
<
393.7
5
°
C
,
0.0277
+
2.788
·
1
0
-
4
u
if
393.7
5
∘
C
≤
u
<
797
.
2
°
C
,
-
0.91323
+
1.459
·
1
0
-
3
u
if
797
.
2
°
C
≤
u
.Figure 13
Furnace temperature versus time.We haveλ
(
u
)
=
ρ
(
u
)
C
p
(
u
), where the density is defined by
(36)
ρ
(
u
)
=
{
648
if
u
<
12
0
°
C
,
534.6
if
12
0
°
C
≤
u
≤
65
0
∘
C
,
500.6
if
65
0
°
C
<
u
and the specific heat is defined with a peak between 100 and 17
0
°
C and a value corresponding to the total amount of energy of the two dehydration reactions, 669 kJ:
(37)
C
p
(
u
)
=
{
6.69
·
1
0
4
7
if
10
0
∘
C
≤
u
≤
17
0
°
C
,
950
else
.We supposed that the convective heat transfer coefficienth
0
=
25 and the effective emissivity ϵ
0
=
0.82 are known in x
=
L. We deduce that φ
0
=
σ
ϵ
0 with σ
=
5.6697
·
1
0
-
8
W
/
m
2
K
4, the Stefan-Boltzmann's constant.We want to estimateh
1, supposing that ϵ
1
=
0. The control function η plays here the role of h
1. We choose the following exact expression for h
1: η
ex
(
=
h
1
)
=
5
+
cos
⁡
(
t
/
T
). Then, we can compute the observation u
obs as the solution of the direct problem (3), solve the control problem (9) (with the cost function and its gradient given by (23) or (24) and (25)) by the optimization algorithm (using the CGDYOPT method), and compute error estimations on the control function η at each optimization iteration.For the space stepδ
x
=
0.1cm, the time step δ
t
=
1s, and the initialization η
0
=
η
ex
(
0
), we present in Figure 14(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for both kinds of problems (with distributed observation and boundary observation) and (b) the function η
k
(
t
) for some k (corresponding to the boundary observation). Despite the fact that the nonlinear parameters λ
(
u
) and κ
(
u
) are not regular, we get a good convergence of the method.Simulation on real data withn
l
c
=
2: (a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
k
(
t
) for some k with boundary observation.
(a)
(b)


## 3.1. Numerical Implementation and Outline

As noted in the previous subsection, solving the nonlinear boundary control problem (9) by a gradient method requires, at each iteration of the optimization algorithm, solving the direct problem and its corresponding adjoint problem. In order to solve numerically these two problems, we use first the discretization in space with Lagrange finite elements of order 2 and the derived differential systems are integrated by using variable-order, variable-step-size backward differentiation formulas [21]. Finally, the obtained nonlinear algebraic systems are solved with Newton's method [22] and at each iteration a direct method is used to solve the considered linear system. In order to solve the optimization problems, we have used the gradient method (GM) and the conjugate gradient method (CGM). For the CGM method we have considered the following well-known descent direction methods: Fletcher-Reeves, Polak-Ribière, Hestenes-Stiefel, and the recent method of Dai and Yuan [23]. More precisely, for the iteration index k, we denote by g
(
k
) the numerical approximation of the gradient function g (given by (25)), d
(
k
) the descent direction, ρ
k the descent step, and η
k the numerical approximation of the control function η; at the kth iteration of the algorithm, the considered gradient schemes are as follows.(1)
Initialization:η
0 (given):
(a)
computeu
0
=
ℱ
(
η
0
) by solving the direct problem;
(b)
computeu
~
0 by solving the adjoint problem;
(c)
gradient ofJ at η
0, g
(
0
) is given by (25);
(d)
determine the direction:d
(
0
)
=
-
g
(
0
);
(e)
determineη
1
=
η
0
+
ρ
0
d
(
0
) and initialize k
=
1.
(2)
Computeu
k
=
ℱ
(
η
k
) by solving the direct problem.
(3)
Computeu
~
k by solving the adjoint problem.
(4)
Gradient ofJ at η
k, g
(
k
) is given by (25).
(5)
Determine the directiond
(
k
) by one of the expressions shown in Table 1 (where (
·
;
·
) and ∥
·
∥ are the scalar product and its associated norm).
(6)
Determine the descent stepρ
k by the methods shown in Table 2.
(7)
Determineη
k
+
1
=
η
k
+
ρ
k
d
(
k
).
(8)
If the gradient is sufficiently small then end; else setk
:
=
k
+
1 and go to (
2
). The approximation of the optimal solution η
* is η
k.Table 1
Method
Notation
Descent direction
Parameter
Gradient
Grad
d
(
k
)
=
-
g
(
k
)
—
Fletcher-Reeves CGM
CGFR
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
∥
g
(
k
)
∥
2
∥
g
(
k
-
1
)
∥
2
Polak-Ribieres CGM
CGPR
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
(
g
(
k
)
-
g
(
k
-
1
)
;
g
(
k
)
)
∥
g
(
k
-
1
)
∥
2
Hestenes-Stiefel CGM
CGHS
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
(
g
(
k
)
;
g
(
k
)
-
g
(
k
-
1
)
)
(
d
(
k
-
1
)
;
g
(
k
)
-
g
(
k
-
1
)
)
Dai-Yuan CGM
CGDY
d
(
k
)
=
-
g
(
k
)
+
β
k
-
1
d
(
k
-
1
)
β
k
-
1
=
∥
g
(
k
)
∥
2
(
d
(
k
-
1
)
;
g
(
k
)
-
g
(
k
-
1
)
)Table 2
Method
Notation
Description
Constant step
cst
ρ
ρ
k
=
ρ
Decreasing step
dim
⁡
ρ
ρ
0
=
ρ,  ρ
k
=
ρ
k
-
1
2 if J
(
η
k
)
>
J
(
η
k
-
1
) (not changed elsewhere)
Minimum step
min
⁡
ρ
k
=
min
⁡
⁡
(
1
,
1
∥
g
(
k
)
∥
∞
)
Optimal step
opt
ρ
k is an approximate solution of min
ρ
≥
0
⁡
J
(
η
k
+
ρ
d
(
k
)
)Remark 9.
An approximate value of the optimal step is obtained using the linear approximation ofℱ. For example, in the case of the distributed observations control problem, we have
(26)
ρ
k
≃
-
(
u
k
-
u
obs
,
w
k
)
L
2
(
𝒬
)
∥
w
k
∥
L
2
(
𝒬
)
,
where w
k is solution of (13) for the iteration k.Remark 10.
If the exact control functionη
ex is known, we can measure the efficiency of the method with the following relative error on η
k:
(27)
errel
(
η
k
)
=
∥
η
ex
-
η
k
∥
L
2
(
𝒬
)
∥
η
ex
∥
L
2
(
𝒬
)
.Then, to solve the optimal control problem (9), we have developed a specific software, based on Comsol and Matlab tools, taking into account the following: first the nature of the nonlinearity in the operators λ and κ and the nonlinear radiative term on the boundary of the domain and second the nature of the adjoint problem which is backward in time and coupled with the direct problem. In order to validate our approach, we have studied several examples in different situations and used different descent direction methods for the optimization algorithm. In this paper we present only two application examples and a realistic simulation model. For the first application, we consider the model without the radiative boundary conditions and we assume that the observation u
obs is an analytical given function. In the second one, the observation is computed by solving the model with the radiative boundary conditions, corresponding to given convective heat transfer parameters. For the realistic situation, we consider the real parameters and data given in [4, 16] and we construct the observation model as in the second application. With these examples, we validate the method and discuss two numerical questions: the improvement of the computed solution at final time and the stability of the algorithm for noisy data.The numerical simulations are performed on a computer with a processor 2.8 GHz Core 2 Duo and 4 Gb memory and taken between 1400 and 2200 s for 50 optimization iterations (depending on the descent direction method and on the treated example).


## 3.2. Numerical Examples and Validation

In this section, we denote byη
ex the exact value of the control function η that we want to approximate by η
k, at the kth iteration of the gradient algorithm. Moreover, we assume that T
=
L
=
1 and we fix the time step δ
t and space step δ
x to 1
0
-
3.To simplify the presentation, we detail only the results for the control on the convective heat transfer coefficienth
1 and give some comments for the control on the convective heat transfer coefficient h
0.


### 3.2.1. Example without Radiation Conditions

In this example, the control functionη plays the role of h
1; the observation is the given functions u
obs
=
e
t
+
L
-
x and η
ex
(
=
h
1
)
=
1
+
cos
⁡
(
t
). The other operators and data of the model are given by h
0
=
1
,
φ
0
=
φ
1
=
0
,
λ
(
u
)
=
1
+
u, κ
(
u
)
=
1
+
u
2, u
0
(
x
)
=
e
L
-
x
,
u
f
(
t
)
=
[
2
+
e
2
(
t
+
L
)
]
e
t
+
L
,
u
a
(
t
)
=
[
1
-
(
1
+
e
2
t
)
/
(
1
+
cos
⁡
t
)
]
e
t and we have added the following function F
(
x
,
t
)
=
e
2
(
t
+
L
-
x
)
(
1
-
3
e
t
+
L
-
x
) such that (3) becomes
(28)
λ
(
u
)
∂
u
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
=
F
to ensure that u
obs is the analytical solution of the direct problem (corresponding to η
ex).Remark 11.
We want to emphasize on the fact that the addition of the artificial right-hand memberF does not change the formulation of the adjoint problem, the expression of the gradient of the objective function, and the optimization algorithm.To validate our approach, we have tested all the gradient methods described above. With all these gradient methods we get computed valuesη
k which converge towards the exact value η
ex: the relative error (27) decreases when k increases (see Figures 2, 3, and 4). But the convergence speed depends on the method and on the initial value. So, we are going to discuss these two points for the distributed control problem (23) and then compare with the results of the boundary control problems (24).Figure 2
Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k for different initializations with the CGDYOPT method.Example of Section3.2.1—η
k
(
t
) for some k with different initializations: (a) η
0
=
η
ex
(
0
)
=
2 and (b) η
0
=
η
ex
(
T
)
=
1
+
cos
⁡
(
T
).
(a)
(b)Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k with (a) CGPR method and different optimization steps and (b) an optimal step and different gradient methods.
(a)
(b)First, to use the optimization algorithm we have chosen an initial valueη
0 for the control function η. We present on Figure 2 the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for different initializations and on Figure 3 the function η
k
(
t
) for some k (computed with the CGDYOPT method: Dai-Yuan conjugate gradient with optimal step) and for these different initializations: (a) η
0
=
2.2, (b) η
0
=
η
ex
(
0
), (c) η
0
=
1.7, (d) η
0
=
η
ex
(
T
), (e) η
0
=
1, and (f) η
0
=
[
η
ex
(
T
)
-
η
ex
(
0
)
]
(
t
/
T
)
+
η
ex
(
0
). As expected, the method is convergent for each initialization but the convergence speed depends on the proximity of the initial value η
0 from the exact value η
ex but the best accuracy that can be obtained with the method depends on the gap between η
0
(
T
) and η
ex
(
T
).Indeed, for each numerical simulation, we note thatη
k
(
T
)
=
η
0 for each k (see Figure 3 for the function η
k
(
t
) with (a) η
0
=
η
ex
(
0
)
=
2 and (b) η
0
=
η
ex
(
T
)
=
1
+
cos
⁡
(
T
)): the last value of η
k is never modified. For this kind of optimization problem, there is no final observation and the gradient (like the descent direction) is always zero for t
=
T. Then, the parameter cannot be well approximated near the final time T. Two methods will be proposed in the next subsection to solve this problem. Of course, this has no effect when the initialization is a continuous function with η
0
(
T
)
=
η
ex
(
T
) (see Figure 3(b)). From a physical point of view, the good choice for the initial value can be η
0 such that η
0
(
0
)
=
η
ex
(
0
) because the parameters are assumed known at ambient temperature (corresponding to t
=
0). It is why we use this initialization for the next examples.To compare the different gradient methods, we fix the initialization toη
0
=
η
ex
(
T
) and we present on Figure 4(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for the CGPR method with different optimization steps and on Figure 4(b) the convergence curves for the different gradient methods with optimal step. We conclude that the choice of the optimal step gives a significant improvement compared to other possibilities. As expected the conjugate gradient methods have a better convergence speed than the gradient method. Among the conjugate gradient methods, the Dai-Yuan method seems to be the better choice because it gives the better result for almost all iterations. This is in agreement with the literature results in optimization. It is why we present mainly the results for this method for the next examples.Now, we want to compare the distributed observations control problem and the boundary observations control problem. We present on Figure5 the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for the CGDY method and η
0
=
η
ex
(
0
). In this case, the convergence is a little bit better for the boundary control problem but it depends on the initialization. The important point is that the results are close and it is sufficient to measure the values of u only in x
=
0 and x
=
L.Figure 5
Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k for the distributed and boundary control problems.Remark 12.
To explore the difference between the control on the cold side of the domain (x
=
L) and the control on the hot side of the domain (x
=
0), we have tested the control on the convective heat transfer coefficient h
0 (the function η plays here the role of h
0). We have considered the same data as in the previous example except for the values of h
0 and h
1 that are exchanged η
ex
(
=
h
0
)
=
1
+
cos
⁡
(
t
)
,
h
1
=
1, and the related data u
f
(
t
)
=
[
1
+
(
1
+
e
2
(
t
+
L
)
)
/
(
1
+
cos
⁡
t
)
]
e
t
+
L
,
u
a
(
t
)
=
-
e
3
t, and we obtained similar results to the previous ones.


### 3.2.2. How to Improve the Result at Final Time?

We propose here two methods to overcome the difficulty to compute the solution of the inverse problem at the final time, the gradient of the functional being null for each step of the algorithm. The first possibility is to add a term with final observation to the functional and the second to deduce the value of the solution at final time by a linear approximation of previous values as presented in [24]. We present these methods for the distributed control problem.The first method consists of replacing the expression ofJ in (23) by
(29)
J
(
η
)
=
1
2
∥
u
-
u
obs
∥
L
2
(
𝒬
)
2
+
β
2
∥
u
(
T
)
-
u
ob
s
T
∥
L
2
(
Ω
)
2
,
where β is a positive constant. More precisely, 𝒟 is the identity function multiplied by β in (9). The expression of the gradient (25) is unchanged, but the definition of the adjoint problem (21) is modified because of the new value of 𝒟.The choice of the value ofβ is important. If it is too large the algorithm does not converge to the expected solution for the same initialization because, with only a final observation (problem of exact controllability [25, 26]), different values of η give solutions u of the direct problem with the same final value. If it is too small, there is no change with the new functional. As we can see in Figure 6, a good choice consists to give the final observation the same importance than each of the other observations with β being a value close to the time step δ
t
=
1
0
-
3. We use here the data of the example of Section 3.2.1 defined in the previous section.Example of Section3.2.1 with final observation: (a) log
⁡
10
(
errel
(
η
k
)
) versus k for different β and (b) η
k
(
t
) for some k with β
=
0.001.
(a)
(b)The second method consists of a linear approximation of the solution at final time by previous values of the solution. In [24], is proposed the following formula:
(30)
η
k
(
T
)
=
2
η
k
(
T
-
δ
t
)
-
η
k
(
T
-
2
δ
t
)
which is efficient if the time step is not too small. For δ
t
=
1
0
-
3, the value of η
k
(
T
) converges very slowly to the exact value; see Figure 7. This convergence can be accelerated with the following generalization of the formula:
(31)
η
k
(
T
)
=
2
η
k
(
T
-
n
l
c
δ
t
)
-
η
k
(
T
-
2
n
l
c
δ
t
)
with a small integer n
l
c which should not be too large to avoid instability. We present on Figure 8(a) the convergence curves for different values of n
l
c and on Figure 8(b)  η
k computed with n
l
c
=
2. We will use this value in the sequel.Figure 7
Example of Section3.2.1—linear approximation—η
k
(
t
) for some k with n
l
c
=
1.Example of Section3.2.1—linear approximation (a) log
⁡
10
(
errel
(
η
k
)
) versus k for different n
l
c and (b) η
k
(
t
) for some k with n
l
c
=
2.
(a)
(b)


### 3.2.3. Example with Radiation Conditions

For this second kind of tests, there is no more analytical expression for the observationu
obs which is now computed as the solution of the direct problem (3)-(4)-(5) with the exact value η
ex of the control function η. Then, solving the control problem from an initialization η
0 of η, we compute approximations η
k of η
ex using the CGDYOPT method.As in the previous section we detail the results for the control on the convective heat transfer coefficienth
1 and give a remark for the control on the convective heat transfer coefficient h
0.With the datah
0
=
1
,
h
1
=
1
+
cos
⁡
(
t
)
,
φ
0
=
φ
1
=
1
0
-
4, λ
(
u
)
=
1
+
u
,
κ
(
u
)
=
1
+
u
2
,
u
a
=
10
,
u
f
(
t
)
=
20
+
ln
⁡
(
t
)
,
u
0
=
u
a, we compute the observation u
obs as numerical solution of (3)-(4)-(5). It is a problem with radiation boundary conditions.The observation being known on the boundary and because of the time dependence ofh
0 and h
1, we can consider an optimization problem without radiation condition on the boundary. Indeed, the boundary conditions of the direct problem become
(32)
-
κ
(
u
(
0
)
)
∂
u
∂
x
(
0
)
=
H
0
(
u
f
-
u
(
0
)
)
,
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
∂
x
(
L
)
=
H
1
(
u
a
-
u
(
L
)
)
,
in
(
0
,
T
)
,
with
(33)
H
0
=
h
0
+
φ
0
(
u
f
+
u
obs
)
(
u
f
2
+
u
obs
2
)
,
H
1
=
h
1
+
φ
1
(
u
a
+
u
obs
)
(
u
a
2
+
u
obs
2
)
.The functionsh
0 and H
0, h
1 and H
1 are compared in Figure 9.Example of Section3.2.3—comparing data (a) h
0 and H
0 and (b) h
1 and H
1.
(a)
(b)In this example, the control functionη plays the role of H
1 with η
ex
=
H
1 and we solve the same optimization problem as in Example of Section 3.2.1 with F
=
0 and (
h
0
,
h
1
) replaced by (
H
0
,
H
1
). The functions u
obs and u
a being known, we can deduce φ
1 from H
1 and h
1 or h
1 from H
1 and φ
1.For the initializationη
0
=
η
ex
(
0
), we present on Figure 10(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for both kinds of problems (with distributed observations and boundary observations) and (b) the function η
k
(
t
) for some k (corresponding to the boundary observation). The results are good but contrary to what is expected, we note again a better result with the boundary observations (it depends on the initialization, e.g., with η
0
=
η
ex
(
T
) the method with distributed observations is a little bit better).Example of Section3.2.3 with n
l
c
=
2: (a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
k
(
t
) for some k with boundary observation.
(a)
(b)Remark 13.
We observed that the result is worse if we initialize our algorithm withη
0
=
η
ex
(
T
) (instead of η
0
=
η
ex
(
0
)). It can be explained by the fact that for η
0
=
η
ex
(
T
), η
k is not changed for t small because without a source term in (3), the value of H
1, used in the boundary condition at the point x
=
L, has no significant effect immediately on u. This comes from the fact that the source term u
f is only on the boundary at the point x
=
0, and there is a delay depending on the parameters of the equation. We checked that there is no more delay when the source term and the control are at the same point x
=
0.


### 3.2.4. Numerical Stability of the Algorithm

In real problems, the observationu
obs comes from experimental data with measurement errors and it is interesting to study the influence of these errors on the result given by the algorithm. As in the stability study presented in [24] for a linear problem, the noisy data are generated at each point number i of the mesh by the formula
(34)
u
~
obs
,
i
=
u
obs
,
i
(
1
+
δ
ζ
i
)
,
where δ gives the relative noise level and ζ
i is a random number obtained by a normal distribution of mean zero and unit standard deviation (Matlab function randn).On Figure11, we present the results for the example of Section 3.2.1 with boundary observation and final observation (β
=
10
-
2) (the results are similar with n
l
c
=
2 instead of final observation) and initialization η
0
(
t
)
=
η
ex
(
0
), δ
t
=
1
0
-
3
,
δ
x
=
1
0
-
2. For the 3 noise levels (1
%, 3
%, and 5
%), the computed solution remains close to the exact solution with a gap comparable to the noise level.Example of Section3.2.1 with β
=
0.01 and boundary observation—(a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
20
(
t
), for different values of noise level δ.
(a)
(b)The results are completely different for the example of Section3.2.3 (see Figure 12). For noise level δ
≥
0.0005, the computed solution does not remain close to the exact solution. The radiation condition on the boundary, with term (
|
u
a
|
3
u
a
-
|
u
|
3
u
), is very sensitive to the perturbations of the data.Figure 12
Example of Section3.2.3 with β
=
0.01 and boundary observation—η
7
(
t
) for different values of noise level δ.


## 3.2.1. Example without Radiation Conditions

In this example, the control functionη plays the role of h
1; the observation is the given functions u
obs
=
e
t
+
L
-
x and η
ex
(
=
h
1
)
=
1
+
cos
⁡
(
t
). The other operators and data of the model are given by h
0
=
1
,
φ
0
=
φ
1
=
0
,
λ
(
u
)
=
1
+
u, κ
(
u
)
=
1
+
u
2, u
0
(
x
)
=
e
L
-
x
,
u
f
(
t
)
=
[
2
+
e
2
(
t
+
L
)
]
e
t
+
L
,
u
a
(
t
)
=
[
1
-
(
1
+
e
2
t
)
/
(
1
+
cos
⁡
t
)
]
e
t and we have added the following function F
(
x
,
t
)
=
e
2
(
t
+
L
-
x
)
(
1
-
3
e
t
+
L
-
x
) such that (3) becomes
(28)
λ
(
u
)
∂
u
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
∂
x
)
=
F
to ensure that u
obs is the analytical solution of the direct problem (corresponding to η
ex).Remark 11.
We want to emphasize on the fact that the addition of the artificial right-hand memberF does not change the formulation of the adjoint problem, the expression of the gradient of the objective function, and the optimization algorithm.To validate our approach, we have tested all the gradient methods described above. With all these gradient methods we get computed valuesη
k which converge towards the exact value η
ex: the relative error (27) decreases when k increases (see Figures 2, 3, and 4). But the convergence speed depends on the method and on the initial value. So, we are going to discuss these two points for the distributed control problem (23) and then compare with the results of the boundary control problems (24).Figure 2
Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k for different initializations with the CGDYOPT method.Example of Section3.2.1—η
k
(
t
) for some k with different initializations: (a) η
0
=
η
ex
(
0
)
=
2 and (b) η
0
=
η
ex
(
T
)
=
1
+
cos
⁡
(
T
).
(a)
(b)Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k with (a) CGPR method and different optimization steps and (b) an optimal step and different gradient methods.
(a)
(b)First, to use the optimization algorithm we have chosen an initial valueη
0 for the control function η. We present on Figure 2 the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for different initializations and on Figure 3 the function η
k
(
t
) for some k (computed with the CGDYOPT method: Dai-Yuan conjugate gradient with optimal step) and for these different initializations: (a) η
0
=
2.2, (b) η
0
=
η
ex
(
0
), (c) η
0
=
1.7, (d) η
0
=
η
ex
(
T
), (e) η
0
=
1, and (f) η
0
=
[
η
ex
(
T
)
-
η
ex
(
0
)
]
(
t
/
T
)
+
η
ex
(
0
). As expected, the method is convergent for each initialization but the convergence speed depends on the proximity of the initial value η
0 from the exact value η
ex but the best accuracy that can be obtained with the method depends on the gap between η
0
(
T
) and η
ex
(
T
).Indeed, for each numerical simulation, we note thatη
k
(
T
)
=
η
0 for each k (see Figure 3 for the function η
k
(
t
) with (a) η
0
=
η
ex
(
0
)
=
2 and (b) η
0
=
η
ex
(
T
)
=
1
+
cos
⁡
(
T
)): the last value of η
k is never modified. For this kind of optimization problem, there is no final observation and the gradient (like the descent direction) is always zero for t
=
T. Then, the parameter cannot be well approximated near the final time T. Two methods will be proposed in the next subsection to solve this problem. Of course, this has no effect when the initialization is a continuous function with η
0
(
T
)
=
η
ex
(
T
) (see Figure 3(b)). From a physical point of view, the good choice for the initial value can be η
0 such that η
0
(
0
)
=
η
ex
(
0
) because the parameters are assumed known at ambient temperature (corresponding to t
=
0). It is why we use this initialization for the next examples.To compare the different gradient methods, we fix the initialization toη
0
=
η
ex
(
T
) and we present on Figure 4(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for the CGPR method with different optimization steps and on Figure 4(b) the convergence curves for the different gradient methods with optimal step. We conclude that the choice of the optimal step gives a significant improvement compared to other possibilities. As expected the conjugate gradient methods have a better convergence speed than the gradient method. Among the conjugate gradient methods, the Dai-Yuan method seems to be the better choice because it gives the better result for almost all iterations. This is in agreement with the literature results in optimization. It is why we present mainly the results for this method for the next examples.Now, we want to compare the distributed observations control problem and the boundary observations control problem. We present on Figure5 the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for the CGDY method and η
0
=
η
ex
(
0
). In this case, the convergence is a little bit better for the boundary control problem but it depends on the initialization. The important point is that the results are close and it is sufficient to measure the values of u only in x
=
0 and x
=
L.Figure 5
Example of Section3.2.1—log
⁡
10
(
errel
(
η
k
)
) versus k for the distributed and boundary control problems.Remark 12.
To explore the difference between the control on the cold side of the domain (x
=
L) and the control on the hot side of the domain (x
=
0), we have tested the control on the convective heat transfer coefficient h
0 (the function η plays here the role of h
0). We have considered the same data as in the previous example except for the values of h
0 and h
1 that are exchanged η
ex
(
=
h
0
)
=
1
+
cos
⁡
(
t
)
,
h
1
=
1, and the related data u
f
(
t
)
=
[
1
+
(
1
+
e
2
(
t
+
L
)
)
/
(
1
+
cos
⁡
t
)
]
e
t
+
L
,
u
a
(
t
)
=
-
e
3
t, and we obtained similar results to the previous ones.


## 3.2.2. How to Improve the Result at Final Time?

We propose here two methods to overcome the difficulty to compute the solution of the inverse problem at the final time, the gradient of the functional being null for each step of the algorithm. The first possibility is to add a term with final observation to the functional and the second to deduce the value of the solution at final time by a linear approximation of previous values as presented in [24]. We present these methods for the distributed control problem.The first method consists of replacing the expression ofJ in (23) by
(29)
J
(
η
)
=
1
2
∥
u
-
u
obs
∥
L
2
(
𝒬
)
2
+
β
2
∥
u
(
T
)
-
u
ob
s
T
∥
L
2
(
Ω
)
2
,
where β is a positive constant. More precisely, 𝒟 is the identity function multiplied by β in (9). The expression of the gradient (25) is unchanged, but the definition of the adjoint problem (21) is modified because of the new value of 𝒟.The choice of the value ofβ is important. If it is too large the algorithm does not converge to the expected solution for the same initialization because, with only a final observation (problem of exact controllability [25, 26]), different values of η give solutions u of the direct problem with the same final value. If it is too small, there is no change with the new functional. As we can see in Figure 6, a good choice consists to give the final observation the same importance than each of the other observations with β being a value close to the time step δ
t
=
1
0
-
3. We use here the data of the example of Section 3.2.1 defined in the previous section.Example of Section3.2.1 with final observation: (a) log
⁡
10
(
errel
(
η
k
)
) versus k for different β and (b) η
k
(
t
) for some k with β
=
0.001.
(a)
(b)The second method consists of a linear approximation of the solution at final time by previous values of the solution. In [24], is proposed the following formula:
(30)
η
k
(
T
)
=
2
η
k
(
T
-
δ
t
)
-
η
k
(
T
-
2
δ
t
)
which is efficient if the time step is not too small. For δ
t
=
1
0
-
3, the value of η
k
(
T
) converges very slowly to the exact value; see Figure 7. This convergence can be accelerated with the following generalization of the formula:
(31)
η
k
(
T
)
=
2
η
k
(
T
-
n
l
c
δ
t
)
-
η
k
(
T
-
2
n
l
c
δ
t
)
with a small integer n
l
c which should not be too large to avoid instability. We present on Figure 8(a) the convergence curves for different values of n
l
c and on Figure 8(b)  η
k computed with n
l
c
=
2. We will use this value in the sequel.Figure 7
Example of Section3.2.1—linear approximation—η
k
(
t
) for some k with n
l
c
=
1.Example of Section3.2.1—linear approximation (a) log
⁡
10
(
errel
(
η
k
)
) versus k for different n
l
c and (b) η
k
(
t
) for some k with n
l
c
=
2.
(a)
(b)


## 3.2.3. Example with Radiation Conditions

For this second kind of tests, there is no more analytical expression for the observationu
obs which is now computed as the solution of the direct problem (3)-(4)-(5) with the exact value η
ex of the control function η. Then, solving the control problem from an initialization η
0 of η, we compute approximations η
k of η
ex using the CGDYOPT method.As in the previous section we detail the results for the control on the convective heat transfer coefficienth
1 and give a remark for the control on the convective heat transfer coefficient h
0.With the datah
0
=
1
,
h
1
=
1
+
cos
⁡
(
t
)
,
φ
0
=
φ
1
=
1
0
-
4, λ
(
u
)
=
1
+
u
,
κ
(
u
)
=
1
+
u
2
,
u
a
=
10
,
u
f
(
t
)
=
20
+
ln
⁡
(
t
)
,
u
0
=
u
a, we compute the observation u
obs as numerical solution of (3)-(4)-(5). It is a problem with radiation boundary conditions.The observation being known on the boundary and because of the time dependence ofh
0 and h
1, we can consider an optimization problem without radiation condition on the boundary. Indeed, the boundary conditions of the direct problem become
(32)
-
κ
(
u
(
0
)
)
∂
u
∂
x
(
0
)
=
H
0
(
u
f
-
u
(
0
)
)
,
in
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
∂
x
(
L
)
=
H
1
(
u
a
-
u
(
L
)
)
,
in
(
0
,
T
)
,
with
(33)
H
0
=
h
0
+
φ
0
(
u
f
+
u
obs
)
(
u
f
2
+
u
obs
2
)
,
H
1
=
h
1
+
φ
1
(
u
a
+
u
obs
)
(
u
a
2
+
u
obs
2
)
.The functionsh
0 and H
0, h
1 and H
1 are compared in Figure 9.Example of Section3.2.3—comparing data (a) h
0 and H
0 and (b) h
1 and H
1.
(a)
(b)In this example, the control functionη plays the role of H
1 with η
ex
=
H
1 and we solve the same optimization problem as in Example of Section 3.2.1 with F
=
0 and (
h
0
,
h
1
) replaced by (
H
0
,
H
1
). The functions u
obs and u
a being known, we can deduce φ
1 from H
1 and h
1 or h
1 from H
1 and φ
1.For the initializationη
0
=
η
ex
(
0
), we present on Figure 10(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for both kinds of problems (with distributed observations and boundary observations) and (b) the function η
k
(
t
) for some k (corresponding to the boundary observation). The results are good but contrary to what is expected, we note again a better result with the boundary observations (it depends on the initialization, e.g., with η
0
=
η
ex
(
T
) the method with distributed observations is a little bit better).Example of Section3.2.3 with n
l
c
=
2: (a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
k
(
t
) for some k with boundary observation.
(a)
(b)Remark 13.
We observed that the result is worse if we initialize our algorithm withη
0
=
η
ex
(
T
) (instead of η
0
=
η
ex
(
0
)). It can be explained by the fact that for η
0
=
η
ex
(
T
), η
k is not changed for t small because without a source term in (3), the value of H
1, used in the boundary condition at the point x
=
L, has no significant effect immediately on u. This comes from the fact that the source term u
f is only on the boundary at the point x
=
0, and there is a delay depending on the parameters of the equation. We checked that there is no more delay when the source term and the control are at the same point x
=
0.


## 3.2.4. Numerical Stability of the Algorithm

In real problems, the observationu
obs comes from experimental data with measurement errors and it is interesting to study the influence of these errors on the result given by the algorithm. As in the stability study presented in [24] for a linear problem, the noisy data are generated at each point number i of the mesh by the formula
(34)
u
~
obs
,
i
=
u
obs
,
i
(
1
+
δ
ζ
i
)
,
where δ gives the relative noise level and ζ
i is a random number obtained by a normal distribution of mean zero and unit standard deviation (Matlab function randn).On Figure11, we present the results for the example of Section 3.2.1 with boundary observation and final observation (β
=
10
-
2) (the results are similar with n
l
c
=
2 instead of final observation) and initialization η
0
(
t
)
=
η
ex
(
0
), δ
t
=
1
0
-
3
,
δ
x
=
1
0
-
2. For the 3 noise levels (1
%, 3
%, and 5
%), the computed solution remains close to the exact solution with a gap comparable to the noise level.Example of Section3.2.1 with β
=
0.01 and boundary observation—(a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
20
(
t
), for different values of noise level δ.
(a)
(b)The results are completely different for the example of Section3.2.3 (see Figure 12). For noise level δ
≥
0.0005, the computed solution does not remain close to the exact solution. The radiation condition on the boundary, with term (
|
u
a
|
3
u
a
-
|
u
|
3
u
), is very sensitive to the perturbations of the data.Figure 12
Example of Section3.2.3 with β
=
0.01 and boundary observation—η
7
(
t
) for different values of noise level δ.


## 3.3. Simulation on Real Data

We now apply our method to a realistic problem. To that purpose, we take similarly data as in [4, 16]. We consider a gypsum plasterboard, with thickness L
=
4cm, exposed to fire in x
=
0, during T
=
1200 s. The surroundings temperature is u
a
=
2
0
∘C and the furnace temperature is u
f
=
u
a
+
345
log
⁡
10
(
1
+
8
t
/
60
) (see Figure 13). The initial temperature is u
0
=
u
a. The thermal conductivity is given by
(35)
κ
(
u
)
=
{
0.25
if
u
<
112
.
5
°
C
,
0.1116
+
6.578
·
1
0
-
5
u
if
112
.
5
∘
C
≤
u
<
393.7
5
°
C
,
0.0277
+
2.788
·
1
0
-
4
u
if
393.7
5
∘
C
≤
u
<
797
.
2
°
C
,
-
0.91323
+
1.459
·
1
0
-
3
u
if
797
.
2
°
C
≤
u
.Figure 13
Furnace temperature versus time.We haveλ
(
u
)
=
ρ
(
u
)
C
p
(
u
), where the density is defined by
(36)
ρ
(
u
)
=
{
648
if
u
<
12
0
°
C
,
534.6
if
12
0
°
C
≤
u
≤
65
0
∘
C
,
500.6
if
65
0
°
C
<
u
and the specific heat is defined with a peak between 100 and 17
0
°
C and a value corresponding to the total amount of energy of the two dehydration reactions, 669 kJ:
(37)
C
p
(
u
)
=
{
6.69
·
1
0
4
7
if
10
0
∘
C
≤
u
≤
17
0
°
C
,
950
else
.We supposed that the convective heat transfer coefficienth
0
=
25 and the effective emissivity ϵ
0
=
0.82 are known in x
=
L. We deduce that φ
0
=
σ
ϵ
0 with σ
=
5.6697
·
1
0
-
8
W
/
m
2
K
4, the Stefan-Boltzmann's constant.We want to estimateh
1, supposing that ϵ
1
=
0. The control function η plays here the role of h
1. We choose the following exact expression for h
1: η
ex
(
=
h
1
)
=
5
+
cos
⁡
(
t
/
T
). Then, we can compute the observation u
obs as the solution of the direct problem (3), solve the control problem (9) (with the cost function and its gradient given by (23) or (24) and (25)) by the optimization algorithm (using the CGDYOPT method), and compute error estimations on the control function η at each optimization iteration.For the space stepδ
x
=
0.1cm, the time step δ
t
=
1s, and the initialization η
0
=
η
ex
(
0
), we present in Figure 14(a) the convergence curves, log
⁡
10
(
errel
(
η
k
)
) versus k, for both kinds of problems (with distributed observation and boundary observation) and (b) the function η
k
(
t
) for some k (corresponding to the boundary observation). Despite the fact that the nonlinear parameters λ
(
u
) and κ
(
u
) are not regular, we get a good convergence of the method.Simulation on real data withn
l
c
=
2: (a) log
⁡
10
(
errel
(
η
k
)
) versus k and (b) η
k
(
t
) for some k with boundary observation.
(a)
(b)


## 4. Commentary and Conclusion

Optimal boundary heating control strategies for time-dependent thermal problems in spatially 1-dimensional domains, related to dehydration of gypsum plasterboards exposed to fire, are developed. Optimal heating strategies and calibration of process models (parameter identification) are obtained as solutions of certain minimization problems and are computed from conjugate gradient method by considering the following well-known descent direction methods: Fletcher-Reeves, Polak-Ribière, Hestenes-Stiefel, and Dai and Yuan. Other choices of control variables can be envisaged. The numerical results show the efficiency of the developed method.This method can be applied without additional argument to estimate emissivities. Indeed, we can introduceϕ
=
(
φ
0
,
φ
1
) and the following optimal control problem:Find(
u
,
h
,
ϕ
) such that the following objective function:
(38)
J
(
h
,
ϕ
)
=
1
2
∫
‍
∫
𝒬
‍
|
𝒞
(
u
-
u
obs
)
|
2
d
x
d
t
+
1
2
∫
Ω
‍
|
𝒟
(
u
(
T
)
-
u
obs
T
)
|
2
d
x
+
1
2
∫
0
T
‍
γ
0
|
u
(
0
)
-
u
obs
0
|
2
+
γ
1
|
u
(
L
)
-
u
obs
L
|
2
d
t
+
1
2
∫
0
T
‍
|
ℳ
h
|
2
d
t
+
1
2
∫
0
T
‍
|
𝒩
ϕ
|
2
d
tis minimized with respect to(
h
,
ϕ
)
∈
𝒴
ad
=
V
ad
×
W
ad subject to (3)–(5), where 𝒴
ad is the set of admissible controls with
(39)
W
ad
=
{
ϕ
=
(
φ
0
,
φ
1
)
∈
(
L
2
(
0
,
T
)
)
2
such
that
(
L
2
(
0
,
T
)
)
2
0
≤
a
i
≤
φ
i
≤
b
i
a
.
e
.
in
(
0
,
T
)
,
for
i
=
0,1
}
,
V
ad
=
{
h
=
(
h
0
,
h
1
)
∈
(
L
2
(
0
,
T
)
)
2
such
that
(
L
2
(
0
,
T
)
)
2
0
≤
c
i
≤
h
i
≤
d
i
a
.
e
.
in
(
0
,
T
)
,
for
i
=
0,1
}
,
where u
=
ℱ
(
h
,
ϕ
), γ
0 and γ
1 are predefined nonnegative weights, ℳ
=
diag
⁡
(
α
0
,
α
1
), and 𝒩
=
diag
⁡
(
β
0
,
β
1
) are predefined nonnegative weights such that ∑
i
=
0,1
‍
(
α
i
+
β
i
)
≠
0. The operators 𝒞 and 𝒟 are unbounded operators on L
2
(
Ω
) satisfying (∀
v
∈
H
1
(
Ω
)):
(40)
∥
𝒞
v
∥
L
∞
(
Ω
)
2
≤
δ
1
∥
v
∥
L
2
(
Ω
)
2
+
δ
2
∥
v
∥
H
1
(
Ω
)
2
,
∥
𝒟
v
∥
L
2
(
Ω
)
2
≤
δ
1
∥
v
∥
L
2
(
Ω
)
2
+
δ
2
∥
v
∥
H
1
(
Ω
)
2
,
with δ
1
+
δ
2
>
0 and δ
i
≥
0, for i
=
1,2.With the same kind of calculus as previously shown, we obtain a more general form of Theorem8.Theorem 14.
If J attains a (local) minimum at a point(
h
*
,
ϕ
*
)
∈
𝒴
a
d, then
(41)
h
0
*
=
Π
[
c
0
,
d
0
]
(
(
-
u
f
+
u
*
(
0
)
)
u
~
*
(
0
)
α
0
)
,
h
1
*
=
Π
[
c
1
,
d
1
]
(
(
-
u
a
+
u
*
(
L
)
)
u
~
*
(
L
)
α
1
)
,
φ
0
*
=
Π
[
a
0
,
b
0
]
(
(
-
|
u
f
|
3
u
f
+
|
u
*
(
0
)
|
3
u
*
(
0
)
)
u
~
*
(
0
)
β
0
)
,
φ
1
*
=
Π
[
a
1
,
b
1
]
(
(
-
|
u
a
|
3
u
a
+
|
u
*
(
L
)
|
3
u
*
(
L
)
)
u
~
*
(
L
)
β
1
)
,
or in the variational inequality formulation (for all (
h
,
ϕ
)
∈
𝒴
a
d)
(42)
∫
0
T
‍
(
h
1
-
h
1
*
)
(
(
u
a
-
u
*
(
L
)
)
u
~
*
(
L
)
+
α
1
h
1
*
)
d
t
+
∫
0
T
‍
(
φ
1
-
φ
1
*
)
×
(
(
|
u
a
|
3
u
a
-
|
u
*
(
L
)
|
3
u
*
(
L
)
)
u
~
*
(
L
)
+
β
1
φ
1
*
)
d
t
+
∫
0
T
‍
(
h
0
-
h
0
*
)
(
(
u
f
-
u
*
(
0
)
)
u
~
*
(
0
)
+
α
0
h
0
*
)
d
t
+
∫
0
T
‍
(
φ
0
-
φ
0
*
)
(
(
|
u
f
|
3
u
f
-
|
u
*
(
0
)
|
3
u
*
(
0
)
)
u
~
*
(
0
)
llllllllllllllllllllllllllllllllllllll
+
β
0
φ
0
*
(
|
u
f
|
3
u
f
-
|
u
*
(
0
)
|
3
u
*
(
0
)
)
)
d
t
≥
0
,
where u
*
=
ℱ
(
h
*
,
ϕ
*
) and   u
~
*
=
ℱ
*
(
h
*
,
ϕ
*
) with the function u
~
=
ℱ
*
(
h
,
ϕ
) which is the unique solution of the adjoint (costate) problem (with initial value given at final time T) given by
(43)
-
λ
(
u
)
∂
u
~
∂
t
-
∂
∂
x
(
κ
(
u
)
∂
u
~
∂
x
)
+
G
(
u
)
∂
u
~
∂
x
=
𝒞
*
𝒞
(
u
-
u
o
b
s
)
,
hhhhhhhhhh
i
n
𝒬
,
-
κ
(
u
(
0
)
)
∂
u
~
∂
x
(
0
)
=
-
(
h
0
+
4
φ
0
|
u
(
0
)
|
3
)
u
~
+
γ
0
(
u
(
0
)
-
u
o
b
s
0
)
,
i
n
(
0
,
T
)
,
κ
(
u
(
L
)
)
∂
u
~
∂
x
(
L
)
=
-
(
h
1
+
4
φ
1
|
u
(
L
)
|
3
)
u
~
+
γ
1
(
u
(
L
)
-
u
o
b
s
L
)
,
i
n
(
0
,
T
)
,
λ
(
u
(
·
,
T
)
)
u
~
(
·
,
T
)
=
𝒟
*
𝒟
(
u
(
·
,
T
)
-
u
o
b
s
T
)
,
i
n
Ω
.
Moreover, the gradient ofJ at any element (
h
,
ϕ
) of 𝒴
a
d can be given by
(44)
∂
J
∂
h
(
h
,
ϕ
)
=
(
(
u
f
-
u
(
0
)
)
u
~
(
0
)
+
α
0
h
0
,
(
u
f
-
u
(
0
)
)
(
u
a
-
u
(
L
)
)
u
~
(
L
)
+
α
1
h
1
)
,
∂
J
∂
ϕ
(
h
,
ϕ
)
=
(
(
|
u
f
|
3
u
f
-
|
u
(
0
)
|
3
u
(
0
)
)
u
~
(
0
)
+
β
0
φ
0
,
(
|
u
f
|
3
u
f
-
|
u
(
0
)
|
3
u
(
0
)
)
(
|
u
a
|
3
u
a
-
|
u
(
L
)
|
3
u
(
L
)
)
u
~
(
L
)
+
β
1
φ
1
)
.It would be interesting to use the here developed method with observations coming from experimental data and a more complete description of the dehydration process including reactions for high temperatures. It is clear that to get even closer to a more realistic calculation it is necessary to study in the future the mathematical analysis and numerical validation of the controllability of the system, such as what is presented, for example, for linear parabolic problems in [27] and the references therein. In order to take into account the influence of noises and fluctuations, we can also study robust control problems by using the approach developed in Belmiloudi's Book [19].



---

*Source: 634712-2014-04-13.xml*