text:
\section{Non-autonomous Linear-Quadratic Regulator}\label{sec:LQR}
In this section, we formulate the Riccati-feedback approach for the Stefan problem.
We focus on the non-autonomous character of this problem, which is induced
by the moving interface and the consequently changing sub-domains.
This results in time-dependent coefficients in \cref{eq:sss}.

In order to derive a feedback-based stabilization of the Stefan problem,
we use the LQR approach (e.g.,~\cite{Son98}).
We use this approach because it is well studied for related types of problems, e.g., convection diffusion equations~\cite{Wei16}, and demonstrates promising performance for these.

To formulate the control problem,
we define a quadratic cost functional $J(\output\ensuremath{^h},\ensuremath{\control_\Delta}\ensuremath{^h})$,
tracking the deviation of the output
from the desired output as well as penalizing the control costs with
a weight factor $ 0 <\ensuremath{\lambda} \in \mathbb{R}$. The cost functional is thus defined as
\begin{equation}\label{eq:cost}
  J(\output\ensuremath{^h},\ensuremath{\control_\Delta}\ensuremath{^h}) = \frac{1}{2}\int_0^{\ensuremath{\ttime_\text{end}}} \nrm{\output\ensuremath{^h} - \output\ensuremath{^h}_{d}}^2 + \ensuremath{\lambda} \nrm{\ensuremath{\control_\Delta}\ensuremath{^h}}^2\;\diff t.
\end{equation}
We minimize this cost functional subject to
the linear time-varying system~\labelcref{eq:sss}.
With the matrices from \cref{eq:StefanMats} the LQR problem reads
\begin{equation}\begin{aligned}\label{eq:LQR}
	&\min\limits_{\ensuremath{\control_\Delta}\ensuremath{^h}} J(\output\ensuremath{^h}, \ensuremath{\control_\Delta}\ensuremath{^h}) \\
	\text{subject to} &\\
	&\ensuremath{\mathcal{M}}\ensuremath{\dot{\state}_\Delta}\ensuremath{^h} 	= \ensuremath{\mathcal{A}}\ensuremath{\state_\Delta}\ensuremath{^h} + \ensuremath{\hat{\cB}}\ensuremath{\control_\Delta}\ensuremath{^h}, \\
  	&\phantom{\ensuremath{\mathcal{M}}\ensuremath{\dot{\state}_\Delta}\ensuremath{^h}}\makebox[0pt][r]{\output\ensuremath{^h}} 			= \ensuremath{\mathcal{C}}\ensuremath{\state_\Delta}\ensuremath{^h}.
\end{aligned}\end{equation}
The unique solution to the LQR problem \labelcref{eq:LQR} is (see, e.g., \cite{Meh91})
\begin{equation}\label{eq:u}
  \ensuremath{\control_\Delta}\ensuremath{^h} = -\ensuremath{\mathcal{K}} \ensuremath{\state_\Delta}\ensuremath{^h},
\end{equation}
where the feedback gain matrix
\begin{equation}\label{eq:K}
  \ensuremath{\mathcal{K}} = \frac{1}{\ensuremath{\lambda}}\ensuremath{\hat{\cB}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}} = \frac{1}{\sqrt{\ensuremath{\lambda}}}\ensuremath{\mathcal{B}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}
\end{equation}
requires the solution $\ensuremath{\mathbf{X}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{n}}$
of a differential Riccati equation (DRE).
For the Stefan problem, this is the large-scale matrix-valued non-autonomous generalized DRE
\begin{equation}\label{eq:DREpre}
      - \ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}}{(\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}})} = \ensuremath{\mathcal{C}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathcal{C}}+\ensuremath{\mathcal{A}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}+\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}}\ensuremath{\mathcal{A}}
      -\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{B}}\cB\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}.
\end{equation}
All coefficients of \cref{eq:DREpre} can be non-autonomous (but we skip the $(\ensuremath{t})$-dependency for better readability).
The coefficients of the DRE, at each time instance \ensuremath{t}{}, are the matrices
$\ensuremath{\mathcal{A}}\time, \ensuremath{\mathcal{M}}\time$, $\dot\ensuremath{\mathcal{M}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{n}}$, $\ensuremath{\mathcal{B}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{m}}$,
and $\ensuremath{\mathcal{C}}\time \in \mathbb{R}^{\ensuremath{r} \times \ensuremath{n}}$ from \cref{eq:StefanMats},
where the input matrix $\ensuremath{\mathcal{B}}\time = \frac{1}{\sqrt{\ensuremath{\lambda}}}\ensuremath{\hat{\cB}}\time$ is scaled.
In order to solve the DRE, the time-derivative requires special treatment due to the time-dependent mass matrix \ensuremath{\mathcal{M}}\time{}.
The left-hand side in \cref{eq:DREpre} can be computed applying the chain rule,
\begin{equation*}
 	- \ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}}{(\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}})} = - \dot\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}} - \ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \dot\ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}} - \ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \dot\ensuremath{\mathcal{M}}.
\end{equation*}
We subtract the two terms containing $\dot\ensuremath{\mathcal{M}}\time$ from \cref{eq:DREpre} to obtain a DRE
with the time-derivative of \ensuremath{\mathcal{M}}\time{} moved to the right-hand side:
\begin{equation}\label{eq:DRE}
      - \ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \dot\ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}} = \ensuremath{\mathcal{C}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathcal{C}}+(\dot\ensuremath{\mathcal{M}} + \ensuremath{\mathcal{A}})\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}+\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}}(\dot\ensuremath{\mathcal{M}} + \ensuremath{\mathcal{A}})
      -\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{B}}\cB\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}.
\end{equation}
Besides the challenges arising from solving a large-scale matrix-valued DRE, the time-dependent matrices, especially the presence
of $\dot\ensuremath{\mathcal{M}}\time$, impose additional difficulties, which we address in this section.

To numerically compute the feedback gain matrix $\ensuremath{\mathcal{K}}\ensuremath{_\kiter} = \ensuremath{\mathcal{K}}\ensuremath{(\tk)}$
and the feedback control $\ensuremath{u}\ensuremath{_\kiter}=\ensuremath{\control_\Delta}\ensuremath{^h}\ensuremath{(\tk)}$ for $\ensuremath{\ttime\kiterind} \in \ensuremath{\discreteset_{\text{bwd}}}$,
the solution of the DRE \labelcref{eq:DRE} is required.
We use efficient low-rank methods for the computation of the numerical solution of the DRE,
$\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} = \ensuremath{\mathbf{X}}\ensuremath{(\tk)}$, since we assume that
\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}{} has a low (numerical) rank motivated by~\cite{Sti18}.
For this, we use only a small number of inputs and outputs in our experiments,
i.e.\ $\ensuremath{m},\ensuremath{r} \ll \ensuremath{n}$.
As introduced in~\cite{LanMS15}, which is motivated by~\cite{BenLT09}, \ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}{} can, thus, be approximated to high accuracy by the decomposition
\begin{equation*}
  \ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} \approx L\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} D\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} L\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}},
\end{equation*}
where the matrices $L\ensuremath{_\kiter} \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{s}}$ and $D\ensuremath{_\kiter} \in \mathbb{R}^{\ensuremath{s} \times \ensuremath{s}}$
have rank $\ensuremath{s} \ll \ensuremath{n}$.
We use the low-rank non-autonomous backward differentiation formula (BDF) (see \cref{algo:BDF}) to solve the non-autonomous
DRE \labelcref{eq:DRE}. This method is described in greater detail
in~\cite[Section 2]{BarBSS21a}.
\input{bdf_algo/BDF_algo.tex}
Compared to an open-loop control problem, the DRE replaces the adjoint equations in the LQR setting.
Consequently, we solve the DRE backwards in time and have the time-steps \ensuremath{\discreteset_{\text{bwd}}}{} (see \cref{eq:timestepsbwd}) in reversed order as an input to \cref{algo:BDF}.
In each step of the BDF method, we solve an algebraic Riccati equation (ARE)
\begin{equation}\begin{aligned}\label{eq:ARE}
      0 =\
      &\ensuremath{\mathcal{C}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}\ensuremath{\mathcal{S}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{C}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}
      +(\dot\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} + \ensuremath{\mathcal{A}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind})\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}
      +\ensuremath{\mathcal{M}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}(\dot\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} + \ensuremath{\mathcal{A}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind})\\
      &-\ensuremath{\mathcal{M}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{B}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{B}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}
\end{aligned}\end{equation}
for the low-rank solution
$\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} \approx L\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} D\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} L\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}$.
The matrices in \cref{eq:ARE} are constructed in \cref{line:BDF3,line:BDF4,line:BDF5,line:BDF6,line:BDF7}.
The coefficients $\alpha_j$ and $\beta$ are taken from the literature,
e.g.~\cite{AscP98}.
The BDF method is a multistep method.
Thus, to start a BDF method of order \(\ensuremath{\wp} > 1\), the terminal values $\ensuremath{\mathbf{X}}_0,\ldots,\ensuremath{\mathbf{X}}_{\ensuremath{\wp}-1}$ are required with
sufficient accuracy to obtain the desired order of convergence.
These values can be computed with sufficiently small time-steps of
the order \(\ensuremath{\wp}-1\) method.
This procedure is repeated recursively for the order \(\ensuremath{\wp}-1\) method to compute terminal values for its start.
For details see~\cite[Algorithm 2]{BarBSS21a}, or the implementation in~\cite[mess\_bdf\_dre.m]{SaaKB21-mmess-2.1}.
This method requires an additional input parameter that we set to $\ensuremath{n_{\ord}} = 10$ as in~\cite{BarBSS21a}.

\cref{algo:BDF} is embedded in the open source software package \mbox{\textsf{M-}\mess{}}{}~2.1~\cite{SaaKB21-mmess-2.1}, where it
benefits from efficient solvers for the ARE \labelcref{eq:ARE}.

To have an overview of the three steps for the Riccati-feedback stabilization of the Stefan problem,
we collect the overall procedure in \cref{algo:LQR}.
\IncMargin{1em}
\begin{algorithm}[t]
	\SetAlgoNoLine
	\caption{LQR for the Stefan problem}
	\label{algo:LQR}
	\BlankLine
	solve open-loop control problem to get the reference trajectory \newline
	and \ensuremath{\mathcal{A}}\time, \ensuremath{\mathcal{M}}\time, $\dot\ensuremath{\mathcal{M}}\time$, \ensuremath{\mathcal{B}}\time, \ensuremath{\mathcal{C}}\time{} of the linearized problem\label{line:LQR1}\;
	solve DRE \labelcref{eq:DRE} with \cref{algo:BDF} to get $\ensuremath{\mathcal{K}}\ensuremath{_\kiter}, \ensuremath{k} = 1,\ldots,\ensuremath{n_t}$\label{line:LQR2}\;
	apply \ensuremath{\mathcal{K}}\ensuremath{_\kiter}{} in a forward simulation of the Stefan problem\label{line:LQR3}\;
\end{algorithm}
\DecMargin{1em}
The method we choose for the forward simulations of the Stefan problem in \cref{line:LQR1,line:LQR3}
is a fractional-step-theta scheme and for solving the DRE in \cref{line:LQR2} the non-autonomous BDF methods, which we described above.
The time discretization \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}}{}, defined in \cref{eq:timesteps}, is used in all three steps.
In \cref{line:LQR2}, it is used backwards in time (\ensuremath{\discreteset_{\text{bwd}}}, \cref{eq:timestepsbwd}) and,
for the second forward solve in \cref{line:LQR3}, additional time steps can be added adaptively.
We use the time-adaptivity to prevent numerical instabilities
that can occur with feedback controls that have very large variation, as demonstrated in~\cite{BarBSS21a}.
The specific time-adaptivity is tailored to include the relative change of the feedback control.
A detailed description of this time-adaptive fractional-step-theta scheme with a relative control-based indicator function
is presented in~\cite[Section 3]{BarBSS21a}.

We provide the references for the codes of the methods above in \cref{sec:code}
and showcase the behavior of \cref{algo:BDF,algo:LQR} by means of several numerical experiments.

\section{Conclusions}\label{sec:conclusions}
In this work, we address the challenging task to derive, compute, and apply a feedback control for the two-dimensional two-phase Stefan problem.
The particular challenges for the Stefan problem lie in the non-linearities, the DAE structure, as well as
the moving inner boundary.
Additionally, the applied feedback control approach results in a generalized DRE with time-dependent large-scale matrices.
We address this task with a sharp interface representation and mesh movement techniques.
Resulting from this coupling of the Stefan problem with mesh movement, we have a detailed description of the linearization, discretization and matrix assembly.
Regarding this, we particularly elaborate on how Dirichlet boundary conditions can be treated
and display the DAE structure in the resulting matrices.

To obtain a feedback control for the Stefan problem, we apply the LQR approach and treat large-scale non-autonomous DREs with a corresponding low-rank BDF method.
More specifically, we also include the time-dependent mass matrix and its derivative.
The feedback control resulting from this approach is applied in a
forward simulation of the closed-loop Stefan problem
where we handle the numerical difficulties that arise with
a time-adaptive fractional-step-theta scheme specifically adapted to this process.

Through several numerical experiments, we demonstrate how effectively our methods and the resulting feedback controls perform.
The performance of the feedback controls strongly depends on the choice of control parameters like the weight factor in the cost functional
and the selected inputs and outputs.
The outputs are particularly important since in our formulation of the Stefan problem, the interface position to be controlled,
is not explicitly available as an output.
Thus, the outputs need to indicate the interface deviation reliably.

Our methods include several stages of approximations, like the linearization, implicit index reduction,
and low-rank representation of the DRE solution in an iterative solver.
Like this, we are able to make the numerical solution of this large-scale problem feasible
while the feedback controls still perform as known from problem types, that are already well studied in this regard.

Future research could investigate other promising methods for the non-autonomous DRE like, e.g., splitting schemes.
Additionally, a quantitative error analysis is still to be done in future work.
\section{Discretization and Linearization}\label{sec:discretization}
In order to apply the LQR approach in \cref{sec:LQR},
we need to formulate the Stefan problem in a standard state-space format,
which is linear and semi-discretized in space.
Thus, we describe how to transform the coupled DAE system of \cref{eq:heat,eq:V}
into
\begin{equation}\label{eq:sss}\begin{aligned}
  \ensuremath{\mathcal{M}}\dot\ensuremath{x}\ensuremath{^h} 	&= \ensuremath{\mathcal{A}}\ensuremath{x}\ensuremath{^h} + \ensuremath{\hat{\cB}}\ensuremath{u}\ensuremath{^h}, \\
  \output\ensuremath{^h} 			&= \ensuremath{\mathcal{C}}\ensuremath{x}\ensuremath{^h},
\end{aligned}\end{equation}
in this section.
To generate the square matrices $\ensuremath{\mathcal{A}}\time, \ensuremath{\mathcal{M}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{n}}$,
the input matrix $\ensuremath{\hat{\cB}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{m}}$, and
output matrix $\ensuremath{\mathcal{C}}\time \in \mathbb{R}^{\ensuremath{r} \times \ensuremath{n}}$,
we spatially discretize and linearize the Stefan problem.
Further, we take special care on how the boundary conditions \labelcref{eq:heat4,eq:V2}
are treated in the definition of the matrices for \cref{eq:sss}
since they are of particular importance for our feedback control problem.

For the spatial discretization, we use the finite element method (FEM)
on a mesh of triangular cells $\ensuremath{\mathscr{Q}{\ensuremath{^h}}}\time = \{\ensuremath{Q}\time\}$ that changes over time, driven by the movement of the interface.
The interface \ensuremath{\Gamma_{\text{int}}}\time{} itself is represented explicitly and sharply
through facets that are aligned with \ensuremath{\Gamma_{\text{int}}}\time.
In order to track \ensuremath{\Gamma_{\text{int}}}\time{} with the mesh, we use the semi-discrete
extended interface movement \ensuremath{\Upsilon}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}{} to adapt the mesh inside
the whole domain, as described in~\cite{BaePS10,BaePS13}.
In this way, we prevent mesh tangling and too strong deformations.
See~\cite{BarBHetal18a} for a detailed description of our implementation.

In order to define the linearization of the Stefan problem, we require a trajectory generated
by applying an open-loop control approach to the nonlinear problem. For this, we take the open-loop control approach from~\cite{BarBHetal18a}.
This (desired) reference trajectory contains the semi-discrete reference solutions \ensuremath{\tilde \temp}\ensuremath{^h}\time, \ensuremath{\tilde{\V}}\ensuremath{^h}\ensuremath{(\nabla \tempRef\ensuremath{^h}\time)}, \ensuremath{{\Gamma}_{\text{int,ref}}(\ensuremath{t})}, and
the heat conductivities \ensuremath{\tilde \alphac}\time, depending on the reference trajectory.

Combining these two techniques, we derive the linearized, semi-discrete version of \cref{eq:heat1}
\begin{equation}\label{eq:heatRef}
  \dot\ensuremath{\varTheta}\ensuremath{^h} - \ensuremath{\Upsilon}\ensuremath{^h} \cdot \nabla \ensuremath{\tilde \temp}\ensuremath{^h} - \ensuremath{\tilde \alphac} \Delta \ensuremath{\varTheta}\ensuremath{^h} = 0, \qquad\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega},
\end{equation}
by using the semi-discrete states ($\ensuremath{\varTheta}\ensuremath{^h}\time, \ensuremath{\Upsilon}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}$) and, in particular,
by replacing \ensuremath{\varTheta}\ensuremath{^h}\time{} in the convection term by \ensuremath{\tilde \temp}\ensuremath{^h}\time{} and \ensuremath{\alpha}\time{}, which
depends on \ensuremath{\varTheta}\ensuremath{^h}\time{}, by the reference heat conductivity \ensuremath{\tilde \alphac}\time{}:

