# Extremum Seeking for Stefan PDE with Moving Boundary

Maurício Linhares Galvão<EMAIL_ADDRESS>Tiago Roux Oliveira
Dept. of Electronics and Telecommunication Engineering, State University of
Rio de Janeiro (UERJ), Rio de Janeiro, RJ 20550-900, Brazil.
<EMAIL_ADDRESS>Miroslav Krstic  Dept. of Mechanical and Aerospace
Engineering, University of California - San Diego (UCSD), La Jolla, CA
92093-0411, USA<EMAIL_ADDRESS>
###### Abstract

This paper presents the design and analysis of the extremum seeking for static
maps with input passed through a partial differential equation (PDE) of the
diffusion type defined on a time-varying spatial domain whose boundary
position is governed by an ordinary differential equation (ODE). This is the
first effort to pursue an extension of extremum seeking from the heat PDE to
the Stefan PDE. We compensate the average-based actuation dynamics by a
controller via backstepping transformation for the moving boundary, which is
utilized to transform the original coupled PDE-ODE into a target system whose
exponential stability of the average equilibrium of the average system is
proved. The discussion for the delay-compensated extremum seeking control of
the Stefan problem is also presented and illustrated with numerical
simulations.

###### keywords:

Adaptive Control, Extremum Seeking, Partial Differential Equations, Stefan
Problem, Delays, Averaging Theory, Backstepping in Infinite Dimensions.

††journal: arxiv.org

##  1 Introduction

Extremum seeking (ES) is a non-model based approach in the field of adaptive
control which searches in real time extremum points (maximum or minimum) of a
performance index of a system. This method has received great attention in the
control community by the means of facing control problems when the plant has
imperfections in its model or uncertainties dynamics [1].

ES was first introduced in [2] for maximizing power transfer to a tram car.
Along the history, the number of publications concerning ES remained low until
it’s first general stability proof for stable dynamic systems with unknown
output maps was carried out [1]. Since then, important studies in theory and
applications were developed, such as [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

Reference [14] was the first one to handle partial differential equations
(PDEs) with ES scheme, addressing the design and analysis of multi-variable
static maps subject to arbitrarily long time delays. The delays pointed out by
the authors can be modelled as first-order hyperbolic transport PDEs [15].
This idea has enabled the development of extension to other classes of PDEs
with fixed domains [16, 17].

On the other hand, a large number of applications in various areas appear as
moving boundary or phase change problems, such as [18] and [19]. Usually,
these kind of problems arise in heat conduction situations and need to be
solved in a time-dependent space domain with a moving boundary condition. For
this reason diffusion PDEs with moving boundaries, known as “Stefan problem”
have been studied actively for the last few decades [20], [21].

The dynamics of the the position of the moving boundary in the Stefan problem
is governed by an ODE which depends on the PDE state, generating a nonlinear
coupling of the PDE and ODE dynamics, increasing the complexity of the problem
when compared to conventional analyses for PDEs of fixed domains (not
depending on time or states) and ODEs.

In this paper, we develop an ES controller for the Stefan problem. The
integrator that is usually employed in the ES scheme can be leveraged as part
of the Stefan model, just like the one proposed in [22]. The objective of the
ES will be to find the maximizer interface $s^{*}$ of some unknown map
$Q(s^{*})$ aiming to regulate the phase change interface position to a value
that attains the extremum. For this purpose, we design a compensator of the
heat PDE with moving boundary and the probing signal, which is the result of
solving the problem of generating a sinusoid at the distal end of a boundary-
actuated heat equation. As a further contribution, we also study the effect of
delays in the PDE actuation dynamics by compensating it via predictor feedback
[23].

An important discussion is about the validation of the Stefan model. Although
the usual sinusoidal movement provoked by the ES algorithm may violate the
phase maintenance when the extremum is achieved or during the transient
either, we can keep the phase maintenance at least for the average system,
thus, preserving the convergence analysis. The main theoretical contribution
with respect to the previous conference paper [24] is the inclusion of the
proof for the main theorem not presented before due to space limitations.

This paper is organized as follows. The problem statement is presented in
Sections 2 and 3, containing the relation of the ES scheme with the one-phase
Stefan problem. In Section 4, the control related problem is designed. Section
5 describes the stability analysis of the closed-loop average system with
exponential convergence in $\mathscr{H}_{1}$-norm of the distributed
temperature. The convergence of the moving boundary to the desired equilibrium
is discussed in Section 6. Section 7 shows the design of a compensator of the
transport-heat PDE cascade in case of delays in the actuation dynamics.
Numerical simulations are provided in Section 8, followed by conclusions in
Section 9.

Norms and Notations: We denote the partial derivatives of a function $u(x,t)$
as $\partial_{x}u(x,t)=\partial u(x,t)/\partial x$,
$\partial_{t}u(x,t)=\partial u(x,t)/\partial t$. We conveniently use the
compact form $u_{x}(x,t)$ and $u_{t}(x,t)$ for the former and the latter,
respectively. The two-norm of a finite-dimensional ODE state vector $X(t)$ is
denoted by single bars, $|X(t)|$. In contrast, norms of functions (of $x$) are
denoted by double bars. We denote the spatial $\mathcal{L}_{2}[0,D]$ norm of
the PDE state $u(x,t)$ as
$\|u(t)\|^{2}_{\mathcal{L}_{2}([0,D])}\coloneqq\int_{0}^{D}u^{2}(x,t)\,dx$,
where we drop the index $\mathcal{L}_{2}([0,D])$ in the following, hence
$\|\cdot\|=\|\cdot\|_{\mathcal{L}_{2}([0,D])}$, if not otherwise specified.
Moreover, the $\mathscr{H}_{1}$-norm is given by
$\|u(t)\|^{2}_{\mathscr{H}_{1}}=\|u(t)\|^{2}_{\mathcal{L}_{2}}+\|u_{x}(t)\|^{2}_{\mathcal{L}_{2}}$.
As defined in [25], a vector function $f(t,\epsilon)\in\mathbb{R}^{n}$ is said
to be of order $\mathcal{O}(\epsilon)$ over an interval $[t_{1},t_{2}]$, if
$\exists k,\bar{\epsilon}$ : $|f(t,\epsilon)|\leq k\epsilon,\forall
t\in[t_{1},t_{2}]$. In most cases, we provide no precise estimates for the
constants $k$ and $\bar{\epsilon}$, and we use $\mathcal{O}(\epsilon)$ to be
interpreted as an order of magnitude relation for sufficiently small
$\epsilon$.

##  2 One-phase Stefan Problem

Figure 1: Schematic of one-phase Stefan problem [26, 27]. The temperature
profile in the solid phase is assumed to be a uniform melting temperature.

The physical model which describes the 1-D Stefan problem in a pure one-
component material of length $L$ is described in Figure 1. The domain $[0,L]$
is divided in two sub-domains $[0,s(t)]$ and $[s(t),L]$ which represents the
liquid phase and the solid phase, respectively. The system is controlled by
the heat flux $q_{c}(t)$ at $x=0$, because we are dealing with a Neumann
boundary actuation as shown below:

$\displaystyle T_{t}(x,t)$ $\displaystyle=\alpha T_{xx}(x,t),\quad
x\in(0,s(t)),\quad\alpha=\dfrac{k}{\rho C_{p}}$ (1) $\displaystyle-
kT_{x}(0,t)$ $\displaystyle=q_{c}(t)$ (2) $\displaystyle T(s(t),t)$
$\displaystyle=T_{m}$ (3) $\displaystyle\dot{s}(t)$ $\displaystyle=-\beta
T_{x}(s(t),t),\quad\beta=\dfrac{k}{\rho\Delta H^{*}},$ (4)

where $T(x,t)$, $T_{m}$, $q_{c}(t)$, $k$, $\rho$, $C_{p}$ and $\Delta H^{*}$
are the distributed temperature of the liquid phase, melting temperature,
manipulated heat flux, liquid heat conductivity, liquid density, liquid heat
capacity and latent heat of fusion, respectively. Equations (2) and (3) are
the boundary conditions of the system and (4) is the Stefan condition, which
describes the dynamic of the moving boundary. Figure 2 shows the block diagram
of the PDE-ODE cascade represented by equations (1)-(4).

Figure 2: The cascade of the PDE dynamics and the ODE system.

##  3 Problem Statement

For the sake of simplicity, we consider actuation dynamics which are described
by a heat equation with $\alpha,\beta,k=1$, $\theta(t)\in\mathbb{R}$ and the
propagated actuator $\Theta(t)\in\mathbb{R}$ given by

$\displaystyle\dot{\Theta}(t)=\dot{s}(t)$
$\displaystyle=-\alpha_{x}(s(t),t)),\quad x\in(0,s(t))$ (5)
$\displaystyle\partial_{t}\alpha(x,t)$
$\displaystyle=\partial_{xx}\alpha(x,t)$ (6) $\displaystyle\alpha(s(t),t)$
$\displaystyle=0$ (7) $\displaystyle-\partial_{x}\alpha(0,t)$
$\displaystyle=\theta(t),$ (8)

