# Novel Methodologies for Solving the Inverse Unsteady Heat Transfer Problem
of Estimating the Boundary Heat Flux in Continuous Casting Molds

Umberto Emil Morelli1,2,3*<EMAIL_ADDRESS>,  Patricia Barral1,2 ,
Peregrina Quintela1,2 ,  Gianluigi Rozza3 and  Giovanni Stabile3 1Universidade
de Santiago de Compostela, Santiago de Compostela, Spain  2 Technological
Institute for Industrial Mathematics (ITMATI), Santiago de Compostela, Spain
3Scuola Internazionale Superiore di Studi Avanzati (SISSA), Trieste, Italy

###### Abstract.

In this work, we investigate the estimation of the transient mold-slab heat
flux in continuous casting molds given some thermocouples measurements in the
mold plates. Mathematically, we can see this problem as the estimation of a
Neumann boundary condition given pointwise state observations in the interior
of the domain. We formulate it in a deterministic inverse problem setting.
After introducing the industrial problem, we present the mold thermal model
and related assumptions. Then, we formulate the boundary heat flux estimation
problem in a deterministic inverse problem setting using a sequential approach
according to the sequentiality of the temperature measurements. We consider
different formulations of the inverse problem. For each one, we develop novel
direct methodologies exploiting a space parameterization of the heat flux and
the linearity of the mold model. We construct these methods to be divided into
a computationally expensive offline phase that can be computed before the
process starts, and a cheaper online phase to be performed during the casting
process. To conclude, we test the performance of the proposed methods in two
benchmark cases.

###### Key words and phrases:

Inverse Problem, Heat Transfer, Continuous Casting, Optimal Control, Data
Assimilation, Boundary Condition Estimation

Funded by the European Union’s Horizon 2020 research and innovation programme
under the Marie Skaodowska-Curie Grant Agreement No. 765374. It also was
partially supported by the Ministry of Economy, Industry and Competitiveness
through the Plan Nacional de I+D+i (MTM2015-68275-R), by the Agencia Estatal
de Investigacion through project [PID2019-105615RB-I00/ AEI /
10.13039/501100011033], by the European Union Funding for Research and
Innovation - Horizon 2020 Program - in the framework of European Research
Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD
project 681447 ”Advanced Reduced Order Methods with Applications in
Computational Fluid Dynamics” and INDAM-GNCS project ”Advanced intrusive and
non-intrusive model order reduction techniques and applications”, 2019.

##  1\. Introduction

Most of the steel produced everyday worldwide is made by Continuous Casting
(CC).[1] Continuous casters have been around for many decades now and a long
sequence of improvements have increased through the years their productivity
(i.e. the casting speed) and the quality of the casted products.

To motivate and contextualize this research, we provide a brief overview on
the CC process. It starts by tapping the liquid metal from the ladle into the
tundish. In the tundish, the metal flow is regulated and smoothed. Through the
Submerged Entry Nozzle (SEN), the metal is drained into a mold. The role of
the mold is to cool down the steel until it has a solid skin which is thick
and cool enough to be supported by rollers in the secondary cooling region.

At the outlet of the mold, the metal is still molten in its inner region.
Supported by rollers, it is cooled until complete solidification by directly
spraying water over it. At the end of this secondary cooling region, the
casting is completed. This is just a brief overview on the CC process. We
refer the interested reader to Irving’s monograph on the subject.[2]

In this work, we focus on CC of thin slabs, i.e. slabs with rectangular cross
section with thickness smaller than 70 mm and width between 1 and 1.5 m.
Thanks to the small thickness, the solidification in the slab is relatively
fast, consequently the casting speed is generally high, between 7 and 14
meters per minute.

Thin slab molds are made of four different plates: two wide plates and two
lateral plates, all made of copper (see Figure 1). In general, lateral plates
can be moved or changed to modify the slab section dimensions. The geometry of
these plates is more complex than one can expect: they have drilled channels
where the cooling water flows, slots in the outside face for thermal
expansion, thermocouples, and fastening bolts. To compensate the shrinkage of
the slab with the cooling and minimize the gap, the molds are tapered.
Moreover, the upper portion of the mold forms a funnel to accommodate the SEN.

Figure 1. Schematic of a horizontal section of the mold (the casting direction
is perpendicular to the image).

Due to the high casting speed and the related strong thermal gradient, several
complex and coupled phenomena related to steel flow, solidification,
mechanics, and heat transfer appear in the mold region. This complexity makes
the mold the most critical part of the CC process. Here, safety and
productivity issues must be addressed.

For example, a common issue is the sticking of the steel to the mold. In this
case, it is essential to quickly detect the problem and reduce the casting
speed, otherwise it can lead to dangerous events that could force the shutdown
of the caster. Less frequent but more catastrophic events are the liquid
break-out and the excessive increase of the mold temperature. The former is
due to a non-uniform cooling of the metal with the skin being so thin to
break. The latter is generally considered the most dangerous event in a
casting plant. In fact, if the mold temperature is high enough to cause the
boiling of the cooling water, we have a dramatic decrease in the heat
extraction. Then, the temperature in the mold quickly rises, that could cause
the melting of the mold itself. Both these incidents are very dangerous and
costly. In fact, they generally require the shutdown of the caster, the
substitution of expensive components and an extended turnaround.

For all these reasons, the early detection of problems in the mold is crucial
for a safe and productive operation of continuous casters. Their detection
becoming more difficult as casting speed (thus productivity) of the casters
increases.

Until now, operators faced all these problems by equipping the molds with
sensors. Among other parameters, they measure the pointwise temperature of the
mold by thermocouples (see Figure 1) and the cooling water temperature as well
as its flow at the inlet and outlet of the cooling system. On one hand,
thermocouples temperatures are used to have insight of the mold temperature
field. On the other, the water temperature rise is used to approximate the
heat extracted from the steel.

This approach allowed to run continuous casters for decades. Nevertheless, it
has several drawbacks: it relies on the experience of operators, gives very
limited information about the heat flux at the mold-slab interface, and is
customized for each geometry so it requires new effort to be applied to new
designs. So, with the always increasing casting speed of modern casters, a new
and more reliable tool for analyzing the mold behavior is necessary.

According to CC operators and designers, knowing the local heat flux between
mold and slab is the most important information in monitoring the mold.
Moreover, we should estimate it in real-time for the early detection of issues
and a proper monitoring. By considering the mold plates to be our domain and
focusing our interest on its thermal behavior, the mold-slab heat flux can be
seen as a Neumann Boundary Condition (BC) in the model. To compute its value,
we pose the following inverse problem: given the temperature measurements
provided by the thermocouples, estimate the boundary heat flux at the mold-
slab interface. In a previous publication [3], we developed a novel
methodology for the solution of this problem using a steady-state mold model.
The present work is an extension of the previous one considering the more
challenging unsteady-state case.

After deriving the mold heat transfer model in Section 2, we discuss in
Section 3 the steel-mold heat flux estimation problem and propose novel
methodologies for its solution. Finally, we design in Section 4 some numerical
benchmark test cases that we use to study the performance of the proposed
inverse solvers.

##  2\. Mathematical Model

A detailed description of the physical phenomena that occur in the mold region
of a caster can be found in our previous work[3]. Here, we only mention that
the physical phenomena happening in the interior of the mold are extremely
complex and tightly coupled (thermodynamic reactions, multiphase flow, free
liquid surfaces and interfaces, solidification, etc.). Then, monitoring the
casting by simulating all of them from the SEN to the secondary cooling region
would be extremely complex and computationally expensive to deal with,
especially for real-time applications.

According to CC operators, to monitor the mold behavior it is sufficient to
know the mold-slab heat flux. Then, given the mold plate physical properties,
its geometry and the cooling water temperature, our approach is to solve an
inverse problem having as input data the temperature measurements made by the
thermocouples that are buried inside the mold plates.

As mentioned, the mold solid plates are the computational domain while the
mold-slab heat flux is a Neumann BC on a portion of its boundary to be
determined as solution of an inverse problem. Then, the direct problem
corresponds to a model for the heat transfer in the mold plates. In the rest
of this section, we describe the mold thermal model that we use in the present
investigation and the related assumptions.

In modeling the thermal behavior of the mold, we consider the following well
established assumptions[3]:

  * •

The copper mold is assumed a homogeneous and isotropic solid material.

  * •

The thermal expansion of the mold and its mechanical distortion are
negligible.

  * •

The material properties are assumed constant.

  * •

The boundaries in contact with air are assumed adiabatic.

  * •

The heat transmitted by radiation is neglected.

  * •

The cooling water temperature is known at the inlet and outlet of the cooling
system. Moreover, it is assumed to be constant in time and linear with respect
to the $z$ coordinate (see Figure 2).

  * •

No boiling in the water is assumed.

We refer to our previous work[3] for the motivations related to these
assumptions.

According to these assumptions, we consider in the following an unsteady-state
three-dimensional heat conduction model posed on the (solid) copper mold, with
a convective BC in the portion of the boundary in contact with the cooling
water, a Neumann BC in the portion of the boundary in contact with the steel,
and adiabatic BC in the portion of the boundary in contact with air.

After introducing some notation and the computational domain, we devote the
present section to the formulation of the mold model. As common when dealing
with inverse problems, we refer to it as the direct problem. We conclude this
section by discussing its numerical discretization.

###  2.1. Computational Domain and Notation

Consider a solid domain, $\Omega$, which is assumed to be an open Lipschitz
bounded subset of ${\rm I\\!R}^{3}$, with smooth boundary $\Gamma$ (see Figure
2). Let $\Gamma=\Gamma_{s_{in}}\cup\Gamma_{s_{ex}}\cup\Gamma_{sf}$ where
$\mathring{\Gamma}_{s_{in}}$, $\mathring{\Gamma}_{s_{ex}}$ and
$\mathring{\Gamma}_{sf}$ are disjoint sets. Moreover, given $t_{f}\in{\rm
I\\!R}^{+}$, we consider the time domain $(0,t_{f}]$. The Eulerian Cartesian
coordinate vector is denoted by $\mathbf{x}\in\Omega$ and
$\mathbf{n}(\mathbf{x})$ represents the unit normal vector that is directed
outwards from $\Omega$ at point $\mathbf{x}\in\Gamma$.

Figure 2. Schematic of the mold domain, $\Omega$, and its boundaries (images
taken from Morelli et al.[3]).

In this setting, $\Omega$ corresponds to the region of the space occupied by
the mold. The interface between the mold and the cooling system is denoted by
$\Gamma_{sf}$. While $\Gamma_{s_{in}}$ is the portion of the mold boundary in
contact with the solidifying steel. Finally, we denote the remaining part of
the mold boundary with $\Gamma_{s_{ex}}$.

###  2.2. Direct Problem

We shall assume all along the following assumptions on the data:

  1. (H1)

The thermal conductivity is constant and strictly positive: $k_{s}\in{\rm
I\\!R}^{+}$.

  2. (H2)

There is no heat source inside the mold domain.

  3. (H3)

The density and specific heat are constant and strictly positive: $\rho\in{\rm
I\\!R}^{+}$, $C_{p}\in{\rm I\\!R}^{+}$.

  4. (H4)

The heat transfer coefficient on $\Gamma_{sf}$ is constant and strictly
positive: $h\in{\rm I\\!R}^{+}$.

  5. (H5)

The cooling water temperature, $T_{f}$, is known, constant in time, and
belongs to $L^{q}(\Gamma_{sf})$.

  6. (H6)

The initial temperature, $T_{0}$, is known and belongs to $L^{2}(\Omega)$.

  7. (H7)

The steel-mold heat flux, $g$, belongs to
$L^{r}(0,t_{f};L^{q}(\Gamma_{s_{in}}))$.

In (H5), (H7) we assume that $r,q\in(2,+\infty)$ and

(2.1) $\frac{1}{r}+\frac{1}{q}<\frac{1}{2}.$

Notice that it implies $r,q>2$.

Under the assumptions (H1)-(H7), we propose the following three-dimensional,
unsteady-state, heat conduction model

######  Problem 2.1.

Find $T$ such that

(2.2) $\rho C_{p}\frac{\partial T}{\partial t}-k_{s}\Delta T=0,\quad\text{in
}\Omega\times(0,t_{f}],$

with BCs and Initial Condition (IC)

(2.3) $\begin{cases}-k_{s}\nabla T\cdot\mathbf{n}=g&\text{on
}\Gamma_{s_{in}}\times(0,t_{f}]\\\ -k_{s}\nabla T\cdot\mathbf{n}=0&\text{on
}\Gamma_{s_{ex}}\times(0,t_{f}],\\\ -k_{s}\nabla
T\cdot\mathbf{n}=h(T-T_{f})&\text{on }\Gamma_{sf}\times(0,t_{f}],\\\
T(\cdot,0)=T_{0}&\text{in }\Omega.\end{cases}$

A weak solution is now defined by testing against a smooth function and
formally integrating by parts.

######  Definition 2.1.

We say that a function $T\in C([0,t_{f}];L^{2}(\Omega))\cap
L^{2}(0,t_{f};H^{1}(\Omega))$ is a weak solution of Problem 2.3 on $[0,t_{f}]$
for some $t_{f}>0$ if

(2.4) $\displaystyle-\rho
C_{p}\int_{0}^{t_{f}}\int_{\Omega}T(\mathbf{x},t)\frac{\partial\psi(\mathbf{x},t)}{\partial
t}d\mathbf{x}dt+k_{s}\int_{0}^{t_{f}}\int_{\Omega}\nabla
T(\mathbf{x},t)\nabla\psi(\mathbf{x},t)d\mathbf{x}dt+\int_{0}^{t_{f}}\int_{\Gamma_{sf}}hT(\mathbf{x},t)\psi(\mathbf{x},t)d\Gamma
dt=$ $\displaystyle\rho
C_{p}\int_{\Omega}T_{0}(\mathbf{x})\psi(\mathbf{x},0)d\mathbf{x}-\int_{0}^{t_{f}}\int_{\Gamma_{s_{in}}}g(\mathbf{x},t)\psi(\mathbf{x},t)d\Gamma
dt+\int_{0}^{t_{f}}\int_{\Gamma_{sf}}hT_{f}(\mathbf{x})\psi(\mathbf{x},t)d\Gamma
dt,$

for all $\psi\in H^{1}(0,t_{f};H^{1}(\Omega))$ that satisfy
$\psi(\cdot,t_{f})=0$ in $\Omega$.

Thanks to Nittka[4] (its Theorem 2.11 and Corollary 2.13), we have

######  Theorem 2.2.

Let assumptions (H1)-(H7) and (2.1) hold. Then, there exists a unique weak
solution of Problem 2.3 on $[0,t_{f}]$.

Finally, we recall Theorem 3.3 in Nittka[4]

######  Theorem 2.3.

Let assumptions (H1)-(H7) and (2.1) hold. Then, the weak solution $T$ of
Problem 2.3 is in $C([0,t_{f}];C(\overline{\Omega}))$. So, in particular,
$T(\mathbf{x},t)\rightarrow T_{0}(\mathbf{x})$ uniformly on
$\overline{\Omega}$ as $t\rightarrow 0$.

Regarding the numerical solution of Problem 2.3, we use the finite volume
method for its discretization. Given a tessellation $\mathcal{T}$ of the
domain, $\Omega$, we write the discrete unknown $(T_{C}(t))_{C\in\mathcal{T}}$
as the real vector $\mathbf{T}(t)$, belonging to ${\rm I\\!R}^{N_{h}}$ with
$N_{h}=\text{size}(\mathcal{T})$. Then, we write the spatially discretized
problem as

(2.5) $\rho C_{p}M\frac{d\mathbf{T}(t)}{dt}+A\mathbf{T}(t)=\mathbf{b}(t),\quad
t\in(0,t_{f}],$

where $M\in{\rm I\\!M}^{N_{h}\times N_{h}}$ is the mass matrix, $A\in{\rm
I\\!M}^{N_{h}\times N_{h}}$ is the stiffness matrix and $\mathbf{b}\in{\rm
I\\!R}^{N_{h}}$ the source term. The value of each element of $M$, $A$ and
$\mathbf{b}$ depends on the particular finite volume scheme for the
discretization and the mesh used. Since our problem is a classic diffusion
problem, we refer for further details regarding the finite volume
discretization to the Eymard’s monograph.[5]

To discretize (2.5) in time, we divide the time interval of interest into
$N_{T}$ regular steps

(2.6) $t^{0}=0,\quad t^{n+1}=t^{n}+\Delta t,\quad
n=0,\dots,N_{T}-1,\quad\Delta t=\frac{t_{f}}{N_{T}}.$

From now on, we denote by $f^{j}$ an approximation of a given function $f(t)$
at time $t^{j}$.

For the time discretization, we consider the implicit Euler scheme. It is a
first-order implicit scheme. With this discretization, (2.5) becomes

(2.7) $(\rho C_{p}M+\Delta tA)\mathbf{T}^{n+1}=\rho
C_{p}M\mathbf{T}^{n}+\Delta t\mathbf{b}^{n+1},\quad n=0,\dots,N_{T}-1.$

Notice that, thanks to hypotheses (H1)-(H7), matrices $A$ and $M$ are time
independent.

##  3\. Inverse Problem

In this section, we discuss the formulation and solution of the boundary heat
flux estimation problem. We state it in an inverse problem setting using data
assimilation. We begin this section with a literature survey, then we discuss
the mathematical formulation of the problem and, finally, the methodology that
we developed for its solution.