In order to get the matrices for \cref{eq:sss}, we pose the semi-discrete variational formulations
of \cref{eq:heatRef} together with its boundary conditions and \cref{eq:V}:
\begin{align*}
    0 &= \int\limits_{\ensuremath{\Omega}} \dot\ensuremath{\varTheta}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}x
	- \int\limits_{\ensuremath{\Omega}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \nabla \ensuremath{\tilde \temp}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}x
	+ \int\limits_{\ensuremath{\Omega}} \ensuremath{\alpha} \nabla\ensuremath{\varTheta}\ensuremath{^h}\cdot\nabla\ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}x
	- \int\limits_{\ensuremath{\Gamma_{\ensuremath{u}}}} \ensuremath{\kk_l} \ensuremath{u}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}s,\\
    0 &= - \int\limits_{\ensuremath{\Omega}} \nabla\ensuremath{\Upsilon}\ensuremath{^h}\cdot\nabla\ensuremath{\hat{\textsf{v}}\ensuremath{^h}}  \mathrm{d}x
	  + \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s
	  - \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \frac{1}{\ensuremath{\ell}}[\ensuremath{k}(\nabla \ensuremath{\varTheta}\ensuremath{^h})]_l^s \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s.
\end{align*}
We denote the semi-discrete test functions with \ensuremath{\textsf{v}\ensuremath{^h}}\time{} and \ensuremath{\hat{\textsf{v}}\ensuremath{^h}}\time{}
and reformulate these variational formulations into a matrix-based form:
\begin{align*}
    \dot\ensuremath{\varTheta}\ensuremath{^h}		&= \ensuremath{ A_{\temp \temp}}\ensuremath{\varTheta}\ensuremath{^h} + \ensuremath{ A_{\temp \V}}\ensuremath{\Upsilon}\ensuremath{^h} + \ensuremath{B_\temp}\ensuremath{u}\ensuremath{^h}, 	& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega}, \\
    0					&= \ensuremath{ A_{\V \temp}}\ensuremath{\varTheta}\ensuremath{^h} + \ensuremath{ A_{\V \V}}\ensuremath{\Upsilon}\ensuremath{^h}, 						& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega}.
\end{align*}
Here, the coefficient matrices are defined via the inner products
\begin{equation}\label{eq:blockmatrix}\begin{aligned}
    \langle \ensuremath{ M_{\temp}}\ensuremath{\varTheta}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle 	&= \int\limits_{\ensuremath{\Omega}} \ensuremath{\varTheta}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}}\, \mathrm{d}x, \\
    \langle \ensuremath{ A_{\temp \temp}}\ensuremath{\varTheta}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle 	&= - \int\limits_{\ensuremath{\Omega}} \ensuremath{\alpha} \nabla \ensuremath{\varTheta}\ensuremath{^h} \cdot \nabla \ensuremath{\textsf{v}\ensuremath{^h}}\, \mathrm{d}x, \\
    \langle\ensuremath{ A_{\temp \V}}\ensuremath{\Upsilon}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle	&= \int\limits_{\ensuremath{\Omega}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \nabla \ensuremath{\tilde \temp}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}}\, \mathrm{d}x, \\
    \langle \ensuremath{ A_{\V \V}}\ensuremath{\Upsilon}\ensuremath{^h}, \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \rangle	&= \int\limits_{\ensuremath{\Omega}} \nabla\ensuremath{\Upsilon}\ensuremath{^h}\cdot\nabla\ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}x
						- \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s \\
    \langle \ensuremath{ A_{\V \temp}}\ensuremath{\varTheta}\ensuremath{^h}, \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \rangle	&= \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \frac{1}{\ensuremath{\ell}}[\ensuremath{k}(\nabla \ensuremath{\varTheta}\ensuremath{^h})]_l^s \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s, \\
    \langle \ensuremath{B_\temp}\ensuremath{u}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle 	&= \int\limits_{\ensuremath{\Gamma_{\ensuremath{u}}}} \ensuremath{\kk_l} \ensuremath{u}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}s.
\end{aligned}\end{equation}
With these definitions, we can formulate the semi-discrete linearized Stefan problem in the format of \cref{eq:sss}:
\begin{equation}\label{eq:blocksss}\begin{aligned}
\begin{blockmatrix}
  \mblock(0,1)\ensuremath{ M_{\temp}}\vphantom{\Big(}(1.6,1)	\mblock(1.6,1)0(1.6,1)
  \mblock(0,0)0\vphantom{\Big(}(1.6,1)		\mblock(1.6,0)0(1.6,1)
\end{blockmatrix}
\ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}} \begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\varTheta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\Upsilon}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	\end{blockmatrix}
  &= \begin{blockmatrix}
  \mblock(0,1)\ensuremath{ A_{\temp \temp}}\vphantom{\Big(}(2,1)		\mblock(2,1){\ensuremath{ A_{\temp \V}}}(2,1)
  \mblock(0,0)\ensuremath{ A_{\V \temp}}\vphantom{\Big(}(2,1)		\mblock(2,0)\ensuremath{ A_{\V \V}}(2,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\varTheta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\Upsilon}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
+
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{B_\temp}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)0\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
\ensuremath{\mathbf{u}}\ensuremath{^h}, \\ \\
\ensuremath{\mathbf{y}}\ensuremath{^h}\ &=
\begin{blockmatrix}
  \mblock(1,4)\ensuremath{C_\temp}(1,1)	\mblock(2,4)0(1,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,1)\ensuremath{\varTheta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,0)\ensuremath{\Upsilon}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}.
\end{aligned}\end{equation}
Several choices are
plausible for the output matrix $\ensuremath{C_\temp}\time$.
We will introduce some in \cref{sec:experiments}.
With the zero-blocks in \cref{eq:blocksss}, the DAE structure is clearly visible.
Since, at every time instance \ensuremath{t}{}, \ensuremath{ A_{\V \V}}\time{} represents the Poisson operator with a Dirichlet
boundary part,
it is always non-singular.
Further, the mass matrix with respect to the temperature \ensuremath{ M_{\temp}}\time{} is symmetric positive definite
and, in particular, always non-singular as well.
Thus, \cref{eq:blocksss} is a DAE in semi-explicit form of differential index 1 (see e.g.~\cite{KunM06}) and we can
apply the implicit index-reduction techniques from~\cite{morFreRM08} for this type of DAEs.

Next, we discuss the treatment of the boundary conditions in the FEM matrices in \cref{eq:blockmatrix,eq:blocksss}.
This is of special importance for our objective to use these matrices to compute a
feedback control for the Stefan problem in \cref{sec:LQR}.
The subject of attention in our feedback control problem is the position of \ensuremath{\Gamma_{\text{int}}}\time{},
which is directly linked to the boundary conditions \labelcref{eq:heat4,eq:V2}
since they are defined on \ensuremath{\Gamma_{\text{int}}}\time{}.
\cref{eq:heat4,eq:V2} are also important for the coupling of \cref{eq:heat,eq:V}.
This coupling is of particular importance for the feedback control problem
since the interface movement is defined in the Stefan condition \labelcref{eq:V2}
and the control is applied in \cref{eq:heat2}.
One common approach to handle Dirichlet boundary conditions, such as \cref{eq:heat4,eq:V2}, is to
remove the rows and columns corresponding to the degrees of freedom (DOFs)
at the related boundary regions from the FEM matrices.
However, these DOFs are related to the coupling of \ensuremath{\varTheta}\ensuremath{^h}\time{} and \ensuremath{\Upsilon}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}, and,
thus, the coupling of the temperature and the interface movement.
By removing these DOFs, the matrices would lose important information, which is
necessary for the computation of a feedback control.
Thus, the corresponding DOFs on \ensuremath{\Gamma_{\text{int}}}\time{} are to be present in the FEM matrices.
In order to treat the boundary conditions \labelcref{eq:heat4,eq:V2} appropriately
for our feedback control problem
we, from now on, consider the difference states
\begin{align*}
  \ensuremath{\temp_\Delta}\ensuremath{^h} 	&= \ensuremath{\varTheta}\ensuremath{^h} - \ensuremath{\tilde \temp}\ensuremath{^h}, \\
  \ensuremath{\V_\Delta}\ensuremath{^h} 		&= \ensuremath{\Upsilon}\ensuremath{^h} - \ensuremath{\tilde{\V}}\ensuremath{^h}.
\end{align*}
We formulate our feedback control problem in terms of $\ensuremath{\temp_\Delta}\ensuremath{^h}\time$ and $\ensuremath{\V_\Delta}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}$.
As a result, the desired state, which we stabilize with a feedback control, is the all-zero-state.
The semi-discrete version of the Dirichlet condition in \cref{eq:heat4} in terms of the difference state $\ensuremath{\temp_\Delta}\ensuremath{^h}\time$ is
equivalent to
\begin{equation}\label{eq:BCequiv}
  0 = \ensuremath{\ensuremath{\varTheta}_{\Delta,\text{melt}}} - \ensuremath{\temp_\Delta}\ensuremath{^h}, \qquad\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{int}}}.
\end{equation}
Further, both terms $\ensuremath{\ensuremath{\varTheta}_{\Delta,\text{melt}}}\time = \ensuremath{\ensuremath{\varTheta}_\text{melt}}\time - \ensuremath{\tilde \ensuremath{\varTheta}_\text{melt}}\time = 0$ and $\restr{\ensuremath{\temp_\Delta}\ensuremath{^h}\time}{\ensuremath{\Gamma_{\text{int}}}\time} = 0$ in
\cref{eq:BCequiv} are constant.
Thus, also the time derivative equals zero, $\ensuremath{\dot{\temp}_\Delta}\ensuremath{^h}\time = 0$,
and we add this equation to \cref{eq:BCequiv}.
As a result, we can formulate the modified condition \labelcref{eq:BCDiff} to replace \cref{eq:heat4} in terms of the difference state
for the FEM matrices, yielding
\begin{equation}\label{eq:BCDiff}
  \ensuremath{\dot{\temp}_\Delta}\ensuremath{^h} = \ensuremath{\ensuremath{\varTheta}_{\Delta,\text{melt}}} - \ensuremath{\temp_\Delta}\ensuremath{^h}, \qquad\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{int}}}.
\end{equation}
To incorporate \cref{eq:BCDiff} into the matrices, we denote $\ensuremath{\temp_{\Delta,\Gamma}}\ensuremath{^h}\time = \restr{\ensuremath{\temp_\Delta}\ensuremath{^h}\time}{\ensuremath{\Gamma_{\text{int}}}\time}$ and assume we can split
\begin{align*}
\ensuremath{\temp_\Delta}\ensuremath{^h} =
\begin{blockmatrix}
	  \mblock(0,1)\ensuremath{\bar{\temp}_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,0)\ensuremath{\temp_{\Delta,\Gamma}}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}.