where $\alpha:[0,s(t)]\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ is
$\alpha(x,t)=T(x,t)-T_{m}$ and $s(t)=\Theta(t)$ is the unknown interface
represented as the moving boundary. The output is measured by the unknown
static map with input (5):

$y(t)=Q(\Theta(t)).$ (9)

The ES goal is to optimize an unknown static map Q($\cdot$) using a real-time
optimization control with optimal unknown output $y^{*}$ and optimizer
$\Theta^{*}$ as well as measurable output $y$ and input $\theta$.
Consequently, the control objectives of the Stefan problem are achieved, i.e.,
$\lim\limits_{t\rightarrow\infty}s(t)=s^{*}$ and
$\lim\limits_{t\rightarrow\infty}T(x,t)=T_{m}\,,\forall x\in[0,s^{*}]$, as
illustrated in Figure 3.

Figure 3: Control objective of the Stefan problem. We aim to design a heat
flux input $q_{c}(t)$ such that the interface position $s(t)$ is driven to the
setpoint position $s^{*}$.

The unknown nonlinear map is locally quadratic, such that

$Q(\Theta)=y^{*}+\dfrac{H}{2}(\Theta-\Theta^{*})^{2},$ (10)

where $\Theta^{*},y^{*}\in\mathbb{R}$ and $H<0$ is the Hessian. Hence, the
output of the static map is given by

$y(t)=y^{*}+\dfrac{H}{2}(\Theta(t)-\Theta^{*})^{2}.$ (11)

###  3.1 Demodulation Signals

The demodulation signal $N(t)$ which is used to estimate the Hessian of the
static map by multiplying it with the output $y(t)$ of the static map is
defined in [6] as

$\hat{H}(t)=N(t)y(t)\;\;\text{with}\;\;N(t)=-\dfrac{8}{a^{2}}\cos{(2\omega
t)}$ (12)

whereas the signal $M(t)$ is used to estimate the gradient of the static map
as follows:

$G(t)=M(t)y(t)\;\;\text{with}\;\;M(t)=\frac{2}{a}\sin{(\omega t)}.$ (13)

###  3.2 Additive Probing Signal

The perturbation $S(t)$ is adapted from the basic ES to the case of PDE
actuation dynamics. The trajectory generation problem as in [28] is described
as follows:

$\displaystyle S(t)$ $\displaystyle\coloneqq-\partial_{x}\beta(0,t),\quad
x\in(0,s(t))$ (14) $\displaystyle\partial_{t}\beta(x,t)$
$\displaystyle=\partial_{xx}\beta(x,t)$ (15) $\displaystyle\beta(s(t),t)$
$\displaystyle=0$ (16) $\displaystyle\beta_{x}(s(t),t)$
$\displaystyle=-a\omega\cos{(\omega t)},$ (17)

where $\beta:[0,s(t)]\times\mathbb{R}_{+}\rightarrow\mathbb{R}$. The explicit
solution of (14) is found respectively for the reference trajectory and the
reference solution postulated by a power series [29]:

$\displaystyle\xi(t)=a\sin{(\omega t)}$ (18)
$\displaystyle\beta(x,t)=\sum_{i=0}^{\infty}\dfrac{a_{i}(t)}{i!}[x-\xi(t)]^{i}.$
(19)

We can calculate the first coefficients of the power series replacing the
boundary conditions (16) and (17) at (19), such that

$a_{0}(t)=0,\;\;\;\;\;a_{1}=-\dot{\xi}(t).$ (20)