###  3.1. State of the Art

The literature on inverse heat transfer problems is vast[6, 7, 8, 9]. We refer
to Alifanov’s[10], Orlande’s[11], Beck and Clair’s[12], and Chang’s[13] works
for a detailed review. In the literature, other researchers also investigated
the particular problem of computing the mold-slab heat flux from temperature
measurements in the mold[14, 15, 16, 17]. From a mathematical point of view,
the present problem fits in the framework of estimating a Neumann BC (the heat
flux) having as data pointwise measurements of the state (the temperature)
inside the domain. Such problems were also addressed in investigations not
related to heat transfer[18, 19, 20].

The story of inverse heat transfer problems started in the 50s when aerospace
engineers were interested in knowing the thermal properties of heat shields
and heat fluxes on the surface of space vehicles during re-entry. The first
approach was purely heuristic, then in the 60s and 70s, researchers moved to a
more mathematically formal approach. In fact, most of the regularization
theory that we use nowadays to treat ill-posed problems was developed during
these years[10, 21, 22, 23, 24].

The first approach for estimating the boundary heat flux in CC molds was to
select a heat flux profile, and then by trial and error adapt it to match the
measured temperatures[17]. Pinhero et al.[25] were the first to use an optimal
control framework and regularization methods. They used a steady-state version
of the 2D mold model proposed by Samarasekera and Brimacombe[26] and
parameterized the heat flux with a piecewise constant function. Finally, they
used Tikhonov’s regularization for solving the inverse problem and validated
the results with experimental measurements. A similar approach was used more
recently by Rauter et al.[15, 27, 28]. They estimated the heat flux
transferred from the solidifying steel to the mold wall both in a 2D and 3D
domain. They used a steady-state heat conduction model for the mold and
parameterized the heat flux with a piecewise linear profile in 2D and
symmetric cosine profile in 3D. For the solution of the inverse problem, they
used the Conjugate Gradient Method (CGM) and a mixed GA-SIMPLEX algorithm[29]
in 2D while in 3D they only used the GA-SIMPLEX algorithm. Their results were
also tested with experimental data.

Using a 3D unsteady-state heat conduction model in the strand and the mold
with a Robin condition at the mold-strand interface, Hebi et al.[30, 31]
attempted to estimate the solidification in CC round billets. Similarly to the
present work, they looked for the heat transfer coefficient that minimizes a
distance between measured and computed temperatures at the thermocouples’
points. Assuming the heat transfer coefficient to be piecewise constant, they
iteratively adapted each piece to match the measured temperature. However, in
the validation with plant measurements, they did not obtain good agreement. A
similar approach was used by Gonzalez et al.[32] and Wang et al.[33, 34, 35,
36], the latter using a Neumann condition at the mold-strand interface.

Udayraj et al.[16] applied the conjugate gradient method with adjoint problem
for the solution of the steady-state 2D mold-slab heat flux estimation
problem. This methodology was first proposed by Alifanov[10] for the
regularization of boundary inverse heat transfer problems without the need of
parameterizing the heat flux. However, as we also proved in our previous work
[3], this method underestimates the heat flux away from the measurements. To
overcome this issue, Udayraj et al. proposed to average the computed heat flux
at each step and use the uniform averaged value as initial estimation for the
following step. However, the obtained results were not satisfying.

Since the real-time requirement is common in industrial applications, real-
time methodologies for the solution of these problems have already been
investigated in the literature. In particular, Videcoq et al.[37] used a
Branch Eigenmodes Reduced Model[38] for the real-time identification of the
heat source strength variations in a 3D non-linear inverse heat conduction
problem. Later, for solving the same problem, Girault et al.[39] used the
Modal Identification Method[40] for generating the reduced model. Finally,
Aguado et al.[41] coupled classical harmonic analysis with recent model order
reduction techniques (Proper Generalized Decomposition) to solve in real-time
the transient heat equation at monitored points, also showing the
applicability of their method to inverse problems.

To conclude, also deep learning techniques were investigated. Wang and Yao[42]
used the inverse problem solution technique developed by Hebi et al.[31] and a
set of experimental temperature measurements to train a Neural Network (NN)
for on-line computation. Similarly, Chen et al.[43] used the fuzzy inference
method for estimating the mold heat flux. In both works, they modeled the mold
with a 2D steady-state heat conduction model in the solid and parameterized
the boundary heat flux.

Our contribution to the literature is the development of novel methods for
solving the unsteady-state 3D inverse heat transfer problem in CC molds that
exploits the parameterization of the heat flux. We propose different novel
direct methodologies that exploit an offline-online decomposition. In fact, we
divide them in a computationally expensive offline phase and an online phase
whose computational cost is much smaller. The advantage is that we compute
offline phase once and for all before starting the casting process. Then,
while the machine is running, we only need to solve the cheap online phase.
Moreover, in this work, we design some benchmark cases for this application,
and we use them to test the performances of the proposed methodologies.

###  3.2. Inverse Problem Formulation

Before proceeding with the mathematical formulation of the inverse problem, we
do some technical considerations that will guide us in the process. First, the
thermocouples measure the temperature at the sampling frequency $f_{samp}$.
This sampling frequency is typically of $1$ Hz and we will assume this value
all along this investigation (notice that different values of $f_{samp}$ are
compatible with the following discussion). Second, every sampling period
$T_{samp}=1/f_{samp}=1$ s, the thermocouples provide a new set of
measurements, so we have a regular sequence of measurements in time.

That said, we consider the problem of estimating the heat flux, $g$, on
$\Gamma_{s_{in}}$, in between the last acquired measurement instant and the
previous one. In this way, we follow the sequentiality of the measured data in
our solution procedure according to the real-time purpose of this research.

We introduce the following notation. Let
$\Psi:=\\{\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{P}\\}$ be a
collection of points in $\Omega$ and
$\Upsilon:=\\{\tau^{0},\tau^{1},\dots,\tau^{P_{t}}\\}$ a collection of points
in $[0,t_{f}]$ such that $\tau^{k}=t^{kN_{t}}$ (see Figure 3). According to
the introduced sequential approach, we consider the following restriction of
Problem 2.3 to $(\tau^{k-1},\tau^{k}]$, $1\leq k\leq P_{t}$, as direct problem

######  Problem 3.1.

Let $1\leq k\leq P_{t}$ and $g^{k}(\mathbf{x},t)$ be a given heat flux on
$\Gamma_{s_{in}}\times(\tau^{k-1},\tau^{k}]$. Find $T^{k}$ such that

(3.1) $\rho C_{p}\frac{\partial T^{k}}{\partial t}-k_{s}\Delta
T^{k}=0,\quad\text{in }\Omega\times(\tau^{k-1},\tau^{k}],$

with BCs and IC

(3.2) $\begin{cases}-k_{s}\nabla T^{k}\cdot\mathbf{n}=g^{k}&\text{on
}\Gamma_{s_{in}}\times(\tau^{k-1},\tau^{k}],\\\ -k_{s}\nabla
T^{k}\cdot\mathbf{n}=0&\text{on
}\Gamma_{s_{ex}}\times(\tau^{k-1},\tau^{k}],\\\ -k_{s}\nabla
T^{k}\cdot\mathbf{n}=h(T^{k}-T_{f})&\text{on
}\Gamma_{sf}\times(\tau^{k-1},\tau^{k}],\\\
T^{k}(\cdot,\tau^{k-1})=T^{k-1}(\cdot,\tau^{k-1})&\text{in
}\Omega,\end{cases}$

where $T^{0}(\cdot,\tau^{0})=T_{0}$, being $T_{0}$ the initial temperature.

Figure 3. Time line for the inverse problem.

So basically, we are dividing the time domain into chunks going from one
measurement time to the next one in a way that facilitates the definition of
the inverse problems below. Before formulating it, we introduce some further
notation. We define the application
$(\mathbf{x}_{i},\tau^{k})\in\Psi\times\Upsilon\rightarrow\hat{T}(\mathbf{x}_{i},\tau^{k})\in{\rm
I\\!R}^{+}$, $1\leq i\leq P,1\leq k\leq P_{t}$,
$\hat{T}(\mathbf{x}_{i},\tau^{k})$ being the experimentally measured
temperature at $(\mathbf{x}_{i},\tau^{k})\in\Psi\times\Upsilon$. Moreover, to
simplify the notation, and if there is no room for error, we denote

(3.3) $\hat{T}^{k}(\mathbf{x}_{i}):=\hat{T}(\mathbf{x}_{i},\tau^{k}),\quad
1\leq i\leq P,\ 1\leq k\leq P_{t},$

and we let $T^{k}[g]$ represent the solution of Problem 3.2 corresponding to
heat flux $g$ on $\Gamma_{s_{in}}\times(\tau^{k-1},\tau^{k}]$.

At each measurement interval $k$, $1\leq k\leq P_{t}$, we propose an iterative
procedure, assuming that, for $k\geq 1$, $g^{l}$ and $T^{l}[g^{l}]$, $0\leq
l\leq k-1$, have been computed. Using a least square, deterministic approach,
we state two different inverse problems for Problem 3.2. In the first one, we
consider as functional to be minimized a distance between the measured and
computed temperatures at the thermocouples. Then, we state it as

######  Problem 3.2.

(Inverse) Being $g^{l}$ and $T^{l}[g^{l}]$, $1\leq l\leq k-1$, known, and
given the temperature measurements $\hat{T}^{k}(\mathbf{x}_{i})$, $1\leq i\leq
P$, find $g^{k}\in C(\tau^{k-1},\tau^{k};L^{q}(\Gamma_{s_{in}}))$ which
minimizes the functional

(3.4)
$S_{1}^{k}[g^{k}]=\frac{1}{2}\sum_{i=1}^{P}[T^{k}[g^{k}](\mathbf{x}_{i},\tau^{k})-\hat{T}^{k}(\mathbf{x}_{i})]^{2}.$

Here, we denote $T^{0}[g^{0}]=T_{0}$.

The second inverse problem that we consider includes in the cost functional
the $L^{2}$-norm of the heat flux. Thus, we write it as

######  Problem 3.3.

(Inverse) Being $g^{l}$ and $T^{l}[g^{l}]$, $1\leq l\leq k-1$, known, and
given the temperature measurements $\hat{T}^{k}(\mathbf{x}_{i})$, $1\leq i\leq
P$, find $g^{k}\in C(\tau^{k-1},\tau^{k};L^{q}(\Gamma_{s_{in}}))$ which
minimizes the functional

(3.5)
$S_{2}^{k}[g^{k}]=\frac{1}{2}\sum_{i=1}^{P}[T^{k}[g^{k}](\mathbf{x}_{i},\tau^{k})-\hat{T}^{k}(\mathbf{x}_{i})]^{2}+p_{g}\langle
g^{k}(\tau^{k}),g^{k}(\tau^{k})\rangle_{L^{2}(\Gamma_{s_{in}})},$

where $p_{g}\ \left[\text{K\textsuperscript{2} / W\textsuperscript{2}}\right]$
is a weight applied to the heat flux norm.

###  3.3. Inverse Solver for $S_{1}^{k}$

In this section, we discuss a novel methodology for solving Problem 3.4. In
particular, we mimic the methodology developed by the authors for a steady-
state mold model [3], expanding it to the unsteady case.

We exploit a suitable parameterization of the heat flux, $g^{k}$. To properly
parameterize it, we start by considering that we want to parameterize an
unknown function $g^{k}$ in $L^{r}(\tau^{k-1},\tau^{k};L^{q}(\Omega))$, $1\leq
k\leq P_{t}$. Then, we notice that in thin slab casting molds, the
thermocouples are all located few millimeters inward from $\Gamma_{s_{in}}$.
All together they form a uniform 2D grid on a surface parallel to the
$\Gamma_{s_{in}}$ boundary (see Figure 4(b)). Thus, a possible choice for the
parameterization of $g^{k}$ is to use Radial Basis Functions (RBFs) centered
at the projections of the thermocouples points on $\Gamma_{s_{in}}$[44]. By
using this parameterization, we end up having as many basis functions as
thermocouples. Note that the methodology is very well adapted to the
application in use.

In particular, we parameterize $g^{k}$ by Gaussian RBFs which allow us to
separate the time and space dependence. These are continuous functions with
global support in $\Gamma_{s_{in}}$. However, the following discussion can be
applied to other basis functions.

The parameterization of the boundary heat flux reads (see Prando’s
appendix[45])

(3.6)
$g^{k}(\mathbf{x},t)\approx\mathscr{g}^{k}(\mathbf{x},t)=\sum_{i=1}^{P}g_{i}^{k}(t)\phi_{i}(\mathbf{x}),\quad\text{
for }t\in(\tau^{k-1},\tau^{k}],$

where the $\phi_{i}(\mathbf{x})$ are $P$ known basis functions, and the
$g_{i}^{k}(t)$ are the respective time dependent unknown weights.

To define the RBFs, let $\boldsymbol{\xi}_{i},1\leq i\leq P$, be the
projection of the point $\mathbf{x}_{i}\in\Psi$ on $\Gamma_{s_{in}}$, i.e.
such that

(3.7)
$\boldsymbol{\xi}_{i}=\operatorname*{argmin}_{\boldsymbol{\xi}\in\Gamma_{s_{in}}}\left\lVert\mathbf{x}_{i}-\boldsymbol{\xi}\right\rVert_{2},\quad\mathbf{x}_{i}\in\Psi.$

By centering the RBFs in these points, their expression is

(3.8)
$\phi_{j}(\mathbf{x})=e^{-\left(\eta\left\lVert\mathbf{x}-\boldsymbol{\xi}_{j}\right\rVert_{2}\right)^{2}},\quad\text{
for }j=1,2,\dots,P,$

where $\eta$ is the shape parameter of the Gaussian basis. By increasing
(decreasing) its values, the radial decay of the basis slows down (speeds up).

In this work, we explore two different approaches to the time
parameterization. In the first one, we consider $g_{i}^{k}$ independent of
time

(3.9) $g_{i}^{k}(t)=w_{i}^{k},\quad\text{ for }t\in(\tau^{k-1},\tau^{k}],1\leq
i\leq P,$

being $w_{i}^{k}$ real numbers. In this way, the heat flux is assumed to be
piecewise constant, i.e. constant between consecutive measurement instants.

The second approach is to consider the heat flux to be continuous piecewise
linear in $(0,t_{f}]$, being a polynomial of degree 1 between the sampling
times. Then, we assume the weights $g_{i}^{k}(t)$ to be linear in time in the
interval $(\tau^{k-1},\tau^{k}]$. Moreover, in this second case, the following
continuity is assumed

(3.10)
$g_{i}^{k}(t)|_{t\downarrow\tau^{k-1}}=g_{i}^{k-1}(t)|_{t\uparrow\tau^{k-1}}.$

In turn, we characterize $g_{i}^{k}(t)$ as

(3.11)
$g_{i}^{k}(t)=w_{i}^{k-1}+(t-\tau^{k-1})\frac{w^{k}_{i}-w_{i}^{k-1}}{\tau^{k}-\tau^{k-1}},\quad\text{
for }t\in(\tau^{k-1},\tau^{k}].$

Notice that by doing parameterization (3.6), we change the problem from
estimating a function in an infinite dimensional space at each time interval
$(t^{(k-1)N_{t}},t^{kN_{t}}]=(\tau^{k-1},\tau^{k}]$, to estimating the vector
$\mathbf{w}^{k}=(w_{1}^{k},w_{2}^{k},\dots,w_{P}^{k})^{T}$ in ${\rm
I\\!R}^{P}$, for each $1\leq k\leq P_{t}$.

Now, at each time interval $(\tau^{k-1},\tau^{k}]$, the objective of the
inverse problem is to determine $\mathbf{w}^{k}$ which identifies
$\mathscr{g}^{k}$ once the elements of the basis $\phi_{i}$, $i=1,2,\dots,P$
are fixed. We state the inverse problem as

######  Problem 3.4.

(Inverse) Given the temperature measurements $\hat{T}(\Psi,\tau^{k})$, find
$\hat{\mathbf{w}}^{k}\in{\rm I\\!R}^{P}$, $1\leq k\leq P_{t}$, which minimizes
the functional

(3.12)
$S_{1}^{k}[\mathbf{w}^{k}]=\frac{1}{2}\sum_{i=1}^{P}[T^{k}[\mathbf{w}^{k}](\mathbf{x}_{i},\tau^{k})-\hat{T}^{k}(\mathbf{x}_{i})]^{2},$

where if there is not room for confusion $T^{k}[\mathbf{w}^{k}]$ denotes the
temperature $T^{k}[\mathscr{g}^{k}]$, with $\mathscr{g}^{k}$ defined as in
(3.6) and $g_{i}^{k}(t)$ given by (3.9) or (3.11).

For a later use, we define the general vector $\mathbf{a}^{k}\in{\rm
I\\!R}^{P}$ as the vector of the values of a general field $a(\mathbf{x},t)$
at the measurement points and at the measurement time $\tau^{k}$, such as

(3.13) $(\mathbf{a}^{k})_{i}=a(\mathbf{x}_{i},\tau^{k}).$

Moreover, given $\mathbf{w}^{k}$, we define the residual vector
$\mathbf{R}^{k}[\mathbf{w}^{k}]\in{\rm I\\!R}^{P}$ as