\end{align*}
In order to ensure the Dirichlet condition on \ensuremath{\Gamma_{\text{int}}}\time{}, we add \cref{eq:BCDiff} to the related block matrices, which then read
\begin{equation*}
  \ensuremath{ \tilde{M}_{\temp}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ M_{\temp}}\vphantom{\Big(}(1.2,1)	\mblock(1.2,2)0(1.2,1)
  \mblock(0,1)0\vphantom{\Big(}(1.2,1)		\mblock(1.2,1)I(1.2,1)
\end{blockmatrix},\
  \ensuremath{ \tilde{A}_{\temp \temp}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ A_{\temp \temp}}\vphantom{\Big(}(2,1)		\mblock(2,2)\ensuremath{ A_{\temp \tempGamma}}(2,1)
  \mblock(0,1)0\vphantom{\Big(}(2,1)		\mblock(2,1){-I}(2,1)
\end{blockmatrix},\
  \ensuremath{ \tilde{A}_{\temp \V}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ A_{\temp \V}}\vphantom{\Big(}(2,1)
  \mblock(0,1)0\vphantom{\Big(}(2,1)
\end{blockmatrix},\
  \ensuremath{ \tilde{A}_{\V \temp}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ A_{\V \temp}}\vphantom{\Big(}(2,1)		\mblock(2,2)\ensuremath{ A_{\V \tempGamma}}(2,1)
\end{blockmatrix}.
\end{equation*}
With this, the conditions in \cref{eq:heat2,eq:V2} can be incorporated explicitly into the definition of
\ensuremath{ A_{\V \temp}}\time{} and \ensuremath{B_\temp}\time.
Regarding the remaining boundary conditions, the Neumann condition in \cref{eq:heat5} does not need any extra attention, since it is incorporated
automatically.
The Dirichlet boundary conditions in \cref{eq:heat3,eq:V3,eq:V4} can be handled by the common approach
to remove the rows and columns corresponding to the DOFs at the related boundary regions from the matrices.

\cref{eq:blocksss} with \cref{eq:BCDiff} incorporated reads
\begin{equation}\label{eq:blocksssGamma}\begin{aligned}
\begin{blockmatrix}
  \mblock(0,1)\ensuremath{ \tilde{M}_{\temp}}\vphantom{\Big(}(1.6,1)	\mblock(1.6,1)0(1.6,1)
  \mblock(0,0)0\vphantom{\Big(}(1.6,1)		\mblock(1.6,0)0(1.6,1)
\end{blockmatrix}
\ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}} \begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\temp_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\V_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	\end{blockmatrix}
  &= \begin{blockmatrix}
  \mblock(0,1)\ensuremath{ \tilde{A}_{\temp \temp}}\vphantom{\Big(}(2,1)		\mblock(2,1){\ensuremath{ \tilde{A}_{\temp \V}}}(2,1)
  \mblock(0,0)\ensuremath{ \tilde{A}_{\V \temp}}\vphantom{\Big(}(2,1)		\mblock(2,0)\ensuremath{ A_{\V \V}}(2,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\temp_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\V_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
+
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{B_\temp}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)0\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
\ensuremath{\mathbf{u}}\ensuremath{^h}, \\ \\
\ensuremath{\mathbf{y}}\ensuremath{^h}\ &=
\begin{blockmatrix}
  \mblock(1,4)\ensuremath{C_\temp}(1,1)	\mblock(2,4)0(1,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,1)\ensuremath{\temp_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,0)\ensuremath{\V_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}.
\end{aligned}\end{equation}
Notably, this modification preserves the DAE structure and index of \cref{eq:blocksss}.
Consequently, we follow~\cite{morFreRM08} and apply the Schur complement to remove the algebraic conditions.
This yields the new system matrices for the Stefan problem
\begin{equation}\label{eq:StefanMats}\begin{aligned}
  \ensuremath{\mathcal{M}}\time &= \ensuremath{ \tilde{M}_{\temp}}\time, \\
  \ensuremath{\mathcal{A}}\time &= \ensuremath{ \tilde{A}_{\temp \temp}}\time - \ensuremath{ \tilde{A}_{\temp \V}}\time\ensuremath{ A_{\V \V}}^{-1}\time\ensuremath{ \tilde{A}_{\V \temp}}\time, \\
  \ensuremath{\hat{\cB}}\time &= \ensuremath{B_\temp}\time, \\
  \ensuremath{\mathcal{C}}\time &= \ensuremath{C_\temp}\time.
\end{aligned}\end{equation}
With the matrices from \cref{eq:StefanMats}, we can transform the Stefan problem into the formulation of \cref{eq:sss}
and are able to apply the LQR approach for the computation of a feedback control.
In order to have a computationally efficient method, we keep the sparse structure of the matrices and
never compute the generally dense Schur complement for \ensuremath{\mathcal{A}}\time{} explicitly.
Instead we apply \ensuremath{\mathcal{A}}\time{} implicitly and use the matrices from \cref{eq:blocksssGamma}.
For details, see~\cite{morFreRM08}, or the numerical implementation in~\cite{SaaKB21-mmess-2.1}.

In \cref{sec:LQR}, we describe the related numerical methods.
For this, we further discretize the Stefan problem in time.
For the simulation forward in time, we use the reference time-steps
\begin{equation}\label{eq:timesteps}\begin{aligned}
  0 &= \ensuremath{t}_0 < \ensuremath{t}_1 < \ldots < \ensuremath{t}_{\ensuremath{n_t} - 1} < \ensuremath{t}_{\ensuremath{n_t}} = \ensuremath{\ttime_\text{end}}, \\
  \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}} &= \{\ensuremath{t}_0, \ensuremath{t}_1, \ldots,\ensuremath{t}_{\ensuremath{n_t} - 1}, \ensuremath{t}_{\ensuremath{n_t}}\}.
\end{aligned}\end{equation}
In \cref{{sec:LQR}}, we solve differential Riccati equations backwards in time.
For this we use the same time-steps $\hat \ensuremath{t}_{k} = \ensuremath{t}_{\ensuremath{n_t} - k}$ in reversed order
\begin{equation}\label{eq:timestepsbwd}\begin{aligned}
  \ensuremath{\ttime_\text{end}} &= \hat \ensuremath{t}_0 > \hat \ensuremath{t}_1 > \ldots > \hat \ensuremath{t}_{\ensuremath{n_t} - 1} > \hat \ensuremath{t}_{\ensuremath{n_t}} = \ensuremath{t}_0, \\
  \ensuremath{\discreteset_{\text{bwd}}} &= \{\hat \ensuremath{t}_0, \hat \ensuremath{t}_1, \ldots,\hat \ensuremath{t}_{\ensuremath{n_t} - 1}, \hat \ensuremath{t}_{\ensuremath{n_t}}\}.
\end{aligned}\end{equation}
We use the same time-steps because the matrices from \cref{eq:StefanMats},
which form the coefficients of the
differential Riccati equations, are
assembled during the forward simulation.
When we apply a time-adaptive method in the forward simulation, additional time-steps can be added to \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}}{}.
\section{Numerical Experiments}\label{sec:experiments}
With several numerical experiments, we demonstrate that the LQR approach is
applicable to this type of control problem, i.e., the Stefan problem, regardless
of several numerical challenges and approximations that arise
similar to other types of problems, where this approach is applied successfully.
The Stefan problem is non-linear, which we address with a linearization.
Further, the coefficients of the corresponding DRE are large-scale matrices.
Thus, we use a low-rank representation of the DRE solution to make it feasible
in regard of the computational cost and the memory requirements.
Additionally, the Stefan problem is a DAE as well as, e.g., the Navier-Stokes equations,
for which the LQR approach is applied successfully in \cite{Wei16}.
To handle the DAE structure, we use an implicit index reduction technique.

In addition,
for the Stefan problem, the matrices are time-dependent and
we approximate the time-derivative of the mass matrix by centered differences.
Consequently, we use the non-autonomous BDF method to solve the DRE and compute
the feedback gain matrices.
Then, at intermediate time-steps that are
generated by the time-adaptive fractional-step-theta scheme,
we interpolate the feedback gain matrices.
An alternative is to interpolate the DRE solution with the corresponding low-rank factors and compute the feedback gain matrices at intermediate time-steps as in \cref{line:BDF8} of \cref{algo:BDF}.
However, this is more expensive and not necessary in our setting.
All these techniques can introduce additional approximation errors.
Still, the computed feedback controls can steer the interface back to
the desired trajectory successfully in our experiments, similar to previously studied problem types.

With several numerical experiments, we illustrate the
performance of the feedback stabilization, which is computed with \cref{algo:LQR}.
The open-loop computations in \cref{line:LQR1} are reported in~\cite{BarBHetal18a} and
the assessment of the runtime performance of different methods for the solution of the DRE in
\cref{line:LQR2} is left to~\cite{BarBSS21a}.
Instead, we focus on the behavior of the feedback stabilization in \cref{line:LQR3}.
This behavior is strongly influenced by the weight factor \ensuremath{\lambda}{} in the cost
functional as well as different choices of inputs and outputs, as the experiments demonstrate.
In order to
asses the robustness of our computed feedback controls, we are interested in uncertainties in the
cooling, i.e., we investigate different perturbations \ensuremath{\varphi}\time{} to the Dirichlet
boundary condition at \ensuremath{\Gamma_{\text{cool}}}\time{}.

\begin{figure}[t]
  \centering
  \vspace{0em}\includegraphics{tikz/externalize/domain_in_out.pdf
  \caption{Input, Perturbation (left), and Output Areas (right)}%
  \label{fig:inout}
\end{figure}
We depict several boundary regions for perturbations (\ensuremath{\Gamma_{\text{cool}}}\time{}) as well as
inputs (\ensuremath{\Gamma_{\ensuremath{u}}}\time{}) and outputs (\Gammaoutx\time{}) in \cref{fig:inout}.
As for \ensuremath{\Omega}\time{}, the boundary regions for inputs $\Gammainx[1]\time, \ldots, \Gammainx[4]\time$
are constant in this representation of the domain.
However, we do not restrict ourselves to this case and, thus, keep the time-dependence in our notation.

The particular domain \ensuremath{\Omega}\time{}, we choose for the experiments, is a rectangle $[0, 0.5] \times [0, 1]$
with the initial interface position at height $0.5$ as a horizontal line.
For the discretization, we choose a mesh of triangles with \(3\,899\) vertices and 401 time-steps.
This, using standard P1 elements in \textsf{FEniCS}, results in \(3\,899\) DOFs for \ensuremath{\varTheta}\time{} and \(7\,798\) DOFs for \ensuremath{\V\ofGradTemp}.
After the removal of the DOFs corresponding to Dirichlet boundary conditions (see \cref{sec:discretization}),
the size of the matrices is \(\ensuremath{n}{} = 11\,429\).
The model parameters are listed in \cref{tab:para}.
\begin{table}[t]
\centering
\begin{tabular}{|l|l|l|l|l|l|l|l|}
\hline
\ensuremath{\ttime_\text{end}}	&	\ensuremath{\kk_s}	&	\ensuremath{\kk_l}		&	\ensuremath{\ensuremath{\varTheta}_\text{cool}}	&	\ensuremath{\ensuremath{\varTheta}_\text{melt}}	&	\ensuremath{\ell}	&	\ensuremath{\ensuremath{\varTheta}_0} 									&	\ensuremath{\timestep\kiterind} \\
\hline
$1$			&	$6$	&	$10$	&	$-1$	&	$0$		&	$10$	&	$\ensuremath{\ensuremath{\varTheta}_\text{cool}} \cdot (1 - 2 \cdot \ensuremath{\state_2})$	&	 $\in \left[10^{-4}, 2.5\cdot10^{-3}\right]$\\
\hline
\end{tabular}
\caption{Stefan Problem Model Parameters}
\label{tab:para}
\end{table}

\subsection*{Code Availability}\label{sec:code}
All codes and data to reproduce the presented results are available at~\cite{LQR_Stefan_codes}.
The non-autonomous BDF method is incorporated in the software package \mbox{\textsf{M-}\mess{}}{}~2.1~\cite{SaaKB21-mmess-2.1}

\subsection{Experiment 1}\label{sec:ex1}
For this first experiment, we focus on the influence of weight parameters and outputs
on the performance of the feedback control.
The desired interface trajectory, which we intend to stabilize with the feedback control,
is a flat horizontal line moving from its
initial height at \(0.5\) downward by \(0.004\).
Since we assume the Stefan problem to be asymptotically stable,
the interface returns to the desired trajectory after a perturbation without a feedback control as well.
However, this is not achieved within the time horizon.
With the feedback control, we intend to prevent the interface deviation to a certain extend and to steer the interface back in a much shorter time after a perturbation occurs.

The interface deviation is caused by the perturbation \ensuremath{\varphi}\time{}, which is  generated randomly in the form of three scalar values in the range of $[-\ensuremath{\ensuremath{\varTheta}_\text{cool}}, \ensuremath{\ensuremath{\varTheta}_\text{cool}}]$.
These are applied for a period of four time-steps of \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}}{} at the times $0.1$, $0.3$, and $0.5$ on $\ensuremath{\Gamma_{\text{cool}}}\time =
\Gammacoolx[1]\time\cup\Gammacoolx[2]\time$:
\begin{equation*}
  \ensuremath{\varTheta}	= \ensuremath{\ensuremath{\varTheta}_\text{cool}} + \ensuremath{\varphi}, \qquad \text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{cool}}}.
\end{equation*}
The trajectory of \ensuremath{\varphi}\time{} is pictured in the top part of \cref{fig:weights}.
The scaling of the y-axis is relative to \ensuremath{\ensuremath{\varTheta}_\text{cool}}.

Here, for the LQR problem, we use a single input \ensuremath{\control_\ensuremath{\mathcal{K}}}\time{} with the same value on
$\Gammainx[1]\time\cup\ldots\cup\Gammainx[4]\time$, i.e.\ $\ensuremath{\mathcal{B}}\time \in \mathbb{R}^{\ensuremath{n} \times 1}$.
The top part of \cref{fig:weights} shows the control for three different LQR
designs, resulting from three different combinations of weights and outputs:
\begin{align*}
  &\controlx[1]: &&\ensuremath{\lambda} = 10^{-4}, 		&&\text{2 outputs:}\ \Gammaoutx[3],\Gammaoutx[4],\\
  &\controlx[2]: &&\ensuremath{\lambda} = 10^{-6}, 		&&\text{2 outputs:}\ \Gammaoutx[3],\Gammaoutx[4],\\
  &\controlx[3]: &&\ensuremath{\lambda} = 1.6\cdot10^{-2}, 	&&\text{7 outputs:}\ \Gammaoutx[1],\ldots,\Gammaoutx[6],\ensuremath{\Gamma_{\text{int}}}.
\end{align*}
The outputs on \Gammaoutx[1]\time{}, \Gammaoutx[2]\time{}, \Gammaoutx[5]\time{},  and \Gammaoutx[6]\time{}
measure averaged temperatures on the corresponding interval, while
\Gammaoutx[3]\time{} and \Gammaoutx[4]\time{} represent point measurements of the
temperature.
The output at \ensuremath{\Gamma_{\text{int}}}\time{} monitors the difference of the interface movement
to the desired movement
\begin{equation*}
  \langle \ensuremath{\mathcal{C_{\text{int}}}} , \ensuremath{\textsf{v}\ensuremath{^h}} \rangle
    = \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \Big(\frac{1}{\ensuremath{\ell}}[\ensuremath{\kk_s}(\nabla \ensuremath{\temp_\Delta})_s
    - \ensuremath{\kk_l}(\nabla \ensuremath{\temp_\Delta})_l]\Big) \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}s
\end{equation*}
and is defined via the jump term of the Stefan condition (\cref{eq:V2}).
This output alone is not suited for an effective control-design.
It would generate a large output, and thus an active feedback-response, only while the perturbation is actively
driving the interface away from the desired trajectory.
However, it can not detect a difference in the position of the interface.
Consequently, it would not steer the interface back but would keep it on a ``parallel trajectory".

Thus, in the first two LQR designs, the controls \controlx[1]\time{} and \controlx[2]\time{}
are based on two outputs that measure the temperature at the desired interface position on the boundary, such that $\ensuremath{\mathcal{C}}\time{} \in \mathbb{R}^{2 \times \ensuremath{n}}$, while
\controlx[3]\time{}, in the third setting, uses seven outputs, i.e.\ $\ensuremath{\mathcal{C}}\time{} \in \mathbb{R}^{7 \times \ensuremath{n}}$.
With this, we compute the feedback gain matrices according to \cref{sec:LQR} with \cref{algo:BDF} and set $\ensuremath{\wp} = 1$.
Then, we simulate the closed-loop system with these feedback gain matrices together with the perturbations.
\begin{figure}[tp]
  \centering
  \includegraphics{tikz/externalize/Stefan_feedback_adaptive.pdf
  \caption{Perturbation and feedback (top), relative interface position (middle),
  			and time-step sizes (bottom) for different weights and outputs
	    (\nameref{sec:ex1})}%
  \label{fig:weights}
\end{figure}
The resulting feedback controls \controlx[1]\time{}, \controlx[2]\time{},
and \controlx[3]\time{} are displayed in the top part of \cref{fig:weights} and
the interface positions in the middle part of \cref{fig:weights}.
The interface positions are relative to the desired interface position
at the point $\ensuremath{\ensuremath{x}^*}\time$ with the largest deviation on the interface:
\begin{align*}
  \ensuremath{\Gamma_{\text{int},\Delta}(\ensuremath{t},\statexstar)} &= \ensuremath{{\Gamma}_{\text{int,ref}}(\ensuremath{t}, \statexstar)} - \ensuremath{\Gamma_{\text{int}}}(\ensuremath{t}, \ensuremath{\ensuremath{x}^*}), \\
  \ensuremath{\ensuremath{x}^*} &= \argmax_{\ensuremath{\state_1} \in [0,0.5]}\ \left|\ensuremath{{\Gamma}_{\text{int,ref}}(\ensuremath{t}, \ensuremath{\state_1})} - \ensuremath{\Gamma_{\text{int}}}(\ensuremath{t}, \ensuremath{\state_1})\right|.
\end{align*}

The feedback control \controlx[1]\time{} is most active shortly after the
perturbation starts and it steers the interface back to the desired position in
much reduced time, as expected from theory.
Again as expected, a smaller weight factor \ensuremath{\lambda}{} allows the control
\controlx[2]\time{} to react with a larger magnitude input while the perturbation is active.
It is, therefore, able to stop the interface from deviating earlier and drives
it back even faster.

It should be mentioned that smaller weights reduce the contribution of the control cost term
in the cost functional \labelcref{eq:cost}.
This term can also be interpreted as a regularization term.
Consequently, smaller weights decrease the regularity of the feedback control problem.
We note that this results in larger computational cost for the solution of the DRE due to slower convergence
of the internal iterative solvers used in the BDF method, but refer to~\cite[Section~4.2]{BarBSS21a}
for a more detailed discussion of this issue.

The weights for the LQR settings of \controlx[1]\time{} and \controlx[2]\time{}
are not directly comparable to the setting for \controlx[3]\time, since this is based on different outputs.
The outputs at \Gammaoutx[5]\time{} and \Gammaoutx[6]\time{} allow it to detect the temperature perturbation earlier
and the output that monitors the movement of \ensuremath{\Gamma_{\text{int}}}\time{} can observe the deviation of the interface earlier.
Thus, \controlx[3]\time{} is most active immediately after the perturbation starts and moves the interface back much faster.

Note, that we use a time-adaptive fractional-step-theta scheme
with a relative feedback control-based indicator function
to simulate the Stefan problem together with a feedback control.
For this example, the bottom part of \cref{fig:weights}
shows the adaptive time-step sizes.
While the feedback control is inactive, like at the beginning or end of the simulation,
the method chooses the largest time-step size of $\ensuremath{\timestep\kiterind} = 0.0025$.
As soon as the feedback control is active, the method chooses the minimal time-step size of $\ensuremath{\timestep\kiterind} = 10^{-4}$.
This prevents numerical instabilities that can occur with very small \ensuremath{\lambda}{} and thus very active feedback controls
as is demonstrated in~\cite[Section~4.2]{BarBSS21a}.

\nameref{sec:ex1} demonstrates the influence of the chosen weight parameter \ensuremath{\lambda}{} in the cost functional (\cref{eq:cost}).
Smaller weight parameters increase the impact of the deviation term in the cost functional indirectly by
decreasing the relevance of the control costs which are measured by the second term and scaled by \ensuremath{\lambda}.
Further, the choice and number of outputs has a significant impact on the performance of the feedback control.

The next experiment demonstrates that our method works as well with increased time-dependence of the data
in the non-autonomous DRE.

\subsection{Experiment 2}\label{sec:ex2}
In our second experiment, the desired interface trajectory is, again, a flat horizontal line moving from its
initial height at \(0.5\), this time, upward by \(0.1\).
This is a longer distance interface movement than in \nameref{sec:ex1} on the same time horizon and, in turn,
results in a stronger time-dependence of the matrices for the DRE.\@
The investigated LQR design is chosen as in the second case in \nameref{sec:ex1}.
\begin{figure}[t]
  \includegraphics{tikz/externalize/plot_moveup.pdf
  \caption{Perturbation and feedback (top) and relative interface position (bottom) for an interface moving upwards
	    (\nameref{sec:ex2})}%
  \label{fig:up}
\end{figure}

That means, we use the same weight ($\ensuremath{\lambda} = 10^{-6}$) and outputs as for \controlx[2]\time{}
above and call  the corresponding feedback control \controltilde[2]\time.
This time we use four randomly generated perturbations applied at the times $0.1, 0.3, 0.5$, and $0.7$, which, together with
\controltilde[2]\time, can be found in the top part of \cref{fig:up}.
Again, the perturbations cause the interface to deviate from the desired trajectory.
The feedback control \controltilde[2]\time{} behaves similarly to the previous experiment.
It stops the interface from deviating and drives it back to the desired position,
as expected from the theory.
Additionally, the time-adaptivity behaves analogously to \nameref{sec:ex1}.

The experiment showcases the performance of our method with strongly time-varying coefficients in the DRE.
We implemented more experiments with desired interface trajectories that move
considerably.
Our feedback control showed the same performance for all of them. This demonstrates
that our solver can also cope with the stronger time dependence without difficulties.

\subsection{Experiment 3}\label{sec:ex3}
Our third experiment revisits the basic task as in \nameref{sec:ex1}.
This time, we generate two different perturbations $\pertx[1]\time,
\pertx[2]\time \in [-\ensuremath{\ensuremath{\varTheta}_\text{cool}}, \ensuremath{\ensuremath{\varTheta}_\text{cool}}]$, acting on the left (\Gammacoolx[1]\time{}) and right (\Gammacoolx[2]\time{}) part
of the Dirichlet boundary at the bottom of \ensuremath{\Omega}\time{} (see \cref{fig:inout}),
with three random values each, which are applied at the times $0.1, 0.4325$, and $0.765$.
They are displayed in the top part of \cref{fig:curv}.
We apply these two perturbation functions to
\begin{align*}
  \ensuremath{\varTheta}	&= \ensuremath{\ensuremath{\varTheta}_\text{cool}} + \pertx[1], \qquad \text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \Gammacoolx[1],\\
  \ensuremath{\varTheta}	&= \ensuremath{\ensuremath{\varTheta}_\text{cool}} + \pertx[2], \qquad \text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \Gammacoolx[2].
\end{align*}
These perturbations will not only move the interface away from the desired position
but also add some curvature to the deviated interface.
To stabilize this interface position, we use six individual inputs at $\Gammainx[1]\time,\ldots,\Gammainx[6]\time$
($\ensuremath{\mathcal{B}}\time{} \in \mathbb{R}^{\ensuremath{n} \times 6}$) and the same two outputs at \Gammaoutx[3]\time{}, \Gammaoutx[4]\time{}
as in the previous experiments ($\ensuremath{\mathcal{C}}\time{} \in \mathbb{R}^{2 \times \ensuremath{n}}$).
For the weight factor in the cost functional we use $\ensuremath{\lambda} = 10^{-9}$.
Again, this weight factor is not directly comparable to the previous experiments
since we use a different combination of inputs and outputs.

\begin{figure}[tp]
  \includegraphics{tikz/externalize/plot_curvature.pdf
  \caption{Two different perturbations (top), and perturbed and feedback controlled relative interface positions (bottom)
	    (\nameref{sec:ex3})}%
  \label{fig:curv}
\end{figure}
The resulting perturbed and feedback controlled interfaces
are displayed in the bottom part of \cref{fig:curv} at the
time points $\ensuremath{t} \in [0.465, 0.56, 0.798, 1]$.
With the first two perturbations, the interface is pushed downwards and assumes a distinct curvature.
Right after the second perturbation ($\ensuremath{t} = 0.465$), the feedback controlled interface is already essentially
back to the desired interface, but still has undesired curvature.
Well before the third perturbation, at $\ensuremath{t} = 0.56$, the feedback controlled interface
is almost flat again, as the desired interface.
At the two endpoints of the interface, where the outputs are located,
the controlled interface is directly at the desired position.
The remaining curvature can not be measured by the cost functional and is, thus, not
actively removed.
For the given setup, this is the expected behavior of the feedback control.
It is similar after the first perturbation.
The third perturbation moves the interface upwards above the desired position ($\ensuremath{t} = 0.798$).
At the end of the time interval ($\ensuremath{t} = 1$), the feedback controlled interface is again almost flat and back
to the desired position at the ends while the uncontrolled interface is still far away.

This experiment shows that the feedback control can move the perturbed interface back to the desired trajectory
and additionally control the curvature of the interface, to a certain extent.
With the chosen input and output setting, the position of the interface is corrected as fast as in the previous experiments.
The curvature correction requires some more time but is still performed during the considered time interval
up to the level that the cost functional can measure with the given two outputs.

\medskip
In all three experiments, we are able to stabilize the desired interface position
with the proposed feedback control approach
despite the different challenges that are posed by the Stefan problem.

\section{Introduction}
The solidification and melting of materials is an active and intensively studied field with numerous applications.
This phase-change problem can be modeled by a non-linear PDE and is often called Stefan problem after J.~Stefan
who describes it in his works~\cite{Ste89a,Ste89b,Ste90}.
In a certain domain, the temperature of the material is either below, above, or equal to
the specific melting temperature of the material.
Accordingly, the domain is split into a solid and a liquid phase, which are separated by an interface or inner boundary.
J.~Stefan formulates what is now called Stefan condition in~\cite{Ste89a}, which couples the time-derivative of the interface position
with the jump of the temperature gradient along the interface.
On the other hand, J.~Stefan was not the first to consider this type of problem.
G.~Lam{\'e} and B.~P.~Clapeyron were concerned with this in an earlier work as well~\cite{LamC31},
such that it is also called Lam{\'e}--Clapeyron problem.

Several books address the Stefan problem, e.g.,~\cite{Rub71,NieCM11,Gup18,KogK20a}.
An extensive historical survey of the Stefan problem can be found in the book by L.~I.~Ruben\v{s}te\u{\i}n~\cite[Introduction: \S 1]{Rub71}.
Early works on the Stefan problem usually consider the one-dimensional case, while higher dimensional cases where only studied decades after the original publication.

The aim of this manuscript is to compute and apply feedback control to the two-dimensional two-phase Stefan problem.
While this problem has been studied in combination with open-loop controls several times, see
e.g.~\cite{Zie08, Ber10, AntNS14, AntNS15, BarBHetal18a} and the reference therein, only
recently, also closed-loop, i.e. feedback control, for the Stefan problem has been discussed in~\cite{KogK19,KogK20,KogK20a}.
However, these works address only the one-dimensional case.
On the one hand, the novelty of our work lies in the consideration of feedback control for the two-phase Stefan problem in two spatial dimensions
in contrast to the one-phase or one-dimensional Stefan problem.
On the other hand, we apply the linear-quadratic regulator (LQR) approach to the Stefan problem, which goes beyond the types of problems
that have been studied in connection with this approach.

In particular, the application of LQR requires the treatment of the non-linearities, the differential-algebraic nature,
and the time-varying character of the Stefan problem.
In this manuscript, we address the details on how to transform it into a linear ordinary differential equation,
apply the LQR approach to compute feedback controls, and then use these feedback controls
to stabilize the interface position in the non-linear differential-algebraic problem.
In order to do this, we linearize the Stefan problem around a reference trajectory and assemble a generalized differential Riccati equation (DRE).
For the derivation of the LQR problem and the resulting DRE we refer to~\cite{Rei72, BitG91}.
An extensive study on DREs and their numerous applications can be found in~\cite{AboFIetal03}.
Since the Stefan condition is an algebraic equation, which is coupled to the Stefan problem,
the DRE we are considering results from a differential-algebraic equation (DAE).
Details on generalized DREs can be found in~\cite{KunM90} and for DREs resulting from DAEs, see~\cite{KunM89}.

To solve the Stefan problem numerically and assemble the DRE, we discretize in space using the
finite element method (FEM) by applying the software \textsf{FEniCS}~\cite{dolfin}.
This results in a large-scale matrix-valued DRE with time-dependent coefficients, denoted as a non-autonomous DRE.\
While we observe in numerical experiments that, usually, its solution can be well
approximated by a low-rank factorization, this is proven theoretically
only for the autonomous DRE in~\cite{Sti18}.
Well known methods, which use the low-rank structure of the numerical solution,
are splitting schemes~\cite{Sti15,Sti15a,Sti18a,OstPW18,MenOPetal18},
Rosenbrock and Peer methods~\cite{Men12,LanMS15,Lan17,BenL18}
as well as the backward differentiation formulas (BDF)~\cite{BenM04, Men12, LanMS15, BenM18}.
Krylov subspace methods~\cite{BehBH21,KosM20,KirSi19,GueHJetal18} and exponential integrators~\cite{LiZL20}
for DREs have been developed recently as well.
In~\cite{Men12, BenM18} Rosenbrock and BDF methods and in~\cite{LanMS15,Lan17} also Peer methods are studied for a non-autonomous DRE
where the mass matrix is constant.
In extension of this, splitting schemes and BDF methods are developed in~\cite{BarBSS21a} for non-autonomous DREs with a time-varying mass matrix.
In contrast to the BDF methods, the splitting schemes require that
the coefficients can be decomposed into a time-dependent scalar function times a constant matrix.
The non-autonomous DREs we consider go beyond this case.
We use mesh movement techniques to track the interface, and, as a consequence, the single matrix entries change very differently, possibly for all matrices.
Thus, the BDF methods are
the most promising method for the non-autonomous DRE resulting from the Stefan problem and we apply the non-autonomous BDF method from~\cite{BarBSS21a}.

\paragraph{Structure of the Manuscript}
In this manuscript, we derive and apply several feedback controls to the two-dimensional two-phase Stefan problem.
For this, we state the equations that model the Stefan problem and the mesh movement that we use to track the interface position (\cref{sec:stefan}).
Then, we linearize the resulting system and semi-discretize it in space (\cref{sec:discretization}).
With the resulting time-dependent matrices we formulate a non-autonomous DRE resulting from the LQR approach (\cref{sec:LQR}).
We solve this DRE with the non-autonomous BDF method and use the resulting feedback gain matrices to compute a feedback control (also \cref{sec:LQR}).
We apply the described methods in several numerical experiments for different parameter settings and specify where all our codes
and data are available (\cref{sec:experiments}).

\paragraph{Notation}
In most equations, we omit the time-dependence \time{}, the spatial dependence $(\ensuremath{x})$, or the combination of both $(\ensuremath{t}, \ensuremath{x})$.
This is supposed to improve the readability of the equations.
\section{Non-autonomous Linear-Quadratic Regulator}\label{sec:LQR}
In this section, we formulate the Riccati-feedback approach for the Stefan problem.
We focus on the non-autonomous character of this problem, which is induced
by the moving interface and the consequently changing sub-domains.
This results in time-dependent coefficients in \cref{eq:sss}.

In order to derive a feedback-based stabilization of the Stefan problem,
we use the LQR approach (e.g.,~\cite{Son98}).
We use this approach because it is well studied for related types of problems, e.g., convection diffusion equations~\cite{Wei16}, and demonstrates promising performance for these.

To formulate the control problem,
we define a quadratic cost functional $J(\output\ensuremath{^h},\ensuremath{\control_\Delta}\ensuremath{^h})$,
tracking the deviation of the output
from the desired output as well as penalizing the control costs with
a weight factor $ 0 <\ensuremath{\lambda} \in \mathbb{R}$. The cost functional is thus defined as
\begin{equation}\label{eq:cost}
  J(\output\ensuremath{^h},\ensuremath{\control_\Delta}\ensuremath{^h}) = \frac{1}{2}\int_0^{\ensuremath{\ttime_\text{end}}} \nrm{\output\ensuremath{^h} - \output\ensuremath{^h}_{d}}^2 + \ensuremath{\lambda} \nrm{\ensuremath{\control_\Delta}\ensuremath{^h}}^2\;\diff t.
\end{equation}
We minimize this cost functional subject to
the linear time-varying system~\labelcref{eq:sss}.
With the matrices from \cref{eq:StefanMats} the LQR problem reads
\begin{equation}\begin{aligned}\label{eq:LQR}
	&\min\limits_{\ensuremath{\control_\Delta}\ensuremath{^h}} J(\output\ensuremath{^h}, \ensuremath{\control_\Delta}\ensuremath{^h}) \\
	\text{subject to} &\\
	&\ensuremath{\mathcal{M}}\ensuremath{\dot{\state}_\Delta}\ensuremath{^h} 	= \ensuremath{\mathcal{A}}\ensuremath{\state_\Delta}\ensuremath{^h} + \ensuremath{\hat{\cB}}\ensuremath{\control_\Delta}\ensuremath{^h}, \\
  	&\phantom{\ensuremath{\mathcal{M}}\ensuremath{\dot{\state}_\Delta}\ensuremath{^h}}\makebox[0pt][r]{\output\ensuremath{^h}} 			= \ensuremath{\mathcal{C}}\ensuremath{\state_\Delta}\ensuremath{^h}.
\end{aligned}\end{equation}
The unique solution to the LQR problem \labelcref{eq:LQR} is (see, e.g., \cite{Meh91})
\begin{equation}\label{eq:u}
  \ensuremath{\control_\Delta}\ensuremath{^h} = -\ensuremath{\mathcal{K}} \ensuremath{\state_\Delta}\ensuremath{^h},
\end{equation}
where the feedback gain matrix
\begin{equation}\label{eq:K}
  \ensuremath{\mathcal{K}} = \frac{1}{\ensuremath{\lambda}}\ensuremath{\hat{\cB}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}} = \frac{1}{\sqrt{\ensuremath{\lambda}}}\ensuremath{\mathcal{B}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}
\end{equation}
requires the solution $\ensuremath{\mathbf{X}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{n}}$
of a differential Riccati equation (DRE).
For the Stefan problem, this is the large-scale matrix-valued non-autonomous generalized DRE
\begin{equation}\label{eq:DREpre}
      - \ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}}{(\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}})} = \ensuremath{\mathcal{C}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathcal{C}}+\ensuremath{\mathcal{A}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}+\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}}\ensuremath{\mathcal{A}}
      -\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{B}}\cB\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}.