The general expression $a_{i}(t)=\dot{a}_{i-2}(t)-a_{i-1}(t)\dot{\xi}(t)$ is
obtained by substituting (19) at (15). We give the analytic expression of the
first four coefficients of the series (19) so that one can see how the
successive derivatives of $\xi(t)$ appear

$\displaystyle a_{2}(t)=\dot{\xi}(t)^{2}$ (21) $\displaystyle
a_{3}(t)=\ddot{\xi}(t)-\dot{\xi}(t)^{3}$ (22) $\displaystyle
a_{4}(t)=\ddot{\xi}(t)^{2}+\ddot{\xi}(t)\dot{\xi}(t)+\dot{\xi}(t)^{4}.$ (23)

The trajectory generation solution which provide all terms of the power series
(19) is given by [30]

$\beta(x,t)=\sum_{i=0}^{\infty}\dfrac{1}{(2i)!}\dfrac{\partial^{i}}{\partial
t^{i}}[x-\xi(t)]^{2i}.$ (24)

Although (24) is not an explicit expression, choosing suitable values for $a$
and $\omega$ in (18), the series converges with few iterations of the infinite
sum, getting the desirable sinusoid signal $\xi(t)$ in the output of the
integrator.

According to (14), we take the spatial derivative of (24) and substitute
$x=0$, therefore, we arrive at the final expression of

$S(t)=-\sum_{i=0}^{\infty}\dfrac{1}{(2i-1)!}\dfrac{\partial^{i}}{\partial
t^{i}}[-a\sin{(\omega t)}]^{2i-1}.$ (25)

###  3.3 Estimation Errors and PDE-Error Dynamics

Since our objective is to find $\Theta^{*}$, which corresponds to the optimal
unknown actuator $\theta(t)$, we introduce the following estimation errors

$\hat{\theta}(t)=\theta(t)-S(t),\;\;\;\hat{\Theta}(t)=\Theta(t)-a\sin{(\omega
t)},$ (26)
$\tilde{\theta}(t)\coloneqq\hat{\theta}(t)-\Theta^{*},\;\;\;\vartheta(t)\coloneqq\hat{\Theta}(t)-\Theta^{*},$
(27)

reminding that $\Theta(t)\coloneqq s(t)$. Combining $\hat{\Theta}(t)$ in (26)
and (27) we get the relation between the propagated estimation error
$\vartheta(t)$, the propagated input $\Theta(t)$ and the optimizer of the
static map $\Theta^{*}$

$\Theta(t)-\Theta^{*}=\vartheta(t)+a\sin{(\omega t)}.$ (28)

Let us define

$u(x,t)=\alpha(x,t)-\beta(x,t),$ (29) $\hat{\theta}(t)=U(t).$ (30)

By (5)-(8) and (14)-(17) with the help of (26) and (27), we have our original
system:

$\displaystyle\dot{\vartheta}(t)$ $\displaystyle=-u_{x}(s(t),t),\quad
x\in(0,s(t))$ (31) $\displaystyle u_{t}(x,t)$ $\displaystyle=u_{xx}(x,t)$ (32)
$\displaystyle u(s(t),t)$ $\displaystyle=0$ (33) $\displaystyle-u_{x}(0,t)$
$\displaystyle=U(t).$ (34)

##  4 Control Design

###  4.1 Stefan Compensation

We consider the PDE-ODE cascade (31)-(34) and use the backstepping
transformation

$\begin{split}w(x,t)&=u(x,t)-\bar{K}\int_{x}^{s(t)}(x-\sigma)u(\sigma,t)\,dy\\\
&-\bar{K}(x-s(t))\vartheta(t)\end{split}$ (35)

with $\bar{K}>0$ is an arbitrary controller gain. Equation (35) transforms
(31)-(34) into the target system:

$\displaystyle\dot{\vartheta}(t)$
$\displaystyle=-\bar{K}\vartheta(t)-w_{x}(s(t),t),\,\,\,\,x\in(0,s(t))$ (36)
$\displaystyle w_{t}(x,t)$
$\displaystyle=w_{xx}(x,t)+\bar{K}\dot{s}(t)\vartheta(t)$ (37) $\displaystyle
w_{x}(0,t)$ $\displaystyle=0$ (38) $\displaystyle w(s(t),t)$
$\displaystyle=0.$ (39)

The compensation controller can be obtained by taking the derivative of (35)
with respect to $t$ and $x$ respectively along of the solution (31)-(34) and
substituting $x=0$:

$U(t)=-\bar{K}\left(\vartheta(t)+\int_{0}^{s(t)}u(x,t)\,dx\right).$ (40)

###  4.2 Implementable Extremum Seeking Control Law

Since we have no measurement on $\vartheta(t)$, (40) is not applicable
directly. Thus, introducing a result of [6], the average version of the
gradient and Hessian estimates are calculated by

$G_{\rm av}(t)=H\vartheta_{\rm av}(t),\;\;\;\hat{H}_{\rm av}(t)=H.$ (41)

Averaging (40), choosing $\bar{K}=KH$ with $K<0$ and plugging the average
gradient and Hessian estimates (41), we obtain

$U_{\rm av}(t)=-KG_{\rm av}(t)-KH\int_{0}^{s_{\rm av}(t)}u_{\rm av}(x,t)\,dx.$
(42)

We introduce a low-pass filter to the controller with the purpose of applying
the average theorem for infinite-dimensional systems [31] in the following
stability proof, such that

$U(t)=\dfrac{c}{s+c}\Bigg{\\{}K\Bigg{[}G(t)+\hat{H}(t)\int_{0}^{s(t)}u(x,t)\,dx\Bigg{]}\Bigg{\\}},$
(43)

for $c>0$ sufficiently large.

Adapting the original scheme in [14] and combining (5)-(8) with the proposed
boundary control law (43), the closed-loop ES with actuation Stefan PDE
dynamics is shown in Figure 4.

Figure 4: ES control loop applied to the one-phase Stefan problem.

##  5 Stability Analysis

The following theorem summarizes the stability properties for the average
version of the error-dynamics (31)-(34).

###### Theorem 1