(3.14)
$(\mathbf{R}^{k}[\mathbf{w}^{k}])_{i}:=(\mathbf{T}^{k}[\mathbf{w}^{k}])_{i}-(\hat{\mathbf{T}}^{k})_{i},\quad
i=1,2,\dots,P.$

Thanks to (3.14), we rewrite the cost functional (3.12) as

(3.15)
$S_{1}^{k}[\mathbf{w}^{k}]=\frac{1}{2}\mathbf{R}^{k}[\mathbf{w}^{k}]^{T}\mathbf{R}^{k}[\mathbf{w}^{k}].$

To minimize it, we write the critical point equation

(3.16) $\frac{\partial S_{1}^{k}[\hat{\mathbf{w}}^{k}]}{\partial
w_{j}^{k}}=\sum_{i=1}^{P}(R^{k}[\hat{\mathbf{w}}^{k}])_{i}\frac{\partial(\mathbf{T}^{k}[\hat{\mathbf{w}}^{k}])_{i}}{\partial
w_{j}^{k}}=0,\text{ for }j=1,2,\dots,P.$

Thus, for each $k$, $1\leq k\leq P_{t}$, the solution of this equation will
provide the weights vector $\hat{\mathbf{w}}^{k}$ corresponding to a critical
point of $S_{1}^{k}$.

To explicitly obtain from (3.16) an equation for the weights that minimize our
functional $S^{k}_{1}$, we exploit the linearity of Problem 3.2. To derive it,
we consider separately the piecewise constant (3.9) and the piecewise linear
(3.11) cases.

####  3.3.1. Piecewise Constant Approximation of the Heat Flux

Suppose to have the solutions to the following auxiliary problems

######  Problem 3.5.

For each $i$, $1\leq i\leq P$, find $T_{\phi_{i}}$ such that

(3.17) $\rho C_{p}\frac{\partial T_{\phi_{i}}}{\partial t}-k_{s}\Delta
T_{\phi_{i}}=0,\quad\text{in }\Omega\times(\tau^{0},\tau^{1}],$

with BCs and IC

(3.18) $\begin{cases}-k_{s}\nabla
T_{\phi_{i}}\cdot\mathbf{n}=\phi_{i}&\text{on
}\Gamma_{s_{in}}\times(\tau^{0},\tau^{1}],\\\ -k_{s}\nabla
T_{\phi_{i}}\cdot\mathbf{n}=0&\text{on
}\Gamma_{s_{ex}}\times(\tau^{0},\tau^{1}],\\\ -k_{s}\nabla
T_{\phi_{i}}\cdot\mathbf{n}=hT_{\phi_{i}}&\text{on
}\Gamma_{sf}\times(\tau^{0},\tau^{1}],\\\
T_{\phi_{i}}(\cdot,\tau^{0})=0&\text{in }\Omega.\end{cases}$

######  Problem 3.6.

For each $k$, $1\leq k\leq P_{t}$, find $T_{IC}^{k}$ such that

(3.19) $\rho C_{p}\frac{\partial T_{IC}^{k}}{\partial t}-k_{s}\Delta
T_{IC}^{k}=0,\quad\text{in }\Omega\times(\tau^{k-1},\tau^{k}],$

with BCs and IC

(3.20) $\begin{cases}-k_{s}\nabla T_{IC}^{k}\cdot\mathbf{n}=0&\text{on
}(\Gamma_{s_{in}}\cup\Gamma_{s_{ex}})\times(\tau^{k-1},\tau^{k}],\\\
-k_{s}\nabla T_{IC}^{k}\cdot\mathbf{n}=h\left(T_{IC}^{k}-T_{f}\right)&\text{on
}\Gamma_{sf}\times(\tau^{k-1},\tau^{k}],\\\
T_{IC}^{k}(\cdot,\tau^{k-1})=T^{k-1}(\cdot,\tau^{k-1})&\text{in
}\Omega,\end{cases}$

with

(3.21) $T_{IC}^{1}(\cdot,\tau^{0})=T_{0}.$

Notice that Problem 3.18 does not depend on the measurement instants index
$k$. Then, we define it only in the first interval $(\tau^{0},\tau^{1}]$ and,
if needed, translate it such as

(3.22)
$T_{\phi_{i}}^{k}(\mathbf{x},t)=T_{\phi_{i}}(\mathbf{x},t-\tau^{k-1}),\text{
for }t\in(\tau^{k-1},\tau^{k}].$

We can now state

######  Theorem 3.1.

Given $T_{IC}^{k}$ and $\mathbf{w}^{k},1\leq k\leq P_{t}$ (and so
$\mathscr{g}^{k}$ defined as (3.6)) and $T_{\phi_{i}},1\leq i\leq P$, the
function defined as

(3.23)
$T^{k}[\mathscr{g}^{k}]=\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k},$

is the solution to Problem 3.2 associated with the heat flux
$\mathscr{g}(x,t)$ which in each measurement subinterval
$(\tau^{k-1},\tau^{k}]$ coincides with $\mathscr{g}^{k}(x,t)$ given by (3.6)
with $g_{i}^{k}(t)$ as in (3.9).

###### Proof.

Substituting (3.23) in (3.1), taking into account (3.9) so $g_{i}^{k}$ is
constant in $(\tau^{k-1},\tau^{k}]$, and considering (3.17) and (3.19), we
have

(3.24) $\displaystyle\rho
C_{p}\frac{\partial\left(\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k}\right)}{\partial
t}-k_{s}\Delta\left(\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k}\right)=\sum_{i=1}^{P}g_{i}^{k}\left(\rho
C_{p}\frac{\partial T^{k}_{\phi_{i}}}{\partial t}-k_{s}\Delta
T^{k}_{\phi_{i}}\right)+$ $\displaystyle\rho C_{p}\frac{\partial
T_{IC}^{k}}{\partial t}-k_{s}\Delta T_{IC}^{k}=0\text{ in
}\Omega\times(\tau^{k-1},\tau^{k}].$

Similarly, for the BCs we have

(3.25)
$-k_{s}\nabla\left(\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k}\right)\cdot\mathbf{n}=\sum_{i=1}^{P}g_{i}^{k}\left(-k_{s}\nabla
T^{k}_{\phi_{i}}\cdot\mathbf{n}\right)=\sum_{i=1}^{P}g_{i}^{k}\phi_{i}=\mathscr{g}^{k},\text{
on }\Gamma_{s_{in}}\times(\tau^{k-1},\tau^{k}],$ (3.26)
$-k_{s}\nabla\left(\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k}\right)\cdot\mathbf{n}=0,\text{
on }\Gamma_{s_{ex}}\times(\tau^{k-1},\tau^{k}],$

and

(3.27)
$-k_{s}\nabla\left(\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k}\right)\cdot\mathbf{n}=h\left[\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}+T_{IC}^{k}-T_{f}\right]=h\left(T^{k}[\mathscr{g}^{k}]-T_{f}\right),\text{
on }\Gamma_{sf}\times(\tau^{k-1},\tau^{k}].$

With respect to the IC, at each interval we must proceed by induction. For
$k=1$, thanks to (3.21),

(3.28) $T^{1}[g^{1}](\cdot,\tau^{0})=T^{0}(\cdot).$

For $k>1$, taking into account (3.18) and (3.20)

(3.29)
$T^{k}[g^{k}](\cdot,\tau^{k-1})=\sum_{i=1}^{P}g_{i}^{k}T^{k}_{\phi_{i}}(\cdot,\tau^{k-1})+T_{IC}^{k}(\cdot,\tau^{k-1})=T^{k-1}(\cdot,\tau^{k-1}),\text{
in }\Omega.$

This ends the proof. ∎

##### Solving the minimization problem for $S_{1}^{k}$ with piecewise constant
approximation of the heat flux

Thanks to (3.22) and (3.23), (3.16) can be written as

(3.30)
$\displaystyle\sum_{i=1}^{P}(R^{k}[\hat{\mathbf{w}}^{k}])_{i}\frac{\partial\left(\sum_{l=1}^{P}g_{l}^{k}\mathbf{T}_{\phi_{l}}+\mathbf{T}_{IC}^{k}\right)_{i}}{\partial
w_{j}^{k}}=\sum_{i=1}^{P}(R^{k}[\hat{\mathbf{w}}^{k}])_{i}\frac{\partial\left(\sum_{l=1}^{P}w_{l}^{k}\mathbf{T}_{\phi_{l}}+\mathbf{T}_{IC}^{k}\right)_{i}}{\partial
w_{j}^{k}}=$
$\displaystyle\mathbf{R}^{k}[\hat{\mathbf{w}}^{k}]^{T}\left(\mathbf{T}_{\phi_{j}}\right)=$
$\displaystyle\ 0,\quad\text{ for }j=1,2,\dots,P,$

where $(\mathbf{T}_{\phi_{l}})_{i}$ is the vector containing the values of the
field $T_{\phi_{l}}(\mathbf{x}_{i},\tau_{1})$ at the measurement points. We
recall that $\mathbf{T}_{\phi_{l}}$ is independent of $k$.

Let us define the matrix $\Theta$ in ${\rm I\\!M}^{P\times P}$ such that

(3.31) $\Theta_{i,j}=T_{\phi_{j}}(\mathbf{x}_{i},\tau^{1}).$

Equation (3.30) can now be written as

(3.32) $\Theta^{T}\mathbf{R}^{k}[\hat{\mathbf{w}}^{k}]=\mathbf{0}.$

Using (3.23), the vector associated to the solution of the direct problem at
the measurement points, $\mathbf{T}^{k}[\mathbf{w}^{k}]\in{\rm I\\!R}^{P}$,
for each $k$, can be written as

(3.33)
$\mathbf{T}^{k}[\hat{\mathbf{w}}^{k}]=\sum_{j=1}^{P}w^{k}_{j}\mathbf{T}_{\phi_{j}}+\mathbf{T}^{k}_{IC}=\Theta\hat{\mathbf{w}}^{k}+\mathbf{T}^{k}_{IC}.$

Recalling the definition of $\mathbf{R}^{k}$ and (3.33), we have

(3.34)
$\Theta^{T}\mathbf{R}^{k}[\hat{\mathbf{w}}^{k}]=\Theta^{T}(\Theta\hat{\mathbf{w}}^{k}+\mathbf{T}^{k}_{IC}-\hat{\mathbf{T}}^{k})=\mathbf{0}.$

Therefore, for each $k$, $1\leq k\leq P_{t}$, a solution of the inverse
problem, $\hat{\mathbf{w}}^{k}$, is obtained by solving the linear system

(3.35)
$\Theta^{T}\Theta\hat{\mathbf{w}}^{k}=\Theta^{T}(\hat{\mathbf{T}}^{k}-\mathbf{T}^{k}_{IC}).$

It is important to notice that the system matrix, $\Theta^{T}\Theta$, is
$k$-independent.

Equation (3.35) is generally called the normal equation. By solving this
linear system, we obtain the weights, $\hat{\mathbf{w}}^{k}$, that correspond
to a critical point of the functional $S_{1}^{k}$, defined by (3.12). As
mentioned, $\Theta$ is constant. So, we can compute it once and for all in an
offline phase.

The proposed methodology for the solution of the inverse Problem 3.4 is
summarized in Algorithm 1. It is important to notice that, for each time
interval $(\tau^{k-1},\tau^{k}]$, $T_{\phi_{i}}$ and the related vector do not
change but $T_{IC}^{k}$ does because its IC depends on the temperature field
at time $\tau^{k-1}$.

Algorithm 1 Inverse solver for the solution of Problem 3.4 with piecewise
constant parameterization in time of the heat flux, $g$.

OFFLINE
Input RBF shape parameter, $\eta$; thermocouples measurement points and times,
$\Psi,\Upsilon$

1:Setup RBF parameterization by (3.8)

2:Compute $T_{\phi_{i}}$ for $i=1,2,\dots,P$ by solving Problem 3.18

3:Assemble matrix $\Theta$ by (3.31)

ONLINE
Input Initial condition, $T_{0}$

1:Set $k=1$

2:while $k\leq P_{t}$ do

3: Read the thermocouples measurements, $\hat{\mathbf{T}}^{k}$

4: Compute $T_{IC}^{k}$ by solving Problem 3.20

5: Assemble $\mathbf{T}_{IC}^{k}$

6: Compute $\hat{\mathbf{w}}^{k}$ by solving (3.35)

7: Compute $g_{i}^{k}(t)$ by (3.9)

8: Compute the heat flux $\mathscr{g}(\mathbf{x},t)$ for
$t\in(\tau^{k-1},\tau^{k}]$ by (3.6)

9: Use (3.23) to compute $T^{k}[g^{k}]$

10: $k=k+1$

11:end while

Notice that, in this setting, (3.35) is an affine map from the observations,
$\hat{\mathbf{T}}^{k}$, to the heat flux weights, $\mathbf{w}^{k}$.
Consequently, we have that the existence and uniqueness of the solution of the
inverse problem depends on the invertibility of the matrix $\Theta^{T}\Theta$.
We can easily see that the matrix is symmetric and positive semi-definite. In
general, however, we cannot ensure that it is invertible. In fact, the
invertibility depends on the choice of the basis function, the computational
domain, and the BCs.

Before moving to the piecewise linear case, we recall the offline-online
decomposition of Algorithm 1. In the offline phase, we compute $T_{\phi_{i}}$
for $i=1,2,\dots,P$ by solving Problem 3.18 and assemble the related matrix
$\Theta$. Then, in the online phase, we input the measurements
$\hat{\mathbf{T}}$, solve Problem 3.20 and the linear system (3.35).

For the choice made when selecting the basis functions, this linear system has
the dimensions of the number of thermocouples (quite small, in general).
However, solution of Problem 3.20 involves the solution of a full order model
whose computational cost depends on the discretization size. Consequently,
this method is not suitable for real-time as it is. To achieve real-time
performances, we need to apply model order reduction techniques. This will be
the subject of our future work.

As a final remark, we notice, that for the application of this method,
linearity of the direct problem is essential. In fact, it is a necessary
condition for Theorem 3.1.

####  3.3.2. Piecewise Linear Heat Flux

Suppose to have the solution to

######  Problem 3.7.

For each $i$, $1\leq i\leq P$, find $T_{d_{i}}$ such that

(3.36) $\rho C_{p}\frac{\partial T_{d_{i}}}{\partial t}-k_{s}\Delta
T_{d_{i}}=-\rho C_{p}T_{\phi_{i}},\quad\text{in
}\Omega\times(\tau^{0},\tau^{1}],$

with BCs and IC

(3.37) $\begin{cases}-k_{s}\nabla T_{d_{i}}\cdot\mathbf{n}=0&\text{on
}(\Gamma_{s_{in}}\cup\Gamma_{s_{ex}})\times(\tau^{0},\tau^{1}],\\\
-k_{s}\nabla T_{d_{i}}\cdot\mathbf{n}=hT_{d_{i}}&\text{on
}\Gamma_{sf}\times(\tau^{0},\tau^{1}],\\\ T_{d_{i}}(\cdot,\tau^{0})=0&\text{in
}\Omega.\end{cases}$

From the solution of this problem, we define

(3.38)
$T_{d_{i}}^{k}(\mathbf{x},t)=T_{d_{i}}(\mathbf{x},t-\tau^{k-1}),\quad\text{for
}t\in(\tau^{k-1},\tau^{k}].$

Then, we state

######  Theorem 3.2.

Given $T_{IC}^{k}$ and $\mathbf{w}^{k},1\leq k\leq P_{t}$ (and so
$\mathscr{g}^{k}$ defined as (3.6)), $T_{d_{i}}$ and $T_{\phi_{i}}$, $1\leq
i\leq P$, then the function defined as

(3.39)
$T^{k}[\mathscr{g}^{k}]=\sum_{i=1}^{P}\left(g_{i}^{k}T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k},$

is the solution to Problem 3.2 associated with the heat flux
$\mathscr{g}(x,t)$ which in each measurement subinterval
$(\tau^{k-1},\tau^{k}]$ coincides with $\mathscr{g}^{k}(x,t)$ given by (3.6),
with $g_{i}^{k}(t)$ as in (3.11).

###### Proof.

Substituting (3.39) in (3.1) and considering (3.17) and (3.19), we have

(3.40) $\displaystyle\rho
C_{p}\frac{\partial\left[\sum_{i=1}^{P}\left(g_{i}^{k}(t)T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k}\right]}{\partial
t}-k_{s}\Delta\left[\sum_{i=1}^{P}\left(g_{i}^{k}(t)T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k}\right]=$
$\displaystyle\rho
C_{p}\left[\sum_{i=1}^{P}{g_{i}^{k}}^{\prime}T^{k}_{\phi_{i}}+\sum_{i=1}^{P}g_{i}^{k}(t)\left(\frac{\partial
T^{k}_{\phi_{i}}}{\partial t}-k_{s}\Delta
T^{k}_{\phi_{i}}\right)\right]+\sum_{i=1}^{P}{g_{i}^{k}}^{\prime}\left(\rho
C_{p}\frac{\partial T^{k}_{d_{i}}}{\partial t}-k_{s}\Delta
T^{k}_{d_{i}}\right)+\rho C_{p}\frac{\partial T_{IC}^{k}}{\partial
t}-k_{s}\Delta T_{IC}^{k}=$ $\displaystyle\rho
C_{p}\sum_{i=1}^{P}{g_{i}^{k}}^{\prime}T^{k}_{\phi_{i}}+\sum_{i=1}^{P}{g_{i}^{k}}^{\prime}\left(\rho
C_{p}\frac{\partial T^{k}_{d_{i}}}{\partial t}-k_{s}\Delta
T^{k}_{d_{i}}\right)=0,\text{ in }\Omega\times(\tau^{k-1},\tau^{k}].$