\end{equation}
All coefficients of \cref{eq:DREpre} can be non-autonomous (but we skip the $(\ensuremath{t})$-dependency for better readability).
The coefficients of the DRE, at each time instance \ensuremath{t}{}, are the matrices
$\ensuremath{\mathcal{A}}\time, \ensuremath{\mathcal{M}}\time$, $\dot\ensuremath{\mathcal{M}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{n}}$, $\ensuremath{\mathcal{B}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{m}}$,
and $\ensuremath{\mathcal{C}}\time \in \mathbb{R}^{\ensuremath{r} \times \ensuremath{n}}$ from \cref{eq:StefanMats},
where the input matrix $\ensuremath{\mathcal{B}}\time = \frac{1}{\sqrt{\ensuremath{\lambda}}}\ensuremath{\hat{\cB}}\time$ is scaled.
In order to solve the DRE, the time-derivative requires special treatment due to the time-dependent mass matrix \ensuremath{\mathcal{M}}\time{}.
The left-hand side in \cref{eq:DREpre} can be computed applying the chain rule,
\begin{equation*}
 	- \ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}}{(\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}})} = - \dot\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}} - \ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \dot\ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}} - \ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}} \dot\ensuremath{\mathcal{M}}.
\end{equation*}
We subtract the two terms containing $\dot\ensuremath{\mathcal{M}}\time$ from \cref{eq:DREpre} to obtain a DRE
with the time-derivative of \ensuremath{\mathcal{M}}\time{} moved to the right-hand side:
\begin{equation}\label{eq:DRE}
      - \ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \dot\ensuremath{\mathbf{X}} \ensuremath{\mathcal{M}} = \ensuremath{\mathcal{C}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathcal{C}}+(\dot\ensuremath{\mathcal{M}} + \ensuremath{\mathcal{A}})\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}+\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}}(\dot\ensuremath{\mathcal{M}} + \ensuremath{\mathcal{A}})
      -\ensuremath{\mathcal{M}}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{B}}\cB\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{\mathcal{M}}.
\end{equation}
Besides the challenges arising from solving a large-scale matrix-valued DRE, the time-dependent matrices, especially the presence
of $\dot\ensuremath{\mathcal{M}}\time$, impose additional difficulties, which we address in this section.

To numerically compute the feedback gain matrix $\ensuremath{\mathcal{K}}\ensuremath{_\kiter} = \ensuremath{\mathcal{K}}\ensuremath{(\tk)}$
and the feedback control $\ensuremath{u}\ensuremath{_\kiter}=\ensuremath{\control_\Delta}\ensuremath{^h}\ensuremath{(\tk)}$ for $\ensuremath{\ttime\kiterind} \in \ensuremath{\discreteset_{\text{bwd}}}$,
the solution of the DRE \labelcref{eq:DRE} is required.
We use efficient low-rank methods for the computation of the numerical solution of the DRE,
$\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} = \ensuremath{\mathbf{X}}\ensuremath{(\tk)}$, since we assume that
\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}{} has a low (numerical) rank motivated by~\cite{Sti18}.
For this, we use only a small number of inputs and outputs in our experiments,
i.e.\ $\ensuremath{m},\ensuremath{r} \ll \ensuremath{n}$.
As introduced in~\cite{LanMS15}, which is motivated by~\cite{BenLT09}, \ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}{} can, thus, be approximated to high accuracy by the decomposition
\begin{equation*}
  \ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} \approx L\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} D\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} L\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}},
\end{equation*}
where the matrices $L\ensuremath{_\kiter} \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{s}}$ and $D\ensuremath{_\kiter} \in \mathbb{R}^{\ensuremath{s} \times \ensuremath{s}}$
have rank $\ensuremath{s} \ll \ensuremath{n}$.
We use the low-rank non-autonomous backward differentiation formula (BDF) (see \cref{algo:BDF}) to solve the non-autonomous
DRE \labelcref{eq:DRE}. This method is described in greater detail
in~\cite[Section 2]{BarBSS21a}.
\input{bdf_algo/BDF_algo.tex}
Compared to an open-loop control problem, the DRE replaces the adjoint equations in the LQR setting.
Consequently, we solve the DRE backwards in time and have the time-steps \ensuremath{\discreteset_{\text{bwd}}}{} (see \cref{eq:timestepsbwd}) in reversed order as an input to \cref{algo:BDF}.
In each step of the BDF method, we solve an algebraic Riccati equation (ARE)
\begin{equation}\begin{aligned}\label{eq:ARE}
      0 =\
      &\ensuremath{\mathcal{C}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}\ensuremath{\mathcal{S}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{C}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}
      +(\dot\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} + \ensuremath{\mathcal{A}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind})\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}
      +\ensuremath{\mathcal{M}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}} \ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}(\dot\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} + \ensuremath{\mathcal{A}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind})\\
      &-\ensuremath{\mathcal{M}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{B}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{B}}\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}\ensuremath{\mathcal{M}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind}
\end{aligned}\end{equation}
for the low-rank solution
$\ensuremath{\mathbf{X}}\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} \approx L\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} D\ensuremath{^{\phantom{\ensuremath{^{\textsf{T}}}}}\kiterind} L\ensuremath{_\kiter}\ensuremath{^{\textsf{T}}}$.
The matrices in \cref{eq:ARE} are constructed in \cref{line:BDF3,line:BDF4,line:BDF5,line:BDF6,line:BDF7}.
The coefficients $\alpha_j$ and $\beta$ are taken from the literature,
e.g.~\cite{AscP98}.
The BDF method is a multistep method.
Thus, to start a BDF method of order \(\ensuremath{\wp} > 1\), the terminal values $\ensuremath{\mathbf{X}}_0,\ldots,\ensuremath{\mathbf{X}}_{\ensuremath{\wp}-1}$ are required with
sufficient accuracy to obtain the desired order of convergence.
These values can be computed with sufficiently small time-steps of
the order \(\ensuremath{\wp}-1\) method.
This procedure is repeated recursively for the order \(\ensuremath{\wp}-1\) method to compute terminal values for its start.
For details see~\cite[Algorithm 2]{BarBSS21a}, or the implementation in~\cite[mess\_bdf\_dre.m]{SaaKB21-mmess-2.1}.
This method requires an additional input parameter that we set to $\ensuremath{n_{\ord}} = 10$ as in~\cite{BarBSS21a}.