Assume the model validity conditions $T_{\rm av}(x,t)>T_{m}$, $\dot{s}_{\rm
av}(t)>0$, $x\in(0,s_{\rm av}(t))$ and $s_{0}<s_{\rm av}(t)<s^{*}$ are
satisfied at least in the average sense, $\forall t\geq 0$, $s^{*}=\Theta^{*}$
and for initial conditions $(T_{\rm av}(x,0),s_{0})$ compatible with the
control law $U(t)$ in (43). Then, for a sufficiently large $c>0$, the average
version of the closed-loop system (31)-(34) is exponentially stable in the
sense of the norm
$\|u_{\rm{av}}(t)\|^{2}_{\mathcal{H}_{1}}+\left|\vartheta_{\rm{av}}(t)\right|^{2}$,
i.e.,

$\|u_{\rm av}(t)\|_{\mathcal{H}_{1}}^{2}+\left|\vartheta_{\rm
av}(t)\right|^{2}\leq M(\|u_{\rm av}(0)\|_{\mathcal{H}_{1}}^{2}+\vartheta_{\rm
av}(0)^{2})e^{-nt}\,,\quad\forall t\geq 0\,,$ (44)

and appropriate constants $M,n>0$.

Proof. The proof is carried out in Steps 1 to 3 below.

####  Step 1: Average closed-loop system

The average version of the system (31)-(34) is

$\displaystyle\dot{\vartheta}_{\rm av}(t)=-(u_{\rm av})_{x}(s_{\rm
av}(t),t),\,\,\,\,x\in(0,s_{\rm av}(t))$ (45) $\displaystyle(u_{\rm
av})_{t}(x,t)=(u_{\rm av})_{xx}(x,t)$ (46) $\displaystyle u_{\rm av}(s_{\rm
av}(t),t)=0$ (47) $\displaystyle\begin{split}&\dfrac{d}{dt}(u_{\rm
av})_{x}(0,t)=-c(u_{\rm av})_{x}(0,t)\\\ &-cKH\left[\vartheta_{\rm
av}(t)+\int_{0}^{s_{\rm av}(t)}u_{\rm av}(x,t)\,dx\right],\end{split}$ (48)

where the low-pass filter is represented in the state-space form. To derive
(48), we plug the relationships $\vartheta(t)+a\sin(\omega
t)=\Theta(t)-\Theta^{*}$, $G(t)=M(t)y(t)$ and (12) into (43). With the help of
the identities $2\sin^{2}(\omega t)=1-\cos({2\omega t})$, $2\sin({\omega
t})\cos({2\omega t})=\sin({3\omega t})-\sin({\omega t})$, $4\sin^{3}({\omega
t})=3\sin({\omega t})-\sin({3\omega t})$ and $4\sin^{2}(\omega t)\cos({2\omega
t})=2\cos({2\omega t})-\cos({4\omega t})-1$ and applying averaging, we arrive
at (48).

The backstepping transformation

$\begin{split}w(x,t)&=u_{\rm av}(x,t)-KH\int_{x}^{s_{\rm
av}(t)}(x-\sigma)u_{\rm av}(\sigma,t)\,d\sigma\\\ &-KH(x-s_{\rm
av}(t))\vartheta_{\rm av}(t)\end{split}$ (49)

maps the average error-dynamics (45)-(48) into the exponentially stable target
system after assuming $c$ $\rightarrow$ +$\infty$ for the sake of simplicity.
Consequently,

$\displaystyle\dot{\vartheta}_{\rm av}(t)=-KH\vartheta_{\rm
av}(t)-w_{x}(s_{\rm av}(t),t),\,\,\,\,x\in(0,s_{\rm av}(t))$ (50)
$\displaystyle w_{t}(x,t)=w_{xx}(x,t)+KH\dot{s}_{\rm av}(t)\vartheta_{\rm
av}(t)$ (51) $\displaystyle w_{x}(0,t)=0$ (52) $\displaystyle w(s_{\rm
av}(t),t)=0.$ (53)

####  Step 2: Inverse transformation

To ensure the equivalent stability property between the target system and the
original system, the invertibility of the transformation (49) needs to be
guaranteed. Assume the inverse transformation that maps (50)-(53) into
(45)-(48):

$\begin{split}u_{\rm av}(x,t)&=w(x,t)+\int_{x}^{s_{\rm
av}(t)}k(x-\sigma)w(\sigma,t)\,d\sigma\\\ &+\phi(x-s_{\rm
av}(t))\vartheta_{\rm av}(t),\end{split}$ (54)

where $k(x-\sigma)$ and $\phi(x-s_{\rm av}(t))$ are the kernel functions.
Taking the derivatives with respect to $t$ and $x$, respectively, along the
solution of (50)-(53), the functions $\phi(x)$ and $k(x-\sigma)$ must satisfy:

$\phi^{{}^{\prime\prime}}(x)=-KH\phi(x),\,\,\,\phi(0)=0,\,\,\,\phi^{{}^{\prime}}=KH\,,$
(55) $k(x-s_{\rm av}(t))=\phi(x-s_{\rm av}(t))\,,$ (56)
$\phi^{{}^{\prime}}(x-s_{\rm av}(t))=KH\left(1+\int_{x}^{s_{\rm
av}(t)}k(x-\sigma)\,d\sigma\right).$ (57)

The solutions of the gain kernel can be deduced from (55)-(57), such that

$\phi(x)=KH\sqrt{\dfrac{1}{KH}}\sin({\sqrt{KH}x})\,,$ (58)
$k(x-\sigma)=\phi(x-\sigma).$ (59)

Hence, replacing (58) and (59) into (54), we have the following inverse
transformation:

$\begin{split}&u_{\rm av}(x,t)=w(x,t)\\\ &+\int_{x}^{s_{\rm
av}(t)}KH\sqrt{\dfrac{1}{KH}}\sin({\sqrt{KH}(x-\sigma)})w(\sigma,t)\,d\sigma\\\
&+KH\sqrt{\dfrac{1}{KH}}\sin({\sqrt{KH}(x-s_{\rm av}(t))})\vartheta_{\rm
av}(t).\end{split}$ (60)

####  Step 3: Exponential stability

We prove the exponential stability of the average closed-loop system based on
the target system (45)-(48) using the Lyapunov method. We consider the
following Lyapunov functional:

$V=V_{1}+V_{2}+V_{3}\,,$ (61) $V_{1}=\dfrac{1}{2}\int_{0}^{s_{\rm
av}(t)}w(x,t)^{2}\,dx\,,$ (62) $V_{2}=\dfrac{1}{2}\int_{0}^{s_{\rm
av}(t)}w_{x}(x,t)^{2}\,dx\,,$ (63) $V_{3}=\rho\dfrac{1}{2}\vartheta_{\rm
av}(t)^{2}\,.$ (64)

Taking the derivative of (62) with respect to $t$:

$\begin{split}\dot{V_{1}}&=-\int_{0}^{s_{\rm av}(t)}w_{x}(x,t)^{2}\,dx\\\
&+KH\dot{s}_{\rm av}(t)\vartheta_{\rm av}(t)\int_{0}^{s_{\rm
av}(t)}w(x,t)\,dx.\end{split}$ (65)

Taking the derivative of (63) with respect to $t$:

$\begin{split}&\dot{V_{2}}=w_{x}(s_{\rm
av}(t),t)w_{t}(x,t)+\dfrac{1}{2}\dot{s}_{\rm av}(t)w_{x}(s_{\rm
av}(t),t)^{2}\\\ &\\!-\\!KH\dot{s}_{\rm av}(t)\vartheta_{\rm
av}(t)w_{x}(s_{\rm av}(t),t)-\int\displaylimits_{0}^{s_{\rm
av}(t)}w_{xx}(x,t)^{2}\,dx.\end{split}$ (66)

Using the relationship $w_{t}(s_{\rm av}(t),t)=-\dot{s}_{\rm
av}(t)w_{x}(s_{\rm av}(t),t)$ and replacing it into (66), we obtain

$\begin{split}\dot{V_{2}}&=-\int_{0}^{s_{\rm
av}(t)}w_{xx}(x,t)^{2}\,dx-\dfrac{1}{2}\dot{s}_{\rm av}(t)w_{x}(s_{\rm
av}(t),t)^{2}\\\ &-KH\dot{s}_{\rm av}(t)\vartheta_{\rm av}(t)w_{x}(s_{\rm
av}(t),t).\end{split}$ (67)

Taking the derivative of (64) with respect to $t$, lead us to

$\dot{V_{3}}=-\rho KH\vartheta_{\rm av}(t)^{2}-\rho\vartheta_{\rm
av}(t)w_{x}(s_{\rm av}(t),t).$ (68)

Substituting the terms (65), (67) and (68) into the time derivative of (61)
and using the Young’s inequality in $-\rho\vartheta_{\rm av}(t)w_{x}(s_{\rm
av}(t),t)$, $KH\dot{s}_{\rm av}(t)\vartheta_{\rm av}(t)\int_{0}^{s_{\rm
av}(t)}w(x,t)\,dx$ and $-KH\dot{s}_{\rm av}(t)\vartheta_{\rm
av}(t)w_{x}(s_{\rm av}(t),t)$, we have

$\begin{split}\dot{V}&\leq-\int_{0}^{s_{\rm
av}(t)}w_{xx}(x,t)^{2}\,dx-\int_{0}^{s_{\rm av}(t)}w_{x}(x,t)^{2}\,dx\\\
&-\dfrac{\rho KH}{2}\vartheta_{\rm av}(t)^{2}+\dfrac{\rho}{2KH}w_{x}(s_{\rm
av}(t),t)^{2}\\\ &+\dot{s}_{\rm av}(t)\left(\dfrac{s^{*}}{2}\int_{0}^{s_{\rm
av}(t)}w(x,t)^{2}+(KH)^{2}\vartheta_{\rm av}(t)^{2}\right).\end{split}$ (69)

By choosing $\rho=\dfrac{KH}{4s^{*}}$ and applying Poincaré and Agmon’s
inequalities at $\int_{0}^{s_{\rm av}(t)}w(x,t)^{2}\,dx$ and $w_{x}(s_{\rm
av}(t),t)^{2}$, respectively, we obtain

$\begin{split}\dot{V}\leq&-\dfrac{1}{8{s^{*}}^{2}}\int\displaylimits_{0}^{s_{\rm
av}(t)}w_{x}(x,t)^{2}dx-\dfrac{1}{4{s^{*}}^{2}}\int\displaylimits_{0}^{s_{\rm
av}(t)}w(x,t)^{2}dx\\\ &+\dot{s}_{\rm
av}(t)\left(\dfrac{s^{*}}{2}\int_{0}^{s_{\rm
av}(t)}\\!\\!\\!w(x,t)^{2}\,dx\\!+\\!(KH)^{2}\vartheta_{\rm
av}(t)^{2}\right)-\dfrac{\rho KH}{2}\vartheta_{\rm av}(t)^{2}\\\
\leq&-mV+n\dot{s}_{\rm av}(t)V,\end{split}$ (70)

where

$n=\max\,\,\\{1,\,8s^{*}KH\\},\,\,\,\,\,m=\min\,\\{1/4{s^{*}}^{2},\,KH\\}.$
(71)

The term $n\dot{s}_{\rm av}(t)V$ on the right-hand side of (70) does not let
us to directly conclude exponential stability. To deal with it, a new Lyapunov
function candidate $W$ is defined according to

$W(t)=V(t)e^{-ns_{\rm av}(t)}.$ (72)

The time derivative of (72) can be calculated using (70) such that

$\dot{W}(t)=(\dot{V}(t)-n\dot{s}_{\rm av}(t)V(t))e^{-ns_{\rm av}(t)}\leq-
mW(t).$ (73)

Taking into account (61), we can establish the following relationship:

$\|w(t)\|_{\mathcal{H}_{1}}^{2}+\rho\vartheta_{\rm
av}(t)^{2}\\!\leq\\!e^{ns^{*}}(\|w_{0}\|_{\mathcal{H}_{1}}^{2}+\rho\vartheta_{\rm
av}(0)^{2})e^{-mt}\,,\quad w(x,0)\\!=\\!w_{0}\,.$ (74)

Hence, we can conclude the existence of a positive constant $M>0$ using the
inverse transformation (49) combined with Young’s and Cauchy-Schwarz
inequalities, such that

$\|u_{\rm av}(t)\|_{\mathcal{H}_{1}}^{2}+\vartheta_{\rm
av}(t)^{2}\\!\leq\\!M(\|u_{\rm
av0}\|_{\mathcal{H}_{1}}^{2}\\!+\\!\vartheta_{\rm
av}(0)^{2})e^{-mt}\,,\quad\\!u_{\rm{av}}(x,0)\\!=\\!u_{\rm av0}\,,$ (75)