Similarly, for the BCs we have

(3.41) $\displaystyle-
k_{s}\nabla\left(\sum_{i=1}^{P}\left(g_{i}^{k}T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k}\right)\cdot\mathbf{n}=\sum_{i=1}^{P}g_{i}^{k}\left(-k_{s}\nabla
T^{k}_{\phi_{i}}\cdot\mathbf{n}\right)=\sum_{i=1}^{P}g_{i}^{k}\phi_{i}=\mathscr{g}^{k},\text{
on }\Gamma_{s_{in}}\times(\tau^{k-1},\tau^{k}],$ (3.42)
$-k_{s}\nabla\left(\sum_{i=1}^{P}\left(g_{i}^{k}T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k}\right)\cdot\mathbf{n}=0\text{
on }\Gamma_{s_{ex}}\times(\tau^{k-1},\tau^{k}],$

and

(3.43) $\displaystyle-
k_{s}\nabla\left(\sum_{i=1}^{P}\left(g_{i}^{k}T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k}\right)\cdot\mathbf{n}=h\left[\sum_{i=1}^{P}\left(g_{i}^{k}T^{k}_{\phi_{i}}+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}\right)+T_{IC}^{k}-T_{f}\right]=h\left(T^{k}[\mathscr{g}^{k}]-T_{f}\right),$
$\displaystyle\text{ on }\Gamma_{sf}\times(\tau^{k-1},\tau^{k}].$

With respect to the IC, at each interval we must proceed by induction. For
$k=1$,

(3.44) $T^{1}[g^{1}](\cdot,\tau^{0})=T^{0}(\cdot)\text{ in }\Omega.$

For $k>1$, thanks to (3.18), (3.20) and (3.37), we deduce,

(3.45)
$T^{k}[g^{k}](\cdot,\tau^{k-1})=\sum_{i=1}^{P}\left(g_{i}^{k}(\tau^{k-1})T^{k}_{\phi_{i}}(\cdot,\tau^{k-1})+{g_{i}^{k}}^{\prime}T^{k}_{d_{i}}(\cdot,\tau^{k-1})\right)+T_{IC}^{k}(\cdot,\tau^{k-1})=T^{k-1}(\cdot,\tau^{k-1}),\text{
in }\Omega.$

∎

##### Solving the minimization problem for $S_{1}^{k}$ with piecewise linear
approximation of the heat flux

Thanks to Theorem 3.2, we have

(3.46) $\displaystyle\frac{\partial\mathbf{T}^{k}[\mathbf{w}^{k}]}{\partial
w_{j}^{k}}=\frac{\partial\left[\sum_{l=1}^{P}\left(g_{l}^{k}(\tau^{k})\mathbf{T}^{k}_{\phi_{l}}(\tau^{k})+{g_{l}^{k}}^{\prime}\mathbf{T}^{k}_{d_{l}}(\tau^{k})\right)+\mathbf{T}_{IC}^{k}(\tau^{k})\right]}{\partial
w_{j}^{k}}=$
$\displaystyle\frac{\partial\left\\{\sum_{l=1}^{P}\left[w_{l}^{k}\mathbf{T}^{k}_{\phi_{l}}(\tau^{k})+(w_{l}^{k}-w_{l}^{k-1})f_{samp}\mathbf{T}^{k}_{d_{l}}(\tau^{k})\right]+\mathbf{T}_{IC}^{k}(\tau^{k})\right\\}}{\partial
w_{j}^{k}}=\mathbf{T}^{k}_{\phi_{j}}(\tau^{k})+f_{samp}\mathbf{T}^{k}_{d_{j}}(\tau^{k})=\mathbf{T}_{\phi_{j}}(\tau^{1})+f_{samp}\mathbf{T}_{d_{j}}(\tau^{1}).$

Thanks to (3.46), (3.16) in the linear case is rewritten as

(3.47)
$\mathbf{R}^{k}[\hat{\mathbf{w}}^{k}]^{T}(\mathbf{T}_{\phi_{j}}+f_{samp}\mathbf{T}_{d_{j}})=0,\text{
for all }j=1,2,\dots,P.$

Let us define the matrices $\tilde{\Theta},\Theta_{d}$ in ${\rm
I\\!M}^{P\times P}$ such that

(3.48)
$\left(\Theta_{d}\right)_{i,j}:=f_{samp}T_{d_{j}}(\mathbf{x}_{i},\tau^{1}),\quad\left(\tilde{\Theta}\right)_{i,j}:=\left(\Theta\right)_{i,j}+\left(\Theta_{d}\right)_{i,j},$

both independent of the index $k$.

Using (3.48) and (3.39), we have

(3.49)
$\displaystyle\mathbf{T}^{k}[\hat{\mathbf{w}}^{k}]=\sum_{l=1}^{P}\left(g_{l}^{k}(\tau^{k})\mathbf{T}_{\phi_{l}}+{g_{l}^{k}}^{\prime}\mathbf{T}_{d_{l}}\right)+\mathbf{T}_{IC}^{k}=\sum_{l=1}^{P}\left[w_{l}^{k}\mathbf{T}_{\phi_{l}}+(w_{l}^{k}-w_{l}^{k-1})f_{samp}\mathbf{T}_{d_{l}}\right]+\mathbf{T}_{IC}^{k}=\tilde{\Theta}\hat{\mathbf{w}}^{k}-\Theta_{d}\hat{\mathbf{w}}^{k-1}+\mathbf{T}^{k}_{IC}.$

Recalling the definition of $\mathbf{R}^{k}$, taking into account (3.48) and
(3.49), system (3.47) can be written as

(3.50)
${(\tilde{\Theta})}^{T}\mathbf{R}^{k}[\hat{\mathbf{w}}^{k}]={(\tilde{\Theta})}^{T}\left(\tilde{\Theta}\hat{\mathbf{w}}^{k}-\Theta_{d}\hat{\mathbf{w}}^{k-1}+\mathbf{T}^{k}_{IC}-\hat{\mathbf{T}}^{k}\right)=\mathbf{0}.$

Therefore, a solution of the inverse problem $S_{1}^{k}$, considering
$g^{k}(t)$ piecewise linear, is obtained by solving the linear system

(3.51)
${(\tilde{\Theta})}^{T}\tilde{\Theta}\hat{\mathbf{w}}^{k}={(\tilde{\Theta})}^{T}\left(\hat{\mathbf{T}}^{k}+\Theta_{d}\hat{\mathbf{w}}^{k-1}-\mathbf{T}^{k}_{IC}\right).$

Notice that as in the piecewise constant case, the system matrix,
$(\tilde{\Theta})^{T}\tilde{\Theta}$, is $k$-independent.

We summarize the proposed methodology for the solution of the inverse Problem
3.4 in Algorithm 2. Similarly to the piecewise constant case, for each time
interval $(\tau^{k-1},\tau^{k}]$, $T_{\phi_{i}}$, $T_{d_{i}}$ and the related
vectors do not change but $T_{IC}^{k}$ does because its IC depends on the
temperature field at time $\tau^{k-1}$.

Algorithm 2 Inverse solver for the solution of Problem 3.4 with piecewise
linear parameterization in time of the heat flux, $g$.

OFFLINE
Input RBF shape parameter, $\eta$; thermocouples measurement points and times,
$\Psi,\Upsilon$

1:Setup RBF parameterization by (3.8)

2:Compute $T_{\phi_{i}}$ for $i=1,2,\dots,P$ by solving Problem 3.18

3:Compute $T_{d_{i}}$ for $i=1,2,\dots,P$ by solving Problem 3.37

4:Assemble matrices $\tilde{\Theta}$ and $\Theta_{d}$

ONLINE
Input Initial condition, $T_{0}$

1:Set $k=1$

2:while $k\leq P_{t}$ do

3: Read the thermocouples measurements, $\hat{\mathbf{T}}^{k}$

4: Compute $T_{IC}^{k}$ by solving Problem 3.20

5: Assemble $\mathbf{T}_{IC}^{k}$

6: Compute $\hat{\mathbf{w}}^{k}$ by solving (3.51)

7: Compute $g_{i}^{k}(t)$ by (3.11)

8: Compute the heat flux $\mathscr{g}(\mathbf{x},t)$ for
$t\in(\tau^{k-1},\tau^{k}]$ by (3.6)

9: Use (3.49) to compute $T^{k}[g^{k}]$

10: $k=k+1$

11:end while

Also in this setting, (3.51) is an affine map from the observations,
$\hat{\mathbf{T}}^{k}$, to the heat flux weights, $\hat{\mathbf{w}}^{k}$.
Consequently, we have that the existence and uniqueness of the solution of the
inverse problem depends on the invertibility of the matrix
${(\tilde{\Theta})}^{T}\tilde{\Theta}$. It is symmetric and positive semi-
definite. However, we cannot ensure that it is invertible. In fact, the
invertibility depends on the choice of the basis functions, the computational
domain, and the BCs.

We notice that, also in the piecewise linear case, the offline-online
decomposition holds. Moreover, the linear system dimensions are the same of
the piecewise constant case and the linearity of the direct problem is still a
necessary condition also for Theorem 3.2.

###  3.4. Inverse Solver for $S_{2}^{k}$

In this section, we discuss the solution of the inverse Problem 3.5. In
particular, we extend the previously developed methodologies adapting them to
the cost function $S_{2}^{k}$ as defined in (3.5).

As shown in detail in Section 4.1, the piecewise linear inverse solver of
Section 3.3.2 presents instability issues, under certain conditions. Then, we
stated this second inverse problem with the purpose of stabilizing the
solution. To do this, we designed Problem 3.5 from Problem 3.4 by adding to
the cost function (3.5) a term that penalizes the heat flux norm $\langle
g^{k}(\tau^{k}),g^{k}(\tau^{k})\rangle_{L^{2}(\Gamma_{s_{in}})}$.

Also in this case, we exploit the parameterization of the heat flux (3.6). As
a consequence, we introduce the inverse problem in terms of $\mathbf{w}^{k}$
as

######  Problem 3.8.

(Inverse) Given the temperature measurements $\hat{T}(\Psi,\Upsilon)$, find
$\hat{\mathbf{w}}^{k}\in{\rm I\\!R}^{P}$, $1\leq k\leq P_{t}$, which minimizes
the functional

(3.52)
$S_{2}^{k}[\mathbf{w}^{k}]=\frac{1}{2}\sum_{i=1}^{P}[T^{k}[\mathbf{w}^{k}](\mathbf{x}_{i},\tau^{k})-\hat{T}^{k}(\mathbf{x}_{i})]^{2}+p_{g}\int_{\Gamma_{s_{in}}}\left(\sum_{j=1}^{P}w_{j}^{k}\phi_{j}(\mathbf{x})\right)\left(\sum_{q=1}^{P}w_{q}^{k}\phi_{q}(\mathbf{x})\right)d\Gamma.$

Notice that (3.52) holds true for both the piecewise constant and linear
parameterization of the heat flux since

(3.53) $g_{i}^{k}(\tau^{k})=w_{i}^{k},$

for (3.9) as well as for (3.11).

Considering the second term of the right hand side of (3.52), we can write

(3.54)
$\int_{\Gamma_{s_{in}}}\left(\sum_{l=1}^{P}w_{l}^{k}\phi_{l}(\mathbf{x})\right)\left(\sum_{q=1}^{P}w_{q}^{k}\phi_{q}(\mathbf{x})\right)d\Gamma=\sum_{l=1}^{P}\sum_{q=1}^{P}w_{l}^{k}w_{q}^{k}\int_{\Gamma_{s_{in}}}\phi_{l}(\mathbf{x})\phi_{q}(\mathbf{x})d\Gamma.$

Let us define the vectors of ${\rm I\\!R}^{P^{2}}$

(3.55)
$\boldsymbol{\phi}_{\phi}=\begin{bmatrix}\int_{\Gamma_{s_{in}}}\phi_{1}(\mathbf{x})\phi_{1}(\mathbf{x})d\Gamma&\int_{\Gamma_{s_{in}}}\phi_{1}(\mathbf{x})\phi_{2}(\mathbf{x})d\Gamma&\cdots&\int_{\Gamma_{s_{in}}}\phi_{P}(\mathbf{x})\phi_{P}(\mathbf{x})d\Gamma\\\
\end{bmatrix}^{T},$

and

(3.56)
$\mathbf{a}_{w}^{k}=\begin{bmatrix}{w_{1}^{k}}^{2}&w_{1}^{k}w_{2}^{k}&\cdots&{w_{P}^{k}}^{2}\end{bmatrix}^{T}.$

We can now rewrite (3.54) as

(3.57)
$\sum_{l=1}^{P}\sum_{q=1}^{P}w_{l}^{k}w_{q}^{k}\int_{\Gamma_{s_{in}}}\phi_{l}(\mathbf{x})\phi_{q}(\mathbf{x})d\Gamma=\boldsymbol{\phi}_{\phi}^{T}\mathbf{a}_{w}^{k}.$

Furthermore, deriving (3.57) with respect to the weights, we obtain

(3.58)
$\frac{\partial\left[\int_{\Gamma_{s_{in}}}\left(\sum_{l=1}^{P}w_{l}^{k}\phi_{l}(\mathbf{x})\right)\left(\sum_{q=1}^{P}w_{q}^{k}\phi_{q}(\mathbf{x})\right)d\Gamma\right]}{\partial
w_{j}^{k}}=\frac{\partial\boldsymbol{\phi}_{\phi}^{T}\mathbf{a}_{w}^{k}}{\partial
w_{j}^{k}}=\boldsymbol{\phi}_{\phi}^{T}\frac{\partial\mathbf{a}_{w}^{k}}{\partial
w_{j}^{k}},\quad\text{for }j=1,2,\dots,P.$

Considering the case $j=1$, we have

(3.59) $\boldsymbol{\phi}_{\phi}^{T}\frac{\partial\mathbf{a}_{w}^{k}}{\partial
w_{1}^{k}}=\boldsymbol{\phi}_{\phi}^{T}\setcounter{MaxMatrixCols}{15}\begin{bmatrix}2w_{1}^{k}&w_{2}^{k}&\cdots&w_{P}^{k}&w_{2}^{k}&0&\cdots&0&w_{3}^{k}&0&\cdots&0&w_{P}^{k}&0&\cdots\end{bmatrix}^{T},$

that we can rewrite as

(3.60) $\boldsymbol{\phi}_{\phi}^{T}\frac{\partial\mathbf{a}_{w}^{k}}{\partial
w_{1}^{k}}=2\boldsymbol{\phi}_{\phi_{1}}w_{1}^{k}+\boldsymbol{\phi}_{\phi_{2}}w_{2}^{k}+\boldsymbol{\phi}_{\phi_{3}}w_{3}^{k}+\cdots+\boldsymbol{\phi}_{\phi_{P}}w_{P}^{k}+\boldsymbol{\phi}_{\phi_{P+1}}w_{2}^{k}+\boldsymbol{\phi}_{\phi_{2P+1}}w_{3}^{k}+\cdots+\boldsymbol{\phi}_{\phi_{(P-1)P+1}}w_{P}^{k}.$

Now, by noticing that

(3.61)
$\boldsymbol{\phi}_{\phi_{(r-1)P+s}}=\int_{\Gamma_{s_{in}}}\phi_{r}(\mathbf{x})\phi_{s}(\mathbf{x})d\Gamma=\boldsymbol{\phi}_{\phi_{(s-1)P+r}},\quad
1\leq r,s\leq P,$

we obtain

(3.62) $\boldsymbol{\phi}_{\phi}^{T}\frac{\partial\mathbf{a}_{w}^{k}}{\partial
w_{1}^{k}}=2\boldsymbol{\phi}_{\phi_{1}}w_{1}^{k}+2\boldsymbol{\phi}_{\phi_{2}}w_{2}^{k}+2\boldsymbol{\phi}_{\phi_{3}}w_{3}^{k}+\cdots+2\boldsymbol{\phi}_{\phi_{P}}w_{P}^{k}=2\boldsymbol{\phi}_{\phi_{1:P}}^{T}\mathbf{w}^{k}.$

Similarly, if we consider the general case, we have

(3.63)
$\displaystyle\boldsymbol{\phi}_{\phi}^{T}\frac{\partial\mathbf{a}_{w}^{k}}{\partial
w_{j}^{k}}=2\boldsymbol{\phi}_{\phi_{(j-1)P+1}}w_{1}^{k}+2\boldsymbol{\phi}_{\phi_{(j-1)P+2}}w_{2}^{k}+2\boldsymbol{\phi}_{\phi_{(j-1)P+3}}w_{3}^{k}+\cdots+2\boldsymbol{\phi}_{\phi_{(j-1)P+P}}w_{P}^{k}=2\boldsymbol{\phi}_{\phi_{(j-1)P+1:(j-1)P+P}}^{T}\mathbf{w}^{k},\quad\text{for
}j=1,2,\dots,P.$

Therefore, thanks to (3.58) and (3.63), we can write