\cref{algo:BDF} is embedded in the open source software package \mbox{\textsf{M-}\mess{}}{}~2.1~\cite{SaaKB21-mmess-2.1}, where it
benefits from efficient solvers for the ARE \labelcref{eq:ARE}.

To have an overview of the three steps for the Riccati-feedback stabilization of the Stefan problem,
we collect the overall procedure in \cref{algo:LQR}.
\IncMargin{1em}
\begin{algorithm}[t]
	\SetAlgoNoLine
	\caption{LQR for the Stefan problem}
	\label{algo:LQR}
	\BlankLine
	solve open-loop control problem to get the reference trajectory \newline
	and \ensuremath{\mathcal{A}}\time, \ensuremath{\mathcal{M}}\time, $\dot\ensuremath{\mathcal{M}}\time$, \ensuremath{\mathcal{B}}\time, \ensuremath{\mathcal{C}}\time{} of the linearized problem\label{line:LQR1}\;
	solve DRE \labelcref{eq:DRE} with \cref{algo:BDF} to get $\ensuremath{\mathcal{K}}\ensuremath{_\kiter}, \ensuremath{k} = 1,\ldots,\ensuremath{n_t}$\label{line:LQR2}\;
	apply \ensuremath{\mathcal{K}}\ensuremath{_\kiter}{} in a forward simulation of the Stefan problem\label{line:LQR3}\;
\end{algorithm}
\DecMargin{1em}
The method we choose for the forward simulations of the Stefan problem in \cref{line:LQR1,line:LQR3}
is a fractional-step-theta scheme and for solving the DRE in \cref{line:LQR2} the non-autonomous BDF methods, which we described above.
The time discretization \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}}{}, defined in \cref{eq:timesteps}, is used in all three steps.
In \cref{line:LQR2}, it is used backwards in time (\ensuremath{\discreteset_{\text{bwd}}}, \cref{eq:timestepsbwd}) and,
for the second forward solve in \cref{line:LQR3}, additional time steps can be added adaptively.
We use the time-adaptivity to prevent numerical instabilities
that can occur with feedback controls that have very large variation, as demonstrated in~\cite{BarBSS21a}.
The specific time-adaptivity is tailored to include the relative change of the feedback control.
A detailed description of this time-adaptive fractional-step-theta scheme with a relative control-based indicator function
is presented in~\cite[Section 3]{BarBSS21a}.

We provide the references for the codes of the methods above in \cref{sec:code}
and showcase the behavior of \cref{algo:BDF,algo:LQR} by means of several numerical experiments.

\section{Two-dimensional Two-phase Stefan Problem}\label{sec:stefan}
In this section, we define the equations describing the Stefan problem, its boundary conditions and initial values.
This involves equations characterizing the temperature and interface movement.
Our goal is to
use a Riccati feedback control approach
to control the interface position.
For this,
we choose a sharp interface representation in the formulation of the Stefan problem,
opposed to, e.g., the level set representation of the interface in \cite{Ber10}.
As in~\cite{BaePS10,BaePS13}, we extend the interface movement to the whole domain.
This ensures the mesh regularity for the semi-discretized Stefan problem that is described in \cref{sec:discretization}.

At each \(\ensuremath{t} \in [0, \ensuremath{\ttime_\text{end}}]\), the domain is $\ensuremath{\Omega}\time \subset
\mathbb{R}^2$.
One instance is illustrated in \cref{fig:domain}.
We split the domain \ensuremath{\Omega}\time{} into the two regions corresponding to the two phases.
These are the region where the material is in its solid phase \ensuremath{\domain_s}\time{} and accordingly, the region \ensuremath{\domain_l}\time{} related to the liquid phase.
The two phases are separated by the interface \ensuremath{\Gamma_{\text{int}}}\time.
This inner phase-boundary can move such that its position is time-dependent.
Thus, also the two phases \ensuremath{\domain_s}\time{} and \ensuremath{\domain_l}\time{} are time-dependent
and, as a consequence, so is the whole domain \ensuremath{\Omega}\time{} and its boundary regions.
The boundary of \ensuremath{\Omega}\time{} is separated into \ensuremath{\Gamma_{\ensuremath{u}}}\time{}, \ensuremath{\Gamma_{\text{cool}}}\time{} and \ensuremath{\Gamma_\text{N}}\time{}
as depicted in \cref{fig:domain}.
Note that the outer shape of \ensuremath{\Omega}\time{} is constant for the realization chosen in this manuscript.
Thus, the time-dependence of \ensuremath{\Omega}\time{} is not absolutely necessary even though
its sub-domains are time-dependent.
However, we keep the time-dependence in our notation in order to not restrict our methods
to this case.

\includetikz{domain}

We follow the definition of the Stefan problem from~\cite{BarBHetal18a}. However, here we use
it in a more compact form, i.e.\ omitting the couplings with the Navier-Stokes equations and the interface graph formulation.
While the Navier-Stokes equations alone add additional algebraic constraints, making the DAE harder to classify,
both the Navier-Stokes equations and the interface graph formulation add more nonlinearities to the problem.
In order to develop the general numerical strategy,
we first want to study the feedback control problem for this simplified setting without these couplings,
before delving into the additional technical challenges of the full problem formulation.

We denote the temperature as \ensuremath{\varTheta}\time{} and model it with the partial
differential \cref{eq:heat}:
\begin{subequations}\label{eq:heat}\begin{align}
    \label{eq:heat1}
    \dot\ensuremath{\varTheta} - \ensuremath{\Upsilon} \cdot \nabla \ensuremath{\varTheta} - \ensuremath{\alpha} \Delta \ensuremath{\varTheta}	&= 0,  				& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega}, \\
    \label{eq:heat2}
    \partial_{\ensuremath{\boldsymbol{n}}}\ensuremath{\varTheta}						&= \ensuremath{u},			& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\ensuremath{u}}}, \\
    \label{eq:heat3}
    \ensuremath{\varTheta}								&= \ensuremath{\ensuremath{\varTheta}_\text{cool}},			& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{cool}}}, \\
    \label{eq:heat4}
    \ensuremath{\varTheta}								&= \ensuremath{\ensuremath{\varTheta}_\text{melt}},			& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{int}}}, \\
    \label{eq:heat5}
    \partial_{\ensuremath{\boldsymbol{n}}} \ensuremath{\varTheta}						&= 0,				& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_\text{N}}, \\
    \label{eq:heat6}
    \ensuremath{\varTheta}(0)								&= \ensuremath{\ensuremath{\varTheta}_0},			& &\text{on}\ \ensuremath{\Omega}.
\end{align}\end{subequations}

In \cref{eq:heat2}, we apply the control \ensuremath{u}\time{} as a Neumann condition on the control boundary \ensuremath{\Gamma_{\ensuremath{u}}}\time.
The \cref{eq:heat3,eq:heat4} describe the Dirichlet conditions on the cooling boundary \ensuremath{\Gamma_{\text{cool}}}\time{} and the interface \ensuremath{\Gamma_{\text{int}}}\time{}
with the constants \ensuremath{\ensuremath{\varTheta}_\text{cool}}{} and \ensuremath{\ensuremath{\varTheta}_\text{melt}}{}, respectively.
\cref{eq:heat6} presents the initial condition with the initial temperature distribution \ensuremath{\ensuremath{\varTheta}_0}.
The heat conductivities in the solid phase \ensuremath{\kk_s}{} and in the liquid phase \ensuremath{\kk_l}{} are collected in \ensuremath{\alpha}\ensuremath{(\temp\time)}{}:
\begin{equation*}
    \ensuremath{\alpha} = \begin{cases}
		      \ensuremath{\kk_s}, & \text{on}\ \ensuremath{\domain_s},\\
		      \ensuremath{\kk_l}, & \text{on}\ \ensuremath{\domain_l}.\\
               \end{cases}
\end{equation*}
In \cref{eq:heat1}, the temperature is coupled with the extended interface movement \ensuremath{\V\ofGradTemp}{}.
For every $\ensuremath{t} \in (0, \ensuremath{\ttime_\text{end}}]$, we model \ensuremath{\V\ofGradTemp}{} with a system of algebraic equations,
in the sense that they do not contain any time-derivatives:

\begin{subequations}\label{eq:V}\begin{align}
    \label{eq:V1}
    \Delta \ensuremath{\Upsilon}		&= 0, 									& &\text{on}\ \ensuremath{\Omega}, \\
    \label{eq:V2}
    \ensuremath{\Upsilon} 			&= \Big(\frac{1}{\ensuremath{\ell}}[\ensuremath{k}(\nabla \ensuremath{\varTheta}\time)]_l^s\Big) \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}},	& &\text{on}\ \ensuremath{\Gamma_{\text{int}}},\\
    \label{eq:V3}
    \ensuremath{\Upsilon}			&= 0, 									& &\text{on}\ \ensuremath{\Gamma_{\text{cool}}} \cup \ensuremath{\Gamma_{\ensuremath{u}}}, \\
    \label{eq:V4}
    \ensuremath{\Upsilon} \cdot \ensuremath{\boldsymbol{n}}	&= 0, 									& &\text{on}\ \ensuremath{\Gamma_\text{N}}, \\
    \label{eq:V5}
    \ensuremath{\Upsilon}(0)			&= 0, 									& &\text{on}\ \ensuremath{\Omega},
\end{align}\end{subequations}

On the interface, $\restr{\ensuremath{\Upsilon}}{\ensuremath{\Gamma_{\text{int}}}}\ensuremath{(\nabla \temp\time)} = \ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}$ is the interface movement in normal direction,
where \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}}\time{} is the unit normal vector pointing from \ensuremath{\domain_s}\time{} to \ensuremath{\domain_l}\time.
\ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}{} is coupled to \ensuremath{\varTheta}\time{} through the Stefan condition~\eqref{eq:V2}.
Here, \ensuremath{\ell}{} is the latent heat constant and
\begin{equation}\label{stefan}
   [k(\nabla \ensuremath{\varTheta})]_l^s = \ensuremath{\kk_s}\partial_{\ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}}}\restr{\ensuremath{\varTheta}}{\ensuremath{\domain_s}} - \ensuremath{\kk_l}\partial_{-\ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}}}\restr{\ensuremath{\varTheta}}{\ensuremath{\domain_l}}
\end{equation}
is the jump of the temperature gradient along \ensuremath{\Gamma_{\text{int}}}\time.
\cref{eq:V3,eq:V4} ensure that the outer boundaries of \ensuremath{\Omega}\time{} do not move
such that the outer shape of the domain does not change.

We use \ensuremath{\V\ofGradTemp}{} in \cref{sec:discretization} for the mesh movement.
The Stefan problem, which is described by the system of \cref{eq:heat,eq:V},
is a non-linear system of DAEs on a
time-varying domain.
Due to the coupling of \ensuremath{\varTheta}\time{} and \ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}{} in~\eqref{eq:V2}
and the two temperature-dependent phases in the definition of \ensuremath{\alpha}\ensuremath{(\temp\time)}{},
both, \ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}{} and \ensuremath{\alpha}\ensuremath{(\temp\time)}{}, depend on \ensuremath{\varTheta}\time{}.
Thus, the terms $\ensuremath{\V\ofGradTemp} \cdot \nabla \ensuremath{\varTheta}\time$ and $\ensuremath{\alpha}\ensuremath{(\temp\time)} \Delta \ensuremath{\varTheta}\time$ in \cref{eq:heat1} are nonlinear.

\section{Conclusions}\label{sec:conclusions}
In this work, we address the challenging task to derive, compute, and apply a feedback control for the two-dimensional two-phase Stefan problem.
The particular challenges for the Stefan problem lie in the non-linearities, the DAE structure, as well as
the moving inner boundary.
Additionally, the applied feedback control approach results in a generalized DRE with time-dependent large-scale matrices.
We address this task with a sharp interface representation and mesh movement techniques.
Resulting from this coupling of the Stefan problem with mesh movement, we have a detailed description of the linearization, discretization and matrix assembly.
Regarding this, we particularly elaborate on how Dirichlet boundary conditions can be treated
and display the DAE structure in the resulting matrices.

To obtain a feedback control for the Stefan problem, we apply the LQR approach and treat large-scale non-autonomous DREs with a corresponding low-rank BDF method.
More specifically, we also include the time-dependent mass matrix and its derivative.
The feedback control resulting from this approach is applied in a
forward simulation of the closed-loop Stefan problem
where we handle the numerical difficulties that arise with
a time-adaptive fractional-step-theta scheme specifically adapted to this process.

Through several numerical experiments, we demonstrate how effectively our methods and the resulting feedback controls perform.
The performance of the feedback controls strongly depends on the choice of control parameters like the weight factor in the cost functional
and the selected inputs and outputs.
The outputs are particularly important since in our formulation of the Stefan problem, the interface position to be controlled,
is not explicitly available as an output.
Thus, the outputs need to indicate the interface deviation reliably.

Our methods include several stages of approximations, like the linearization, implicit index reduction,
and low-rank representation of the DRE solution in an iterative solver.
Like this, we are able to make the numerical solution of this large-scale problem feasible
while the feedback controls still perform as known from problem types, that are already well studied in this regard.

Future research could investigate other promising methods for the non-autonomous DRE like, e.g., splitting schemes.
Additionally, a quantitative error analysis is still to be done in future work.
\section{Discretization and Linearization}\label{sec:discretization}
In order to apply the LQR approach in \cref{sec:LQR},
we need to formulate the Stefan problem in a standard state-space format,
which is linear and semi-discretized in space.
Thus, we describe how to transform the coupled DAE system of \cref{eq:heat,eq:V}
into
\begin{equation}\label{eq:sss}\begin{aligned}
  \ensuremath{\mathcal{M}}\dot\ensuremath{x}\ensuremath{^h} 	&= \ensuremath{\mathcal{A}}\ensuremath{x}\ensuremath{^h} + \ensuremath{\hat{\cB}}\ensuremath{u}\ensuremath{^h}, \\
  \output\ensuremath{^h} 			&= \ensuremath{\mathcal{C}}\ensuremath{x}\ensuremath{^h},
\end{aligned}\end{equation}
in this section.
To generate the square matrices $\ensuremath{\mathcal{A}}\time, \ensuremath{\mathcal{M}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{n}}$,
the input matrix $\ensuremath{\hat{\cB}}\time \in \mathbb{R}^{\ensuremath{n} \times \ensuremath{m}}$, and
output matrix $\ensuremath{\mathcal{C}}\time \in \mathbb{R}^{\ensuremath{r} \times \ensuremath{n}}$,
we spatially discretize and linearize the Stefan problem.
Further, we take special care on how the boundary conditions \labelcref{eq:heat4,eq:V2}
are treated in the definition of the matrices for \cref{eq:sss}
since they are of particular importance for our feedback control problem.

For the spatial discretization, we use the finite element method (FEM)
on a mesh of triangular cells $\ensuremath{\mathscr{Q}{\ensuremath{^h}}}\time = \{\ensuremath{Q}\time\}$ that changes over time, driven by the movement of the interface.
The interface \ensuremath{\Gamma_{\text{int}}}\time{} itself is represented explicitly and sharply
through facets that are aligned with \ensuremath{\Gamma_{\text{int}}}\time.
In order to track \ensuremath{\Gamma_{\text{int}}}\time{} with the mesh, we use the semi-discrete
extended interface movement \ensuremath{\Upsilon}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}{} to adapt the mesh inside
the whole domain, as described in~\cite{BaePS10,BaePS13}.
In this way, we prevent mesh tangling and too strong deformations.
See~\cite{BarBHetal18a} for a detailed description of our implementation.

In order to define the linearization of the Stefan problem, we require a trajectory generated
by applying an open-loop control approach to the nonlinear problem. For this, we take the open-loop control approach from~\cite{BarBHetal18a}.
This (desired) reference trajectory contains the semi-discrete reference solutions \ensuremath{\tilde \temp}\ensuremath{^h}\time, \ensuremath{\tilde{\V}}\ensuremath{^h}\ensuremath{(\nabla \tempRef\ensuremath{^h}\time)}, \ensuremath{{\Gamma}_{\text{int,ref}}(\ensuremath{t})}, and
the heat conductivities \ensuremath{\tilde \alphac}\time, depending on the reference trajectory.

Combining these two techniques, we derive the linearized, semi-discrete version of \cref{eq:heat1}
\begin{equation}\label{eq:heatRef}
  \dot\ensuremath{\varTheta}\ensuremath{^h} - \ensuremath{\Upsilon}\ensuremath{^h} \cdot \nabla \ensuremath{\tilde \temp}\ensuremath{^h} - \ensuremath{\tilde \alphac} \Delta \ensuremath{\varTheta}\ensuremath{^h} = 0, \qquad\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega},
\end{equation}
by using the semi-discrete states ($\ensuremath{\varTheta}\ensuremath{^h}\time, \ensuremath{\Upsilon}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}$) and, in particular,
by replacing \ensuremath{\varTheta}\ensuremath{^h}\time{} in the convection term by \ensuremath{\tilde \temp}\ensuremath{^h}\time{} and \ensuremath{\alpha}\time{}, which
depends on \ensuremath{\varTheta}\ensuremath{^h}\time{}, by the reference heat conductivity \ensuremath{\tilde \alphac}\time{}:

In order to get the matrices for \cref{eq:sss}, we pose the semi-discrete variational formulations
of \cref{eq:heatRef} together with its boundary conditions and \cref{eq:V}:
\begin{align*}
    0 &= \int\limits_{\ensuremath{\Omega}} \dot\ensuremath{\varTheta}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}x
	- \int\limits_{\ensuremath{\Omega}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \nabla \ensuremath{\tilde \temp}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}x
	+ \int\limits_{\ensuremath{\Omega}} \ensuremath{\alpha} \nabla\ensuremath{\varTheta}\ensuremath{^h}\cdot\nabla\ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}x
	- \int\limits_{\ensuremath{\Gamma_{\ensuremath{u}}}} \ensuremath{\kk_l} \ensuremath{u}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}s,\\
    0 &= - \int\limits_{\ensuremath{\Omega}} \nabla\ensuremath{\Upsilon}\ensuremath{^h}\cdot\nabla\ensuremath{\hat{\textsf{v}}\ensuremath{^h}}  \mathrm{d}x
	  + \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s
	  - \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \frac{1}{\ensuremath{\ell}}[\ensuremath{k}(\nabla \ensuremath{\varTheta}\ensuremath{^h})]_l^s \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s.
\end{align*}
We denote the semi-discrete test functions with \ensuremath{\textsf{v}\ensuremath{^h}}\time{} and \ensuremath{\hat{\textsf{v}}\ensuremath{^h}}\time{}
and reformulate these variational formulations into a matrix-based form:
\begin{align*}
    \dot\ensuremath{\varTheta}\ensuremath{^h}		&= \ensuremath{ A_{\temp \temp}}\ensuremath{\varTheta}\ensuremath{^h} + \ensuremath{ A_{\temp \V}}\ensuremath{\Upsilon}\ensuremath{^h} + \ensuremath{B_\temp}\ensuremath{u}\ensuremath{^h}, 	& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega}, \\
    0					&= \ensuremath{ A_{\V \temp}}\ensuremath{\varTheta}\ensuremath{^h} + \ensuremath{ A_{\V \V}}\ensuremath{\Upsilon}\ensuremath{^h}, 						& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega}.
\end{align*}
Here, the coefficient matrices are defined via the inner products
\begin{equation}\label{eq:blockmatrix}\begin{aligned}
    \langle \ensuremath{ M_{\temp}}\ensuremath{\varTheta}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle 	&= \int\limits_{\ensuremath{\Omega}} \ensuremath{\varTheta}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}}\, \mathrm{d}x, \\
    \langle \ensuremath{ A_{\temp \temp}}\ensuremath{\varTheta}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle 	&= - \int\limits_{\ensuremath{\Omega}} \ensuremath{\alpha} \nabla \ensuremath{\varTheta}\ensuremath{^h} \cdot \nabla \ensuremath{\textsf{v}\ensuremath{^h}}\, \mathrm{d}x, \\
    \langle\ensuremath{ A_{\temp \V}}\ensuremath{\Upsilon}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle	&= \int\limits_{\ensuremath{\Omega}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \nabla \ensuremath{\tilde \temp}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}}\, \mathrm{d}x, \\
    \langle \ensuremath{ A_{\V \V}}\ensuremath{\Upsilon}\ensuremath{^h}, \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \rangle	&= \int\limits_{\ensuremath{\Omega}} \nabla\ensuremath{\Upsilon}\ensuremath{^h}\cdot\nabla\ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}x
						- \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \ensuremath{\Upsilon}\ensuremath{^h} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s \\
    \langle \ensuremath{ A_{\V \temp}}\ensuremath{\varTheta}\ensuremath{^h}, \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \rangle	&= \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \frac{1}{\ensuremath{\ell}}[\ensuremath{k}(\nabla \ensuremath{\varTheta}\ensuremath{^h})]_l^s \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}} \cdot \ensuremath{\hat{\textsf{v}}\ensuremath{^h}} \mathrm{d}s, \\
    \langle \ensuremath{B_\temp}\ensuremath{u}\ensuremath{^h}, \ensuremath{\textsf{v}\ensuremath{^h}} \rangle 	&= \int\limits_{\ensuremath{\Gamma_{\ensuremath{u}}}} \ensuremath{\kk_l} \ensuremath{u}\ensuremath{^h} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}s.
\end{aligned}\end{equation}
With these definitions, we can formulate the semi-discrete linearized Stefan problem in the format of \cref{eq:sss}:
\begin{equation}\label{eq:blocksss}\begin{aligned}
\begin{blockmatrix}
  \mblock(0,1)\ensuremath{ M_{\temp}}\vphantom{\Big(}(1.6,1)	\mblock(1.6,1)0(1.6,1)
  \mblock(0,0)0\vphantom{\Big(}(1.6,1)		\mblock(1.6,0)0(1.6,1)
\end{blockmatrix}
\ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}} \begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\varTheta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\Upsilon}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	\end{blockmatrix}
  &= \begin{blockmatrix}
  \mblock(0,1)\ensuremath{ A_{\temp \temp}}\vphantom{\Big(}(2,1)		\mblock(2,1){\ensuremath{ A_{\temp \V}}}(2,1)
  \mblock(0,0)\ensuremath{ A_{\V \temp}}\vphantom{\Big(}(2,1)		\mblock(2,0)\ensuremath{ A_{\V \V}}(2,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\varTheta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\Upsilon}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
+
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{B_\temp}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)0\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
\ensuremath{\mathbf{u}}\ensuremath{^h}, \\ \\
\ensuremath{\mathbf{y}}\ensuremath{^h}\ &=
\begin{blockmatrix}
  \mblock(1,4)\ensuremath{C_\temp}(1,1)	\mblock(2,4)0(1,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,1)\ensuremath{\varTheta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,0)\ensuremath{\Upsilon}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}.
\end{aligned}\end{equation}
Several choices are
plausible for the output matrix $\ensuremath{C_\temp}\time$.
We will introduce some in \cref{sec:experiments}.
With the zero-blocks in \cref{eq:blocksss}, the DAE structure is clearly visible.
Since, at every time instance \ensuremath{t}{}, \ensuremath{ A_{\V \V}}\time{} represents the Poisson operator with a Dirichlet
boundary part,
it is always non-singular.
Further, the mass matrix with respect to the temperature \ensuremath{ M_{\temp}}\time{} is symmetric positive definite
and, in particular, always non-singular as well.
Thus, \cref{eq:blocksss} is a DAE in semi-explicit form of differential index 1 (see e.g.~\cite{KunM06}) and we can
apply the implicit index-reduction techniques from~\cite{morFreRM08} for this type of DAEs.

Next, we discuss the treatment of the boundary conditions in the FEM matrices in \cref{eq:blockmatrix,eq:blocksss}.
This is of special importance for our objective to use these matrices to compute a
feedback control for the Stefan problem in \cref{sec:LQR}.
The subject of attention in our feedback control problem is the position of \ensuremath{\Gamma_{\text{int}}}\time{},
which is directly linked to the boundary conditions \labelcref{eq:heat4,eq:V2}
since they are defined on \ensuremath{\Gamma_{\text{int}}}\time{}.
\cref{eq:heat4,eq:V2} are also important for the coupling of \cref{eq:heat,eq:V}.
This coupling is of particular importance for the feedback control problem
since the interface movement is defined in the Stefan condition \labelcref{eq:V2}
and the control is applied in \cref{eq:heat2}.
One common approach to handle Dirichlet boundary conditions, such as \cref{eq:heat4,eq:V2}, is to
remove the rows and columns corresponding to the degrees of freedom (DOFs)
at the related boundary regions from the FEM matrices.
However, these DOFs are related to the coupling of \ensuremath{\varTheta}\ensuremath{^h}\time{} and \ensuremath{\Upsilon}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}, and,
thus, the coupling of the temperature and the interface movement.
By removing these DOFs, the matrices would lose important information, which is
necessary for the computation of a feedback control.
Thus, the corresponding DOFs on \ensuremath{\Gamma_{\text{int}}}\time{} are to be present in the FEM matrices.
In order to treat the boundary conditions \labelcref{eq:heat4,eq:V2} appropriately
for our feedback control problem
we, from now on, consider the difference states
\begin{align*}
  \ensuremath{\temp_\Delta}\ensuremath{^h} 	&= \ensuremath{\varTheta}\ensuremath{^h} - \ensuremath{\tilde \temp}\ensuremath{^h}, \\
  \ensuremath{\V_\Delta}\ensuremath{^h} 		&= \ensuremath{\Upsilon}\ensuremath{^h} - \ensuremath{\tilde{\V}}\ensuremath{^h}.
\end{align*}
We formulate our feedback control problem in terms of $\ensuremath{\temp_\Delta}\ensuremath{^h}\time$ and $\ensuremath{\V_\Delta}\ensuremath{^h}\ensuremath{(\nabla \temp\ensuremath{^h}\time)}$.
As a result, the desired state, which we stabilize with a feedback control, is the all-zero-state.
The semi-discrete version of the Dirichlet condition in \cref{eq:heat4} in terms of the difference state $\ensuremath{\temp_\Delta}\ensuremath{^h}\time$ is
equivalent to
\begin{equation}\label{eq:BCequiv}
  0 = \ensuremath{\ensuremath{\varTheta}_{\Delta,\text{melt}}} - \ensuremath{\temp_\Delta}\ensuremath{^h}, \qquad\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{int}}}.
\end{equation}
Further, both terms $\ensuremath{\ensuremath{\varTheta}_{\Delta,\text{melt}}}\time = \ensuremath{\ensuremath{\varTheta}_\text{melt}}\time - \ensuremath{\tilde \ensuremath{\varTheta}_\text{melt}}\time = 0$ and $\restr{\ensuremath{\temp_\Delta}\ensuremath{^h}\time}{\ensuremath{\Gamma_{\text{int}}}\time} = 0$ in
\cref{eq:BCequiv} are constant.
Thus, also the time derivative equals zero, $\ensuremath{\dot{\temp}_\Delta}\ensuremath{^h}\time = 0$,
and we add this equation to \cref{eq:BCequiv}.
As a result, we can formulate the modified condition \labelcref{eq:BCDiff} to replace \cref{eq:heat4} in terms of the difference state
for the FEM matrices, yielding
\begin{equation}\label{eq:BCDiff}
  \ensuremath{\dot{\temp}_\Delta}\ensuremath{^h} = \ensuremath{\ensuremath{\varTheta}_{\Delta,\text{melt}}} - \ensuremath{\temp_\Delta}\ensuremath{^h}, \qquad\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{int}}}.
\end{equation}
To incorporate \cref{eq:BCDiff} into the matrices, we denote $\ensuremath{\temp_{\Delta,\Gamma}}\ensuremath{^h}\time = \restr{\ensuremath{\temp_\Delta}\ensuremath{^h}\time}{\ensuremath{\Gamma_{\text{int}}}\time}$ and assume we can split
\begin{align*}
\ensuremath{\temp_\Delta}\ensuremath{^h} =
\begin{blockmatrix}
	  \mblock(0,1)\ensuremath{\bar{\temp}_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,0)\ensuremath{\temp_{\Delta,\Gamma}}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}.