which completes the proof. $\square$

##  6 Asymptotic Convergence to a Neighborhood of the
Extremum Point

In Theorem 1, we prove that the average closed-loop system (45)-(48) is
exponentially stable. However, there is no suitable averaging theorem for
moving-boundary PDE systems and it remains as an open problem in the
literature.

If this theorem existed such that employed for PDEs of fixed domains [31],
then we would apply it to (31)-(34) and conclude for the non-average system
the existence of a unique exponentially stable periodic solution in $t$ of
period $\Pi\coloneqq 2\pi/\omega$, denoted by $\vartheta^{\Pi}(t)$,
$u^{\Pi}(x,t)$, satisfying

$\displaystyle{}\left(|\vartheta^{\Pi}(t)|^{2}+\|u^{\Pi}(t)\|^{2}+\|u_{x}^{\Pi}(t)\|^{2}\right)^{1/2}\leq\mathcal{O}(1/\omega)\,,\quad\forall
t\geq 0\,.$ (76)

On the other hand, the asymptotic convergence to a neighborhood of the
extremum point would be proved taking the absolute value of (28):

$|\Theta(t)-\Theta^{*}|=|\vartheta(t)+a\sin{(\omega t)}|,$ (77)

and writing (77) in terms of the periodic solution $\vartheta^{\Pi}(t)$:
$|\Theta(t)\\!-\\!\Theta^{*}|\\!=\\!|\vartheta(t)\\!-\\!\vartheta^{\Pi}(t)\\!+\\!\vartheta^{\Pi}(t)\\!+\\!a\sin{(\omega
t)}|$. By applying again the appropriate averaging theorem, one would have
$\vartheta(t)\\!-\\!\vartheta^{\Pi}(t)\\!\to\\!0$ exponentially and,
consequently,

$\limsup_{t\to\infty}|\Theta(t)\\!-\\!\Theta^{*}|=\limsup_{t\to\infty}|\vartheta^{\Pi}(t)+a\sin{(\omega
t)}|.$ (78)

Finally, with (76) we would arrive to

$\displaystyle\limsup_{t\to\infty}|\Theta(t)-\Theta^{*}|=\mathcal{O}\left(|a|+1/\omega\right).$
(79)

In order to show the convergence of the output $y(t)$, we could follow the
same steps employed for $\Theta(t)$ by plugging (28) into (11), such that

$\begin{split}\limsup_{t\to\infty}|y(t)-y^{*}|=\limsup_{t\to\infty}|H\vartheta^{2}(t)+Ha^{2}\sin{(\omega
t)}^{2}|.\end{split}$ (80)

Hence, rewriting (80) in terms of $\vartheta^{\Pi}(t)$ and again with the help
of (76), we finally get

$\displaystyle\limsup_{t\to\infty}|y(t)-y^{*}|=\mathcal{O}\left(|a|^{2}+1/\omega^{2}\right).$
(81)

##  7 Delay-Compensated Control

In this section, we formulate the Stefan problem for ES with delays in the
actuation dynamics. We have omitted it before since the formulation of the
Stefan problem for ES without delays itself has not been done before in the
literature yet either. Moreover, the delay would complicate the understanding
of the new contribution for the Stefan model. In this delayed scenario, we
show in the following the relation between the designed control law and a
state prediction used for delay compensation.

The diffusion equation of the temperature in the liquid-phase with actuator
delay represented by Figure 5 is described by:

$\displaystyle T_{t}(x,t)$ $\displaystyle=\alpha T_{xx}(x,t),\quad
x\in(0,s(t)),\quad\alpha=\dfrac{k}{\rho C_{p}}$ (82) $\displaystyle-
kT_{x}(0,t)$ $\displaystyle=q_{c}(t-D)$ (83) $\displaystyle T(s(t),t)$
$\displaystyle=T_{m}$ (84) $\displaystyle\dot{s}(t)$ $\displaystyle=-\beta
T_{x}(s(t),t),\quad\beta=\dfrac{k}{\rho\Delta H^{*}},$ (85)

where $D$ is the input delay.

Figure 5: Schematic of one-phase Stefan problem with actuator delay [32].

As assumed in Section 4.2, $\alpha=\beta=k=1$. Based on reference [32], the
control law $U(t)$ for our extremum seeking problem of the closed-loop system
consisting of the plant (82)-(85) would be:

$\begin{split}U(t)=-\bar{K}\bigg{(}\int_{t-D}^{t}U(\psi)\,d\psi+\int_{0}^{s(t)}u(x,t)\,dx+\vartheta(t)\bigg{)}.\end{split}$
(86)

In the next, we are going to represent equation (40) as $U_{d}(t)$ to differ
from (86). The objective of this section is to prove the equivalence between
the delay-compensated control (86) and the predictor-based feedback. In other
words, $U(t)\equiv U_{d}(t+D)$, which can be rewritten as

$U_{d}(t+D)=-\bar{K}\left(\int\displaylimits_{0}^{s(t+D)}u(x,t+D)\,dx+\vartheta(t+D)\right).$
(87)

Integrating $\dot{s}(t)=-u_{x}(s(t),t)$, from $t$ to $t+D$, yields

$s(t+D)=s(t)-\int_{t}^{t+D}u_{x}(s(\tau),\tau)\,d\tau.$ (88)

Now, integrating $u_{t}(x,t)=u_{tt}(x,t)$ in time, from $t$ to $t+D$, and in
space, from $0$ to $s(t+D)$, respectively, we get

$\begin{split}\int_{0}^{s(t+D)}&\int_{t}^{t+D}u_{t}(x,t)\,dt\,dx=\\\
&\int_{0}^{s(t+D)}\int_{t}^{t+D}u_{xx}(x,t)\,dt\,dx\end{split}$ (89)
$\begin{split}&\int_{0}^{s(t+D)}u(x,t+D)\,dx=\int_{0}^{s(t+D)}u(x,t)\,dx\\\
&+\int_{t}^{t+D}u_{x}(s(t+D),\tau)\,d\tau+\int_{t-D}^{t}U_{d}(\varphi)\,d\varphi.\end{split}$
(90)

Substituting (88) and (90) in (87), we obtain