(3.64)
$\frac{\partial\left[\int_{\Gamma_{s_{in}}}\left(\sum_{l=1}^{P}w_{l}^{k}\phi_{l}(\mathbf{x})\right)\left(\sum_{q=1}^{P}w_{q}^{k}\phi_{q}(\mathbf{x})\right)d\Gamma\right]}{\partial
w_{j}^{k}}=\boldsymbol{\phi}_{\phi}^{T}\frac{\partial\mathbf{a}_{w}^{k}}{\partial
w_{j}^{k}}=2\boldsymbol{\phi}_{\phi_{(j-1)P+1:(j-1)P+P}}^{T}\mathbf{w}^{k},\quad\text{for
}j=1,2,\dots,P.$

Let us define the matrix $\Phi\in{\rm I\\!M}^{P\times P}$ such that

(3.65)
$\Phi_{r,s}:=\int_{\Gamma_{s_{in}}}\phi_{r}(\mathbf{x})\phi_{s}(\mathbf{x})d\Gamma.$

If we now consider the minimization of $S_{2}^{k}$ with respect to the
weights, $\mathbf{w}^{k}$, as in (3.16), we have

(3.66) $\frac{\partial S_{2}^{k}[\hat{\mathbf{w}}^{k}]}{\partial
w_{j}^{k}}=\sum_{i=1}^{P}(R^{k}[\hat{\mathbf{w}}^{k}])_{i}\frac{\partial(\mathbf{T}^{k}[\hat{\mathbf{w}}^{k}])_{i}}{\partial
w_{j}^{k}}+2p_{g}\boldsymbol{\phi}_{\phi_{(j-1)P+1:(j-1)P+P}}^{T}\mathbf{w}^{k}=0,\text{
for }j=1,2,\dots,P.$

Considering the piecewise constant case, thanks to (3.33) and (3.58), we
rewrite (3.66) as

(3.67)
$\Theta^{T}(\Theta\hat{\mathbf{w}}^{k}+\mathbf{T}^{k}_{IC}-\hat{\mathbf{T}}^{k})+2p_{g}\Phi\hat{\mathbf{w}}^{k}=\mathbf{0},$

being $\Theta$ the matrix defined in (3.31). Therefore, for each $k$, $1\leq
k\leq P_{t}$, the solution of the inverse problem, $\hat{\mathbf{w}}^{k}$, is
obtained by solving the linear system

(3.68)
$\left(\Theta^{T}\Theta+2p_{g}\Phi\right)\hat{\mathbf{w}}^{k}=\Theta^{T}(\hat{\mathbf{T}}^{k}-\mathbf{T}^{k}_{IC}).$

Similarly, for the piecewise linear case, thanks to (3.49) and (3.58), we
rewrite (3.66) as

(3.69)
${(\tilde{\Theta})}^{T}\left(\tilde{\Theta}\hat{\mathbf{w}}^{k}-\Theta_{d}\hat{\mathbf{w}}^{k-1}+\mathbf{T}^{k}_{IC}-\hat{\mathbf{T}}^{k}\right)+2p_{g}\Phi\hat{\mathbf{w}}^{k}=\mathbf{0},$

being $\tilde{\Theta}$ and $\Theta_{d}$ the matrices defined in (3.48).
Therefore, for each $k$, $1\leq k\leq P_{t}$, a solution of the inverse
problem, $\hat{\mathbf{w}}^{k}$, is obtained by solving the linear system

(3.70)
$\left(\tilde{\Theta}^{T}\tilde{\Theta}+2p_{g}\Phi\right)\hat{\mathbf{w}}^{k}=\tilde{\Theta}^{T}(\hat{\mathbf{T}}^{k}+\Theta_{d}\hat{\mathbf{w}}^{k-1}-\mathbf{T}^{k}_{IC}).$

The resulting inverse solvers are straightforward modifications of Algorithm 1
and 2. Then, we show them in the following Algorithm 3 and 4.

Algorithm 3 Inverse solver for the solution of Problem 3.5 with piecewise
constant parameterization in time of the heat flux, $g$.

OFFLINE
Input RBF shape parameter, $\eta$; thermocouples measurement points and times,
$\Psi,\Upsilon$; cost functional parameter, $p_{g}$

1:Setup RBF parameterization by (3.8)

2:Compute $T_{\phi_{i}}$ for $i=1,2,\dots,P$ by solving Problem 3.18

3:Assemble matrix $\Theta$ by (3.31)

4:Assemble matrix $\Phi$ by (3.65)

ONLINE
Input Initial condition, $T_{0}$

1:Set $k=1$

2:while $k\leq P_{t}$ do

3: Read the thermocouples measurements, $\hat{\mathbf{T}}^{k}$

4: Compute $T_{IC}^{k}$ by solving Problem 3.20

5: Assemble $\mathbf{T}_{IC}^{k}$

6: Compute $\hat{\mathbf{w}}^{k}$ by solving (3.68)

7: Compute $g_{i}^{k}(t)$ by (3.9)

8: Compute the heat flux $\mathscr{g}(\mathbf{x},t)$ for
$t\in(\tau^{k-1},\tau^{k}]$ by (3.6) and (3.9)

9: Use (3.23) to compute $T^{k}[\mathscr{g}^{k}]$

10: $k=k+1$

11:end while

Algorithm 4 Inverse solver for the solution of Problem 3.5 with piecewise
linear parameterization in time of the heat flux, $g$.

OFFLINE
Input RBF shape parameter, $\eta$; thermocouples measurement points and times,
$\Psi,\Upsilon$; cost functional parameter, $p_{g}$

1:Setup RBF parameterization by (3.8)

2:Compute $T_{\phi_{i}}$ for $i=1,2,\dots,P$ by solving Problem 3.18

3:Compute $T_{d_{i}}$ for $i=1,2,\dots,P$ by solving Problem 3.37

4:Assemble matrices $\tilde{\Theta}$ and $\Theta_{d}$

5:Assemble matrix $\Phi$ by (3.65)

ONLINE
Input Initial condition, $T_{0}$

1:Set $k=1$

2:while $k\leq P_{t}$ do

3: Read the thermocouples measurements, $\hat{\mathbf{T}}^{k}$

4: Compute $T_{IC}^{k}$ by solving Problem 3.20

5: Assemble $\mathbf{T}_{IC}^{k}$

6: Compute $\hat{\mathbf{w}}^{k}$ by solving (3.70)

7: Compute $g_{i}^{k}(t)$ by (3.11)

8: Compute the heat flux $\mathscr{g}(\mathbf{x},t)$ for
$t\in(\tau^{k-1},\tau^{k}]$ by (3.6) and (3.11)

9: Use (3.39) to compute $T^{k}[\mathscr{g}^{k}]$

10: $k=k+1$

11:end while

As a final remark, notice that for $p_{g}=0\frac{K^{2}}{W^{2}}$, we end up
with the same solution as for $S_{1}^{k}$.

###  3.5. Regularization

After the development of novel inverse solvers, we provide a brief discussion
about regularization. It is well known that inverse problems as the ones here
considered are ill-posed. This means that for our problem at least one of the
following properties does not hold: for all admissible data, a solution
exists; for all admissible data, the solution is unique; the solution depends
continuously on the data [46]. In our discussion, we turned the infinite
dimensional inverse Problem 3.4 into the solution of the discrete linear
systems (3.35) and (3.51) by making some assumptions on the heat flux (i.e.
parameterizing it). In this new setting, if the matrices $\Theta^{T}\Theta$
and $\tilde{\Theta}^{T}\tilde{\Theta}$ are invertible, we have the existence
of a unique solution for our inverse problem.

As we will see in the numerical tests section, it turns out that these
matrices are very ill-conditioned. This can cause the matrix to be numerically
rank deficient, losing the uniqueness of a solution. However, this is not the
only concern. We still have the problem of a continuous dependence of the
solution on the data. The ill-conditioning of the linear system causes that,
if we have some noise in the data vector (as usual in an industrial
measurement equipment), the solution of the linear system diverges from the
correct value.

To address both these problems, we require regularization. There are several
techniques available for regularizing a discrete ill-posed problem as the
present one. In general, they are divided into direct methods like Truncated
Singular Values Decomposition (TSVD) and Tikhonov regularization, and
iterative methods such as the conjugate gradient method. For a deep discussion
of all regularization methods, we refer the interested reader to Hansen’s
monograph on the subject[47].

In the present investigation, we use TSVD. To briefly describe this
regularization technique, we denote the Singular Values Decomposition (SVD) of
a matrix $K$ by

(3.71) $K=U\Sigma
V^{T}=\sum_{i=1}^{r}\mathbf{u}_{i}\sigma_{i}\mathbf{v}_{i}^{T},$

where $\sigma_{i}$ denotes the i-th singular value of $K$ (numbered according
to their decreasing value), $r$ denotes the last no null singular value (i.e.
the rank of $K$), $\mathbf{u}_{i}$ and $\mathbf{v}_{i}$ are the i-th column of
the semi-unitary matrices $U$ and $V$ respectively (both belonging to ${\rm
I\\!M}^{P\times r}$), and $\Sigma$ is the square matrix of ${\rm
I\\!M}^{r\times r}$ such that $\Sigma_{ii}=\sigma_{i}$ and $\Sigma_{ij}=0$ if
$i\neq j$. Then, given $\alpha_{TSVD}\leq r$, the TSVD regularized solution of
the general linear system $K\mathbf{z}=\mathbf{c}$ is

(3.72)
$\mathbf{z}=\sum_{i=1}^{\alpha_{TSVD}}\left(\frac{\mathbf{u}_{i}^{T}\mathbf{c}}{\sigma_{i}}\right)\mathbf{v}_{i}.$

This solution differs from the least square solution only in that the sum is
truncated at $i=\alpha_{TSVD}$ instead of $i=r$. In this way, we cut off the
smallest singular values that are responsible of the errors propagation. For a
detailed discussion on the solution of discrete ill-posed inverse problems, we
refer the reader to Hansen’s monograph on the subject [47].

Together with the classical aforementioned regularization method, we
investigate also the regularization by discretization[48, 49]. Using this
method, we exploit the regularizing properties of coarsening the time and/or
space discretization to improve the heat flux estimation. In the next section,
we will test the performance of these regularization methods also by adding
noise to the thermocouples measurements.

###  3.6. Discretization Selection Algorithm

To conclude this section, we propose an algorithm for the automated selection
of some of the parameters required by Algorithm 4. As will be shown in Section
4, the numerical tests highlight that this algorithm is very sensitive to the
mesh and time discretization refinement as well as to the parameter $p_{g}$.
We anticipate here that this inverse solver shows severe instabilities for
fine discretizations. However, these instabilities are effectively eliminated
for values of $p_{g}$ that are above a threshold that depends on the
discretization refinement. In for these values of $p_{g}$, we notice a drastic
decrease of the dependency of the algorithm from the discretization.

However, the uncontrolled increase in $p_{g}$ does not lead to a monotonic
improvement of the inverse solver performances. As can be observed in the
numerical results of Section 4 (see Figures 15, 17, 19, 31, 33), the
dependency of the algorithm from $p_{g}$ is such that it is unstable for low
values of $p_{g}$ then, increasing further $p_{g}$, it sharply achieves an
optimum of performance before reaching a plateau at which the algorithm is
stable but the term $\langle
g^{k}(\tau^{k}),g^{k}(\tau^{k})\rangle_{L^{2}(\Gamma_{s_{in}})}$ in (3.5)
overcomes the measurements distance one, $S_{1}^{k}$ defined in (3.4). Thus,
for too high values of $p_{g}$, we have a stable algorithm that is almost
independent from the discretization refinement but that provides poor heat
flux estimations.

To allow an industrial use of the proposed inverse solver, the objective of
this section is to develop a method for automatically selecting the $\Delta
t$, the mesh and the value of $p_{g}$ such that the algorithm is stable and
accurately estimates the mold-steel heat flux.

In developing such method, we assume to have available a reliable dataset of
thermocouples measurements, $\hat{T}_{train}(\Psi,\Upsilon_{train})$, that we
can use to perform this tuning offline. Moreover, we assume that,
independently from the mold physical parameters and the heat flux values, this
inverse solver always shows the previously described behavior with respect to
$p_{g}$. In particular, we assume that, for values of $p_{g}$ higher than a
problem specific threshold, the algorithm is stable for all the
discretizations and independent from them (i.e. we obtain similar solutions
for any given mesh and $\Delta t$).

We recall, that in the real industrial case, we do not have any information
about the true heat flux that we want to estimate. Thus, this selection
methodology cannot be based on the heat flux estimation error. However,
Figures 21 and 35 show that the measurement discrepancy functional $S_{1}^{k}$
and the heat flux estimation error have a similar behavior as functions of
$p_{g}$ and we will use this quantity to determine the quality of the heat
flux estimation.

All that said, we begin by selecting an ordered set of meshes $(\Delta x_{1},$
$\Delta x_{2},$ $\dots,$ $\Delta x_{n_{M}})$ and an ordered set of timestep
sizes $(\Delta t_{1},$$\Delta t_{2},$$\dots,$$\Delta t_{n_{t}})$. We order
them from the finest to the coarsest discretization (i.e. $\Delta x_{1}<$
$\Delta x_{2}<$ $\dots<$ $\Delta x_{n_{M}}$ and $\Delta t_{1}<$ $\Delta
t_{2}<$ $\dots<$ $\Delta t_{n_{t}}$). Then, our first objective is to identity
a $p_{g}^{0}$ within the aforementioned stability region.

To do it, we start with a tentative $p_{g}^{0}$. For this value of the
parameter, we solve the inverse problem on the training measurement dataset
for all $\Delta x$ and $\Delta t$. Let us denote by $T[\Delta x,\Delta t]$ the
corresponding solution. Having done so, we compute

(3.73) $\Delta T:=\max_{i,j,q,p}{\left\lVert T[\Delta x_{i},\Delta
t_{j}]-T[\Delta x_{q},\Delta
t_{p}]\right\rVert_{L^{\infty}\left((0,t_{f}];L^{2}(\Omega_{s})\right)}}.$

If we have

(3.74) $\Delta T>\Delta x_{n_{M}}+\Delta t_{n_{t}},$

we consider that the solution is too dependent on the discretization
refinement. Then, we increase the value of $p_{g}^{0}$ and redo the
calculations until

(3.75) $\Delta T\leq\Delta x_{n_{M}}+\Delta t_{n_{t}},$

is satisfied.

Once we find a value of $p_{g}^{0}$ within the stability region, we choose the
discretization setup ($\Delta x^{1}$, $\Delta t^{1}$) that corresponds to the
minimum for $p_{g}^{0}$ of

(3.76) $m_{S}[\Delta x,\Delta t,p_{g}]:=mean_{k}\left(S_{1}^{k}[\Delta
x,\Delta t,p_{g}]\right).$

Once we select $\Delta x^{1}$ and $\Delta t^{1}$, we choose $p_{g}^{1}$ as the
value of $p_{g}$ that minimizes $m_{S}[\Delta x^{1},\Delta t^{1},p_{g}]$.
Then, we fix $p_{g}=p_{g}^{1}$ and we solve again the inverse problem for all
the considered discretization setups. If the previously selected
discretization is the one that corresponds to the lowest value of
$m_{S}[\Delta x,\Delta t,p_{g}^{1}]$, we choose $\Delta x^{1}$, $\Delta
t^{1}$, and $p_{g}=p_{g}^{1}$, and we stop the process. Otherwise, we continue
iterating by selecting $\Delta x^{2}$ and $\Delta t^{2}$ as the ones
corresponding to the smallest $m_{S}[\Delta x,\Delta t,p_{g}^{1}]$ and looking
for the $p_{g}^{2}$ that minimizes $m_{S}[\Delta x^{2},\Delta t^{2},p_{g}]$,
and so on. We summarize all this process in Algorithm 5.

Algorithm 5 Offline selection of the mesh, the timestep size and $p_{g}$ for
the inverse solver in Algorithm 4.

Input Ordered set of meshes, $(\Delta x_{1},$$\Delta x_{2},$$\dots,$$\Delta
x_{n_{M}})$; ordered set of timestep sizes, $(\Delta t_{1},$$\Delta
t_{2},$$\dots,$$\Delta t_{n_{t}})$; $p_{g}^{0}$; training set,
$\hat{T}_{train}(\Psi,\Upsilon_{train})$

1:while $\Delta T>\Delta x_{n_{M}}+\Delta t_{n_{t}}$ do$\triangleright$
Identify stability region

2: for $i=1$ to $n_{M}$ do

3: for $j=1$ to $n_{t}$ do

4: Solve the inverse problem on the training set, $\hat{T}_{train}$, by using
Algorithm 4 with $\Delta x_{i}$, $\Delta t_{j}$ and $p_{g}=p_{g}^{0}$

5: Compute $m_{S}[\Delta x_{i},\Delta t_{j},p_{g}^{0}]$ by (3.76)

6: end for

7: end for

8: Compute $\Delta T$ by (3.75)

9: if $\Delta T>\Delta x_{n_{M}}+\Delta t_{n_{t}}$ then

10: $p_{g}^{0}=10p_{g}^{0}$

11: end if

12:end while

13:Choose $\Delta x^{0}$ and $\Delta t^{0}$ corresponding to
$\min_{i,j}m_{S}[\Delta x_{i},\Delta t_{j},p_{g}^{0}]$ for $1\leq i\leq
n_{M}$, $1\leq j\leq n_{t}$