\end{align*}
In order to ensure the Dirichlet condition on \ensuremath{\Gamma_{\text{int}}}\time{}, we add \cref{eq:BCDiff} to the related block matrices, which then read
\begin{equation*}
  \ensuremath{ \tilde{M}_{\temp}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ M_{\temp}}\vphantom{\Big(}(1.2,1)	\mblock(1.2,2)0(1.2,1)
  \mblock(0,1)0\vphantom{\Big(}(1.2,1)		\mblock(1.2,1)I(1.2,1)
\end{blockmatrix},\
  \ensuremath{ \tilde{A}_{\temp \temp}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ A_{\temp \temp}}\vphantom{\Big(}(2,1)		\mblock(2,2)\ensuremath{ A_{\temp \tempGamma}}(2,1)
  \mblock(0,1)0\vphantom{\Big(}(2,1)		\mblock(2,1){-I}(2,1)
\end{blockmatrix},\
  \ensuremath{ \tilde{A}_{\temp \V}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ A_{\temp \V}}\vphantom{\Big(}(2,1)
  \mblock(0,1)0\vphantom{\Big(}(2,1)
\end{blockmatrix},\
  \ensuremath{ \tilde{A}_{\V \temp}} = \begin{blockmatrix}
  \mblock(0,2)\ensuremath{ A_{\V \temp}}\vphantom{\Big(}(2,1)		\mblock(2,2)\ensuremath{ A_{\V \tempGamma}}(2,1)
\end{blockmatrix}.
\end{equation*}
With this, the conditions in \cref{eq:heat2,eq:V2} can be incorporated explicitly into the definition of
\ensuremath{ A_{\V \temp}}\time{} and \ensuremath{B_\temp}\time.
Regarding the remaining boundary conditions, the Neumann condition in \cref{eq:heat5} does not need any extra attention, since it is incorporated
automatically.
The Dirichlet boundary conditions in \cref{eq:heat3,eq:V3,eq:V4} can be handled by the common approach
to remove the rows and columns corresponding to the DOFs at the related boundary regions from the matrices.

\cref{eq:blocksss} with \cref{eq:BCDiff} incorporated reads
\begin{equation}\label{eq:blocksssGamma}\begin{aligned}
\begin{blockmatrix}
  \mblock(0,1)\ensuremath{ \tilde{M}_{\temp}}\vphantom{\Big(}(1.6,1)	\mblock(1.6,1)0(1.6,1)
  \mblock(0,0)0\vphantom{\Big(}(1.6,1)		\mblock(1.6,0)0(1.6,1)
\end{blockmatrix}
\ensuremath{\frac{\mathrm{d}}{\mathrm{d}t}} \begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\temp_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\V_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	\end{blockmatrix}
  &= \begin{blockmatrix}
  \mblock(0,1)\ensuremath{ \tilde{A}_{\temp \temp}}\vphantom{\Big(}(2,1)		\mblock(2,1){\ensuremath{ \tilde{A}_{\temp \V}}}(2,1)
  \mblock(0,0)\ensuremath{ \tilde{A}_{\V \temp}}\vphantom{\Big(}(2,1)		\mblock(2,0)\ensuremath{ A_{\V \V}}(2,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{\temp_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)\ensuremath{\V_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
+
\begin{blockmatrix}
	  \mblock(0,2)\ensuremath{B_\temp}\vphantom{\Big(}(1.2,1)
	  \mblock(0,1)0\vphantom{\Big(}(1.2,1)
\end{blockmatrix}
\ensuremath{\mathbf{u}}\ensuremath{^h}, \\ \\
\ensuremath{\mathbf{y}}\ensuremath{^h}\ &=
\begin{blockmatrix}
  \mblock(1,4)\ensuremath{C_\temp}(1,1)	\mblock(2,4)0(1,1)
\end{blockmatrix}
\begin{blockmatrix}
	  \mblock(0,1)\ensuremath{\temp_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
	  \mblock(0,0)\ensuremath{\V_\Delta}\ensuremath{^h}\vphantom{\Big(}(1.2,1)
\end{blockmatrix}.
\end{aligned}\end{equation}
Notably, this modification preserves the DAE structure and index of \cref{eq:blocksss}.
Consequently, we follow~\cite{morFreRM08} and apply the Schur complement to remove the algebraic conditions.
This yields the new system matrices for the Stefan problem
\begin{equation}\label{eq:StefanMats}\begin{aligned}
  \ensuremath{\mathcal{M}}\time &= \ensuremath{ \tilde{M}_{\temp}}\time, \\
  \ensuremath{\mathcal{A}}\time &= \ensuremath{ \tilde{A}_{\temp \temp}}\time - \ensuremath{ \tilde{A}_{\temp \V}}\time\ensuremath{ A_{\V \V}}^{-1}\time\ensuremath{ \tilde{A}_{\V \temp}}\time, \\
  \ensuremath{\hat{\cB}}\time &= \ensuremath{B_\temp}\time, \\
  \ensuremath{\mathcal{C}}\time &= \ensuremath{C_\temp}\time.
\end{aligned}\end{equation}
With the matrices from \cref{eq:StefanMats}, we can transform the Stefan problem into the formulation of \cref{eq:sss}
and are able to apply the LQR approach for the computation of a feedback control.
In order to have a computationally efficient method, we keep the sparse structure of the matrices and
never compute the generally dense Schur complement for \ensuremath{\mathcal{A}}\time{} explicitly.
Instead we apply \ensuremath{\mathcal{A}}\time{} implicitly and use the matrices from \cref{eq:blocksssGamma}.
For details, see~\cite{morFreRM08}, or the numerical implementation in~\cite{SaaKB21-mmess-2.1}.

In \cref{sec:LQR}, we describe the related numerical methods.
For this, we further discretize the Stefan problem in time.
For the simulation forward in time, we use the reference time-steps
\begin{equation}\label{eq:timesteps}\begin{aligned}
  0 &= \ensuremath{t}_0 < \ensuremath{t}_1 < \ldots < \ensuremath{t}_{\ensuremath{n_t} - 1} < \ensuremath{t}_{\ensuremath{n_t}} = \ensuremath{\ttime_\text{end}}, \\
  \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}} &= \{\ensuremath{t}_0, \ensuremath{t}_1, \ldots,\ensuremath{t}_{\ensuremath{n_t} - 1}, \ensuremath{t}_{\ensuremath{n_t}}\}.
\end{aligned}\end{equation}
In \cref{{sec:LQR}}, we solve differential Riccati equations backwards in time.
For this we use the same time-steps $\hat \ensuremath{t}_{k} = \ensuremath{t}_{\ensuremath{n_t} - k}$ in reversed order
\begin{equation}\label{eq:timestepsbwd}\begin{aligned}
  \ensuremath{\ttime_\text{end}} &= \hat \ensuremath{t}_0 > \hat \ensuremath{t}_1 > \ldots > \hat \ensuremath{t}_{\ensuremath{n_t} - 1} > \hat \ensuremath{t}_{\ensuremath{n_t}} = \ensuremath{t}_0, \\
  \ensuremath{\discreteset_{\text{bwd}}} &= \{\hat \ensuremath{t}_0, \hat \ensuremath{t}_1, \ldots,\hat \ensuremath{t}_{\ensuremath{n_t} - 1}, \hat \ensuremath{t}_{\ensuremath{n_t}}\}.
\end{aligned}\end{equation}
We use the same time-steps because the matrices from \cref{eq:StefanMats},
which form the coefficients of the
differential Riccati equations, are
assembled during the forward simulation.
When we apply a time-adaptive method in the forward simulation, additional time-steps can be added to \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}}{}.
\section{Numerical Experiments}\label{sec:experiments}
With several numerical experiments, we demonstrate that the LQR approach is
applicable to this type of control problem, i.e., the Stefan problem, regardless
of several numerical challenges and approximations that arise
similar to other types of problems, where this approach is applied successfully.
The Stefan problem is non-linear, which we address with a linearization.
Further, the coefficients of the corresponding DRE are large-scale matrices.
Thus, we use a low-rank representation of the DRE solution to make it feasible
in regard of the computational cost and the memory requirements.
Additionally, the Stefan problem is a DAE as well as, e.g., the Navier-Stokes equations,
for which the LQR approach is applied successfully in \cite{Wei16}.
To handle the DAE structure, we use an implicit index reduction technique.

In addition,
for the Stefan problem, the matrices are time-dependent and
we approximate the time-derivative of the mass matrix by centered differences.
Consequently, we use the non-autonomous BDF method to solve the DRE and compute
the feedback gain matrices.
Then, at intermediate time-steps that are
generated by the time-adaptive fractional-step-theta scheme,
we interpolate the feedback gain matrices.
An alternative is to interpolate the DRE solution with the corresponding low-rank factors and compute the feedback gain matrices at intermediate time-steps as in \cref{line:BDF8} of \cref{algo:BDF}.
However, this is more expensive and not necessary in our setting.
All these techniques can introduce additional approximation errors.
Still, the computed feedback controls can steer the interface back to
the desired trajectory successfully in our experiments, similar to previously studied problem types.

With several numerical experiments, we illustrate the
performance of the feedback stabilization, which is computed with \cref{algo:LQR}.
The open-loop computations in \cref{line:LQR1} are reported in~\cite{BarBHetal18a} and
the assessment of the runtime performance of different methods for the solution of the DRE in
\cref{line:LQR2} is left to~\cite{BarBSS21a}.
Instead, we focus on the behavior of the feedback stabilization in \cref{line:LQR3}.
This behavior is strongly influenced by the weight factor \ensuremath{\lambda}{} in the cost
functional as well as different choices of inputs and outputs, as the experiments demonstrate.
In order to
asses the robustness of our computed feedback controls, we are interested in uncertainties in the
cooling, i.e., we investigate different perturbations \ensuremath{\varphi}\time{} to the Dirichlet
boundary condition at \ensuremath{\Gamma_{\text{cool}}}\time{}.

\begin{figure}[t]
  \centering
  \vspace{0em}\includegraphics{tikz/externalize/domain_in_out.pdf
  \caption{Input, Perturbation (left), and Output Areas (right)}%
  \label{fig:inout}
\end{figure}
We depict several boundary regions for perturbations (\ensuremath{\Gamma_{\text{cool}}}\time{}) as well as
inputs (\ensuremath{\Gamma_{\ensuremath{u}}}\time{}) and outputs (\Gammaoutx\time{}) in \cref{fig:inout}.
As for \ensuremath{\Omega}\time{}, the boundary regions for inputs $\Gammainx[1]\time, \ldots, \Gammainx[4]\time$
are constant in this representation of the domain.
However, we do not restrict ourselves to this case and, thus, keep the time-dependence in our notation.

The particular domain \ensuremath{\Omega}\time{}, we choose for the experiments, is a rectangle $[0, 0.5] \times [0, 1]$
with the initial interface position at height $0.5$ as a horizontal line.
For the discretization, we choose a mesh of triangles with \(3\,899\) vertices and 401 time-steps.
This, using standard P1 elements in \textsf{FEniCS}, results in \(3\,899\) DOFs for \ensuremath{\varTheta}\time{} and \(7\,798\) DOFs for \ensuremath{\V\ofGradTemp}.
After the removal of the DOFs corresponding to Dirichlet boundary conditions (see \cref{sec:discretization}),
the size of the matrices is \(\ensuremath{n}{} = 11\,429\).
The model parameters are listed in \cref{tab:para}.
\begin{table}[t]
\centering
\begin{tabular}{|l|l|l|l|l|l|l|l|}
\hline
\ensuremath{\ttime_\text{end}}	&	\ensuremath{\kk_s}	&	\ensuremath{\kk_l}		&	\ensuremath{\ensuremath{\varTheta}_\text{cool}}	&	\ensuremath{\ensuremath{\varTheta}_\text{melt}}	&	\ensuremath{\ell}	&	\ensuremath{\ensuremath{\varTheta}_0} 									&	\ensuremath{\timestep\kiterind} \\
\hline
$1$			&	$6$	&	$10$	&	$-1$	&	$0$		&	$10$	&	$\ensuremath{\ensuremath{\varTheta}_\text{cool}} \cdot (1 - 2 \cdot \ensuremath{\state_2})$	&	 $\in \left[10^{-4}, 2.5\cdot10^{-3}\right]$\\
\hline
\end{tabular}
\caption{Stefan Problem Model Parameters}
\label{tab:para}
\end{table}

\subsection*{Code Availability}\label{sec:code}
All codes and data to reproduce the presented results are available at~\cite{LQR_Stefan_codes}.
The non-autonomous BDF method is incorporated in the software package \mbox{\textsf{M-}\mess{}}{}~2.1~\cite{SaaKB21-mmess-2.1}

\subsection{Experiment 1}\label{sec:ex1}
For this first experiment, we focus on the influence of weight parameters and outputs
on the performance of the feedback control.
The desired interface trajectory, which we intend to stabilize with the feedback control,
is a flat horizontal line moving from its
initial height at \(0.5\) downward by \(0.004\).
Since we assume the Stefan problem to be asymptotically stable,
the interface returns to the desired trajectory after a perturbation without a feedback control as well.
However, this is not achieved within the time horizon.
With the feedback control, we intend to prevent the interface deviation to a certain extend and to steer the interface back in a much shorter time after a perturbation occurs.

The interface deviation is caused by the perturbation \ensuremath{\varphi}\time{}, which is  generated randomly in the form of three scalar values in the range of $[-\ensuremath{\ensuremath{\varTheta}_\text{cool}}, \ensuremath{\ensuremath{\varTheta}_\text{cool}}]$.
These are applied for a period of four time-steps of \ensuremath{\discreteset_{\text{fwd}}^{\text{ref}}}{} at the times $0.1$, $0.3$, and $0.5$ on $\ensuremath{\Gamma_{\text{cool}}}\time =
\Gammacoolx[1]\time\cup\Gammacoolx[2]\time$:
\begin{equation*}
  \ensuremath{\varTheta}	= \ensuremath{\ensuremath{\varTheta}_\text{cool}} + \ensuremath{\varphi}, \qquad \text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{cool}}}.
\end{equation*}
The trajectory of \ensuremath{\varphi}\time{} is pictured in the top part of \cref{fig:weights}.
The scaling of the y-axis is relative to \ensuremath{\ensuremath{\varTheta}_\text{cool}}.

Here, for the LQR problem, we use a single input \ensuremath{\control_\ensuremath{\mathcal{K}}}\time{} with the same value on
$\Gammainx[1]\time\cup\ldots\cup\Gammainx[4]\time$, i.e.\ $\ensuremath{\mathcal{B}}\time \in \mathbb{R}^{\ensuremath{n} \times 1}$.
The top part of \cref{fig:weights} shows the control for three different LQR
designs, resulting from three different combinations of weights and outputs:
\begin{align*}
  &\controlx[1]: &&\ensuremath{\lambda} = 10^{-4}, 		&&\text{2 outputs:}\ \Gammaoutx[3],\Gammaoutx[4],\\
  &\controlx[2]: &&\ensuremath{\lambda} = 10^{-6}, 		&&\text{2 outputs:}\ \Gammaoutx[3],\Gammaoutx[4],\\
  &\controlx[3]: &&\ensuremath{\lambda} = 1.6\cdot10^{-2}, 	&&\text{7 outputs:}\ \Gammaoutx[1],\ldots,\Gammaoutx[6],\ensuremath{\Gamma_{\text{int}}}.
\end{align*}
The outputs on \Gammaoutx[1]\time{}, \Gammaoutx[2]\time{}, \Gammaoutx[5]\time{},  and \Gammaoutx[6]\time{}
measure averaged temperatures on the corresponding interval, while
\Gammaoutx[3]\time{} and \Gammaoutx[4]\time{} represent point measurements of the
temperature.
The output at \ensuremath{\Gamma_{\text{int}}}\time{} monitors the difference of the interface movement
to the desired movement
\begin{equation*}
  \langle \ensuremath{\mathcal{C_{\text{int}}}} , \ensuremath{\textsf{v}\ensuremath{^h}} \rangle
    = \int\limits_{\ensuremath{\Gamma_{\text{int}}}} \Big(\frac{1}{\ensuremath{\ell}}[\ensuremath{\kk_s}(\nabla \ensuremath{\temp_\Delta})_s
    - \ensuremath{\kk_l}(\nabla \ensuremath{\temp_\Delta})_l]\Big) \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}} \cdot \ensuremath{\textsf{v}\ensuremath{^h}} \mathrm{d}s
\end{equation*}
and is defined via the jump term of the Stefan condition (\cref{eq:V2}).
This output alone is not suited for an effective control-design.
It would generate a large output, and thus an active feedback-response, only while the perturbation is actively
driving the interface away from the desired trajectory.
However, it can not detect a difference in the position of the interface.
Consequently, it would not steer the interface back but would keep it on a ``parallel trajectory".

Thus, in the first two LQR designs, the controls \controlx[1]\time{} and \controlx[2]\time{}
are based on two outputs that measure the temperature at the desired interface position on the boundary, such that $\ensuremath{\mathcal{C}}\time{} \in \mathbb{R}^{2 \times \ensuremath{n}}$, while
\controlx[3]\time{}, in the third setting, uses seven outputs, i.e.\ $\ensuremath{\mathcal{C}}\time{} \in \mathbb{R}^{7 \times \ensuremath{n}}$.
With this, we compute the feedback gain matrices according to \cref{sec:LQR} with \cref{algo:BDF} and set $\ensuremath{\wp} = 1$.
Then, we simulate the closed-loop system with these feedback gain matrices together with the perturbations.
\begin{figure}[tp]
  \centering
  \includegraphics{tikz/externalize/Stefan_feedback_adaptive.pdf
  \caption{Perturbation and feedback (top), relative interface position (middle),
  			and time-step sizes (bottom) for different weights and outputs
	    (\nameref{sec:ex1})}%
  \label{fig:weights}
\end{figure}
The resulting feedback controls \controlx[1]\time{}, \controlx[2]\time{},
and \controlx[3]\time{} are displayed in the top part of \cref{fig:weights} and
the interface positions in the middle part of \cref{fig:weights}.
The interface positions are relative to the desired interface position
at the point $\ensuremath{\ensuremath{x}^*}\time$ with the largest deviation on the interface:
\begin{align*}
  \ensuremath{\Gamma_{\text{int},\Delta}(\ensuremath{t},\statexstar)} &= \ensuremath{{\Gamma}_{\text{int,ref}}(\ensuremath{t}, \statexstar)} - \ensuremath{\Gamma_{\text{int}}}(\ensuremath{t}, \ensuremath{\ensuremath{x}^*}), \\
  \ensuremath{\ensuremath{x}^*} &= \argmax_{\ensuremath{\state_1} \in [0,0.5]}\ \left|\ensuremath{{\Gamma}_{\text{int,ref}}(\ensuremath{t}, \ensuremath{\state_1})} - \ensuremath{\Gamma_{\text{int}}}(\ensuremath{t}, \ensuremath{\state_1})\right|.
\end{align*}

The feedback control \controlx[1]\time{} is most active shortly after the
perturbation starts and it steers the interface back to the desired position in
much reduced time, as expected from theory.
Again as expected, a smaller weight factor \ensuremath{\lambda}{} allows the control
\controlx[2]\time{} to react with a larger magnitude input while the perturbation is active.
It is, therefore, able to stop the interface from deviating earlier and drives
it back even faster.

It should be mentioned that smaller weights reduce the contribution of the control cost term
in the cost functional \labelcref{eq:cost}.
This term can also be interpreted as a regularization term.
Consequently, smaller weights decrease the regularity of the feedback control problem.
We note that this results in larger computational cost for the solution of the DRE due to slower convergence
of the internal iterative solvers used in the BDF method, but refer to~\cite[Section~4.2]{BarBSS21a}
for a more detailed discussion of this issue.

The weights for the LQR settings of \controlx[1]\time{} and \controlx[2]\time{}
are not directly comparable to the setting for \controlx[3]\time, since this is based on different outputs.
The outputs at \Gammaoutx[5]\time{} and \Gammaoutx[6]\time{} allow it to detect the temperature perturbation earlier
and the output that monitors the movement of \ensuremath{\Gamma_{\text{int}}}\time{} can observe the deviation of the interface earlier.
Thus, \controlx[3]\time{} is most active immediately after the perturbation starts and moves the interface back much faster.

Note, that we use a time-adaptive fractional-step-theta scheme
with a relative feedback control-based indicator function
to simulate the Stefan problem together with a feedback control.
For this example, the bottom part of \cref{fig:weights}
shows the adaptive time-step sizes.
While the feedback control is inactive, like at the beginning or end of the simulation,
the method chooses the largest time-step size of $\ensuremath{\timestep\kiterind} = 0.0025$.
As soon as the feedback control is active, the method chooses the minimal time-step size of $\ensuremath{\timestep\kiterind} = 10^{-4}$.
This prevents numerical instabilities that can occur with very small \ensuremath{\lambda}{} and thus very active feedback controls
as is demonstrated in~\cite[Section~4.2]{BarBSS21a}.

\nameref{sec:ex1} demonstrates the influence of the chosen weight parameter \ensuremath{\lambda}{} in the cost functional (\cref{eq:cost}).
Smaller weight parameters increase the impact of the deviation term in the cost functional indirectly by
decreasing the relevance of the control costs which are measured by the second term and scaled by \ensuremath{\lambda}.
Further, the choice and number of outputs has a significant impact on the performance of the feedback control.

The next experiment demonstrates that our method works as well with increased time-dependence of the data
in the non-autonomous DRE.

\subsection{Experiment 2}\label{sec:ex2}
In our second experiment, the desired interface trajectory is, again, a flat horizontal line moving from its
initial height at \(0.5\), this time, upward by \(0.1\).
This is a longer distance interface movement than in \nameref{sec:ex1} on the same time horizon and, in turn,
results in a stronger time-dependence of the matrices for the DRE.\@
The investigated LQR design is chosen as in the second case in \nameref{sec:ex1}.
\begin{figure}[t]
  \includegraphics{tikz/externalize/plot_moveup.pdf
  \caption{Perturbation and feedback (top) and relative interface position (bottom) for an interface moving upwards
	    (\nameref{sec:ex2})}%
  \label{fig:up}
\end{figure}

That means, we use the same weight ($\ensuremath{\lambda} = 10^{-6}$) and outputs as for \controlx[2]\time{}
above and call  the corresponding feedback control \controltilde[2]\time.
This time we use four randomly generated perturbations applied at the times $0.1, 0.3, 0.5$, and $0.7$, which, together with
\controltilde[2]\time, can be found in the top part of \cref{fig:up}.
Again, the perturbations cause the interface to deviate from the desired trajectory.
The feedback control \controltilde[2]\time{} behaves similarly to the previous experiment.
It stops the interface from deviating and drives it back to the desired position,
as expected from the theory.
Additionally, the time-adaptivity behaves analogously to \nameref{sec:ex1}.

The experiment showcases the performance of our method with strongly time-varying coefficients in the DRE.
We implemented more experiments with desired interface trajectories that move
considerably.
Our feedback control showed the same performance for all of them. This demonstrates
that our solver can also cope with the stronger time dependence without difficulties.

\subsection{Experiment 3}\label{sec:ex3}
Our third experiment revisits the basic task as in \nameref{sec:ex1}.
This time, we generate two different perturbations $\pertx[1]\time,
\pertx[2]\time \in [-\ensuremath{\ensuremath{\varTheta}_\text{cool}}, \ensuremath{\ensuremath{\varTheta}_\text{cool}}]$, acting on the left (\Gammacoolx[1]\time{}) and right (\Gammacoolx[2]\time{}) part
of the Dirichlet boundary at the bottom of \ensuremath{\Omega}\time{} (see \cref{fig:inout}),
with three random values each, which are applied at the times $0.1, 0.4325$, and $0.765$.
They are displayed in the top part of \cref{fig:curv}.
We apply these two perturbation functions to
\begin{align*}
  \ensuremath{\varTheta}	&= \ensuremath{\ensuremath{\varTheta}_\text{cool}} + \pertx[1], \qquad \text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \Gammacoolx[1],\\
  \ensuremath{\varTheta}	&= \ensuremath{\ensuremath{\varTheta}_\text{cool}} + \pertx[2], \qquad \text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \Gammacoolx[2].
\end{align*}
These perturbations will not only move the interface away from the desired position
but also add some curvature to the deviated interface.
To stabilize this interface position, we use six individual inputs at $\Gammainx[1]\time,\ldots,\Gammainx[6]\time$
($\ensuremath{\mathcal{B}}\time{} \in \mathbb{R}^{\ensuremath{n} \times 6}$) and the same two outputs at \Gammaoutx[3]\time{}, \Gammaoutx[4]\time{}
as in the previous experiments ($\ensuremath{\mathcal{C}}\time{} \in \mathbb{R}^{2 \times \ensuremath{n}}$).
For the weight factor in the cost functional we use $\ensuremath{\lambda} = 10^{-9}$.
Again, this weight factor is not directly comparable to the previous experiments
since we use a different combination of inputs and outputs.

\begin{figure}[tp]
  \includegraphics{tikz/externalize/plot_curvature.pdf
  \caption{Two different perturbations (top), and perturbed and feedback controlled relative interface positions (bottom)
	    (\nameref{sec:ex3})}%
  \label{fig:curv}
\end{figure}
The resulting perturbed and feedback controlled interfaces
are displayed in the bottom part of \cref{fig:curv} at the
time points $\ensuremath{t} \in [0.465, 0.56, 0.798, 1]$.
With the first two perturbations, the interface is pushed downwards and assumes a distinct curvature.
Right after the second perturbation ($\ensuremath{t} = 0.465$), the feedback controlled interface is already essentially
back to the desired interface, but still has undesired curvature.
Well before the third perturbation, at $\ensuremath{t} = 0.56$, the feedback controlled interface
is almost flat again, as the desired interface.
At the two endpoints of the interface, where the outputs are located,
the controlled interface is directly at the desired position.
The remaining curvature can not be measured by the cost functional and is, thus, not
actively removed.
For the given setup, this is the expected behavior of the feedback control.
It is similar after the first perturbation.
The third perturbation moves the interface upwards above the desired position ($\ensuremath{t} = 0.798$).
At the end of the time interval ($\ensuremath{t} = 1$), the feedback controlled interface is again almost flat and back
to the desired position at the ends while the uncontrolled interface is still far away.

This experiment shows that the feedback control can move the perturbed interface back to the desired trajectory
and additionally control the curvature of the interface, to a certain extent.
With the chosen input and output setting, the position of the interface is corrected as fast as in the previous experiments.
The curvature correction requires some more time but is still performed during the considered time interval
up to the level that the cost functional can measure with the given two outputs.

\medskip
In all three experiments, we are able to stabilize the desired interface position
with the proposed feedback control approach
despite the different challenges that are posed by the Stefan problem.

\section{Introduction}
The solidification and melting of materials is an active and intensively studied field with numerous applications.
This phase-change problem can be modeled by a non-linear PDE and is often called Stefan problem after J.~Stefan
who describes it in his works~\cite{Ste89a,Ste89b,Ste90}.
In a certain domain, the temperature of the material is either below, above, or equal to
the specific melting temperature of the material.
Accordingly, the domain is split into a solid and a liquid phase, which are separated by an interface or inner boundary.
J.~Stefan formulates what is now called Stefan condition in~\cite{Ste89a}, which couples the time-derivative of the interface position
with the jump of the temperature gradient along the interface.
On the other hand, J.~Stefan was not the first to consider this type of problem.
G.~Lam{\'e} and B.~P.~Clapeyron were concerned with this in an earlier work as well~\cite{LamC31},
such that it is also called Lam{\'e}--Clapeyron problem.

Several books address the Stefan problem, e.g.,~\cite{Rub71,NieCM11,Gup18,KogK20a}.
An extensive historical survey of the Stefan problem can be found in the book by L.~I.~Ruben\v{s}te\u{\i}n~\cite[Introduction: \S 1]{Rub71}.
Early works on the Stefan problem usually consider the one-dimensional case, while higher dimensional cases where only studied decades after the original publication.

The aim of this manuscript is to compute and apply feedback control to the two-dimensional two-phase Stefan problem.
While this problem has been studied in combination with open-loop controls several times, see
e.g.~\cite{Zie08, Ber10, AntNS14, AntNS15, BarBHetal18a} and the reference therein, only
recently, also closed-loop, i.e. feedback control, for the Stefan problem has been discussed in~\cite{KogK19,KogK20,KogK20a}.
However, these works address only the one-dimensional case.
On the one hand, the novelty of our work lies in the consideration of feedback control for the two-phase Stefan problem in two spatial dimensions
in contrast to the one-phase or one-dimensional Stefan problem.
On the other hand, we apply the linear-quadratic regulator (LQR) approach to the Stefan problem, which goes beyond the types of problems
that have been studied in connection with this approach.

In particular, the application of LQR requires the treatment of the non-linearities, the differential-algebraic nature,
and the time-varying character of the Stefan problem.
In this manuscript, we address the details on how to transform it into a linear ordinary differential equation,
apply the LQR approach to compute feedback controls, and then use these feedback controls
to stabilize the interface position in the non-linear differential-algebraic problem.
In order to do this, we linearize the Stefan problem around a reference trajectory and assemble a generalized differential Riccati equation (DRE).
For the derivation of the LQR problem and the resulting DRE we refer to~\cite{Rei72, BitG91}.
An extensive study on DREs and their numerous applications can be found in~\cite{AboFIetal03}.
Since the Stefan condition is an algebraic equation, which is coupled to the Stefan problem,
the DRE we are considering results from a differential-algebraic equation (DAE).
Details on generalized DREs can be found in~\cite{KunM90} and for DREs resulting from DAEs, see~\cite{KunM89}.

To solve the Stefan problem numerically and assemble the DRE, we discretize in space using the
finite element method (FEM) by applying the software \textsf{FEniCS}~\cite{dolfin}.
This results in a large-scale matrix-valued DRE with time-dependent coefficients, denoted as a non-autonomous DRE.\
While we observe in numerical experiments that, usually, its solution can be well
approximated by a low-rank factorization, this is proven theoretically
only for the autonomous DRE in~\cite{Sti18}.
Well known methods, which use the low-rank structure of the numerical solution,
are splitting schemes~\cite{Sti15,Sti15a,Sti18a,OstPW18,MenOPetal18},
Rosenbrock and Peer methods~\cite{Men12,LanMS15,Lan17,BenL18}
as well as the backward differentiation formulas (BDF)~\cite{BenM04, Men12, LanMS15, BenM18}.
Krylov subspace methods~\cite{BehBH21,KosM20,KirSi19,GueHJetal18} and exponential integrators~\cite{LiZL20}
for DREs have been developed recently as well.
In~\cite{Men12, BenM18} Rosenbrock and BDF methods and in~\cite{LanMS15,Lan17} also Peer methods are studied for a non-autonomous DRE
where the mass matrix is constant.
In extension of this, splitting schemes and BDF methods are developed in~\cite{BarBSS21a} for non-autonomous DREs with a time-varying mass matrix.
In contrast to the BDF methods, the splitting schemes require that
the coefficients can be decomposed into a time-dependent scalar function times a constant matrix.
The non-autonomous DREs we consider go beyond this case.
We use mesh movement techniques to track the interface, and, as a consequence, the single matrix entries change very differently, possibly for all matrices.
Thus, the BDF methods are
the most promising method for the non-autonomous DRE resulting from the Stefan problem and we apply the non-autonomous BDF method from~\cite{BarBSS21a}.

\paragraph{Structure of the Manuscript}
In this manuscript, we derive and apply several feedback controls to the two-dimensional two-phase Stefan problem.
For this, we state the equations that model the Stefan problem and the mesh movement that we use to track the interface position (\cref{sec:stefan}).
Then, we linearize the resulting system and semi-discretize it in space (\cref{sec:discretization}).
With the resulting time-dependent matrices we formulate a non-autonomous DRE resulting from the LQR approach (\cref{sec:LQR}).
We solve this DRE with the non-autonomous BDF method and use the resulting feedback gain matrices to compute a feedback control (also \cref{sec:LQR}).
We apply the described methods in several numerical experiments for different parameter settings and specify where all our codes
and data are available (\cref{sec:experiments}).

\paragraph{Notation}
In most equations, we omit the time-dependence \time{}, the spatial dependence $(\ensuremath{x})$, or the combination of both $(\ensuremath{t}, \ensuremath{x})$.
This is supposed to improve the readability of the equations.

\section{Two-dimensional Two-phase Stefan Problem}\label{sec:stefan}
In this section, we define the equations describing the Stefan problem, its boundary conditions and initial values.
This involves equations characterizing the temperature and interface movement.
Our goal is to
use a Riccati feedback control approach
to control the interface position.
For this,
we choose a sharp interface representation in the formulation of the Stefan problem,
opposed to, e.g., the level set representation of the interface in \cite{Ber10}.
As in~\cite{BaePS10,BaePS13}, we extend the interface movement to the whole domain.
This ensures the mesh regularity for the semi-discretized Stefan problem that is described in \cref{sec:discretization}.

At each \(\ensuremath{t} \in [0, \ensuremath{\ttime_\text{end}}]\), the domain is $\ensuremath{\Omega}\time \subset
\mathbb{R}^2$.
One instance is illustrated in \cref{fig:domain}.
We split the domain \ensuremath{\Omega}\time{} into the two regions corresponding to the two phases.
These are the region where the material is in its solid phase \ensuremath{\domain_s}\time{} and accordingly, the region \ensuremath{\domain_l}\time{} related to the liquid phase.
The two phases are separated by the interface \ensuremath{\Gamma_{\text{int}}}\time.
This inner phase-boundary can move such that its position is time-dependent.
Thus, also the two phases \ensuremath{\domain_s}\time{} and \ensuremath{\domain_l}\time{} are time-dependent
and, as a consequence, so is the whole domain \ensuremath{\Omega}\time{} and its boundary regions.
The boundary of \ensuremath{\Omega}\time{} is separated into \ensuremath{\Gamma_{\ensuremath{u}}}\time{}, \ensuremath{\Gamma_{\text{cool}}}\time{} and \ensuremath{\Gamma_\text{N}}\time{}
as depicted in \cref{fig:domain}.
Note that the outer shape of \ensuremath{\Omega}\time{} is constant for the realization chosen in this manuscript.
Thus, the time-dependence of \ensuremath{\Omega}\time{} is not absolutely necessary even though
its sub-domains are time-dependent.
However, we keep the time-dependence in our notation in order to not restrict our methods
to this case.

\includetikz{domain}

We follow the definition of the Stefan problem from~\cite{BarBHetal18a}. However, here we use
it in a more compact form, i.e.\ omitting the couplings with the Navier-Stokes equations and the interface graph formulation.
While the Navier-Stokes equations alone add additional algebraic constraints, making the DAE harder to classify,
both the Navier-Stokes equations and the interface graph formulation add more nonlinearities to the problem.
In order to develop the general numerical strategy,
we first want to study the feedback control problem for this simplified setting without these couplings,
before delving into the additional technical challenges of the full problem formulation.

We denote the temperature as \ensuremath{\varTheta}\time{} and model it with the partial
differential \cref{eq:heat}:
\begin{subequations}\label{eq:heat}\begin{align}
    \label{eq:heat1}
    \dot\ensuremath{\varTheta} - \ensuremath{\Upsilon} \cdot \nabla \ensuremath{\varTheta} - \ensuremath{\alpha} \Delta \ensuremath{\varTheta}	&= 0,  				& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Omega}, \\
    \label{eq:heat2}
    \partial_{\ensuremath{\boldsymbol{n}}}\ensuremath{\varTheta}						&= \ensuremath{u},			& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\ensuremath{u}}}, \\
    \label{eq:heat3}
    \ensuremath{\varTheta}								&= \ensuremath{\ensuremath{\varTheta}_\text{cool}},			& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{cool}}}, \\
    \label{eq:heat4}
    \ensuremath{\varTheta}								&= \ensuremath{\ensuremath{\varTheta}_\text{melt}},			& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_{\text{int}}}, \\
    \label{eq:heat5}
    \partial_{\ensuremath{\boldsymbol{n}}} \ensuremath{\varTheta}						&= 0,				& &\text{on}\ (0, \ensuremath{\ttime_\text{end}}] \times \ensuremath{\Gamma_\text{N}}, \\
    \label{eq:heat6}
    \ensuremath{\varTheta}(0)								&= \ensuremath{\ensuremath{\varTheta}_0},			& &\text{on}\ \ensuremath{\Omega}.
\end{align}\end{subequations}

In \cref{eq:heat2}, we apply the control \ensuremath{u}\time{} as a Neumann condition on the control boundary \ensuremath{\Gamma_{\ensuremath{u}}}\time.
The \cref{eq:heat3,eq:heat4} describe the Dirichlet conditions on the cooling boundary \ensuremath{\Gamma_{\text{cool}}}\time{} and the interface \ensuremath{\Gamma_{\text{int}}}\time{}
with the constants \ensuremath{\ensuremath{\varTheta}_\text{cool}}{} and \ensuremath{\ensuremath{\varTheta}_\text{melt}}{}, respectively.
\cref{eq:heat6} presents the initial condition with the initial temperature distribution \ensuremath{\ensuremath{\varTheta}_0}.
The heat conductivities in the solid phase \ensuremath{\kk_s}{} and in the liquid phase \ensuremath{\kk_l}{} are collected in \ensuremath{\alpha}\ensuremath{(\temp\time)}{}:
\begin{equation*}
    \ensuremath{\alpha} = \begin{cases}
		      \ensuremath{\kk_s}, & \text{on}\ \ensuremath{\domain_s},\\
		      \ensuremath{\kk_l}, & \text{on}\ \ensuremath{\domain_l}.\\
               \end{cases}
\end{equation*}
In \cref{eq:heat1}, the temperature is coupled with the extended interface movement \ensuremath{\V\ofGradTemp}{}.
For every $\ensuremath{t} \in (0, \ensuremath{\ttime_\text{end}}]$, we model \ensuremath{\V\ofGradTemp}{} with a system of algebraic equations,
in the sense that they do not contain any time-derivatives:

\begin{subequations}\label{eq:V}\begin{align}
    \label{eq:V1}
    \Delta \ensuremath{\Upsilon}		&= 0, 									& &\text{on}\ \ensuremath{\Omega}, \\
    \label{eq:V2}
    \ensuremath{\Upsilon} 			&= \Big(\frac{1}{\ensuremath{\ell}}[\ensuremath{k}(\nabla \ensuremath{\varTheta}\time)]_l^s\Big) \cdot \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}},	& &\text{on}\ \ensuremath{\Gamma_{\text{int}}},\\
    \label{eq:V3}
    \ensuremath{\Upsilon}			&= 0, 									& &\text{on}\ \ensuremath{\Gamma_{\text{cool}}} \cup \ensuremath{\Gamma_{\ensuremath{u}}}, \\
    \label{eq:V4}
    \ensuremath{\Upsilon} \cdot \ensuremath{\boldsymbol{n}}	&= 0, 									& &\text{on}\ \ensuremath{\Gamma_\text{N}}, \\
    \label{eq:V5}
    \ensuremath{\Upsilon}(0)			&= 0, 									& &\text{on}\ \ensuremath{\Omega},
\end{align}\end{subequations}

On the interface, $\restr{\ensuremath{\Upsilon}}{\ensuremath{\Gamma_{\text{int}}}}\ensuremath{(\nabla \temp\time)} = \ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}$ is the interface movement in normal direction,
where \ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}}\time{} is the unit normal vector pointing from \ensuremath{\domain_s}\time{} to \ensuremath{\domain_l}\time.
\ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}{} is coupled to \ensuremath{\varTheta}\time{} through the Stefan condition~\eqref{eq:V2}.
Here, \ensuremath{\ell}{} is the latent heat constant and
\begin{equation}\label{stefan}
   [k(\nabla \ensuremath{\varTheta})]_l^s = \ensuremath{\kk_s}\partial_{\ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}}}\restr{\ensuremath{\varTheta}}{\ensuremath{\domain_s}} - \ensuremath{\kk_l}\partial_{-\ensuremath{\boldsymbol{\ensuremath{\boldsymbol{n}}_{\text{int}}}}}\restr{\ensuremath{\varTheta}}{\ensuremath{\domain_l}}
\end{equation}
is the jump of the temperature gradient along \ensuremath{\Gamma_{\text{int}}}\time.
\cref{eq:V3,eq:V4} ensure that the outer boundaries of \ensuremath{\Omega}\time{} do not move
such that the outer shape of the domain does not change.

We use \ensuremath{\V\ofGradTemp}{} in \cref{sec:discretization} for the mesh movement.
The Stefan problem, which is described by the system of \cref{eq:heat,eq:V},
is a non-linear system of DAEs on a
time-varying domain.
Due to the coupling of \ensuremath{\varTheta}\time{} and \ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}{} in~\eqref{eq:V2}
and the two temperature-dependent phases in the definition of \ensuremath{\alpha}\ensuremath{(\temp\time)}{},
both, \ensuremath{\ensuremath{\Upsilon}_{\text{int}}}\ensuremath{(\nabla \temp\time)}{} and \ensuremath{\alpha}\ensuremath{(\temp\time)}{}, depend on \ensuremath{\varTheta}\time{}.
Thus, the terms $\ensuremath{\V\ofGradTemp} \cdot \nabla \ensuremath{\varTheta}\time$ and $\ensuremath{\alpha}\ensuremath{(\temp\time)} \Delta \ensuremath{\varTheta}\time$ in \cref{eq:heat1} are nonlinear.

redpajama_set_name: RedPajamaArXiv