$\begin{split}&U_{d}(t+D)=-\bar{K}\bigg{(}\int_{0}^{s(t+D)}u(x,t)\,dx\\\
&+\int_{t}^{t+D}(u_{x}(s(t+D),\tau)-u_{x}(s(\tau),\tau))\,d\tau\\\
&+\int_{t-D}^{t}U_{d}(\varphi)\,d\varphi+\vartheta(t)\bigg{)}.\end{split}$
(91)

Analyzing the second integral of the right-hand side in (91):

$\begin{split}&\int_{t}^{t+D}(u_{x}(s(t+D),\tau)-u_{x}(s(\tau),\tau))\,d\tau=\\\
&\int_{t}^{t+D}\int_{s(\tau)}^{s(t+D)}u_{xx}(x,\tau)\,dx\,d\tau=\\\
&\int_{s(t)}^{s(t+D)}(u(x,s^{-1}(x))-u(x,t))\,dx.\end{split}$ (92)

The boundary condition $u(s(t),t)=0,\,\forall t\geq 0$, implies that
$u(x,s^{-1}(x))=0$. By applying this condition to (92), one has

$\begin{split}&\int_{t}^{t+D}(u_{x}(s(t+D),\tau)-u_{x}(s(\tau),\tau))\,d\tau=\\\
&-\int_{s(t)}^{s(t+D)}u(x,t)\,dx.\end{split}$ (93)

Finally, replacing (93) into (91), we arrive at

$U_{d}(t)=-\bar{K}\left(\int\displaylimits_{t-D}^{t}U(\psi)\,d\psi+\int\displaylimits_{0}^{s(t)}u(x,t)\,dx+\vartheta(t)\right),$
(94)

which is the same as the delay-compensated control (86).

Comparing the prediction control law (94) with our original control law
without delay (40), we can verify that the additional expression on the right-
hand side $\int_{t-D}^{t}U(\psi)\,d\psi$ is similar to the integral of [14]
[Eq. 35] and it can be represented as a basic prediction scheme according to
Figure 6.

Figure 6: Extremum seeking control loop with basic prediction scheme applied
to the one-phase Stefan problem under actuation delay.

The dither signal $S(t)$ in Figure 6, using as reference (24), will be simply
an implementable and advanced version of (14), i.e.,

$S(t)=\beta_{x}(0,t+D).$ (95)

Using the same steps of Section 4 for (86), yields to the following
implementable control law with a low-pass filter:

$\begin{split}U(t)&=\dfrac{c}{s+c}\Bigg{\\{}K\Bigg{[}G(t)+\hat{H}(t)\int_{t-D}^{t}U(\psi)\,d\psi\\\
&+\hat{H}(t)\int_{0}^{s(t)}u(x,t)\,dx\Bigg{]}\Bigg{\\}}.\end{split}$ (96)

##  8 Simulations

The numerical simulation employs the quadratic map described in (10) and the
parameters are chosen as stated in Table 1. For the sake of completeness, we
have considered the general case containing a delay $D=0.5$ seconds in the
actuation dynamics. It is worth to mention that this amount of delay is long
enough to destabilized the closed-loop system, if not properly compensated.

Table 1: Simulation parameters | Symbol | Description | Value
---|---|---|---
| $K$ | controller gain | -0.1
Controller | $a$ | perturbation amplitude | 0.1
parameters | $c$ | pole of the low-pass filter [rad/s] | 10
| $\omega$ | perturbation frequency [rad/s] | 10
| $L$ | spatial domain | 1
| $\Theta^{*}$ | optimizer static map | 0.8
System | $y^{*}$ | optimal value static map | 4
parameters | $H$ | Hessian of the static map | -1
| $s_{0}$ | Initial interface [m] | 0.12
| $T_{0}$ | Initial temperature [${}^{\circ}C$] | 110
| $T_{m}$ | Melting temperature [${}^{\circ}C$] | 100
| $D$ | delay [s] | 0.5

Figure 7 corresponds to the numerical plot of the temperature profile for the
closed-loop system converging in a three-dimensional space (taking into
account the domain $L$ and the time $t$) to a close neighborhood of $T_{m}$.

Figure 8 shows the convergence of the moving boundary to the optimizer
$\Theta^{*}$. The sinusoidal movement of $s(t)$ would violate the usual
conditions for the Stefan problem that the temperature remains above or below
the melting temperature on the whole interval $[0,s(t)]$, forming a periodic
chain of liquid and solid. However, it is well known in physics that the phase
transition is not sharp but gradual. For a small perturbation amplitude, the
system will be operating in the “phase transition range” when it reaches the
periodic steady state at the extremum. A more detailed model is needed to
capture the system dynamics in that narrow phase transition range. In
addition, we could redesign the algorithm in order to introduce vanishing
probing signals and tapering off the perturbation after the extremum
neighborhood is achieved, as studied in [33], [34] and [35]. At last, Figures
9 and 10 show the convergence of the output $y(t)$ to $y^{*}$ and $U(t)$ to 0,
respectively.

Figure 7: The initial states $T(x,0)$ (green curve), $T(0,t)$ (red curve) and
the convergence of $T(s(t),t)$ (blue curve) to $T_{m}=100^{\circ}C$ in a
three-dimensional space for the PDE state $T(x,t)$. The blue curve shows the
expansion of the domain of the Stefan PDE.  Figure 8: Convergence of
$s(t):=\Theta(t)$ to a small neighborhood of $\Theta^{*}$. Figure 9:
Convergence of the output $y(t)$ to $y^{*}$. Figure 10: Convergence of the
control signal $U(t)$ to $0$.

##  9 Conclusions

The proposed approach manages to maximize the static map by searching the
extremum point even in the presence of a Stefan PDE with moving boundary,
possibly including actuator delays as well. Although the actuation dynamics
must be known, no information is assumed from the map parameters. The average
boundary control law to compensate the actuation dynamics employed the
backstepping methodology. Local exponential stability of the average system
was guaranteed and convergence to a small neighborhood of the extremum was
also indicated. An illustration of the benefits of the new extremum seeking
scheme for the Stefan PDE is presented using consistent simulation results.