14:$l=1$, $f=0$

15:while $f=0$ do

16: Find $p_{g}^{l}=\operatorname*{argmin}_{p_{g}}(m_{S}[\Delta x^{l-1},\Delta
t^{l-1},p_{g}])$

17: for $i=1$ to $n_{M}$ do

18: for $j=1$ to $n_{t}$ do

19: Solve the inverse problem on the training set, $\hat{T}_{train}$, by using
Algorithm 4 with $\Delta x_{i}$, $\Delta t_{j}$ and $p_{g}=p_{g}^{l}$

20: Compute $m_{S}[\Delta x_{i},\Delta t_{j},p_{g}^{l}]$ by (3.76)

21: end for

22: end for

23: Choose $\Delta x^{l}$ and $\Delta t^{l}$ corresponding to
$\min_{i,j}m_{S}[\Delta x_{i},\Delta t_{j},p_{g}^{l}]$ for $1\leq i\leq
n_{M}$, $1\leq j\leq n_{t}$

24: if $\Delta x^{l}=\Delta x^{l-1}\And\Delta t^{l}=\Delta t^{l-1}$ then

25: $f=1$

26: end if

27: $l=l+1$

28:end while

This method allows a data-driven, automated selection of the discretization
refinement and the $p_{g}$ parameter. This result comes to the cost of
computing $n_{M}\cdot n_{t}$ solutions to the inverse problem at each
iteration. If the available memory allows it, we can keep in the memory the
results of the offline computations related to each discretization. Otherwise,
we have the recompute every time these offline phases. However, we designed
this algorithm to be used offline. Then, even if it is computationally
expensive, we can run it before the caster starts to work and it only requires
the dataset of thermocouples measurements $\hat{T}_{train}$.

##  4\. Numerical Tests

To test the previously developed methodologies, we design different benchmark
cases. Through these tests, we validate and analyze the performances of the
inverse solvers that we proposed in the previous sections. We design two
benchmarks to perform different tests for the inverse problem solvers proposed
in Section 3.

Notice that all the computations are performed in ITHACA-FV[50, 51] which is a
C++ library based on OpenFOAM[52] developed at SISSA Mathlab.

###  4.1. Benchmark 1

In this section, we test the performances of the inverse solvers proposed in
Section 3 in the reconstruction of a linear in time heat flux, which is non-
linear in space.

####  4.1.1. Setup of the Test Case

To design a numerical test case for the inverse problems, we proceed as
follows: we arbitrarily define a boundary heat flux, $g_{tr}(\mathbf{x},t)$,
and the thermocouples positions, $\Psi$, and sampling frequency, $f_{samp}$.
Then, we solve the direct Problem 2.3 associated with $g_{tr}(\mathbf{x},t)$
in the time domain $(0,t_{f}]$, obtaining the related temperature field.
Finally, we use its values at the thermocouples points and sampling times as
input measurements to the inverse problem, $\hat{\mathbf{T}}$. Using this
approach, we are able to analyze the inverse problem performance in the
reconstruction of the boundary heat flux, $g_{tr}(\mathbf{x},t)$.

Table 1 shows the geometrical and physical parameters selected for the present
benchmark case. In the attempt of mimicking the real industrial situation of
estimating the boundary heat flux in a plate of a CC mold, these parameters
are close to real industrial values. We use the computational domain in Figure
4(a) where $L$, $W$ and $H$ are set as in a real mold plate. Finally, Figure
4(b) shows the thermocouple locations.

(a) Domain

(b) Thermocouples positions

Figure 4. Schematic of the domain used in the benchmark test cases (a) and position of the 100 thermocouples at the plane $y=0.02~{}m$ (b) used for the inverse solver tests (images taken from Morelli et al.[3]). Table 1. Geometrical and physical parameters used for the benchmark test cases. Parameter | Value
---|---
Thermal conductivity, $k_{s}$ |  $383$ W/(m K)
Density, $\rho_{s}$ |  $8940$ kg/m3
Specific heat capacity, $C_{p_{s}}$ |  $390$ J/(kg K)
Heat transfer coefficient, $h$ |  $5.66e4$ W/(m2 $\cdot$ K)
Water temperature, $T_{f}$ |  $350$ K
Initial condition, $T_{0}$ |  $350$ K
$L$ |  $2$ m
$W$ |  $0.1$ m
$H$ |  $1.2$ m
Sampling frequency, $f_{samp}$ |  $1$ Hz
$a$ |  $1100$ W/(m2$\cdot$ s)
$b$ |  $1200$ W/(m2$\cdot$ s)
$c$ |  $3000$ W/(m2$\cdot$ s)
Final time, $t_{f}$ |  $50$ s

To test the effect of the regularization by discretization, we use different
space and time discretizations. For the time discretization, we use
homogeneous time discretization with $\Delta t=0.1s$, $0.2s$, $0.25s$ and
$0.5s$. For the space discretization, we use the uniform, structured,
orthogonal, hexahedral meshes presented in Table 2.

Table 2. Summary of the different meshes used in the numerical tests. | Mesh 1 | Mesh 2 | Mesh 3 | Mesh 4 | Mesh 5
---|---|---|---|---|---
Number of elements | $1.7e5$ | $4.5e4$ | $2.1e4$ | $7.5e3$ | $1.5e3$

For this test case, we select the heat flux $g_{tr}$ to be linear in time and
quadratic in space. In particular, given $g_{1}(\mathbf{x})=bz^{2}+c$ with $b$
and $c$ as in Table 1, we select the heat flux

(4.1)
$g_{tr}(\mathbf{x},t)=-k_{s}\left(0.5tg_{1}(\mathbf{x})+g_{1}(\mathbf{x})\right).$

Moreover, to analyze the performance of the inverse solvers, we introduce the
relative error

(4.2)
$e_{rel}(\mathbf{x},t):=\frac{g_{tr}(\mathbf{x},t)-g_{c}(\mathbf{x},t)}{g_{tr}(\mathbf{x},t)},$

where $g_{c}$ is the heat flux computed with the different methodologies
described in Section 3 that will be tested in the following.

####  4.1.2. Effect of Time and Space Discretization Refinement

Now, provided all the details for the first benchmark setup, we can proceed
presenting the results. Firstly, we show the effects of mesh and time
discretization refinement. To do it, we do not add any noise to the
temperature measurements and do not apply any regularization in the solution
of the linear systems solving them by a LU factorization with full pivoting.

We start by analyzing the case in which we minimize the functional $S_{1}^{k}$
(i.e. $p_{g}=0$). We show in Figure 6 the maximum and mean value of the
$L^{2}$\- and $L^{\infty}$-norm of the relative error, $e_{rel}$, in the
interval $(0,t_{f}]$, as the time and space discretization changes for
Algorithm 1 (i.e. piecewise constant approximation in time of the heat flux).

(a) $L^{2}$-norm of the relative error as a function of $\Delta t$

(b) $L^{\infty}$-norm of the relative error as a function of $\Delta t$

(c) $L^{2}$-norm of the relative error as a function of the mesh size

(d) $L^{\infty}$-norm of the relative error as a function of the mesh size

Figure 6. Maximum (circles) and mean (squares) values of the $L^{2}$\- and
$L^{\infty}$-norm of the relative error, $e_{rel}$, in the interval
$(0,t_{f}]$, for Benchmark 1 as the time and space discretization changes for
Algorithm 1 (piecewise constant time approximation of the heat flux and
$p_{g}=0\frac{K^{2}}{W^{2}}$).

From the figures, we appreciate on one side that the time discretization
coarsening has very little effects on Algorithm 1 with a small decrease of the
error as $\Delta t$ increases. On the other, the space discretization does not
have any effect on this inverse solver.

We now perform the same test for the piecewise linear time approximation of
Algorithm 2. Similarly, Figure 8 shows the maximum and mean value of the
$L^{2}$\- and $L^{\infty}$-norm of the relative error, $e_{rel}$, in the
interval $(0,t_{f}]$, as the time and space discretization changes.

(a) $L^{2}$-norm of the relative error as a function of $\Delta t$

(b) $L^{\infty}$-norm of the relative error as a function of $\Delta t$

(c) $L^{2}$-norm of the relative error as a function of the mesh size

(d) $L^{\infty}$-norm of the relative error as a function of the mesh size

Figure 8. Maximum (circles) and mean (squares) values of the $L^{2}$\- and
$L^{\infty}$-norm of the relative error, $e_{rel}$, in the interval
$(0,t_{f}]$, for Benchmark 1 as the time and space discretization changes for
Algorithm 2 (piecewise linear time approximation of the heat flux and
$p_{g}=0\frac{K^{2}}{W^{2}}$).

In this setting, the obtained results are very different from the previous
case. First of all, we notice a massive influence of both the space and time
discretization refinement on the performances of the inverse solver. As
anticipated in Section 3.5, the regularization by discretization plays an
important role as the algorithm performances are improved by several orders of
magnitude by the coarsening of the discretization. Moreover, when comparing
the results for Algorithm 1 and 2, we notice that the piecewise linear solver
is able to outperform the constant one by three orders of magnitude but is
also very unstable depending on the discretization.

To better understand the behavior of this inverse solver, Figure 9 illustrates
the $L^{2}$-norm of the relative error, $e_{rel}$, as a function of time for
mesh 3 with different $\Delta t$. From these results, we see that the high
errors shown in Figure 8 are caused by diverging oscillations in the
algorithm. However, we also notice from the figure that, coarsening the time
discretization, monotonically reduces such instability until achieving a
stable solution, eventually.

Figure 9. $L^{2}$-norm of the relative error, $e_{rel}$, in Benchmark 1 for
Algorithm 2 (piecewise linear time approximation of the heat flux and
$p_{g}=0\frac{K^{2}}{W^{2}}$). The presented results are obtained with Mesh 3.

####  4.1.3. Effect of Cost Functional Parameter, $p_{g}$

In this section, we analyze the role that the cost functional parameter,
$p_{g}$, in (3.5), has on the performance of the proposed inverse solvers. To
do it, we solve several times this benchmark case using the different meshes
of Table 2 and different timestep sizes. Then, we plot the maximum and mean
value of the $L^{2}$-norm of the relative error, $e_{rel}$, over the entire
interval $t=(0,t_{f}]$ as a function of the cost functional parameter,
$p_{g}$.

We start with Algorithm 3 (i.e. piecewise constant approximation in time of
the heat flux). Figure 11 shows the obtained results for different timestep
sizes and a fixed space discretization.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 11. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 1 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 3 (piecewise
constant time approximation of the heat flux). We show the results for Mesh 3
and different $\Delta t$.

From the results, we notice that increasing the value of $p_{g}$ monotonically
decreases the quality of the reconstruction. Moreover, it is true for all
considered $\Delta t$ with a slight improvement of the performances as the
time discretization gets coarser.

Now, we perform a similar test but this time we keep $\Delta t=0.25~{}s$ and
test the different meshes of Table 2. We illustrate in Figure 13 the obtained
results.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 13. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 1 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 3 (piecewise
constant time approximation of the heat flux). We show the results for $\Delta
t=0.25~{}s$ and different meshes.

This figure confirms that Algorithm 3 is badly affected by the implementation
of the second term of (3.5). In fact, its performance dramatically
deteriorates as soon as this term begins to play a role (i.e. $p_{g}\gtrsim
10^{-12}K^{2}/W^{2}$). Moreover, the results are almost independent from the
discretization refinement. This further confirms the insensibility of this
algorithm from the used discretization.

We continue by performing the same kind of tests on Algorithm 3 (i.e.
piecewise linear approximation of the heat flux in time). We start by testing
different timestep sizes while using Mesh 3 for the space discretization. We
present the results in Figure 15.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 15. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 1 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 4 (piecewise
linear time approximation of the heat flux). We show the results for Mesh 3
and different $\Delta t$.

At first, we notice that this algorithm has a very different behavior with
respect to the piecewise constant case. In this case, the timestep size
dramatically affects the results. We can depict two different behaviors as
$p_{g}$ changes for a chosen $\Delta t$. In the first one (i.e. $\Delta
t=0.1~{}s$, $0.2~{}s$ and $0.25~{}s$), the inverse solver is very unstable and
provides completely useless solutions for low values of $p_{g}$ (i.e.
$p_{g}\lesssim 10^{-12}K^{2}/W^{2}$). As $p_{g}$ increases, the quality of the
approximation rapidly rises up until the error reaches a minimum. Here, we
have stable solutions and a good approximation of the heat flux. For higher
values of $p_{g}$, the error monotonically increases until it reaches a
plateau at $100\%$.

On the other hand, we have a different behavior for $\Delta t=0.5~{}s$. In
this case, the inverse solver performs similarly to the piecewise constant
case, but the quality of the estimation is by almost two orders of magnitude
better. Then, we have stable and accurate solutions for low values of $p_{g}$.
For $p_{g}\gtrsim 10^{-12}K^{2}/W^{2}$ we have a monotonic degradation of the
heat flux estimation until we reach the $100\%$ plateau.

It is interesting to notice that the second term in the functional $S_{2}^{k}$
can make the solver insensible to the discretization refinement. In fact,
after a certain value of $p_{g}$, the relative error norms for the different
$\Delta t$ are almost coincident.

We can see a similar behavior in Figure 17 where we show the results obtained
refining the mesh and keeping $\Delta t=0.25~{}s$. Also in this case, we
notice the two previously described behaviors with the coarsest mesh being
always stable and providing the best results for the lowest values of $p_{g}$.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 17. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 1 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 4 (piecewise
linear time approximation of the heat flux). We show the results for $\Delta
t=0.25s$ and different meshes.

These statements are remarked by the results shown in Figure 19, where we show
the results of the same test as in Figure 8 but for
$p_{g}=5e-11\frac{K^{2}}{W^{2}}$. The obtained results confirm that for some
values of $p_{g}$ we can obtain a stable solver with a moderate dependency on
the discretization refinement.

(a) $L^{2}$-norm of the relative error as a function of $\Delta t$

(b) $L^{\infty}$-norm of the relative error as a function of $\Delta t$

(c) $L^{2}$-norm of the relative error as a function of the mesh size

(d) $L^{\infty}$-norm of the relative error as a function of the mesh size

Figure 19. Maximum (circles) and mean (squares) values of the $L^{2}$\- and
$L^{\infty}$-norm of the relative error, $e_{rel}$, in the interval
$(0,t_{f}]$, for Benchmark 1 as the time and space discretization changes for
Algorithm 4 (piecewise linear time approximation of the heat flux and
$p_{g}=5e-11\frac{K^{2}}{W^{2}}$).

To conclude this analysis, we test the discretization and $p_{g}$ selection
method of Algorithm 5 in this benchmark case. We use the virtual thermocouples
measurements as input training dataset for the algorithm, $\hat{T}_{train}$.
In the test, the algorithm has to select a combination of mesh, timestep size
and $p_{g}$ that provides stable and accurate solutions to this inverse
problem. Before presenting the results of Algorithm 5, we show in Figure 21
the mean value of the temperature discrepancy functional $S_{1}^{k}$, defined
in (3.4), as function of $p_{g}$ for different meshes and $\Delta t$.

(a) Mean of $S_{1}^{k}$ for Mesh 3.

(b) Mean of $S_{1}^{k}$ for $\Delta t=0.25s$.

Figure 21. Mean values of $S_{1}^{k}$ for $1\leq k\leq P_{t}$, for Benchmark 1
as the value of the cost function parameter, $p_{g}$, changes. The results are
obtained using Algorithm 4 (piecewise linear time approximation of the heat
flux). We show the results for Mesh 3 and different $\Delta t$ in (a), and for
$\Delta t=0.25$ and different meshes in (b).

These results show that $m_{S}$ in (3.76) behaves similarly to the relative
error (4.2) as function of $p_{g}$. Notice that it presents the same two
behaviors that we previously described for the relative error depending on the
used discretization setup. Moreover, the $m_{S}$ minima are close to those of
the relative error. For these reasons, in Section 3.6, we used this result in
the creation of the selection criteria for the $p_{g}$ as well as for the mesh
and the $\Delta t$.

That said, we are now ready to test Algorithm 5. With respect to its
implementation, in step 16 of the algorithm, we use the Nelder-Mead method to
find the $p_{g}$ that minimizes $S_{1}^{k}$[53]. To start the algorithm, we
set $p_{g}^{0}=1e-7~{}\frac{K^{2}}{W^{2}}$. Table 3 summarizes the algorithm
behavior.

Table 3. Test of Algorithm 5 for Benchmark 1. Iteration | Mesh | $\Delta t~{}[s]$ | $p_{g}~{}\left[\frac{K^{2}}{W^{2}}\right]$ | $mean_{k}\left(S_{1}^{k}\right)~{}\left[K^{2}\right]$
---|---|---|---|---
0 | 5 | $0.5$ | $1e-7$ | $6.9e3$
1 | 5 | $0.5$ | $3.2e-21$ | $3.9e-1$

From the results in the table, we appreciate that the algorithm chooses the
coarsest discretization since the first iteration. Then, it looks for the
$p_{g}$ that minimizes $S_{1}^{k}$ for this discretization and, not finding a
better discretization setup for this value of $p_{g}$, exits the process.
Comparing the obtained results to the relative error plots of Figures 15 and
17, we confirm that the algorithm is selecting the best configuration in
between all the available.

###  4.2. Benchmark 2

In designing this benchmark case, we use the same geometrical and physical
parameters as in Benchmark 1, but we choose a non linear in time true heat
flux, $g_{tr}$, as in Table 4.

Table 4. Parameters used for the unsteady Benchmark 2. Parameter | Value
---|---
$g_{1}(\mathbf{x})$ | $bz^{2}+c$
$g_{2}(\mathbf{x})$ | $\frac{10c}{1+(x-1)^{2}+z^{2}}$
Heat flux, $g_{tr}(\mathbf{x},t)$ |  $-k_{s}\left[g_{1}+\frac{g_{1}}{2}\sin{\left(2\pi f_{max}\frac{t^{2}}{t_{f}}\right)}+g_{2}e^{-0.1t}\right]$ W/m2
Maximum frequency, $f_{max}$ |  $0.1$ Hz

Also for this benchmark, we test both the piecewise constant and linear
approximation algorithms. In particular, we investigate the regularization
properties of the discretization coarsening, the effects of the $p_{g}$
parameter and the robustness of the algorithms to noise in the measurements.

####  4.2.1. Effect of Time and Space Discretization Refinement

The first test we do is related to the space and time discretization
refinement. To test the effect of changing the space and time discretization
sizes, we reconstruct the heat flux $g_{tr}(\mathbf{x},t)$ for all the meshes
in Table 2 and $\Delta t=0.1$, $0.2$, $0.25$ and $0.5~{}s$. In these tests, we
use the cost functional (3.4) (i.e. $p_{g}=0~{}\frac{K^{2}}{W^{2}}$) and do
not add noise to the measurements.

First, we test Algorithm 1. Figure 23 illustrates the mean and maximum values
in $(0,t_{f}]$ of the $L^{2}$\- and $L^{\infty}$-norm of the relative error
(4.2). The results show a entirely similar behavior to that of the previous
benchmark as described in Section 4.1.2.

(a) $L^{2}$-norm of the relative error as a function of $\Delta t$

(b) $L^{\infty}$-norm of the relative error as a function of $\Delta t$

(c) $L^{2}$-norm of the relative error as a function of the mesh size

(d) $L^{\infty}$-norm of the relative error as a function of the mesh size

Figure 23. Maximum (circles) and mean (squares) values of the $L^{2}$\- and
$L^{\infty}$-norm of the relative error, $e_{rel}$, in the interval
$(0,t_{f}]$, for Benchmark 2 as the time and space discretization changes for
Algorithm 1 (piecewise constant time approximation of the heat flux and
$p_{g}=0\frac{K^{2}}{W^{2}}$).

(a) $L^{2}$-norm of the relative error as a function of $\Delta t$

(b) $L^{\infty}$-norm of the relative error as a function of $\Delta t$

(c) $L^{2}$-norm of the relative error as a function of the mesh size

(d) $L^{\infty}$-norm of the relative error as a function of the mesh size

Figure 25. Maximum (circles) and mean (squares) values of the $L^{2}$\- and
$L^{\infty}$-norm of the relative error, $e_{rel}$, in the interval
$(0,t_{f}]$, for Benchmark 2 as the time and space discretization changes for
Algorithm 2 (piecewise linear time approximation of the heat flux and
$p_{g}=0\frac{K^{2}}{W^{2}}$).

####  4.2.2. Effect of Cost Functional Parameter, $p_{g}$

Also for this benchmark case, we test the effects that the parameter $p_{g}$
in (3.5) has on the inverse solvers performances. We do it by performing the
same tests of Section 4.1.3 but for the present test case. In particular, we
solve this inverse problem using all the meshes in Table 2 and $\Delta
t=0.1,0.2,0.25$ and $0.5~{}s$, for $1e-16\leq p_{g}\leq
1e-6\frac{K^{2}}{W^{2}}$. We use both the piecewise constant and linear
approximation of the heat flux in Algorithm 3 and 4, respectively. Notice that
in these tests, we use the full order algorithms.

First, we present in Figures 27 and 29 the results obtained using the
piecewise constant approximation algorithm. These plots show the mean and
maximum values of the $L^{2}$-norm of the relative error, $e_{rel}$, in the
interval $(0,t_{f}]$ as functions of $p_{g}$ changing the mesh and the $\Delta
t$, respectively.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 27. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 2 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 3 (piecewise
constant time approximation of the heat flux). We show the results for Mesh 3
and different $\Delta t$.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 29. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 2 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 3 (piecewise
constant time approximation of the heat flux). We show the results for $\Delta
t=0.25s$ and different meshes.

From the presented results, we notice that Algorithm 3 has the same behavior
shown in the previous benchmark case. In particular, this inverse solver
confirms to be badly affected by $p_{g}>0\frac{K^{2}}{W^{2}}$. The effect of
$p_{g}$ on its performance is very non linear with a first region of no
effects for $p_{g}\lesssim 1e-11\frac{K^{2}}{W^{2}}$ followed by a steep
degradation and a plateau at $100\%$ relative error for $p_{g}\gtrsim
1e-7\frac{K^{2}}{W^{2}}$. Moreover, the different meshes and $\Delta t$ have
the same behavior and similar values of the relative error. These results
confirm once more the almost insensibility of this inverse solver to the
discretization refinement.

Now, we consider the piecewise linear approximation of Algorithm 4. For this
inverse solver, Figures 31 and 33 show the mean and maximum values of the
$L^{2}$-norm of the relative error, $e_{rel}$, in the interval $(0,t_{f}]$ as
functions of $p_{g}$ changing the mesh and the $\Delta t$, respectively.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 31. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 2 as the value of
the cost function parameter, $p_{g}$, changes for Algorithm 4 (piecewise
linear time approximation of the heat flux). We show the results for Mesh 3
and different $\Delta t$.

(a) Mean relative error norm.

(b) Max. relative error norm.

Figure 33. Mean (a) and maximum (b) values of the $L^{2}$-norm of the relative
error, $e_{rel}$, in the interval $(0,t_{f}]$, for Benchmark 2 as the value of
the cost function parameter, $p_{g}$, changes. The results are obtained using
Algorithm 4 (piecewise linear time approximation of the heat flux). We show
the results for $\Delta t=0.25s$ and different meshes.

Again, the results are very similar to those of the previous benchmark. On one
hand, the coarsest discretizations show a similar behavior to the piecewise
constant approximation case with a monotonic degradation of the performance as
$p_{g}$ increases. On the other hand, we have unstable solutions for small
$p_{g}$ that are stabilized by $p_{g}\gtrsim 5e-11\frac{K^{2}}{W^{2}}$.
However, the accuracy of these solution rapidly decreases as we further
increase $p_{g}$ until we reach the $100\%$ relative error plateau.

To conclude, we test the mesh, $\Delta t$ and $p_{g}$ selection method of
Algorithm 5. In this test, we input to the algorithm the virtual thermocouples
measurements that we compute for this benchmark case. Before presenting the
results of this algorithm, we show in Figure 35 the mean value of $S_{1}^{k}$,
$m_{S}$, as function of $p_{g}$ for different meshes and $\Delta t$.

(a) Mean of $S_{1}^{k}$ for Mesh 3.

(b) Mean of $S_{1}^{k}$ for $\Delta t=0.25s$.

Figure 35. Mean values of $S_{1}^{k}$ for $1\leq k\leq P_{t}$, for Benchmark 2
as the value of the cost function parameter, $p_{g}$, changes. The results are
obtained using Algorithm 4 (piecewise linear time approximation of the heat
flux). We show the results for Mesh 3 and different $\Delta t$ in (a), and for
$\Delta t=0.25$ and different meshes in (b).

Notice that $m_{S}$ has a behavior that is very similar to the relative error
norm shown in Figures 31 and 33. Moreover, its minima are almost correspondent
to the relative error ones. As already mentioned, we used these results as a
guideline in developing Algorithm 5. We present its behavior for the present
benchmark case in Table 5.

Table 5. Test of Algorithm 5 for Benchmark 2. Iteration | Mesh | $\Delta t~{}[s]$ | $p_{g}~{}\left[\frac{K^{2}}{W^{2}}\right]$ | $mean_{k}\left(S_{1}^{k}\right)~{}\left[K^{2}\right]$
---|---|---|---|---
0 | 4 | $0.1$ | $8e-7$ | $8.46e2$
1 | 5 | $0.5$ | $3.6e-11$ | $6.2e0$
2 | 5 | $0.5$ | $3.2e-13$ | $5.5e-1$

In this case, the algorithm does not select since the beginning the coarsest
discretization. In the first iteration, it finds the $p_{g}$ that minimizes
$m_{S}$ for this setup. Then, when comparing it to the other discretizations
in step 23, it selects again the coarsest one. Also in this benchmark case,
this algorithm showed to be able to select the discretization setup and the
value of $p_{g}$ corresponding to the best performance of the inverse solver.

####  4.2.3. Effects of Measurements Noise and Regularization

In this section, we test the effect that adding noise to the measurements
vector, $\hat{\mathbf{T}}$, has in the performances of Algorithm 1 and 2. From
the industrial point of view, this analysis is of particular interest for our
application since in the real case, thermocouples measurements are affected by
noise.

We perform this analysis by adding to the measurements vector the Gaussian
random noise $\boldsymbol{\eta}=\mathcal{N}(\boldsymbol{\mu},\Sigma)$, where
$\boldsymbol{\mu}\in{\rm I\\!R}^{M}$ is the mean vector and $\Sigma\in{\rm
I\\!M}^{M\times M}$ is the covariance matrix. Then, we have

(4.3) $\hat{\mathbf{T}}^{k}_{\eta}=\hat{\mathbf{T}}^{k}+\boldsymbol{\eta}.$

In particular, we choose $\boldsymbol{\eta}$ to be a independent and
identically distributed random variable with zero mean, i.e.
$\boldsymbol{\eta}=\mathcal{N}(\mathbf{0},\omega^{2}I)$ , where $\omega$
denotes the noise standard deviation. To study the effect of noise, we perform
several solutions of the inverse problem using $\hat{\mathbf{T}}^{k}_{\eta}$
as thermocouples measurements. For each test, we compute 200 samples.

We show in Figures 37 and 39 the obtained results for the piecewise constant
and linear algorithm, respectively. In particular, we illustrate for each of
them the mean values over the samples of the mean and maximum of the relative
error (4.2) in $(0,t_{f}]$ (with 90% quantile bars) for different values of
the noise standard deviation, $\omega$. The figure compares the results
obtained using LU with full pivoting and TSVD with different values for the
regularization parameter $\alpha_{TSVD}$.

(a) Mean.

(b) Maximum.

Figure 37. Effect of the noise in the temperature measurements for Algorithm 1
in Benchmark 2. In the figures, we show the mean (a) and maximum (b) values of
the relative error (4.2) in $(0,t_{f}]$ for different values of the noise
standard deviation and using both LU with full pivoting and TSVD for the
solution of inverse problem linear system (3.35). For each case, we performed
200 runs. The markers show the mean values while the bars are the 90%
quantiles. In these computations, we considered $p_{g}=0\frac{K^{2}}{W^{2}}$.

(a) Mean.

(b) Maximum.

Figure 39. Effect of the noise in the temperature measurements for Algorithm 2
in Benchmark 2. In the figures, we show the mean (a) and maximum (b) values of
the relative error (4.2) in $(0,t_{f}]$ for different values of the noise
standard deviation and using both LU with full pivoting and TSVD for the
solution of inverse problem linear system (3.51). For each case, we performed
200 runs. The markers show the mean values while the bars are the 90%
quantiles. In these computations, we considered $p_{g}=0\frac{K^{2}}{W^{2}}$.

The results show a very different dependency from the measurement noise in the
two algorithm. The piecewise constant Algorithm 1 shows in Figure 37 to be
quite robust with respect to these levels of noise. Moreover, the TSVD is
effective in reducing the noise propagation and we are able to keep a
reasonable level of accuracy.

On the other hand, the piecewise linear Algorithm 2 is much more affected by
the noise. By using the TSVD regularization, we have an improvement of the
noise robustness. However, the error rate of increase is much higher than for
the piecewise constant solver.

####  4.2.4. Computational Cost

To conclude this analysis, we present in Tables 6 and 7 the numerical cost of
performing one iteration of the proposed inverse solvers. Notice that all the
computations were performed in serial on a Intel® Core™ i7-8550U CPU
processor. As expected the required CPU time increases with the refinement of
the discretization. Since we are using relative coarse meshes due to the
simplified geometry, the computational cost in many cases meets the real-time
requirement for this application (i.e. $1$ s). However, the meshes required
for the discretization of the real mold geometry are such that we cannot
ensure real-time performances in these cases.

Table 6. Average computational cost for one iteration of Algorithm 1 in Benchmark 2. $\Delta t$ Mesh | Mesh 1 | Mesh 2 | Mesh 3 | Mesh 4 | Mesh 5
---|---|---|---|---|---
$0.1$ s |  $5184.0$ ms |  $1402.6$ ms |  $815.6$ ms |  $482.8$ ms |  $343.4$ ms
$0.2$ s |  $2520.1$ ms |  $741.2$ ms |  $427.4$ ms |  $260.2$ ms |  $192.6$ ms
$0.25$ s |  $2052.9$ ms |  $600.6$ ms |  $350.9$ ms |  $215.8$ ms |  $162.1$ ms
$0.5$ s |  $1115.5$ ms |  $333.4$ ms |  $197.6$ ms |  $128.2$ ms |  $101.0$ ms
Table 7. Average computational cost for one iteration of Algorithm 2 in Benchmark 2. $\Delta t$ Mesh | Mesh 1 | Mesh 2 | Mesh 3 | Mesh 4 | Mesh 5
---|---|---|---|---|---
$0.1$ s |  $5966.4$ ms |  $1621.2$ ms |  $1442.5$ ms |  $517.3$ ms |  $400.1$ ms
$0.2$ s |  $2836.2$ ms |  $907.4$ ms |  $518.1$ ms |  $293.4$ ms |  $206.2$ ms
$0.25$ s |  $2258.8$ ms |  $707.6$ ms |  $382.5$ ms |  $232.4$ ms |  $182.9$ ms
$0.5$ s |  $1195.8$ ms |  $435.7$ ms |  $216.8$ ms |  $140.2$ ms |  $108.0$ ms

##  5\. Conclusions and Future Works

The goal of the present investigation was to develop mathematical tools to
monitor the mold behavior in CC machineries. At industrial experts’
suggestion, we identified the mold-steel heat flux as the quantity of interest
for the mold monitoring. Then, our objective has been to investigate a
methodology for the real-time estimation of this heat flux having as data the
physical properties of the mold, its geometry, the cooling water temperature
and some pointwise temperature measurements in the interior of the mold
plates.

We opted for stating the problem in a data assimilation, optimal control
setting in which we look for the heat flux that minimizes a functional that
includes a measure of the distance between the computed and measured
temperature at the measurement points. However, given a mold-steel heat flux,
we required a mold model to compute the temperature at these points.

In deriving the mold model, we had to take into account the real-time
requirement of this mold monitoring task. Then, we avoided to model all the
mold region, including the steel pool. In fact, it would require us to model
several complex and coupled physical phenomena (heat transfer, steel
solidification, steel and cooling water fluid flows, etc.) happening on
different time and space scales. With this level of complexity, the
computational cost of such simulation would have been unsustainable for real-
time computations. Thus, in modeling the mold, we considered as computational
domain the mold plates only. Hence, our first task has been the derivation of
the three dimensional unsteady-state heat conduction mold model.

Once the mold model has been established, we focused on the mold-steel heat
flux estimation. Notice that, in this setting, this flux is a Neumann BC on a
portion of the boundary of our domain. Then, we can generalize this
mathematical inverse problem as the estimation of a Neumann BC given pointwise
state measurements in the interior of the domain.

In this unsteady-state setting, we used a sequential approach to the inverse
problem. In fact, to provide a real-time solution in this setting means to
stay always at the front of the time line as it stretches. Then, since our
measurements come equally spaced in time by one second, we considered the
problem of estimating the heat flux only in between the last measurement and
the previous one, assuming to have already the solution for older times.

In this framework, we stated two different inverse problems. One looking for
the heat flux that minimizes a measure of the distance between computed and
measured temperature only. While, in the other, we want to minimize this
distance plus an heat flux norm.

For both these inverse problems, we developed novel methodologies for their
solution that exploit a RBFs parameterization of the heat flux in space with
time dependent coefficients. With respect to these coefficients, we considered
both the piecewise constant and the piecewise linear case. It means that the
estimated heat flux is constant or linear in between two contiguous
measurement instants.

These novel methodologies are direct methods that benefit from an offline-
online decomposition. Thanks to this decomposition, we have a first
computationally expensive offline phase, in which we solve several direct
problems. This offline phase is computed once and for all and does not require
any measurement. Then, when the caster starts to work, we only have to collect
the thermocouples measurements and run the online phase which is
computationally much cheaper.

To conclude, we tested the proposed inverse solvers on some benchmark cases.
We designed two benchmark cases for the inverse problem. To design an inverse
problem test, we arbitrarily selected a mold-steel heat flux. We solved the
direct problem for this heat flux obtaining the corresponding temperature
field in the mold. Finally, we located some virtual thermocouples and used the
computed temperature at these points as input data for the inverse solvers.
Then, we tested their ability to reconstruct the previously selected heat
flux. The two benchmark cases share the same domain and physical parameters
which were chosen to resemble the industrial setting. The difference is in the
selected heat fluxes.