## References

  * [1] M. Krstić, H.-H. Wang, Stability of extremum seeking feedback for general nonlinear dynamic systems, Automatica 36 (4) (2000) 595–601.
  * [2] M. Leblanc, Sur l’electrication des chemins de fer au moyen de courants alternatifs de frequence elevee, Revue Generale de l’Electricite 12 (8) (1922) 275–277.
  * [3] M. Krstic, Performance improvement and limitations in extremum seeking control, Systems & Control Letters 39 (2000) 313–326.
  * [4] K. Ariyur, M. Krstić, Real-time optimization by extremum-seeking control, Wiley-Interscience, 2003.
  * [5] C. Zhang, D. Arnold, N. Ghods, A. Siranosian, M. Krstic, Source seeking with nonholonomic unicycle without position measurement and with tuning of forward velocity, Systems & Control Letters 56 (2007) 245–252.
  * [6] A. Ghaffari, M. Krstić, D. Nesić, Multivariable Newton-based extremum seeking, Automatica 48 (2012) 1759–1767.
  * [7] C. Manzie, M. Krstić, Extremum seeking with stochastic perturbations, IEEE Transactions on Automatic Control 54 (3) (2009) 580–585.
  * [8] A. Scheinker, M. Krstic, Minimum-seeking for CLFs: universal semiglobally stabilizing feedback under unknown control directions, IEEE Transactions on Automatic Control 58 (2013) 1107–1122.
  * [9] A. Scheinker, M. Krstic, Extremum seeking with bounded update rates, Systems & Control Letters 63 (2014) 25–31.
  * [10] A. Ghaffari, M. Krstic, S. Seshagiri, Power optimization for photovoltaic microconverters using multivariable Newton-based extremum seeking, IEEE Transactions on Control Systems Technology 22 (2014) 2141–2149.
  * [11] P. Paz, T. R. Oliveira, A. V. Pino, A. P. Fontana, Model-free neuromuscular electrical stimulation by stochastic extremum seeking, IEEE Transactions on Control Systems Technology 28 (1) (2020) 238–253.
  * [12] R. W. S. Pessoa, F. Mendes, T. R. Oliveira, K. Oliveira-Esquerre, M. Krstic, Numerical optimization based on generalized extremum seeking for fast methane production by a modified ADM1, Journal of Process Control 84 (2019) 56–69.
  * [13] D. A. Wepener, J. D. le Roux, I. K. Craig, Extremum seeking control to optimize mineral recovery of a flotation circuit using peak air recovery, Journal of Process Control 129 (2023) 1–13.
  * [14] T. R. Oliveira, M. Krstić, D. Tsubakino, Extremum seeking for static maps with delays, IEEE Transactions on Automatic Control 62 (4) (2017) 1911–1926.
  * [15] T. R. Oliveira, M. Krstic, Inverse optimal extremum seeking under delays, Systems & Control Letters 177 (2023) 1–12.
  * [16] T. R. Oliveira, M. Krstic, Extremum Seeking through Delays and PDEs, SIAM, 2022\.
  * [17] T. R. Oliveira, M. Krstic, Extremum seeking boundary control for PDE–PDE cascades, Systems & Control Letters 155 (2021) 1–15.
  * [18] J. Wettlaufer, Heat flux at the ice-ocean interface, Journal of Geophysical Research: Oceans 96 (C4) (1991) 7215–7236.
  * [19] B. Petrus, J. Bentsman, B. G. Thomas, Enthalpy-based feedback control algorithms for the Stefan problem, 2012 IEEE 51st IEEE Conference on Decision and Control (CDC) (2012) 7037–7042. Maui, USA.
  * [20] J. Ng, S. Dubljevic, Optimal control of convection–diffusion process with time-varying spatial domain: Czochralski crystal growth, Journal of Process Control 21 (2011) 1361–1369.
  * [21] A. Maidi, J.-P. Corriou, Boundary geometric control of a linear Stefan problem, Journal of Process Control 24 (2014) 939–946.
  * [22] C. Zhang, A. Siranosian, M. Krstić, Extremum seeking for moderately unstable systems and for autonomous vehicle target tracking without position measurements, Automatica 43 (10) (2007) 1832–1839.
  * [23] D. Bresch-Pietri, M. Krstić, Delay-adaptive predictor feedback for systems with unknown long actuator delay, IEEE Transactions on Automatic Control 55 (9) (2010) 2106–2112.
  * [24] M. L. Galvao, T. R. Oliveira, M. Krstic, Extremum seeking for Stefan PDE with moving boundary and delays, 17th IFAC Workshop on Time Delay Systems (2022) 222–227. Montreal, Canada.
  * [25] H. K. Khalil, Nonlinear Systems, Prentice-Hall, 2002.
  * [26] S. Koga, M. Krstic, Materials Phase Change PDE Control & Estimation: From Additive Manufacturing to Polar Ice, Birkhaüser, 2020.
  * [27] S. Koga, M. Krstic, Control of Stefan system and applications: a tutorial review, Annual Reviews in Control, Robotics, and Autonomous Systems 5 (2022) 547–577.
  * [28] M. Krstić, A. Smyshlyaev, Boundary control of PDEs: A course on backstepping designs, SIAM, 2008.
  * [29] W. B. Dunbar, N. Petit, P. Rouchon, P. Martin, Motion planning for a nonlinear Stefan problem, ESAIM: Control, Optimisation and Calculus of Variations 9 (2003) (2003) 275–296.
  * [30] C. D. Hill, Parabolic equations in one space variable and the non-characteristic Cauchy problem, Communications on Pure and Applied Mathematics 20 (3) (1967) 619–633.
  * [31] J. Hale, S. V. Lunel, Averaging in infinite dimensions, The Journal of integral equations and applications 2 (4) (1990) 463–494.
  * [32] S. Koga, D. Bresch-Pietri, M. Krstić, Delay compensated control of the Stefan problem and robustness to delay mismatch, International Journal of Robust and Nonlinear Control 30 (6) (2020) 2304–2334.
  * [33] H.-B. Dürr, M. S. Stanković, C. Ebenbauer, K. H. Johansson, Lie bracket approximation of extremum seeking systems, Automatica 49 (6) (2013) 1538–1552.
  * [34] A. Scheinker, M. Krstić, Non-$C^{2}$ lie bracket averaging for nonsmooth extremum seekers, Journal of Dynamic Systems, Measurement, and Control 136 (1).
  * [35] L. Wang, S. Chen, K. Ma, On stability and application of extremum seeking control without steady-state oscillation, Automatica 68 (1) (2016) 18–26.