From the obtained results, we noticed a great difference in the behavior of
the piecewise constant and linear inverse solvers. The former showed a very
stable behavior and insensitivity to the time and space discretization used.
The latter, rather, is very much influenced by the discretization used. In
particular, it can be very unstable when using fine discretizations but this
instability is reduced by coarsening the time and/or space discretization. In
fact, for some discretizations, we achieved very stable and accurate
solutions, eventually.

We also tested the effects that adding the heat flux norm to the minimization
functional has on these inverse solvers. We implemented this new term
multiplying it by a parameter. Then, we tested the effect that its value has
on the solvers performance.

We noticed that the piecewise constant algorithm performance monotonically
deteriorates as this parameter increases. The same goes for the piecewise
linear solver when using the coarsest, stable discretizations. However, the
unstable configurations showed to be positively affected by the addition of
this new term and, for some values of this parameter, we were able to obtain
stable and accurate solutions for all the tested discretizations. While, for
too high values of the parameter, the solution is stable but inaccurate for
all the meshes and timestep sizes. Moreover, we showed that, for values of the
parameter above a threshold, the inverse solver performance is almost
independent from the discretization.

Due to this dependency from the discretization and the functional parameter,
we developed an algorithm for the automatic selection of these quantities. In
the numerical tests, it proved to be able of a selection that corresponds to a
stable and accurate inverse solver.

Testing the inverse solvers for several noise levels showed again a different
behavior between the piecewise constant and linear approximations. The former
is much less sensitive to the measurements noise than the latter. For both,
the TSVD regularization proved to be able to mitigate the noise propagation
and we were able to obtain accurate and stable solutions also in the presence
of noise.

To conclude, we recall that the online phases of the proposed algorithms
require the solution of a full order problem whose computational cost depends
on the mesh and timestep size. As shown in the numerical tests section, it
means that we cannot ensure real-time performance for these algorithms as they
are. Then, in our future work, we will develop model order reduction
techniques that will allow us to reduce the online phase computational time
and make it independent from the discretization.[54, 55]

As a final remark, we discuss the application of the new proposed methodology
to other problems. Recalling that the presented continuous casting problem is
a Neumann BC estimation problem in a unsteady linear setting with pointwise
state measurements in the interior of the domain, we can apply the proposed
methodologies to any problem sharing these features. An example can be a
boundary stress estimation problem in linear elasticity with pointwise
deformation measurements.

Other possible future works on the subject could be related to the study of
theoretical results that can ensure a priori the stability and accuracy of
these inverse solvers with respect to the used discretization. It would
increase the potential of the proposed methodologies as well as their
reliability. In particular, it would be useful for the final user to know a
priori the time and space discretization to select as well as the minimization
functional parameter. Notice that it is needed for the piecewise linear
inverse solver because the piecewise constant one is almost insensible to the
discretization refinement.

In the future, it would also be interesting to investigate the use of a
completely different approach in the solution of this inverse problem. In
particular, thinking about a more proper handling of the measurement noise, we
could think of using a Bayesian approach.[56] Techniques such as ensemble
Kalman filter could be suitable for this problem given the sequentiality of
the measurements. Moreover, considering the real-time requirement of the
application, it would probably require an effective use of model order
reduction techniques to reduce the demanding computational cost of these
techniques.

##  6\. Acknowledge

We would like to acknowledge the financial support of the European Union under
the Marie Sklodowska-Curie Grant Agreement No. 765374. We also acknowledge the
partial support by the Ministry of Economy, Industry and Competitiveness
through the Plan Nacional de I+D+i (MTM2015-68275-R), by the Agencia Estatal
de Investigacion through project [PID2019-105615RB-I00/ AEI /
10.13039/501100011033], by the European Union Funding for Research and
Innovation - Horizon 2020 Program - in the framework of European Research
Council Executive Agency: Consolidator Grant H2020 ERC CoG 2015 AROMA-CFD
project 681447 ”Advanced Reduced Order Methods with Applications in
Computational Fluid Dynamics” and INDAM-GNCS project ”Advanced intrusive and
non-intrusive model order reduction techniques and applications”, 2019.
Moreover, we gratefully thank Gianfranco Marconi, Federico Bianco and Riccardo
Conte from Danieli & C.Officine Meccaniche SpA for helping us in better
understanding the industrial problem and for the fruitful cooperation.

## References

  * [1] World Steel Association . World steel in figures 2018. World Steel Association: Brussels, Belgium 2018\.
  * [2] Irving WR. Continuous casting of steel.  The Institute of Materials (UK) .  1993\.
  * [3] Morelli UE, Barral P, Quintela P, Rozza G, Stabile G. A numerical approach for heat flux estimation in thin slabs continuous casting molds using data assimilation. International Journal for Numerical Methods in Engineering 2021; 122(17): 4541-4574.  doi: https://doi.org/10.1002/nme.6713
  * [4] Nittka R. Inhomogeneous parabolic Neumann problems. Czechoslovak Mathematical Journal 2014; 64(3): 703-742.  doi: https://doi.org/10.48550/arXiv.1108.6227
  * [5] Eymard R, Gallouët T, Herbin R. Finite volume methods. In: . 7 of Handbook of Numerical Analysis. Elsevier. 2000 (pp. 713 - 1018).
  * [6] Ling X, Keanini RG, Cherukuri HP. A non-iterative finite element method for inverse heat conduction problems. International Journal for Numerical Methods in Engineering 2003; 56(9): 1315-1334.  doi: https://doi.org/10.1002/nme.614
  * [7] Loulou T, Scott EP. An inverse heat conduction problem with heat flux measurements. International Journal for Numerical Methods in Engineering 2006; 67(11): 1587-1616.  doi: https://doi.org/10.1002/nme.1674
  * [8] Jin B. Conjugate gradient method for the Robin inverse problem associated with the Laplace equation. International Journal for Numerical Methods in Engineering 2007; 71(4): 433-453.  doi: https://doi.org/10.1002/nme.1949
  * [9] Huang CH, Yan JY. An Inverse Problem In Predicting Temperature Dependent Heat Capacity Per Unit Volume Without Internal Measurements. International Journal for Numerical Methods in Engineering 1996; 39(4): 605-618.  doi: https://doi.org/10.1002/(SICI)1097-0207(19960229)39:4¡605::AID-NME872¿3.0.CO;2-H
  * [10] Alifanov OM. Inverse Heat Transfer Problems.  Moscow Izdatel Mashinostroenie.  1 ed. 1988.
  * [11] Orlande HRB. Inverse Problems in Heat Transfer: New Trends on Solution Methodologies and Applications. Journal of Heat Transfer 2012; 134(3): 1-13.  031011doi: 10.1115/1.4005131
  * [12] Beck JV, Blackwell B, Clair Jr. CRS. Inverse heat conduction: Ill-posed problems.  James Beck .  1985\.
  * [13] Chang CW, Liu CH, Wang CC. Review of Computational Schemes in Inverse Heat Conduction Problems. Smart Science 2018; 6(1): 94-103.  doi: 10.1080/23080477.2017.1408987
  * [14] Şahin HM, Kocatepe K, Kayıkcı R, Akar N. Determination of unidirectional heat transfer coefficient during unsteady-state solidification at metal casting–chill interface. Energy Conversion and Management 2006; 47(1): 19–34.  doi: 10.1016/j.enconman.2005.03.021
  * [15] Ranut P. Optimization and Inverse Problems in Heat Transfer. PhD thesis. Universitá degli Studi di Udine, Via delle Scienze, 206, 33100 Udine UD, Italy; 2012.
  * [16] Udayraj , Chakraborty S, Ganguly S, Chacko E, Ajmani S, Talukdar P. Estimation of surface heat flux in continuous casting mould with limited measurement of temperature. International Journal of Thermal Sciences 2017; 118: 435 - 447.  doi: https://doi.org/10.1016/j.ijthermalsci.2017.05.012
  * [17] Mahapatra R, Brimacombe J, Samarasekera I. Mold behavior and its influence on quality in the continuous casting of steel slabs: Part II. Mold heat transfer, mold flux behavior, formation of oscillation marks, longitudinal off-corner depressions, and subsurface cracks. Metallurgical and Materials Transactions B 1991; 22(6): 875–888.  doi: 10.1007/BF02651164
  * [18] Raymond JP. Optimal control of partial differential equations. Université Paul Sabatier, Internet 2013\.
  * [19] Vitale G, Preziosi L, Ambrosi D. Force traction microscopy: An inverse problem with pointwise observations. Journal of Mathematical Analysis and Applications 2012; 395(2): 788 - 801.  doi: https://doi.org/10.1016/j.jmaa.2012.05.074
  * [20] Huang CH, Chen CW. A boundary element-based inverse-problem in estimating transient boundary conditions with conjugate gradient method. International Journal for Numerical Methods in Engineering 1998; 42(5): 943-965.  doi: https://doi.org/10.1002/(SICI)1097-0207(19980715)42:5¡943::AID-NME395¿3.0.CO;2-V
  * [21] Tikhonov AN. Solution of incorrectly formulated problems and the regularization method. Soviet Math. Dokl. 1963; 4: 1035–1038.
  * [22] Beck JV. Surface heat flux determination using an integral method. Nuclear Engineering and Design 1968; 7(2): 170 - 178.
  * [23] Beck JV. Nonlinear estimation applied to the nonlinear inverse heat conduction problem. International Journal of Heat and Mass Transfer 1970; 13(4): 703-716.
  * [24] Chen CJ, Chiou JS. Prediction of surface temperature and heat flux from an interior temperature response. Letters in Heat and Mass Transfer 1976; 3(6): 539 - 548.
  * [25] Pinheiro C, Samarasekera I, Brimacomb J, Walker B. Mould heat transfer and continuously cast billet quality with mould flux lubrication Part 1 Mould heat transfer. Ironmaking & Steelmaking 2000; 27(1): 37-54.  doi: 10.1179/030192300677363
  * [26] Samarasekera I, Brimacombe J. The influence of mold behavior on the production of continuously cast steel billets. Metallurgical Transactions B 1982; 13(1): 105–116.
  * [27] Michelic S, Rauter W, Erker M, Brandl W, Bernhard C. Heat Transfer in a Round CC Mould: Measurement, Modelling and Validation. In: Steelmaking & Plastic Deformation Study Groups. Associazione Italiana di Metallurgia; 2008; Italy: Paper–30.
  * [28] Ranut P, Persi C, Nobile E, Spagnul S. Estimation of Heat Flux Distribution in a Continuous Casting Mould by Inverse Heat Transfer Algorithms. In: . Volume 2: 31st Computers and Information in Engineering Conference, Parts A and B. The American Society of Mechanical Engineers. ; 2011: 389-398
  * [29] Nelder JA, Mead R. A simplex method for function minimization. The computer journal 1965; 7(4): 308–313.
  * [30] Man Y, Hebi Y, Dacheng F. Real-time Analysis on Non-uniform Heat Transfer and Solidification in Mould of Continuous Casting Round Billets. Isij International - ISIJ INT 2004; 44: 1696-1704.  doi: 10.2355/isijinternational.44.1696
  * [31] Hebi Y, Man Y, Dacheng F. 3-D Inverse Problem Continuous Model for Thermal Behavior of Mould Process Based on the Temperature Measurements in Plant Trial. Isij International - ISIJ INT 2006; 46: 539-545.  doi: 10.2355/isijinternational.46.539
  * [32] Gonzalez M, Goldschmit M, Assanelli A, Dvorkin E, Berdaguer E. Modeling of the solidification process in a continuous casting installation for steel slabs. Metallurgical and Materials Transactions B 2003; 34: 455-473.  doi: 10.1007/s11663-003-0072-3
  * [33] Wang X, Kong L, Du F, et al. Mathematical Modeling of Thermal Resistances of Mold Flux and Air Gap in Continuous Casting Mold Based on an Inverse Problem. ISIJ International 2016; 56: 803-811.  doi: 10.2355/isijinternational.ISIJINT-2015-601
  * [34] Zhang H, Wang W. Mold Simulator Study of Heat Transfer Phenomenon During the Initial Solidification in Continuous Casting Mold. Metallurgical and Materials Transactions B 2017; 48: 779-793.  doi: 10.1007/s11663-016-0901-9
  * [35] Hu P, Wang X, Wei J, Yao M, Guo Q. Investigation of Liquid/Solid Slag and Air Gap Behavior inside the Mold during Continuous Slab Casting. ISIJ International 2018; 58(5): 892-898.  doi: 10.2355/isijinternational.ISIJINT-2017-393
  * [36] Tang L, Yao M, Wang X, Zhang X. Non-uniform thermal behavior and shell growth within mould for wide and thick slab continuous casting. Steel Research International 2012; 83(12): 1203-1213.  doi: 10.1002/srin.201200075
  * [37] Videcoq E, Quemener O, Lazard M, Neveu A. Heat source identification and on-line temperature control by a Branch Eigenmodes Reduced Model. International Journal of Heat and Mass Transfer 2008; 51(19): 4743-4752.  doi: https://doi.org/10.1016/j.ijheatmasstransfer.2008.02.029
  * [38] Videcoq E, Neveu A, Quemener O, Girault M, Petit D. Comparison of Two Nonlinear Model Reduction Techniques: The Modal Identification Method and the Branch Eigenmodes Reduction Method. Numerical Heat Transfer, Part B: Fundamentals 2006; 49(6): 537-558.  doi: 10.1080/10407790500344035
  * [39] Girault M, Videcoq E, Petit D. Estimation of time-varying heat sources through inversion of a low order model built with the Modal Identification Method from in-situ temperature measurements. International Journal of Heat and Mass Transfer 2010; 53(1): 206-219.  doi: https://doi.org/10.1016/j.ijheatmasstransfer.2009.09.040
  * [40] Girault M, Petit D. Identification methods in nonlinear heat conduction. Part II: inverse problem using a reduced model. International Journal of Heat and Mass Transfer 2005; 48(1): 119-133.  doi: https://doi.org/10.1016/j.ijheatmasstransfer.2004.06.033
  * [41] Aguado JV, Huerta A, Chinesta F, Cueto E. Real-time monitoring of thermal processes by reduced-order modeling. International Journal for Numerical Methods in Engineering 2015; 102(5): 991-1017.  doi: https://doi.org/10.1002/nme.4784
  * [42] Wang X, Yao M. Neural networks for solving the inverse heat transfer problem of continuous casting mould. In: . 2. IEEE Circuits and Systems Society. ; 2011: 791-794.
  * [43] Chen H, Su L, Wang G, Shibin W, Zhang L, Luo Z. Fuzzy estimation for heat flux distribution at the slab continuous casting mold surface. International Journal of Thermal Sciences 2014; 83: 80-88.  doi: 10.1016/j.ijthermalsci.2014.04.012
  * [44] Buhmann MD. Radial basis functions: theory and implementations. 12.  Cambridge university press .  2003\.
  * [45] Prando G. Non-Parametric Bayesian Methods for Linear System Identification. PhD thesis. Universitá di Padova, Via 8 Febbraio 1848, 2, 35122 Padova PD, Italy; 2016.
  * [46] Engl H. Regularization of inverse problems.  Dordrecht Boston: Kluwer Academic Publishers .  1996\.
  * [47] Hansen P. Discrete inverse problems : insight and algorithms.  Philadelphia: Society for Industrial and Applied Mathematics .  2010\.
  * [48] Aster R. Parameter estimation and inverse problems.  Amsterdam, Netherlands: Elsevier .  2019\.
  * [49] Kirsch A. Regularization by Discretization: 63–119; New York, NY: Springer New York .  2011
  * [50] Stabile G, Hijazi S, Mola A, Lorenzi S, Rozza G. POD-Galerkin reduced order methods for CFD using Finite Volume Discretisation: vortex shedding around a circular cylinder. Communications in Applied and Industrial Mathematics (2017); 8(1): 210-236.  doi: 10.1515/caim-2017-0011
  * [51] ITHACA-FV . https://mathlab.sissa.it/ithaca-fv; .  Accessed: 2020-10-26.
  * [52] Moukalled F, Mangani L, Darwish M. The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM and Matlab.  Springer Publishing Company, Incorporated.  1st ed. 2015.
  * [53] Olsson DM, Nelson LS. The Nelder-Mead Simplex Procedure for Function Minimization. Technometrics 1975; 17(1): 45-51.  doi: 10.1080/00401706.1975.10489269
  * [54] Strazzullo M, Ballarin F, Rozza G. POD-Galerkin model order reduction for parametrized nonlinear time-dependent optimal flow control: an application to shallow water equations. Journal of Numerical Mathematics 2022; 30(1): 63–84.  doi: doi:10.1515/jnma-2020-0098
  * [55] Negri F, Rozza G, Manzoni A, Quarteroni A. Reduced Basis Method for Parametrized Elliptic Optimal Control Problems. SIAM Journal on Scientific Computing 2013; 35(5): A2316-A2340.  doi: 10.1137/120894737
  * [56] Stabile G, Rosic B. Bayesian identification of a projection-based reduced order model for computational fluid dynamics. Computers & Fluids 2020; 201: 104477.  doi: 10.1016/j.compfluid.2020.104477