State Estimation of the Stefan PDE:
A Tutorial on Design and Applications to Polar Ice and Batteries

Shumon Koga, Miroslav Krstic

Abstract

The Stefan PDE system is a representative model for thermal phase change phenomena, such as melting and solidification, arising in numerous science and engineering processes. The mathematical description is given by a Partial Differential Equation (PDE) of the temperature distribution defined on a spatial interval with a moving boundary, where the boundary represents the liquid-solid interface and its dynamics are governed by an Ordinary Differential Equation (ODE). The PDE-ODE coupling at the boundary is nonlinear and creates a significant challenge for state estimation with provable convergence and robustness.

This tutorial article presents a state estimation method based on PDE backstepping for the Stefan system, using measurements only at the moving boundary. PDE backstepping observer design generates an observer gain by employing a Volterra transformation of the observer error state into a desirable target system, solving a Goursat-form PDE for the transformation’s kernel, and performing a Lyapunov analysis of the target observer error system.

The observer is applied to models of problems motivated by climate change and the need for renewable energy storage: a model of polar ice dynamics and a model of charging and discharging in lithium-ion batteries. The numerical results for polar ice demonstrate a robust performance of the designed estimator with respect to the unmodeled salinity effect in sea ice. The results for an electrochemical PDE model of a lithium-ion battery with a phase transition material show the elimination of more than 15 % error in State-of-Charge estimate within 5 minutes even in the presence of sensor noise.

keywords:

Stefan system, state estimation, distributed parameter systems, backstepping, nonlinear observer, sea ice, lithium-ion batteries.

^†^†journal: Annual Reviews in Control

1 Introduction

1.1 Phase transitions

Ice melts into water in a hot environment. Conversely, water freezes into ice in a cold environment. These liquid-solid change phenomena are called ”phase changes.” In addition to various quotidian settings, they also arise, for example, in sea ice in the polar regions Maykut and Untersteiner (1971), solidification of molten metal in casting Meng and Thomas (2003), extrusion in polymer 3D-printing Valkenaers et al. (2013), laser sintering for metal 3D-printing Agarwala et al. (1995), cryosurgery for cancer treatment Rabin and Shitzer (1998), thermal energy storage in buildings Sharma et al. (2009), etc.

The phase change can be regarded as a growth of the domain of one phase, accompanied with a reduction of the domain of the other phase. Apart from the thermal phase change, a growth of a domain, and the consequent movement of the domain’s boundary, can be seen in some non-thermal chemical, biological, and social dynamics such as tumor growth in a patient’s body Spangler et al. (2016), axonal elongation for neurons’ signal transmission Diehl et al. (2014); Demir et al. (2021), lithiated region in electrodes of lithium-ion batteries Thomas et al. (2002), crystallization process Braatz (2002); Nagy and Braatz (2012); Fujiwara et al. (2005), domain walls in ferroelectric nanomaterials McGilly et al. (2015), spreading of invasive species to a new environment Du and Lin (2010), etc.

1.2 Stefan model

A representative mathematical dynamic model that contains moving boundaries is the “Stefan system”, which is a physical model of thermal phase changes. The associated problem of analyzing and finding the solutions to the Stefan model is referred to as the “Stefan problem”. Since the phase changes are caused by the temperature dynamics, the physical model involves the temperature of each phase which is distributed in space and changes in time. Hence, the mathematical formulation incorporates partial differential equations (PDEs) defined on a time-varying spatial domain, whereas the dynamics of the position of the moving boundary is governed by an ordinary differential equation (ODE) whose input depends on the PDE state. This configuration gives rise to nonlinear coupling of the PDE and ODE dynamics. As a result, though seemingly consisting of just a linear PDE and a linear scalar ODE, the Stefan problem is mathematically peculiar and not amendable to conventional analysis for PDEs and ODEs.

The Stefan problem is named after the Slovenian-Austrian physicist Josef (Jožef) Stefan, who was one of the most distinguished and influential physicists of the 19th century, celebrated for his numerous contributions to thermodynamics and heat transfer from the experimental perspective. Perhaps Stefan’s name is more recognized for the Stefan-Bolzman law, which revealed that a material with temperature $T$ in absolute unit emits a radiative heat transfer which is proportional to $T^{4}$ , through Stefan’s experimental work and his student Ludwig Boltzman’s work on the theoretical foundation.

After the publication of the thermal radiation law, Stefan started to focus of the thickness evolution of polar ice caps motivated by observed data of ice growth and air temperature acquired by British and German explorers during their expeditions. A long time before that, the phase change model by moving boundaries had been studied by Joseph Black in 1762. Franz Neumann developed the solution in his lectures around 1860. However, Neumann’s result had not been published until Weber’s paper in 1901. Stefan developed his analysis of the solution of ice growth and studied the correspondence with the empirical data, which was published in 1891 STEFAN (1891). Since then, the model has been known as the “Stefan problem”, and has been studied widely by researchers from the middle 1900s Crepeau (2007).

The mathematical and numerical analysis of the Stefan problem has been widely covered in the literature. The existence, uniqueness, monotone dependence, stability, and the differentiability of the one-phase Stefan problem have been studied by several references, see Friedman (1959); Cannon (1984); Fasano and Primicerio (1977); Kolodner (1956); Rubinšteĭn (2000). The existence and uniqueness of the classical solution of the two-phase Stefan problem were proven in Cannon and Primicerio (1971b, a) with the temperature boundary conditions and the flux boundary conditions, respectively. Several numerical methods to solve the Stefan problem were investigated by finite difference and finite element methods, see Kutluay et al. (1997); Bonnerot and Jamet (1977) for instance. The comparison of the numerical methods was studied in Javierre et al. (2006).

1.3 Control of Stefan system

The work on control of the Stefan system begun in the last two decades of the 1900s. For instance, feedback control of the Stefan system is developed in Hoffmann and Sprekels (1982) by an “on-off” switching design, and an inverse Stefan problem is studied in Zabaras et al. (1988) by using an integral method and in Pawlow (1987); Pawłow (1990); Zabaras et al. (1992); Kang and Zabaras (1995) by an optimization approach. An optimal feedback control is developed in Barbu et al. (1996) by convex analysis, and a feedback control for a thermostats modeled by a hyperbolic Stefan problem is given in Colli et al. (1999). Motion planning by boundary control has been solved in Dunbar et al. (2003a, b); Petit (2010), overcoming the challenges of the nonlinearity in the Stefan model. Optimal control for the Stefan problem has been developed in Hinze and Ziegenbalg (2007a, b) using a graph-based approach and in Bernauer and Herzog (2011); Alessandri et al. (2018) using a level-set formulation. Boundary geomertic control has been proposed in Maidi and Corriou (2014), and trajectory tracking control for the Vertical Gradient Freeze crystal growth process has been developed in Ecklebe et al. (2019, 2020). Numerous contributions on the feedback control of the Stefan problem with applications to continuous casting of steels have been achieved by Bentsman and coauthors, see Petrus et al. (2010, 2012, 2014, 2017); Chen et al. (2018, 2019, 2020).

1.4 Backstepping designs for Stefan system

Recently, boundary feedback controllers for the Stefan system have been designed via a “backstepping transformation” which has been used for many other classes of infinite-dimensional systems. The initial article Koga et al. (2018) introduces the designs of a state feedback control law, an observer, and an output feedback law for the one-phase Stefan system by proposing a nonlinear backstepping transformation for the moving-boundary Stefan PDE, and proves the exponential stability of the closed-loop system without imposing a priori assumption that the temperature state respects the phase constraints, but by proving that such constraints are actually maintained under proposed feedback. Extensions have been provided in the following articles: Koga et al. (2020) develops a control design with time-delay in the actuator and proved delay-robustness, Koga and Krstic (2020b) develops a control design for the two-phase Stefan problem, and Koga et al. (2021b) shows stability of the closed-loop system under sampled-data control. The backstepping controller and estimator for the Stefan system have been successfully applied to the screw extruder-based polymer 3-D printing Koga et al. (2020b), the laser sintering-based metal 3-D printing Koga et al. (2020), polar ice in the Arctic Koga and Krstic (2020a), lithium-ion batteries Koga et al. (2021a), and energy storage by paraffin with providing the experimental validation Koga et al. (2020a). The comprehensive materials can be found in our book Krstic and Koga (2020).

The state feedback control presented in aforementioned literature, reviewed in Koga and Krstic (2021), requires the entire profile of the temperature in the liquid phase as a given information. Some imaging-based sensors such as thermographic camera (a.k.a. infrared camera or IR camera) enables to capture the temperature profile, however, they include relatively high noise as compared to single point thermal sensors, such as thermocouples. Thus, estimating the entire temperature profile given a boundary measurement is a significant task for the implementation of the control algorithm.

1.5 State estimators for Stefan model

The problem of estimating variables of interest given some measured value is widely known as ”state estimation”. One of the most popular state estimation methods is the “Kalman filter” Kalman (1960) which is an optimal filter for linear dynamical systems with white Gaussian noise in the model and measurements. Another well known concept is the “Luenberger observer” Luenberger (1971), which stabilizes the estimation error at zero in linear deterministic systems. In finite dimensional systems, the observer gain is designed by means of pole placement or a linear matrix inequality Boyd et al. (1994).

While the control problem for the Stefan system has been developed as mentioned above, a state estimation for the Stefan system has been focused relatively few. In Petrus et al. (2017), an online recalibration method has been proposed for the state estimation of the Stefan problem, where the unknown parameters in the model are calibrated with the discrete-time measurements. In Pernsteiner et al. (2021), an state estimation for the latent heat storage system modeled by the Stefan system is developed by the Extended Kalman Filter (EKF) for the reduced-order model of the Stefan system through truncation.

In this tutorial article, we review state estimation of the Stefan system for the purpose of estimating unknown variables through available measurements. The design is employed by a bacsktepping observer, which is one type of Luenberger observer for the Stefan PDE system, where the observer gain is derived by a PDE backstepping method. After several design methods are introduced, we focus on their applications to polar ice and lithium-ion batteries.

2 Stefan System

Refer to caption — Figure 1: Schematic of one-dimensional one-phase Stefan problem. The temperature profile in the solid phase is assumed to be a uniform melting temperature.

We begin with a consideration of a pure one-component (liquid only) material in one dimension as depicted in Fig. 1, namely, with the assumption that the solid phase happens to be at the freezing temperature (and not cooler than freezing). This assumption is made for simplicity, so that our presentation can focus on controlling just the position of the liquid-solid boundary, and not also the temperature in the “distal” solid compartment.

The essence of the dynamics of the Stefan system is the evolution in time of the moving liquid-solid interface, which we denote by $s(t)$ . The $[0,s(t)]$ portion of the domain is liquid, at a varying temperature $T(x,t)$ , whereas the remainder of the domain $[0,L]$ , namely, the subdomain $[s(t),L]$ is solid, at the melting/freezing temperature $T_{m}$ , i.e., on the verge of becoming liquid.

Considering a material which, in its liquid phase, has density $\rho$ and heat capacity $C_{{\rm p}}$ , the local energy conservation law is given by

\displaystyle\rho C_{{\rm p}}T_{t}(x,t)=-q_{x}(x,t),\quad x\in(0,s(t))

(1)

where $q(x,t)$ is a heat flux profile and $T(x,t)$ is a temperature profile. Moreover, the local energy balance at the position of the liquid-solid interface $x=s(t)$ involved with the latent heat leads to the dynamics of the moving boundary

\displaystyle\rho\Delta H^{*}\dot{s}(t)=q(s(t),t).

(2)

The thermal conduction for a melting component obeys the well known Fourier’s Law

\displaystyle q(x,t)=-kT_{x}(x,t),\quad x\in[0,s(t)]

(3)

where $k$ is the thermal conductivity.

At the boundary $x=0$ , heat flux enters as an external source which can be manipulated as a controlled variable, denoted as $q_{c}(t)$ . The boundary condition at $x=s(t)$ must maintain the melting temperature $T_{\rm m}$ , which is the constant threshold level to cause the phase change from the solid to liquid under a static pressure.

By combining the energy conservation (1) and the thermal condition (3) with these two boundary conditions, the time evolution of the temperature profile in the material’s domain can be obtained by

$\displaystyle T_{t}(x,t)=$	$\displaystyle\alpha T_{xx}(x,t),\hskip 8.53581ptx\in(0,s(t)),$	(4)
$\displaystyle-kT_{x}(0,t)=$	$\displaystyle q_{{\mathrm{c}}}(t),$	(5)
$\displaystyle T(s(t),t)=$	$\displaystyle T_{{\mathrm{m}}},$	(6)

where $\alpha:=\frac{k}{\rho C_{{\rm p}}}$ . The initial condition is defined as an arbitrary spatial function for the temperature profile $T(x,0)=T_{0}(x)>T_{\rm m}$ , along with a positive valued interface position $s(0)=s_{0}$ .

If we disregard the dynamics of the moving boundary $s(t)$ , the linear PDE model (4)–(6) may deceive us into thinking that the Stefan model is linear. However, by combining the latent heat energy balance (2) and the thermal conduction (3), we arrive at the so called ”Stefan condition,” given by the following nonlinear ODE

\displaystyle\dot{s}(t)=-\beta T_{x}(s(t),t),

(7)

where $\beta:=\frac{k}{\rho\Delta H^{*}}$ .

There are two requirements for the validity of the model (4)-(7):

	$\displaystyle T(x,t)\geq$	$\displaystyle T_{{\rm m}},\quad\forall x\in(0,s(t)),\quad\forall t>0,$		(8)
	$\displaystyle 0<s(t)<$	$\displaystyle L,\quad\forall t>0.$		(9)

First, the trivial: the liquid phase is not frozen, i.e., the liquid temperature $T(x,t)$ is greater than the melting temperature $T_{\rm m}$ . Second, equally trivially, the material is not entirely in one phase, i.e., the interface remains inside the material’s domain. These physical conditions are also required for the existence and uniqueness of solutions Alexiades and Solomon (2018). Hence, we assume the following for the initial data.

Assumption 1

$0<s_{0}<L$ , $T_{0}(x)\in C^{1}([0,s_{0}];[T_{\rm m},+\infty))$ with $T_{0}(s_{0})=T_{\rm m}$ .

Lemma 1

With Assumption 1, if $q_{\rm c}(t)$ is a bounded piecewise continuous non-negative heat function, i.e.,

\displaystyle q_{{\rm c}}(t)\geq 0,\quad\forall t\geq 0,

(10)

then there exists a unique classical solution for the Stefan problem (4)–(7), which satisfies (8), and

\displaystyle\dot{s}(t)\geq 0,\quad\forall t\geq 0.

(11)

The definition of the classical solution of the Stefan problem is given in Appendix A of Koga et al. (2018). The proof of Lemma 1 is by maximum principle for parabolic PDEs and Hopf’s lemma, as shown in Gupta (2017).

3 State Estimator Design by PDE Backstepping

In this section, we focus on the state estimation of the Stefan system by a Luenberger-type observer approach with designing the observer gain via the backstepping method.

3.1 Estimation of temperature profile by measuring boundary temperature and interface position

3.1.1 Estimator design and main theorem

First, we develop the temperature profile estimation design under the available measurements of a boundary temperature and the interface position. This setup is practical, as the boundary temperature can be measured by a thermocouple and the interface position can be measured by a camera or some optical sensor. For the Stefan system in (4)–(7), the following observer is designed with the statement on the theorem.

Theorem 1

Consider the plant (4)–(7) with the measurements

\displaystyle Y_{1}(t)=s(t),\quad Y_{2}(t)=T(0,t),

(12)

and the following observer

$\displaystyle\hat{T}_{t}(x,t)=$	$\displaystyle\alpha\hat{T}_{xx}(x,t)+p_{1}(x,Y_{1}(t))\left(Y_{2}(t)-\hat{T}(0,t)\right),$	(13)
$\displaystyle\hat{T}_{x}(0,t)=$	$\displaystyle-\frac{q_{{\mathrm{c}}}(t)}{k}+p_{2}(Y_{1}(t))\left(Y_{2}(t)-\hat{T}(0,t)\right),$	(14)
$\displaystyle\hat{T}(Y_{1}(t),t)=$	$\displaystyle T_{{\mathrm{m}}},$	(15)

where $x\in[0,Y_{1}(t)]$ , and the observer gains are

	$\displaystyle p_{1}(x,Y_{1}(t))=$	$\displaystyle\lambda Y_{1}(t)(Y_{1}(t)-x)\frac{I_{2}\left(\sqrt{\frac{\lambda}{\alpha}\left\{Y_{1}(t)^{2}-(x-Y_{1}(t))^{2}\right\}}\right)}{Y_{1}(t)^{2}-(x-Y_{1}(t))^{2}},$		(16)
	$\displaystyle p_{2}(Y_{1}(t))=$	$\displaystyle-\frac{\lambda}{2\alpha}Y_{1}(t)$		(17)

with a gain parameter $\lambda>0$ . Assume that the model validity condition $T(x,t)\geq T_{m}$ is satisfied. Then, for all $\lambda>0$ , the observer error system is exponentially stable in the sense of the norm $||T-\hat{T}||_{{\cal H}_{1}}^{2}$ .

Let $\tilde{u}(x,t)=T(x,t)-\hat{T}(x,t)$ be an estimation error variable. Then, we have a system for error variable as

$\displaystyle\tilde{u}_{t}(x,t)=$	$\displaystyle\alpha\tilde{u}_{xx}(x,t)-p_{1}(x,s(t))\tilde{u}(0,t),\quad 0<x<s(t)$	(18)
$\displaystyle\tilde{u}_{x}(0,t)=$	$\displaystyle-p_{2}(x,s(t))\tilde{u}(0,t),$	(19)
$\displaystyle\tilde{u}(s(t),t)=$	$\displaystyle 0$	(20)

Due to the moving boundary nature, we cannot establish a good target system using the same form of the transformation. Instead, we first consider the observer design for the analogous system on the fixed domain and develop the observer gain with the associated backstepping transformation. After that, we apply the analogous gain and transformation on the moving boundary to the error system (18)–(20), and prove the stability of the associated target system on the moving boundary.

3.1.2 Fixed domain design

Consider the analogous estimation error system on the fixed domain $x\in(0,D)$ given by

$\displaystyle\tilde{u}_{t}(x,t)=$	$\displaystyle\alpha\tilde{u}_{xx}(x,t)-p_{1}(x,D)\tilde{u}(0,t),\quad 0<x<D$	(21)
$\displaystyle\tilde{u}_{x}(0,t)=$	$\displaystyle-p_{2}(D)\tilde{u}(0,t),$	(22)
$\displaystyle\tilde{u}(D,t)=$	$\displaystyle 0,$	(23)

Introduce the transformation

\displaystyle\tilde{u}(x,t)=w(x,t)+\int_{0}^{x}P(x,y)w(y,t)dy,

(24)

which transforms into

$\displaystyle w_{t}(x,t)=$	$\displaystyle\alpha w_{xx}(x,t)-\lambda w(x,t),$	(25)
$\displaystyle w_{x}(0,t)=$	$\displaystyle 0,$	(26)
$\displaystyle w(D,t)=$	$\displaystyle 0.$	(27)

Taking time and spatial derivatives of (24), the conditions for the gain kernel and the observer gain are obtained by

$\displaystyle P_{xx}(x,y)-P_{yy}(x,y)=$	$\displaystyle-\frac{\lambda}{\alpha}P(x,y),$	(28)
$\displaystyle P(x,x)=$	$\displaystyle-\frac{\lambda}{2\alpha}(x-D),$	(29)
$\displaystyle P(D,y)=$	$\displaystyle 0,$	(30)
$\displaystyle p_{1}(x,D)=$	$\displaystyle-\alpha P_{y}(x,0),$	(31)
$\displaystyle p_{2}(D)=$	$\displaystyle-P(0,0).$	(32)

The solution to the gain kernel PDE is derived as

\displaystyle P(x,y)=\lambda^{\prime}(D-x)\frac{I_{1}\left(\sqrt{\lambda^{\prime}\left\{(D-y)^{2}-(D-x)^{2}\right\}}\right)}{\sqrt{\lambda^{\prime}\left\{(D-y)^{2}-(D-x)^{2}\right\}}}.

(33)

By the differentiation formula for the Bessel functions, we have $\frac{d}{dz}\left(\frac{I_{1}(z)}{z}\right)=\frac{I_{2}(z)}{z}$ . Using this formula and some calculus, the observer gains in (31) and (32) are described by

	$\displaystyle p_{1}(x,D)$	$\displaystyle=\alpha\lambda^{\prime 2}D(D-x)\frac{I_{2}\left(z\right)}{z^{2}},\quad z=\sqrt{\lambda^{\prime}\left\{D^{2}-(D-x)^{2}\right\}},$		(34)
	$\displaystyle p_{2}(D)$	$\displaystyle=-\frac{\lambda}{2\alpha}D.$		(35)

Then, in the similar manner, we have

\displaystyle w(x,t)=\tilde{u}(x,t)+\int_{0}^{x}Q(x,y)\tilde{u}(y,t)dy,

(36)

which leads to the conditions of

$\displaystyle Q_{xx}(x,y)-Q_{yy}(x,y)=$	$\displaystyle\frac{\lambda}{\alpha}Q(x,y),$	(37)
$\displaystyle Q(x,x)=$	$\displaystyle\frac{\lambda}{2\alpha}(x-D),$	(38)
$\displaystyle Q(D,y)=$	$\displaystyle 0.$	(39)

The solution is

\displaystyle Q(x,y)=

\displaystyle P(x,y,-\lambda)=-\lambda^{\prime}(D-x)\frac{J_{1}\left(\sqrt{\lambda^{\prime}\left\{(D-y)^{2}-(D-x)^{2}\right\}}\right)}{\sqrt{\lambda^{\prime}\left\{(D-y)^{2}-(D-x)^{2}\right\}}}.

(40)

3.1.3 Analogous observer design on moving boundary domain

Referring to the result of fixed domain, we apply the backstepping observer design

$\displaystyle\hat{u}_{t}(x,t)=$	$\displaystyle\alpha\hat{u}_{xx}(x,t)+p_{1}(x,s(t))(u(0,t)-\hat{u}(0,t)),\quad 0<x<s(t)$	(41)
$\displaystyle\hat{u}(s(t),t)=$	$\displaystyle 0,$	(42)
$\displaystyle\hat{u}_{x}(0,t)=$	$\displaystyle-q_{c}(t)/k+p_{2}(s(t))(u(0,t)-\hat{u}(0,t)),$	(43)

with gains

	$\displaystyle p_{1}(x,s(t))=$	$\displaystyle\frac{\lambda^{2}}{\alpha}s(t)(x-s(t))\frac{I_{2}\left(\sqrt{\frac{\lambda}{\alpha}\left\{s(t)^{2}-(x-s(t))^{2}\right\}}\right)}{\sqrt{\frac{\lambda}{\alpha}\left\{s(t)^{2}-(x-s(t))^{2}\right\}}},$		(44)
	$\displaystyle p_{2}(s(t))=$	$\displaystyle-\frac{\lambda}{2\alpha}s(t).$		(45)

Now, we look at the original model in moving boundary coordinate. Consider the invertible transformation

	$\displaystyle\tilde{w}(x,t)=\tilde{u}(x,t)+\int_{0}^{x}Q(x-s(t),y-s(t))\tilde{u}(y,t)dy,$		(46)
	$\displaystyle\tilde{u}(x,t)=\tilde{w}(x,t)+\int_{0}^{x}P(x-s(t),y-s(t))\tilde{w}(y,t)dy.$		(47)

Then, the target system has the form of

$\displaystyle\tilde{w}_{t}(x,t)=$	$\displaystyle\alpha\tilde{w}_{xx}(x,t)-\lambda\tilde{w}(x,t)$
	$\displaystyle-\dot{s}(t)\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy,$	(48)
$\displaystyle\tilde{w}(s(t),t)=$	$\displaystyle 0,$	(49)
$\displaystyle\tilde{w}_{x}(0,t)=$	$\displaystyle 0,$	(50)

where $\bar{x}=x-s(t)$ , $\bar{y}=y-s(t)$ , and $q(\bar{x},\bar{y})=Q_{x}(x,y)+Q_{y}(x,y)$ .

We prove that the target $\tilde{w}$ -system in (3.1.3)–(50) is stable under $\dot{s}(t)>0$ and $s(t)<s_{r}$ . Consider the Lyapunov function

\displaystyle\tilde{V}_{1}=\frac{1}{2}||\tilde{w}||^{2}

(51)

The time derivative is given by

	$\displaystyle\dot{\tilde{V}}_{1}=$	$\displaystyle-\alpha\int_{0}^{s(t)}\tilde{w}_{x}(x,t)^{2}dx-\lambda\int_{0}^{s(t)}\tilde{w}(x,t)^{2}dx$
		$\displaystyle-\dot{s}(t)\int_{0}^{s(t)}\tilde{w}(x,t)\left(\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)dx.$		(52)

Define

	$\displaystyle\bar{q}=$	$\displaystyle\max_{(x,y)\in[0,s_{r}]}q(\bar{x},\bar{y})^{2},$		(53)
	$\displaystyle\bar{p}=$	$\displaystyle\max_{(y,z)\in[0,s_{r}]}P(\bar{y},\bar{z}).$		(54)

Applying Young’s Cauchy Schwarz inequalities to the second line of (52) with the help of $\dot{s}(t)>0$ and $s(t)\leq s_{r}$ , the following inequality is derived

\displaystyle\dot{\tilde{V}}_{1}

\displaystyle\leq-\alpha||\tilde{w}_{x}||^{2}-\lambda||\tilde{w}||^{2}+\frac{\dot{s}(t)}{2}\left(1+2\bar{q}s_{r}^{2}(1+\bar{p}^{2}s_{r}^{2})\right)||\tilde{w}||^{2}.

(55)

Consider

\displaystyle\tilde{V}_{2}=\frac{1}{2}\int_{0}^{s(t)}\tilde{w}_{x}(x,t)^{2}dx.

(56)

The time derivative is obtained by

\displaystyle\dot{\tilde{V}}_{2}

\displaystyle=-\alpha||\tilde{w}_{xx}||^{2}-\lambda||\tilde{w}_{x}||^{2}-\frac{\dot{s}(t)}{2}\tilde{w}_{x}(s(t),t)^{2}+\dot{s}(t)\int_{0}^{s(t)}\Phi(w(x,t),s(t),x),

(57)

where

\displaystyle\Phi:=

\displaystyle\int_{0}^{s(t)}\tilde{w}_{xx}(x,t)\left(\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)dx.

(58)

Doing integration by parts twice leads to

$\displaystyle\Phi=$	$\displaystyle\tilde{w}_{x}(s(t),t)\left(\int_{0}^{s(t)}q(0,\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)$
	$\displaystyle+\tilde{w}(0,t)\left(\frac{d}{dx}\left(\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)\right)\|_{x=0}$
	$\displaystyle+\int_{0}^{s(t)}\tilde{w}(x,t)\left(\frac{d^{2}}{dx^{2}}\left(\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)dx\right).$	(59)

We calculate

\displaystyle\frac{d}{dx}

\displaystyle\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy|_{x=0}=q(s(t),s(t))\tilde{w}(0,t)

(60)

Here, we see that $q(s(t),s(t))=Q_{x}(0,0)+Q_{y}(0,0)=\left(\frac{d}{dx}Q(x,x)\right)|_{x=0}=\frac{\lambda}{2\alpha}$ . Moreover,

	$\displaystyle\frac{d^{2}}{dx^{2}}\left(\int_{0}^{x}q(\bar{x},\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)$
$\displaystyle=$	$\displaystyle q(\bar{x},\bar{x})\tilde{w}_{x}(x,t)+(q(\bar{x},\bar{x})P(\bar{x},\bar{x})-2q_{\bar{x}}(\bar{x},\bar{x})-q_{\bar{y}}(\bar{x},\bar{x}))\tilde{w}(x,t)$
	$\displaystyle-\int_{0}^{x}\left((2q_{\bar{x}}(\bar{x},\bar{x})+q_{\bar{y}}(\bar{x},\bar{x}))P(\bar{x},\bar{z})+q(\bar{x},\bar{x})P_{\bar{x}}(\bar{x},\bar{z})+q_{\bar{x}}(\bar{x},\bar{z})\right)\tilde{w}(z,t)dz$
	$\displaystyle-\int_{0}^{x}q_{\bar{x}}(\bar{x},\bar{y})\left(\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy.$	(61)

In addition, we have

\displaystyle\int_{0}^{s(t)}q(\bar{x},\bar{x})\tilde{w}(x,t)\tilde{w}_{x}(x,t)dx=

\displaystyle-\frac{\lambda}{4\alpha}\tilde{w}(0,t)^{2}+\frac{1}{2}(q_{\bar{x}}(\bar{x},\bar{x})+q_{\bar{y}}(\bar{x},\bar{x}))\tilde{w}(x,t)^{2}dx.

(62)

Thus,

$\displaystyle\Phi=$	$\displaystyle\tilde{w}_{x}(s(t),t)\left(\int_{0}^{s(t)}q(0,\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)+\frac{\lambda}{4\alpha}\tilde{w}(0,t)^{2}$
	$\displaystyle+\int_{0}^{s(t)}(q(\bar{x},\bar{x})P(\bar{x},\bar{x})-\frac{3}{2}q_{\bar{x}}(\bar{x},\bar{x})-\frac{1}{2}q_{\bar{y}}(\bar{x},\bar{x}))\tilde{w}(x,t)^{2}dx$
	$\displaystyle+\int_{0}^{s(t)}\tilde{w}(x,t)I(\tilde{w}(x,t),x,s(t))dx,$	(63)

where

	$\displaystyle I(w(x,t),x,s(t))$
$\displaystyle=$	$\displaystyle-\int_{0}^{x}\left((2q_{\bar{x}}(\bar{x},\bar{x})+q_{\bar{y}}(\bar{x},\bar{x}))P(\bar{x},\bar{z})+q(\bar{x},\bar{x})P_{\bar{x}}(\bar{x},\bar{z})+q_{\bar{x}}(\bar{x},\bar{z})\right)\tilde{w}(z,t)dz$
	$\displaystyle-\int_{0}^{x}q_{\bar{x}}(\bar{x},\bar{y})\left(\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy.$	(64)

Applying Young’s and Cauchy-Schwarz inequality, we can show that there exist positive constants $M_{1}$ , $M_{2}$ , $M_{3}$ such that

$\displaystyle\tilde{w}_{x}(s(t),t)\left(\int_{0}^{s(t)}q(0,\bar{y})\left(\tilde{w}(y,t)+\int_{0}^{y}P(\bar{y},\bar{z})\tilde{w}(z,t)dz\right)dy\right)$
	$\displaystyle\leq\frac{1}{2}\tilde{w}_{x}(s(t),t)^{2}+M_{1}\|\|\tilde{w}\|\|^{2},$	(65)
$\displaystyle\int_{0}^{s(t)}(q(\bar{x},\bar{x})P(\bar{x},\bar{x})-\frac{3}{2}q_{\bar{x}}(\bar{x},\bar{x})-\frac{1}{2}q_{\bar{y}}(\bar{x},\bar{x}))\tilde{w}(x,t)^{2}dx\leq$	$\displaystyle M_{2}\|\|\tilde{w}\|\|^{2},$	(66)
$\displaystyle\int_{0}^{s(t)}\tilde{w}(x,t)I(\tilde{w}(x,t),x,s(t))dx\leq$	$\displaystyle M_{3}\|\|\tilde{w}\|\|^{2}$	(67)

Furthermore, by Agmon’s inequality, it holds $w(0,t)^{2}\leq 4s_{r}||\tilde{w}||^{2}$ . Therefore, applying all these inequalities to (63) leads to

\displaystyle\Phi\leq\frac{1}{2}\tilde{w}_{x}(s(t),t)^{2}+\frac{\lambda s_{r}}{\alpha}||\tilde{w}_{x}||^{2}+(M_{1}+M_{2}+M_{3})||\tilde{w}||^{2}.

(68)

Applying (68) to (57), we arrive at

\displaystyle\dot{\tilde{V}}_{2}\leq

\displaystyle-\alpha||\tilde{w}_{xx}||^{2}-\lambda||\tilde{w}_{x}||^{2}+\dot{s}(t)\left(\frac{\lambda s_{r}}{\alpha}||\tilde{w}_{x}||^{2}+(M_{1}+M_{2}+M_{3})||\tilde{w}||^{2}\right)

(69)

Thus, defining $\tilde{V}=\tilde{V}_{1}+\tilde{V}_{2}$ and $b=\frac{\alpha}{4s_{r}^{2}}+\lambda$ , we can see that there exists a positive constant $a>0$ such that the following inequality holds

\displaystyle\dot{\tilde{V}}\leq-b\tilde{V}+a\dot{s}(t)\tilde{V},

(70)

from which we conclude Theorem 1.

3.2 Estimation of Both Temperature Profile and Moving Interface by Measuring Only a Boundary Temperature

The most challenging setup for the state estimation of 1D one-phase Stefan problem is to estimate both temperature profile and moving interface position given a measured value of single boundary temperature. This is not only mathematically challenging, but also practically important in some industrial processes such as steel casting as developed in Petrus et al. (2012). We have not established a mathematically rigorous result for this problem, however, we suggest an observer design and investigate the performance in numerical simulation. We consider the plant (4)–(7) with the measurement

\displaystyle Y(t)=T(0,t),

(71)

and the following PDE observer

$\displaystyle\hat{T}_{t}(x,t)=$	$\displaystyle\alpha\hat{T}_{xx}(x,t)+p_{1}(x,\hat{s}(t))\left(Y(t)-\hat{T}(0,t)\right),$	(72)
$\displaystyle\hat{T}_{x}(0,t)=$	$\displaystyle-\frac{q_{{\mathrm{c}}}(t)}{k}+p_{2}(\hat{s}(t))\left(Y(t)-\hat{T}(0,t)\right),$	(73)
$\displaystyle\hat{T}(\hat{s}(t),t)=$	$\displaystyle T_{{\mathrm{m}}},$	(74)

with the ODE observer

\displaystyle\dot{\hat{s}}(t)=-\beta\hat{T}_{x}(\hat{s}(t),t)+l(Y(t)-\hat{T}(0,t))

(75)

where $x\in[0,Y_{1}(t)]$ , the observer gains $p_{1}$ , $p_{2}$ are given in (16), (17) with a gain parameter $\lambda>0$ , and $l>0$ is a gain parameter.

We study the performance of the observer in numerical simulation using parameters of zinc and tuning the gain parameters $\lambda$ and $l$ . Also, we compare with the observer design suggested in Petrus et al. (2012), which is given by a copy of the plant for PDE observer and the same structure in (75) for ODE observer. Fig. 4 depicts the simulation results and its comparison, as stated in its caption, which illustrates the better performance of the proposed estimation compared to the method in Petrus et al. (2012).

4 Sea Ice

4.1 Importance of the Arctic Sea Ice for Global Climate Modeling

The Arctic sea ice has been studied intensively in the field of climate and geoscience. One of the main reasons is due to ice-albedo feedback which influences climate dynamics through the high reflectivity of sea ice. The other reason is the rapid decline of the Arctic sea ice extent in the recent decade shown in several observations. These observations motivate the investigation of future sea ice amount. Several studies have developed a computational model of the Arctic sea ice and performed numerical simulations of the model with initial sea ice temperature profile. However, the spatially distributed temperature in sea ice is difficult to recover in realtime using a limited number of thermal sensors. Hence, the online estimation of the sea ice temperature profile based on some available measurements is crucial for the prediction of the sea ice thickness.

A thermodynamic model for the Arctic sea ice was firstly developed in Maykut and Untersteiner (1971) (hereafter MU71), in which the authors investigated the correspondence between the annual cycle pattern acquired from the simulation and empirical data of Untersteiner (1961). The model involves a temperature diffusion equation evolving on a spatial domain defined as the sea ice thickness. Due to melting or freezing phenomena, the aforementioned spatial domain is time-varying. Such a model is of the Stefan type Gupta (2017) and involves a PDE with a state-dependent moving boundary driven by a Neumann boundary value.

Refined models of MU71 have been suggested in literature. For instance, Semtner Jr (1976) proposed a numerical model to achieve faster and accurate computation of MU71 by discretizing the temperature profile into some layers and neglecting the salinity effect. An energy-conserving model of MU71 was introduced in Bitz and Lipscomb (1999) by taking into account an internal brine pocket melting on surface ablation and the vertically varying salinity profile. Their thermodynamic model was demonstrated by Bitz et al. (2001) using a global climate model with a Lagrangian ice thickness distribution. Combining these two models, Winton (2000) developed an energy-conserving three-layer model of sea ice by treating the upper half of the ice as a variable heat capacity layer.

Remote sensing techniques have been employed to obtain the Arctic sea ice data in several studies. In Hall et al. (2004), the authors suggested an algorithm to calculate sea ice surface temperature using the satellite measured brightness temperatures, which provided an excellent measurement of the actual surface temperature of the sea ice during the Arctic cold period. The Arctic sea ice thickness data were acquired in Kwok and Rothrock (2009) through a satellite called ”ICESat” during 2003-2008 and compared with the data in Rothrock et al. (2008) observed by a submarine during 1958-2000. More recent data describing the evolution of the sea ice thickness have been collected between 2010 and 2014 from the satellite called ”CryoSat-2” Kwok (2015).

On the other hand, state estimation has been studied as a specific type of data assimilation which utilizes the numerical model along with the measured value. For finite dimensional systems associated with noisy measurements, a well-known approach is the Kalman Filter. Another well-known method is the Luenberger type state observer, which reconstructs the state variable from partially measured variables. For the application to sea ice, Fenty and Heimbach (2013) developed an adjoint-based method as an iterative state and parameter estimation for the coupled sea ice-ocean in the Labrador Sea and Baffin Bay to minimize an uncertainty-weighted model-data misfit in a least-square sense as suggested in Wunsch and Heimbach (2007), using Massachusetts Institute of Technology general circulation model (MITgcm) developed in Marshall et al. (1997). In Fenty et al. (2017), the same methodology was applied to reconstruct the global ocean and ice concentration. Their sea ice model is based on the zero-layer approximation of the numerical model in Semtner Jr (1976), which is a crude model lacking internal heat storage and promoting fast melting.

4.2 Thermodynamic Model of Arctic Sea Ice

The thermodynamic model of MU71 describes the time evolution of the sea ice temperature profile in the vertical axis along with its thickness, which also evolves in time due to accumulation or ablation caused by energy balance.

Fig. 5 provides a schematic of the Arctic sea ice model. During the seasons other than summer (July and August), the sea ice is covered by snow, and the surface position of the snow also evolves in time. Let $T_{{\rm s}}(x,t)$ , $T_{{\rm i}}(x,t)$ denote the temperature profile of snow and sea ice, and $h(t)$ and $H(t)$ denote the thickness of snow and sea ice. The total incoming heat flux from the atmosphere is denoted by $F_{{\rm a}}$ , and the heat flux from the ocean is denoted by $F_{{\rm w}}$ . The Arctic sea ice model suggested by MU71 gives governing equations of a Stefan-type free boundary problem formulated as

$\displaystyle F_{{\rm a}}-I_{0}-\sigma(T_{{\rm s}}(-h(t)$	$\displaystyle,t)+273)^{4}+k_{{\rm s}}\frac{\partial T_{{\rm s}}}{\partial x}(-h(t),t)$
$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{ll}0,&{\rm if}\quad T_{{\rm s}}(-h(t),t)<T_{{\rm m}1},\\ -q\dot{h}(t),&{\rm if}\quad T_{{\rm s}}(-h(t),t)=T_{{\rm m}1},\end{array}\right.$	(78)
$\displaystyle\rho_{{\rm s}}c_{0}\frac{\partial T_{{\rm s}}}{\partial t}(x,t)=$	$\displaystyle k_{{\rm s}}\frac{\partial^{2}T_{{\rm s}}}{\partial x^{2}}(x,t),\quad\forall x\in(-h(t),0),$	(79)
$\displaystyle T_{{\rm s}}(0,t)=$	$\displaystyle T_{{\rm i}}(0,t),$	(80)
$\displaystyle k_{{\rm s}}\frac{\partial T_{{\rm s}}}{\partial x}(0,t)=$	$\displaystyle k_{0}\frac{\partial T_{{\rm i}}}{\partial x}(0,t),$	(81)
$\displaystyle\rho c_{{\rm i}}(T_{{\rm i}},S)\frac{\partial T_{{\rm i}}}{\partial t}(x,t)=$	$\displaystyle k_{{\rm i}}(T_{{\rm i}},S)\frac{\partial^{2}T_{{\rm i}}}{\partial x^{2}}(x,t)+I_{0}\kappa_{{\rm i}}e^{-\kappa_{{\rm i}}x},\quad\forall x\in(0,H(t)),$	(82)
$\displaystyle T_{{\rm i}}(H(t),t)=$	$\displaystyle T_{{\rm m}2},$	(83)
$\displaystyle q\dot{H}(t)=$	$\displaystyle k_{{\rm i}}\frac{\partial T_{{\rm i}}}{\partial x}(H(t),t)-F_{{\rm w}},$	(84)

where $I_{0}$ , $\sigma$ , $k_{{\rm s}}$ , $\rho_{{\rm s}}$ , $c_{0}$ , $k_{0}$ , $\rho$ , $T_{{\rm m}1}$ , $T_{{\rm m}2}$ , and $q$ are solar radiation penetrating the ice, Stefan-Boltzmann constant, thermal conductivity of snow, density of snow, heat capacity of pure ice, thermal conductivity of pure ice, density of pure ice, melting point of surface snow, melting point of bottom sea ice, and latent heat of fusion, respectively. The total heat flux from the air is given by

\displaystyle F_{{\rm a}}=(1-\alpha)F_{{\rm r}}+F_{{\rm L}}+F_{{\rm s}}+F_{{\rm l}},

(85)

where $F_{{\rm r}}$ , $F_{{\rm L}}$ , $F_{{\rm s}}$ , $F_{{\rm l}}$ , and $\alpha$ denote the incoming solar short-wave radiation, the long-wave radiation from the atmosphere and clouds, the flux of sensible heat, the latent heat in the adjacent air, and the surface albedo, respectively. The heat capacity and thermal conductivity of the sea ice are affected by the salinity as

\displaystyle c_{{\rm i}}(T_{{\rm i}},S(x))=c_{0}+\frac{\gamma_{1}S(x)}{T_{{\rm i}}(x,t)^{2}},\quad k_{{\rm i}}(T_{{\rm i}},S(x))=k_{0}+\frac{\gamma_{2}S(x)}{T_{{\rm i}}(x,t)},

(86)

where $S(x)$ denotes the salinity in the sea ice. $\gamma_{1}$ and $\gamma_{2}$ represent the weight parameters. The thermodynamic model (4.2)-(84) allows us to predict the future thickness $(h(t),H(t))$ and the temperature profile $(T_{{\rm s}},T_{{\rm i}})$ given the accurate initial data. However, from a practical point of view, it is not feasible to obtain a complete temperature profile due to a limited number of thermal sensors. To deal with the problem, the estimation algorithm is designed so that the state estimation converges to the actual state starting from an initial estimate.

4.3 Annual Cycle Simulation of Sea Ice Thickness

For the computation, we use boundary immobilization method and finite difference semi-discretization Kutluay et al. (1997) with 100-point mesh in space, and the resulting approximated ODEs are calculated by using MATLAB ode15 solver.

Input Parameters

The input parameters are taken from Maykut and Untersteiner (1971) in SI units and Table 2 shows the monthly averaged values of heat fluxes coming from the atmosphere for each month. Table 2 shows the physical parameters of snow and sea ice. Following Bitz and Lipscomb (1999), the salinity profile is described by

\displaystyle S(x)=A\left[1-{\rm cos}\left\{\pi\left(\frac{x}{H(t)}\right)^{\frac{n}{m+\frac{x}{H(t)}}}\right\}\right],

(87)

where $A=1.6$ , $n=0.407$ , and $m=0.573$ .

Table 1: Average monthly values for the energy fluxes.

Symbol	$F_{{\rm r}}$	$F_{{\rm L}}$	$F_{{\rm s}}$	$F_{{\rm l}}$	$\alpha$
Unit	W/m²	W/m²	W/m²	W/m²
Jan.	0	168	19.0	0	$\cdots$
Feb.	0	166	12.3	-0.323	$\cdots$
Mar.	30.7	166	11.6	-0.484	0.83
Apr.	160	187	4.68	-1.45	0.81
May.	286	244	-7.26	-7.43	0.82
Jun.	310	291	-6.30	-11.3	0.78
Jul.	220	308	-4.84	-10.3	0.64
Aug.	145	302	-6.46	-10.7	0.69
Sep.	59.7	266	-2.74	-6.30	0.84
Oct.	6.46	224	1.61	-3.07	0.85
Nov.	0	181	9.04	-0.161	$\cdots$
Dec.	0	176	12.8	-0.161	$\cdots$

Table 2: Physical parameters of snow and sea ice.

Symbol	Meaning	Unit	Value
$\rho_{{\rm s}}$	density (snow)	kg/m³	330
$k_{s}$	conductivity (snow)	W/m/^∘C	0.31
$\rho$	density (ice)	kg/m³	917
$c_{0}$	heat capacity (ice)	J/kg/^∘C	2110
$k_{0}$	conductivity (ice)	W/m/^∘C	2.034
$\gamma_{1}$	weight of heat capacity	kJ ^∘C/kg	18.0
$\gamma_{2}$	weight of conductivity	W/m	0.117
$I_{0}$	solar radiation	W/m²	1.59
$\kappa_{{\rm i}}$	penetration rate	/m	1.5
$T_{{\rm m}1}$	melting temperature of sea ice at surface	^∘C	-0.1
$T_{{\rm m}2}$	melting temperature of sea ice at bottom	^∘C	-1.8

Simulation Test of MU71

Using the given data, firstly the simulation of (4.2)-(84) is performed and showed in Fig. 6 to recover the evolution of $h(t)$ and $H(t)$ in the annual season as in Maykut and Untersteiner (1971). The dynamic behavior of the snow surface and the bottom of sea ice are shown in Fig. 6 (a), and the time evolution of the temperature profile in sea ice is illustrated in Fig. 6 (b). We can see that both of Fig. 6 (a) and (b) have a good agreement with the simulation results shown in Maykut and Untersteiner (1971).

4.4 Temperature Profile Estimation

In this section, we derive the estimation algorithm utilizing some available measurements and show the exponential convergence of the designed estimation to a simplified sea ice model. The ice thickness and surface temperature are measured in several studies Hall et al. (2004); Kwok and Rothrock (2009); Rothrock et al. (2008). It is indeed typical to check observability before observer design, at least for systems on a constant domain (see Moura et al. (2014) for instance). Here, we start with the observer design that is accompanied by a proof of exponential stability, which ensures the states’ detectability.

Simplification of the Model

For the sake of the design and stability proof, we give a simplification on the system (4.2)-(84). The effect of the salinity profile on the physical parameters is assumed to be sufficiently small so that it can be negligible, i.e. $S(x)=0$ . Therefore, the heat equation of the sea ice temperature (82) is rewritten as

\displaystyle\frac{\partial T_{{\rm i}}}{\partial t}(x,t)=

\displaystyle D_{{\rm i}}\frac{\partial^{2}T_{{\rm i}}}{\partial x^{2}}(x,t)+\bar{I}_{0}\kappa_{{\rm i}}e^{-\kappa_{{\rm i}}x},\hskip 2.84526pt\forall x\in(0,H(t)),

(88)

where the diffusion coefficient is defined as $D_{{\rm i}}=k_{0}/\rho c_{0}$ . Next, we impose the following assumptions.

Assumption 2

The thickness $H(t)$ is positive and upper bounded, i.e. there exists $\bar{H}>0$ such that $0<H(t)<\bar{H}$ , for all $t\geq 0$ .

Assumption 3

$\dot{H}(t)$ is bounded, i.e., there exists $M>0$ such that $|\dot{H}(t)|<M$ , for all $t\geq 0$ .

According to Kwok and Rothrock (2009), the observation data of the sea ice’s thickness from the 1950s to 2008 show that the maximum value including the uncertainty is less than 5[m]. Moreover, the largest variation of the thickness in a snow-covered season of a year essentially happens from December to March as an accumulation, and most of the literature shows at most 20 [cm] accumulation per month. Hence, conservatively it is plausible to set $\bar{H}=10$ [m], and $M=50[\textrm{cm/Month}]=1.9\times 10^{-7}\textrm{[m/s]}$ .

Mathematically, the existence of the classical solution of the simple Stefan problem given by (88) and (83)–(84) has been established in literature. We refer the readers to follow Gupta (2017) for the detailed explanation. The solution of the original sea ice model (4.2)–(84) has not been studied due to its high complexity.

Observer Structure

Suppose that the sea ice thickness and the ice surface temperature are obtained as measurements ${\mathcal{Y}}_{1}(t)$ and ${\mathcal{Y}}_{2}(t)$ , i.e.

\displaystyle{\mathcal{Y}}_{1}(t)=

\displaystyle H(t),\quad{\mathcal{Y}}_{2}(t)=T_{{\rm i}}(0,t).

(89)

The state estimate $\hat{T}_{{\rm i}}$ of the sea ice temperature is governed by a copy of the plant (88) and (83)-(84) plus the error injection of $H(t)$ , namely, as follows:

$\displaystyle\frac{\partial\hat{T}_{{\rm i}}}{\partial t}(x,t)=$	$\displaystyle D_{{\rm i}}\frac{\partial^{2}\hat{T}_{{\rm i}}}{\partial x^{2}}(x,t)+\bar{I}_{0}\kappa_{{\rm i}}e^{-\kappa_{{\rm i}}x}-p_{1}(x,t)\left({\mathcal{Y}}_{1}(t)-\hat{H}(t)\right),\hskip 5.69054pt\forall x\in(0,H(t))$	(90)
$\displaystyle\hat{T}_{{\rm i}}(0,t)=$	$\displaystyle{\mathcal{Y}}_{2}(t)-p_{2}(t)\left({\mathcal{Y}}_{1}(t)-\hat{H}(t)\right),$	(91)
$\displaystyle\hat{T}_{{\rm i}}(H(t),t)=$	$\displaystyle T_{{\rm m}2}-p_{3}(t)\left({\mathcal{Y}}_{1}(t)-\hat{H}(t)\right),$	(92)
$\displaystyle\dot{\hat{H}}(t)=p_{4}(t)$	$\displaystyle\left({\mathcal{Y}}_{1}(t)-\hat{H}(t)\right)+\beta\frac{\partial\hat{T}_{{\rm i}}}{\partial x}({\mathcal{Y}}_{1}(t),t)-\frac{F_{{\rm w}}}{q},$	(93)

where $\beta:=\frac{k_{{\rm i}}}{q}$ . Next, we define the estimation error states as

\displaystyle\tilde{T}(x,t):=-(T_{{\rm i}}(x,t)-\hat{T}_{{\rm i}}(x,t)),\quad\tilde{H}(t):=H(t)-\hat{H}(t),

(94)

where the negative sign is added to be consistent with the description developed in Section 3 for the liquid phase. Subtraction of the observer system (90)-(93) from the system (88) and (83)-(84) yields the estimation error system as

$\displaystyle\frac{\partial\tilde{T}}{\partial t}(x,t)=D_{{\rm i}}$	$\displaystyle\frac{\partial^{2}\tilde{T}}{\partial x^{2}}(x,t)-p_{1}(x,t)\tilde{H}(t),\hskip 2.84526pt\forall x\in(0,H(t))$	(95)
$\displaystyle\tilde{T}(0,t)=$	$\displaystyle-p_{2}(t)\tilde{H}(t),$	(96)
$\displaystyle\tilde{T}(H(t),t)=$	$\displaystyle-p_{3}(t)\tilde{H}(t),$	(97)
$\displaystyle\dot{\tilde{H}}(t)=$	$\displaystyle-p_{4}(t)\tilde{H}(t)-\beta\frac{\partial\tilde{T}}{\partial x}(H(t),t).$	(98)

Our goal is to design the observer gains $p_{1}(x,t)$ , $p_{2}(t)$ , $p_{3}(t)$ , $p_{4}(t)$ so that the temperature error $\tilde{T}$ converges to zero. The main theorem is stated as follows.

Theorem 2

Let Assumptions 2 and 3 hold. Consider the estimation error system (95)-(98) with the design of the observer gains

$\displaystyle p_{1}(x,t)=$	$\displaystyle\frac{c\lambda x}{\beta}\frac{I_{1}\left(z\right)}{z}+\left(\frac{\varepsilon H(t)}{D_{{\rm i}}}-\frac{3}{\beta}\right)\lambda^{2}x\frac{I_{2}\left(z\right)}{z^{2}}+\frac{\lambda^{3}x^{3}}{D_{{\rm i}}\beta}\frac{I_{3}\left(z\right)}{z^{3}},$	(99)
$\displaystyle p_{2}(t)=$	$\displaystyle 0,$	(100)
$\displaystyle p_{3}(t)=$	$\displaystyle-\frac{\lambda}{2\beta}H(t)-\varepsilon,$	(101)
$\displaystyle p_{4}(t)=$	$\displaystyle c-\frac{\lambda}{2}\left(1-\frac{\lambda H(t)^{2}}{8D_{{\rm i}}}\right)+\frac{\beta\lambda}{2D_{{\rm i}}}\varepsilon H(t),$	(102)

where $\lambda>0$ , $c>0$ , and $\varepsilon>0$ are positive free parameters, $z$ is defined by

\displaystyle z:=\sqrt{\bar{\lambda}(H(t)^{2}-x^{2})},

(103)

where $\bar{\lambda}:=\frac{\lambda}{D_{{\rm i}}}$ , and $I_{j}(\cdot)$ denotes the modified Bessel function of the $j$ -th kind. Then, there exist positive constants $c^{*}>0$ and $\tilde{M}>0$ such that, for all $c>c^{*}$ , the norm

\displaystyle\Phi(t):=\int_{0}^{H(t)}\tilde{T}(x,t)^{2}dx+\tilde{H}(t)^{2}

(104)

satisfies the following exponential decay

\displaystyle{\Phi}(t)\leq\tilde{M}{\Phi}(0)e^{-\min\{\lambda,c\}t},

(105)

namely, the origin of the estimation error system is exponentially stable in the spatial $L_{2}$ norm.

Remark 1

The observer gains (99)-(102) include the thickness $H(t)$ , so the gains are not precomputed offline, but are easily calculated online, along with the state estimation. Owing to the slow dynamics of the sea ice model, the computation time is much less than the time step size, which enables the real-time computation of the proposed observer.

Remark 2

The measurements (89) are assumed to be noiseless; however, in practice, the measured data accompany with some noise. Preferably the observer needs pre-filtering to deal with the noisy measurements.

To handle the discrete-time measurements in practice as in Petrus et al. (2017), the designed observer should be discretized in time such as Euler or Runge-Kutta methods so that the estimation can be computed at every sampling of the discrete-time measurements. The free parameters $\lambda$ , $c$ , and $\varepsilon$ have their physical units [1/s], [1/s], and [^∘C/m], respectively. Hence we can see the consistency of the physical units in the estimation error system (95)–(98) together with (99)–(102).

Gain Derivation via State Transformation

For the estimation error system (95)–(98), we apply the following invertible transformations:

	$\displaystyle\tilde{T}(x,t)=$	$\displaystyle w(x,t)-\int_{x}^{H(t)}q(x,y)w(y,t)dy-\psi(x,H(t))\tilde{H}(t),$		(106)
	$\displaystyle w(x,t)=$	$\displaystyle\tilde{T}(x,t)-\int_{x}^{H(t)}r(x,y)\tilde{T}(y,t)dy-\phi(x,H(t))\tilde{H}(t),$		(107)

which map the estimation error system (95)-(98) into the following target system:

$\displaystyle w_{t}(x,t)=$	$\displaystyle D_{{\rm i}}w_{xx}(x,t)-\lambda w(x,t)-\dot{H}(t)f(x,H(t))\tilde{H}(t),\hskip 5.69054pt\forall x\in(0,H(t))$	(108)
$\displaystyle w(0,t)=$	$\displaystyle 0,$	(109)
$\displaystyle w(H(t),t)=$	$\displaystyle\varepsilon\tilde{H}(t),$	(110)
$\displaystyle\dot{\tilde{H}}(t)=$	$\displaystyle-c\tilde{H}(t)-\beta w_{x}(H(t),t),$	(111)

where $f(x,H(t))$ is to be determined. Taking the first and second spatial derivatives of the transformation (106), we get

$\displaystyle\tilde{T}_{x}(x,t)=$	$\displaystyle w_{x}(x,t)+q(x,x)w(x,t)$
	$\displaystyle-\int_{x}^{H(t)}q_{x}(x,y)w(y,t)dy-\psi_{x}(x,H(t))\tilde{H}(t),$	(112)
$\displaystyle\tilde{T}_{xx}(x,t)=$	$\displaystyle w_{xx}(x,t)+q(x,x)w_{x}(x,t)+\left(q_{x}(x,x)+\frac{d}{dx}q(x,x)\right)w(x,t)$
	$\displaystyle-\int_{x}^{H(t)}q_{xx}(x,y)w(y,t)dy-\psi_{xx}(x,H(t))\tilde{H}(t).$	(113)

Next, taking the time derivative of (106) along the solution of the target system (108)–(111), using integration by parts, and substituting the boundary condition (110), we get

$\displaystyle\tilde{T}_{t}(x,t)=$	$\displaystyle D_{{\rm i}}w_{xx}(x,t)+D_{{\rm i}}q(x,x)w_{x}(x,t)-(\lambda+D_{{\rm i}}q_{y}(x,x))w(x,t)$
	$\displaystyle+(\beta\psi(x,H(t))-D_{{\rm i}}q(x,H(t)))w_{x}(H(t),t)$
	$\displaystyle+(D_{{\rm i}}\varepsilon q_{y}(x,H(t))+c\psi(x,H(t)))\tilde{H}(t)$
	$\displaystyle+\int_{x}^{H(t)}(\lambda q(x,y)-D_{{\rm i}}q_{yy}(x,y))w(y,t)dy$
	$\displaystyle-\dot{H}(t)\tilde{H}(t)\left(\varepsilon q(x,H(t))+\psi_{H}(x,H(t))\right.$
	$\displaystyle\left.+f(x,H(t))-\int_{x}^{H(t)}q(x,y)f(y,H(t))dy\right).$	(114)

Thus, by (4) and (114), we have

	$\displaystyle\tilde{T}_{t}(x,t)-D_{{\rm i}}\tilde{T}_{xx}(x,t)+p_{1}(x,t)\tilde{H}(t)$
$\displaystyle=$	$\displaystyle-\left(\lambda+2D_{{\rm i}}\frac{d}{dx}q(x,x)\right)w(x,t)$
	$\displaystyle+(\beta\psi(x,H(t))-D_{{\rm i}}q(x,H(t)))w_{x}(H(t),t)$
	$\displaystyle+\left(D_{{\rm i}}\varepsilon q_{y}(x,H(t))+D_{{\rm i}}\psi_{xx}(x,H(t))+c\psi(x,H(t))+p_{1}(x,t)\right)\tilde{H}(t)$
	$\displaystyle+\int_{x}^{H(t)}(\lambda q(x,y)+D_{{\rm i}}q_{xx}(x,y)-D_{{\rm i}}q_{yy}(x,y))w(y,t)dy$
	$\displaystyle-\dot{H}(t)\tilde{H}(t)\left(\varepsilon q(x,H(t))+\psi_{H}(x,H(t))\right.$
	$\displaystyle\left.+f(x,H(t))-\int_{x}^{H(t)}q(x,y)f(y,H(t))dy\right).$	(115)

Substituting $x=0$ and $x=H(t)$ into (106), we get

$\displaystyle\tilde{T}(0,t)+p_{2}(t)\tilde{H}(t)=$	$\displaystyle-\int_{0}^{H(t)}q(0,y)w(y,t)dy$
	$\displaystyle+(p_{2}(t)-\psi(0,H(t)))\tilde{H}(t),$	(116)
$\displaystyle\tilde{T}(H(t),t)+p_{3}(t)\tilde{H}(t)=$	$\displaystyle(\varepsilon-\psi(H,H)+p_{3}(t))\tilde{H}(t).$	(117)

Moreover, substituting $x=H(t)$ into (4) yields

		$\displaystyle\dot{\tilde{H}}(t)+p_{4}(t)\tilde{H}(t)+\beta\tilde{T}_{x}(H(t),t)$
	$\displaystyle=$	$\displaystyle(p_{4}(t)-c+\beta(\varepsilon q(H(t),H(t))-\psi_{x}(H(t),H(t))))\tilde{H}(t).$		(118)

Therefore, for the equations (95)–(98) to hold, the gain kernel functions must satisfy the following conditions:

$\displaystyle q_{xx}(x,y)-q_{yy}(x,y)=$	$\displaystyle-\bar{\lambda}q(x,y),$	(119)
$\displaystyle\frac{d}{dx}q(x,x)=$	$\displaystyle-\frac{\bar{\lambda}}{2},\quad q(0,y)=0,$	(120)
$\displaystyle\beta\psi(x,H(t))=$	$\displaystyle D_{{\rm i}}q(x,H(t)),$	(121)

and the observer gains must satisfy

$\displaystyle p_{1}(x,t)=$	$\displaystyle-D_{{\rm i}}(\varepsilon q_{y}(x,H(t))+\psi_{xx}(x,H))-c\psi(x,H),$	(122)
$\displaystyle p_{2}(t)=$	$\displaystyle\psi(0,H(t)),$	(123)
$\displaystyle p_{3}(t)=$	$\displaystyle\psi(H(t),H(t))-\varepsilon,$	(124)
$\displaystyle p_{4}(t)=$	$\displaystyle c-\beta(\varepsilon q(H(t),H(t))-\psi_{x}(H(t),H(t))),$	(125)

and the function $f(x,H(t))$ must satisfy

\displaystyle f(x,H)+\varepsilon q(x,H)+\psi_{H}(x,H)=\int_{x}^{H}q(x,y)f(y,H)dy.

(126)

The solutions to (119)–(121) are uniquely given by

	$\displaystyle q(x,y)=$	$\displaystyle-\bar{\lambda}x\frac{I_{1}\left(\sqrt{\bar{\lambda}(y^{2}-x^{2})}\right)}{\sqrt{\bar{\lambda}(y^{2}-x^{2})}},$		(127)
	$\displaystyle\psi(x,H(t))=$	$\displaystyle-\frac{\lambda}{\beta}x\frac{I_{1}\left(z\right)}{z},$		(128)

where $z$ is defined by (103). Then, using (127)–(128), the conditions (122)–(125) are led to the explicit formulations of the observer gains given as (99)–(102). In the similar manner, the conditions for the gain kernel functions of the inverse transformation (107) are given by

$\displaystyle r_{xx}(x,y)-r_{yy}(x,y)=$	$\displaystyle\bar{\lambda}r(x,y),$	(129)
$\displaystyle\frac{d}{dx}r(x,x)=$	$\displaystyle\frac{\bar{\lambda}}{2},\quad r(0,y)=0,$	(130)
$\displaystyle\beta\phi(x,H(t))=$	$\displaystyle D_{{\rm i}}r(x,H(t)),$	(131)

and, the function $f(x,H(t))$ is obtained by

\displaystyle f(x,H(t))=r(x,H(t))p_{3}(H(t))+\phi_{H}(x,H(t)).

(132)

The solutions to (129)–(131) are given by

\displaystyle r(x,y)=

\displaystyle\bar{\lambda}x\frac{J_{1}\left(\sqrt{\bar{\lambda}(y^{2}-x^{2})}\right)}{\sqrt{\bar{\lambda}(y^{2}-x^{2})}},\hskip 2.84526pt\phi(x,H)=\frac{\lambda}{\beta}x\frac{J_{1}\left(z\right)}{z},

(133)

where $J_{1}$ is Bessel function of the first kind. Using the solutions (133), the function $f(x,H(t))$ is obtained explicitly by (132), which also satisfies the condition (126). Hence, the transformation from $(\tilde{T},\tilde{H})$ to $(w,\tilde{H})$ is invertible.

Stability Analysis

We prove the exponential stability of the origin of the estimation error system (95)-(98) in the spatial $L_{2}$ norm. First, we show the exponential stability of the origin of the target system (108)-(111). We consider the following Lyapunov functional

\displaystyle V=\frac{1}{2}||w||^{2}+\frac{\varepsilon}{2\beta}\tilde{H}(t)^{2}.

(134)

Taking the time derivative of (134) together with the solution of (108)-(111) yields

	$\displaystyle\dot{V}=$	$\displaystyle-D_{{\rm i}}\|\|w_{x}\|\|^{2}-\lambda\|\|w\|\|^{2}-\frac{\varepsilon c}{\beta}\tilde{H}(t)^{2}+\frac{\dot{H}(t)}{2}\varepsilon^{2}\tilde{H}(t)^{2}$
		$\displaystyle-\dot{H}(t)\tilde{H}(t)\int_{0}^{H(t)}w(x,t)f(x,H(t))dx.$		(135)

Applying Young’s and Cauchy-Schwarz inequalities to the last term in (135) with the help of Assumption 3, and choosing the gain parameter $c$ to satisfy

\displaystyle c>\frac{\beta M^{2}\bar{f}}{\varepsilon\lambda}+\beta M\varepsilon,

(136)

one can obtain the following inequality:

\displaystyle\dot{V}\leq

\displaystyle-\min\{\lambda,c\}V.

(137)

Applying comparison principle to the differential inequality (137), we get

\displaystyle V(t)\leq V(0)e^{-\min\{\lambda,c\}t}.

(138)

Hence, the target system (108)-(111) is exponentially stable at the origin. Due to the invertibility of the transformations (106) and (107), there exist positive constants $\underline{M}>0$ and $\bar{M}>0$ such that for the norm $\Phi(t)$ defined in (104) the inequalities hold $\underline{M}\Phi(t)\leq V(t)\leq\bar{M}\Phi(t)$ . Hence, we obtain (105) by defining $\tilde{M}=\bar{M}/\underline{M}$ , which completes the proof of Theorem 2. Note that the designed backstepping observer achieves faster convergence with a possibility of causing overshoot since the overshoot coefficient $\bar{M}/\underline{M}$ is a monotonically increasing function in the observer gains’ parameters $(\lambda,c)$ .

While we have focused on the simplified PDE (88) to derive a rigorous proof of the proposed state estimation design (90)-(93) with observer gains given by (99)–(102), simulation studies are performed by applying the estimation design to the original thermodynamic model (4.2)-(84) including salinity.

4.5 Numerical Tests of the Sea Ice Estimation

Initial conditions

The simulation results of temperature estimation $\hat{T}_{{\rm i}}$ computed by (90)-(93) along with the available measurements obtained by the online calculation of (4.2)-(84) are shown in Fig. 7. Here the initial temperature profiles are formulated as

	$\displaystyle T_{{\rm s}}(x,0)=$	$\displaystyle\frac{k_{0}(T_{{\rm m}1}-T_{0})}{k_{s}H_{0}}x+T_{0},$		(139)
	$\displaystyle T_{{\rm i}}(x,0)=$	$\displaystyle\frac{T_{{\rm m}1}-T_{0}}{H_{0}}x+T_{0}+a\sin\left(\frac{4\pi x}{H_{0}}\right),$		(140)

where $T_{0}=T_{{\rm i}}(0,0)$ which is obtained by solving fourth order algebraic equation from (4.2) and the input data, and $a$ is set as $a=1$ [C^∘]. The estimated initial temperature is chosen as

\displaystyle\hat{T}_{{\rm i}}(x,0)=\frac{T_{{\rm m}1}-T_{0}}{H_{0}^{2}(1-2d)}(x^{2}-2dH_{0}x)+T_{0}

(141)

with setting $d=1/4$ . Hence, the initial temperature estimate is lower than the actual temperature. This initial condition satisfies the boundary conditions (91) and (92). The initial state of the estimated ice thickness $\hat{H}(0)$ is set as that of the true thickness, i.e., $\hat{H}(0)=H(0)$ , which is feasible because the thickness is actually measured.

Tuning method for gain parameters

The design parameters $(\lambda,c,\varepsilon)$ are selected as follows:
(i) Choose $\varepsilon\approx\beta$ for the norm (134) to be similarly weighted.
(ii) Select $\lambda$ to be the inverse of a desired time constant (i.e., the time at 63% decay of the estimation error is achieved): here we set as one day, leading to $\lambda\approx\frac{1}{24\times 3600}=1.2\cdot 10^{-5}$ .
(iii) Select $c$ sufficiently larger than $\lambda$ so that the decay rate $\min\{\lambda,c\}$ is not reduced and (136) is satisfied.
Finally, these parameters are varied around these reference values until we observe a smooth and sufficiently fast convergence. Throughout the simulation, we see that the minimum value of the time step size in ode solver is more than 1 minute, while the computation time of each time update is less than 0.1 seconds, which shows its real-time implementability as addressed in Remark 1.

Numerical simulation of state estimation

The contour plot of the simulation results of $T_{{\rm i}}(x,t)$ and $\hat{T}_{{\rm i}}(x,t)$ for open-loop estimation by setting all the observer gain to be zero is depicted in Fig. 7 (a), and those for the proposed estimation are depicted in Fig. 7 (b) and Fig. 8 (a)-(b) with observer gains (99)–(102), respectively, by using input data on January. For the proposed estimation, we fix the parameters of $c=$ 3.0 $\times$ 10^-5 and $\varepsilon=$ 1.0 $\times$ 10^-8, and use the parameter of $\lambda=$ 5.0 $\times$ 10^-6 in Fig. 7 (b), $\lambda=$ 1.0 $\times$ 10^-5 in Fig. 8 (a), and $\lambda=$ 5.0 $\times$ 10^-7 in Fig. 8 (b). The figures show that the backstepping observer gain makes the convergence speed of the estimation to the actual value approximately 5 to 10 times faster at every point in sea ice. As seen in Fig. 8, while the larger value of $\lambda$ makes the convergence speed faster, it causes more overshoot beyond the actual temperature. Hence, the tradeoff between the convergence speed and overshoot can be handled by tuning the gain parameter $\lambda$ appropriately, thereby the parameters used in (b) achieve the desired performance. The overshoot behavior is noted at the end of Section 4 from a theoretical perspective. Consequently, the stability properties stated in Theorem 2 for the simplified model can be observed in numerical results of the proposed estimation applied to the original model (4.2)-(84). To visualize the convergence of the estimated temperature profile used in 7 (b) more clearly, Fig. 9 illustrates the profiles of both true temperature (black solid) and estimated temperature (red dash) on January 1st to 3rd in (a)–(c), respectively. We observe that the estimated temperature profile becomes almost the same as the true temperature profile on January 3rd, which is two days after the estimation algorithm runs. Moreover, Fig. 9 (d) depicts the time evolution of $\tilde{H}(t)$ , which is an estimation error of the ice’s thickness. We observe that the error is “enlarged” from $\tilde{H}(0)=0$ due to the error of temperature profile, and returns to zero after the temperature profiles become almost indistinguishable on January 3rd, from which the necessity of the estimator of the ice’s thickness is ensured while the thickness is actually measured.

Finally, we have studied the robustness of the proposed observer by varying the parameters $D_{\rm i}$ , $\beta$ , and $F_{\rm w}$ in the observer (90)-(93) and the gains (99)–(102) to $D_{\rm i}(1+\delta_{1})$ , $\beta(1+\delta_{2})$ , and $F_{\rm w}(1+\delta_{3})$ with setting $\delta_{1}=0.3$ , $\delta_{2}=-0.3$ , and $\delta_{3}=0.4$ . Fig. 10 (a) shows the contour plots of estimated and true temperature profiles and Fig. 10 (b) shows the evolution of $\tilde{H}(t)$ . From both figures, we can see that the observer states converge and stay around the true states with a modest error after 5 days, which illustrates robust performance of the proposed observer under the parameters’ uncertainties.

sectionLithium-Ion Batteries

4.6 Battery Management Systems

Battery management is crucial for safe and efficient use of numerous kinds of electronics such as smartphones and laptops, and electric vehicles. Among several chemical materials used for electrodes of lithium-ion batteries, Lithium Iron Phosphate (LFP) has several attractive features as an active material in lithium-ion batteries such as thermal safety, high energy, and power density Padhi et al. (1997). LFP and other common active materials show unique charge-discharge characteristics due to an underlying crystallographic solid-solid phase transition. Electrochemical models for lithium-ion batteries with single phase materials do not allow to capture these unique characteristics and thus a mathematical description of phase transitions needs to be added to these models. Electrochemical models are of interest for the design of accurate estimation algorithms in battery management systems. Estimation algorithms based on these models provide visibility into operating regimes that induce degradation enabling a larger domain of operation, therefore, increasing the performance of the battery in terms of energy capacity, power capacity, and fast charge rates Chaturvedi et al. (2010); Perez et al. (2015). Electrochemical model-based estimation is challenging for several reasons. First, measurements of lithium concentrations outside specialized laboratory environments is impractical. Second, the concentration dynamics are governed by coupled and nonlinear partial differential algebraic equations (PDAE) Thomas et al. (2002). Finally, the only measurable quantities (voltage and current) are related to dynamic states through a nonlinear function.

Electrochemical models describe the relevant dynamic phenomena in lithium-ion cells: diffusion, intercalation and electrochemical kinematics (see Figure 11). These models predict accurately the internal states of the battery, however, their complexity renders a challenging problem for estimation algorithms. For this reason, most approaches develop estimation algorithms based on simplified models. Among the various simplified models, the single particle model (SPM) has been broadly used in the observer design problem, see Moura et al. (2014); Di Domenico et al. (2010); Wang et al. (2014); Dey et al. (2015); Perez and Moura (2015); Moura et al. (2016); Tang et al. (2017). The main characteristic of the SPM is the use of a single spherical particle to represent diffusion of lithium ions in the intercalation sites of the porous active materials in the electrodes.

LFP has been extensibility used in lithium ion cells due to its thermal stability, cost effectiveness, non-toxic nature, and long cycle life Padhi et al. (1997). An electrochemical model for LFP batteries was proposed in Srinivasan and Newman (2004b) based on a core-shell model, where the concentration at the core is assumed constant and diffusion is allowed for the phase in the shell. The LFP model with phase transition electrode was revisited in Zhang and White (2007) with a more complete core-shell model, allowing diffusion in both phases of an LiCoO₂ cathode.

The estimation problem for batteries with LFP electrodes has been relatively less studied; a particle filter was derived in Schwunk et al. (2013) and a Sequential Monte Carlo filter was derived in Li et al. (2014). The core-shell model proposed for phase transition electrodes is described by a parabolic PDE with a state-dependent moving boundary.

4.7 Electrochemical Model with Phase Change Electrode

The electrochemical model for lithium-ion cells with a phase transition material in the positive electrode follows Srinivasan and Newman (2004b). We restrict the problem to particular initial conditions of the concentration of lithium ions in the particles (i.e. intercalation sites) of the positive electrode and consider only discharge processes. The initial concentration of lithium ions in the particles of the positive electrode follows a core-shell configuration where the core has a constant distribution of lithium ions in a low concentration phase (the $\alpha$ phase), and the shell has a constant distribution of lithium ions in a high concentration phase (the $\beta$ phase). During discharge, the fluxes of lithium ions at the surface of the particles in the positive electrode are positive, thus, increasing the concentration of lithium ions in the shell and the phase boundary is moving to the center, i.e., a shrinking core process, as depicted in Figure 12.

Single Particle Model

The single particle model is a simple electrochemical model that accounts for some phenomena in lithium-ion cells. The main simplification in this model comes from the assumption that a single diffusion equation in an spherical particle can be used to model the diffusion of lithium ions in all the intercalation sites of the active material of each electrode. In the SPM, the ionic molar fluxes $j_{\mathrm{n},\pm}(t)$ on both electrodes are proportional the current density $I(t)$ applied to the cell

j_{\mathrm{n},\pm}(t)=\mp\frac{I(t)}{a_{\mathrm{s},\pm}FL_{\pm}},

(142)

where $a_{\mathrm{s},\pm}=3\epsilon_{\mathrm{s},\pm}/R_{\mathrm{p},\pm}$ is the interfacial area (per unit volume), $\epsilon_{\mathrm{s},\pm}$ is the volume fraction of active material in each electrode, $R_{\mathrm{p},\pm}$ is the averaged radius of the intercalation sites (particles) in the electrodes, $F$ is the Faraday constant, and $L_{\pm}$ is the thickness of each electrode. Throughout this section, the subscripts $+$ and $-$ indicate that the variable corresponds to the positive or negative particle. The concentration dynamics of lithium ions in the negative electrode (single phase) follow the Fick’s law for diffusion

\frac{\partial c_{\mathrm{s},-}}{\partial t}(r,t)=\frac{D_{\rm s,-}}{r^{2}}\frac{\partial}{\partial r}\left[r^{2}\frac{\partial c_{\mathrm{s},-}}{\partial r}(r,t)\right],

(143)

for $r\in(0,R_{\rm p,-})$ , $t>0$ with boundary conditions

	$\displaystyle\frac{\partial c_{\mathrm{s},-}}{\partial r}(0,t)$	$\displaystyle=0,$		(144)
	$\displaystyle D_{\mathrm{s},-}\frac{\partial c_{\mathrm{s},-}}{\partial r}(R_{\mathrm{p},-},t)$	$\displaystyle=-j_{\mathrm{n},-}(t),$		(145)

and initial condition $c_{\rm 0,-}\in\mathcal{C}(0,R_{\rm p,-})$ . Diffusion in the positive particle follows a core-shell model. In the core of the particle, i.e., for $r\in(0,r_{\rm p}(t))$ , lithium ions are in the $\alpha$ -phase. The concentration in the core is assumed to be constant and equal to the equilibrium value of the $\alpha$ -phase, i.e., $c_{s,+}(r)=c_{\mathrm{s},\alpha}$ for all $r\in(0,r_{\rm p}(t))$ . In the shell of the spherical particle, i.e. for $r\in(r_{\rm p}(t),R_{p,+})$ , the concentration of lithium ions is in $\beta$ -phase. The concentration dynamics of lithium-ions in the shell of the positive particle follows the Fick’s law for diffusion

\frac{\partial c_{\mathrm{s},+}}{\partial t}(r,t)=\frac{D_{\rm s,+}}{r^{2}}\frac{\partial}{\partial r}\left[r^{2}\frac{\partial c_{\mathrm{s},+}}{\partial r}(r,t)\right],

(146)

for $r\in(r_{\rm p}(t),R_{\rm p,+})$ with boundary conditions

	$\displaystyle c_{\mathrm{s},+}(r_{\rm p}(t),t)$	$\displaystyle=c_{\mathrm{s},\beta},$		(147)
	$\displaystyle D_{\mathrm{s},+}\frac{\partial c_{\mathrm{s},+}}{\partial r}(R_{\rm p,+},t)$	$\displaystyle=-j_{\mathrm{n},+}(t),$		(148)

and initial conditions $c_{0,+}\in\mathcal{C}(r_{\rm p}(0),R_{\rm p,+})$ . The time-evolution of the moving interface $r_{p}(t)$ is not given explicitly. Instead, mass balance at the moving interface yields the following state-dependent dynamics:

\displaystyle(c_{\mathrm{s},\beta}-c_{\mathrm{s},\alpha})\frac{dr_{\rm p}(t)}{dt}=-D_{\rm s,+}\frac{\partial c_{\mathrm{s},+}}{\partial r}(r_{\rm p}(t),t).

(149)

Overpotentials $\eta_{\pm}(t)$ are found by solving the nonlinear algebraic equation

	$\displaystyle j_{\mathrm{n},\pm}(t)$	$\displaystyle=\frac{i_{0,\pm}(t)}{F}\left[e^{\frac{\alpha_{\mathrm{a}}F}{RT}\eta_{\pm}(t)}-e^{-\frac{\alpha_{\mathrm{c}}F}{RT}\eta_{\pm}(t)}\right],$		(150)
	$\displaystyle i_{0,\pm}(t)$	$\displaystyle=Fk_{\pm}\left[{c}_{\mathrm{ss},\pm}(t)\right]^{\alpha_{\mathrm{c}}}\left[c_{\mathrm{e},0}\left(c_{\mathrm{s},\max,\pm}-c_{\mathrm{ss},\pm}(t)\right)\right]^{\alpha_{\mathrm{a}}},$		(151)

where ${c}_{\mathrm{ss},\pm}(t):={c}_{\mathrm{s},\pm}(R_{\rm p,\pm},t)$ . The electric potential in each electrode is given by

\phi_{\mathrm{s},\pm}(t)=\eta_{\pm}(t)+U_{\pm}(c_{\mathrm{ss},\pm}(t))+R_{\mathrm{f},\pm}Fj_{\mathrm{n},\pm}(t).

(152)

Finally, output voltage is computed as the difference between the electric potential in each electrode

\displaystyle V(t)=\phi_{\mathrm{s},+}(t)-\phi_{\mathrm{s},-}(t).

(153)

Equations (146) -(153) form a complete description of the single particle model with a phase transition electrode, and it provides the following property on the moving interface during the discharge process.

Remark 3

During the single discharge process, the current density $I(t)$ maintains positive, i.e. $I(t)>0$ for $\forall t>0$ . This current positivity ensures the moving interface being shrinking. Furthermore, the initial interface position is less than the cell radius. Hence,

	$\displaystyle\frac{dr_{\rm p}(t)}{dt}<0,$		(154)
	$\displaystyle 0\leq r_{\rm p}(t)<R_{\rm p,+}.$		(155)

Mass Conservation

In this model, the total amount of lithium ions is conserved. The mathematical description of this property is given in the following lemma.

Lemma 2

The total amount of lithium $n_{\rm Li}$ in solid phase ( moles per unit area ) defined as

n_{\rm Li}(t)=\epsilon_{\rm s,-}L_{-}\overline{c}_{\rm{s},-}(t)+\epsilon_{\rm s,+}L_{+}\overline{c}_{\rm{s},+}(t),

(156)

where $\overline{c}_{\rm{s},-}(t)$ and $\overline{c}_{\rm{s},+}(t)$ are the volumetric averages of the concentrations

	$\displaystyle\overline{c}_{\rm{s},-}(t)=$	$\displaystyle\frac{3}{R_{\rm{p},-}^{3}}\int_{0}^{R_{\rm{p},-}}c_{\rm{s},-}(r,t)r^{2}{\rm d}r,$		(157)
	$\displaystyle\overline{c}_{\rm{s},+}(t)=$	$\displaystyle\frac{3}{R_{\rm{p},+}^{3}}\int_{0}^{R_{\rm{p},+}}c_{\rm{s},+}(r,t)r^{2}{\rm d}r,$		(158)

is conserved, namely $dn_{\rm{Li}}(t)/dt=0$ .

Lemma 2 was derived in Klein et al. (2012) for electrodes with a single phase, and we can show that this result extends to electrodes with phase transition materials.

Proof:
In our problem formulation there is a single phase in the negative particle and there are two phases in the positive particle, i.e., $\alpha$ -phase in the core and $\beta$ -phase in the shell. The concentration in $\alpha$ -phase at the core is assumed to be constant (at its equilibrium value $c_{\mathrm{s},\alpha}$ ). Under these assumptions, the time derivative of (156) is given by

	$\displaystyle\frac{dn_{\rm Li}}{dt}(t)=$	$\displaystyle-a_{\rm{s},-}L_{-}j_{\rm{n},-}(t)-a_{\rm{s},+}L_{+}j_{\rm{n},+}(t)-\frac{3\epsilon_{\rm{s},+}L_{+}}{R_{\rm{p},+}^{3}}r_{\rm{p}}^{2}(t)$
		$\displaystyle\times\left[\frac{dr_{\rm{p}}}{dt}(t)\left[c_{\mathrm{s},\beta}-c_{\mathrm{s},\alpha}\right]+D_{\rm s,+}\frac{\partial c_{\mathrm{s},+}}{\partial r}(r_{\rm p}(t),t)\right].$		(159)

Hence, the molar flux equations in (142) and the dynamics of the moving interface in (149) lead to $dn_{\rm{Li}}(t)/dt=0$ . In a more general formulation introduced in Khandelwal et al. (2014, 2015), i.e. when both electrodes have multiple phase transitions not necessarily at the equilibrium, mass conservation of lithium ions is guaranteed with the following interface dynamics

\displaystyle\frac{dr^{\rm[a,b]}_{i}}{dt}(t)=

\displaystyle\frac{1}{c_{\rm{b}}-c_{\rm{a}}}\left[D_{\rm a}\frac{\partial c}{\partial r}(r_{i}^{\rm[a,b]}(t)^{-},t)-D_{\rm b}\frac{\partial c}{\partial r}(r_{i}^{\rm[a,b]}(t)^{+},t)\right],

(160)

where $r^{\rm[a,b]}_{i}$ is the interface radius between any two phases (phase a and phase b) in any electrode. Each phase has a distinct equilibrium $c_{\rm a}$ , $c_{\rm b}$ and diffusion coefficient $D_{a}$ , $D_{b}$ .

4.8 State-of-Charge Estimation

Now, a state estimation algorithm for concentration of lithium ions, in both negative and positive electrodes, is provided in this section from the single particle model. The state observer for the positive electrode is derived via the backstepping method for moving boundary PDEs, and the observer for the negative electrode is derived from the mass conservation property.

Observer for Phase Transition Positive Electrode

The state observer is a copy of the diffusion system (146)-(148) in the positive electrode together with output error injection

	$\displaystyle\frac{\partial\widehat{c}_{\mathrm{s},+}}{\partial t}(r,t)=$	$\displaystyle\frac{D_{\rm s,+}}{r^{2}}\frac{\partial}{\partial r}\left[r^{2}\frac{\partial\widehat{c}_{\mathrm{s},+}}{\partial r}(r,t)\right]$
		$\displaystyle+P(\widehat{r}_{\rm p}(t),r)\left[c_{\mathrm{ss,+}}(t)-\widehat{c}_{\mathrm{s,+}}(R_{\rm p,+},t)\right],$		(161)

for $r\in(\widehat{r}_{\rm p}(t),R_{\rm p,+})$ with boundary conditions

$\displaystyle\widehat{c}_{\mathrm{s},+}(\widehat{r_{\rm p}}(t),t)=$	$\displaystyle c_{\beta},$	(162)
$\displaystyle D_{\mathrm{s},+}\frac{\partial\widehat{c}_{\mathrm{s},+}}{\partial r}(R_{\rm p,+},t)=$	$\displaystyle-j_{\mathrm{n},+}(t)$
	$\displaystyle+Q(\widehat{r_{\rm p}}(t))\left[c_{\mathrm{ss,+}}(t)-\widehat{c}_{\mathrm{s,+}}(R_{\rm p,+},t)\right],$	(163)

and initial conditions $\widehat{c}_{0,+}\in\mathcal{L}^{2}(\widehat{r_{\rm p}}(0),R_{\rm p,+})$ and $\widehat{r_{\rm p}}(0)\in(0,R_{\rm p,+})$ . Observer gains are given by

	$\displaystyle P(\widehat{r_{\rm p}}(t),r)$	$\displaystyle=D_{s,+}\overline{\lambda}^{2}\frac{R_{\rm p,+}}{r}l(t)s(t)\frac{I_{2}\left(z(t)\right)}{z(t)},$		(164)
	$\displaystyle Q(\widehat{r_{\rm p}}(t))$	$\displaystyle=\frac{D_{s,+}}{R_{\rm p,+}}\left(\frac{\overline{\lambda}}{2}s(t)+1\right),$		(165)

where $I_{2}(\cdot)$ is a modified Bessel function of the second kind and

$\displaystyle\overline{\lambda}$	$\displaystyle=\frac{\lambda}{D_{s,+}},$	(166)
$\displaystyle s(t)$	$\displaystyle=R_{\rm p,+}-\widehat{r_{\rm p}}(t),\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ l(t)=r-\widehat{r_{\rm p}}(t),$	(167)
$\displaystyle z(t)$	$\displaystyle=\sqrt{\overline{\lambda}\left[s(t)^{2}-l(t)^{2}\right]}.$	(168)

The parameter $\lambda>0$ is designed to achieve faster convergence of the estimated concentration to true concentration. Moreover, the estimator for the moving interface position is given by the following dynamics:

	$\displaystyle(c_{\mathrm{s},\beta}-c_{\mathrm{s},\alpha})\frac{d\widehat{r_{\rm p}}(t)}{dt}=-\kappa\left[c_{\mathrm{ss,+}}(t)-\widehat{c}_{\mathrm{s,+}}(R_{\mathrm{p},+},t)\right]$
	$\displaystyle-D_{\rm s,+}\frac{\partial\widehat{c}_{\mathrm{s},+}}{\partial r}(\hat{r}_{\rm p}(t),t),$		(169)

where the parameter $\kappa>0$ is designed to achieve fast convergence of the estimated interface position to the true value.

The stability of the estimation error system is theoretically proven for the PDE observer (4.8)–(4.8) with gains (164), (165) under the assumption $\widehat{r_{\rm p}}(t)\equiv r_{\rm p}(t)$ for all $t\geq 0$ in the next section. As the moving interface position $r_{\rm p}(t)$ is unknown in practice, we construct the estimator (169), and use the estimated interface position $\widehat{r_{\rm p}}(t)$ in the gains (164), (165) of PDE observer.

The sign of the observer gain in (169) (first term in the right hand side) is determined based on the monotonic relation, namely, as the surface concentration $c_{\mathrm{ss,+}}(t)$ is increased the moving interface position $r_{\rm p}(t)$ is decreased. Physically, as the battery is discharged, the domain of the lithium rich $\beta$ -phase in the positive electrode is expanded from the outer region. Hence, the observer (169) is designed so that if the measured surface concentration is larger than the estimated surface concentration, the battery is discharged more than estimated, and the domain of $\beta$ -phase for the estimator is driven to be expanded.

Stability Analysis of the Estimation Error System with Known Interface Position

Let $\widetilde{c}_{\mathrm{s},+}(r,t)$ be an estimation error defined by $\widetilde{c}_{\mathrm{s},+}(r,t):=c_{\mathrm{s},+}(r,t)-\widehat{c}_{\mathrm{s},+}(r,t)$ . The stability analysis of the estimation error system is presented in the following theorem.

Theorem 3

Consider the plant PDE (146)–(148) and the PDE observer (4.8)–(4.8) with observer gains (164) and (165) under the properties of (154), (155), and the assumption $\widehat{r_{\rm p}}(t)\equiv r_{\rm p}(t)$ for all $t\geq 0$ . Then, for any initial estimation error $\widetilde{c_{\mathrm{s},+}}(r,0)$ , the estimation error is exponentially stable at the origin in the sense of the norm

\displaystyle\int_{r_{\rm p}(t)}^{R_{\rm p,+}}r^{2}\widetilde{c}_{\mathrm{s},+}(r,t)^{2}{\rm d}r.

(170)

Note that subtracting (4.8)-(4.8) from (146)-(148) under $\widehat{r_{\rm p}}(t)\equiv r_{\rm p}(t)$ yields the estimation error dynamics

$\displaystyle\frac{\partial\widetilde{c}_{\mathrm{s},+}}{\partial t}(r,t)=$	$\displaystyle\frac{D_{\rm s,+}}{r^{2}}\frac{\partial}{\partial r}\left[r^{2}\frac{\partial\widetilde{c}_{\mathrm{s},+}}{\partial r}(r,t)\right]-P(r_{\rm p}(t),r)\widetilde{c}_{\mathrm{s},+}(R_{\rm p,+},t),$	(171)
$\displaystyle\widetilde{c}_{\mathrm{s},+}(r_{\rm p}(t),t)=$	$\displaystyle 0,$	(172)
$\displaystyle D_{\mathrm{s},+}\frac{\partial\widetilde{c}_{\mathrm{s},+}}{\partial r}(R_{\rm p,+},t)=$	$\displaystyle-Q\left(r_{\mathrm{p}}(t)\right)\widetilde{c}_{\mathrm{s,+}}(R_{\rm p,+},t).$	(173)

Change of coordinate: First, we introduce the following change of coordinate and state variable to simplify the structure of the estimation error dynamics in a cartesian coordinate:

$\displaystyle x$	$\displaystyle=R_{\rm p,+}-r,$	(174)
$\displaystyle\widetilde{u}(x,t)$	$\displaystyle=r\widetilde{c}_{\mathrm{s},+}(r,t),$	(175)
$\displaystyle s(t)$	$\displaystyle=R_{\rm p,+}-\widehat{r_{\rm p}}(t).$	(176)

The estimation error dynamics (171)-(173) is rewritten by the new coordinate and state as

$\displaystyle\frac{\partial\widetilde{u}}{\partial t}(x,t)$	$\displaystyle=D_{\rm s,+}\frac{\partial^{2}\widetilde{u}}{\partial x^{2}}(x,t)-\overline{P}(s(t),x)\widetilde{u}(0,t),$	(177)
$\displaystyle\widetilde{u}(s(t),t)$	$\displaystyle=0,$	(178)
$\displaystyle\frac{\partial\widetilde{u}}{\partial x}(0,t)$	$\displaystyle=-\overline{Q}(s(t))\widetilde{u}(0,t),$	(179)

where

	$\displaystyle\overline{P}(s(t),x)$	$\displaystyle=\frac{r}{R_{\rm p,+}}P(r_{\rm p}(t),r),$		(180)
	$\displaystyle\overline{Q}(s(t))$	$\displaystyle=\frac{1}{R_{\rm p,+}}-\frac{1}{D_{\mathrm{s},+}}Q(r_{\rm p}(t)).$		(181)

With respect to the variable (176), the properties (154) and (155) presented in Remark 3 are equivalent to

	$\displaystyle\dot{s}(t)$	$\displaystyle>0,$		(182)
	$\displaystyle 0<s(t)$	$\displaystyle\leq R_{\rm p,+}.$		(183)

Derivation of observer gains: Consider the following invertible transformation from the estimation error $\widetilde{u}(x,t)$ to the transformed state $\widetilde{w}(x,t)$ :

	$\displaystyle\widetilde{w}(x,t)$	$\displaystyle=\widetilde{u}(x,t)+\int_{0}^{x}q(\overline{x},\overline{y})\widetilde{u}(y,t){\rm d}y,$		(184)
	$\displaystyle\widetilde{u}(x,t)$	$\displaystyle=\widetilde{w}(x,t)+\int_{0}^{x}p(\overline{x},\overline{y})\widetilde{w}(y,t){\rm d}y,$		(185)

where $\overline{x}=s(t)-x$ , $\overline{y}=s(t)-y$ . We can show that if the gain kernel functions and the observer gains satisfy the following conditions:

$\displaystyle\frac{\partial^{2}p}{\partial\bar{x}^{2}}(\bar{x},\bar{y})-\frac{\partial^{2}p}{\partial\bar{y}^{2}}(\bar{x},\bar{y})=$	$\displaystyle-\bar{\lambda}p(\bar{x},\bar{y}),$	(186)
$\displaystyle p(\bar{x},\bar{x})=$	$\displaystyle\frac{\bar{\lambda}}{2}\bar{x},$	(187)
$\displaystyle p(0,\bar{y})=$	$\displaystyle 0,$	(188)
$\displaystyle\frac{\partial^{2}q}{\partial\bar{x}^{2}}(\bar{x},\bar{y})-\frac{\partial^{2}q}{\partial\bar{y}^{2}}(\bar{x},\bar{y})=$	$\displaystyle\bar{\lambda}q(\bar{x},\bar{y}),$	(189)
$\displaystyle q(\bar{x},\bar{x})=$	$\displaystyle-\frac{\bar{\lambda}}{2}\bar{x},$	(190)
$\displaystyle q(0,\bar{y})=$	$\displaystyle 0,$	(191)
$\displaystyle\overline{P}(s(t),x)=$	$\displaystyle D_{\rm s,+}p_{\bar{y}}(\bar{x},s(t)),$	(192)
$\displaystyle\overline{Q}(s(t))=$	$\displaystyle-p(s(t),s(t)),$	(193)

then, the following target $\widetilde{w}$ -system is obtained:

$\displaystyle\frac{\partial\widetilde{w}}{\partial t}(x,t)=$	$\displaystyle\leavevmode\nobreak\ D_{\rm s,+}\frac{\partial^{2}\widetilde{w}}{\partial x^{2}}(x,t)-\lambda\widetilde{w}(x,t)+\dot{s}(t)\int_{0}^{x}q^{\prime}(\overline{x},\overline{y})$
	$\displaystyle\times\left(\widetilde{w}(y,t)+\int_{0}^{y}p(\overline{y},\overline{z})\widetilde{w}(z,t){\rm d}z\right)dy,$	(194)
$\displaystyle\widetilde{w}(s(t),t)=$	$\displaystyle\leavevmode\nobreak\ 0,$	(195)
$\displaystyle\frac{\partial\widetilde{w}}{\partial x}(0,t)=$	$\displaystyle\leavevmode\nobreak\ 0,$	(196)

where $q^{\prime}(\overline{x},\overline{y})=\frac{\partial q}{\partial\overline{x}}(\overline{x},\overline{y})+\frac{\partial q}{\partial\overline{y}}(\overline{x},\overline{y})$ . The equations (186)–(191) lead to the following explicit solutions:

	$\displaystyle p(\overline{x},\overline{y})=$	$\displaystyle\overline{\lambda}\overline{x}\frac{I_{1}\left(\sqrt{\overline{\lambda}\left[\overline{y}^{2}-\overline{x}^{2}\right]}\right)}{\sqrt{\overline{\lambda}\left[\overline{y}^{2}-\overline{x}^{2}\right]}},$		(197)
	$\displaystyle q(\overline{x},\overline{y})=$	$\displaystyle-\overline{\lambda}\overline{x}\frac{J_{1}\left(\sqrt{\overline{\lambda}\left[\overline{y}^{2}-\overline{x}^{2}\right]}\right)}{\sqrt{\overline{\lambda}\left[\overline{y}^{2}-\overline{x}^{2}\right]}},$		(198)

with a modified Bessel function $I_{1}(\cdot)$ and a Bessel function $J_{1}(\cdot)$ of the first kind, respectively. Substituting the solution (197) to the conditions (192), (193) (note that $\frac{dI_{1}(z)}{{\rm d}z}=\frac{I_{2}(z)}{z}$ for all $z$ ), and taking back to the original coordinate and variables, the observer gains are derived as (164) and (165).

Stability proof: We consider the time evolution of the following Lyapunov function:

\displaystyle W(t)=\frac{1}{2}\int_{0}^{s(t)}\widetilde{w}(x,t)^{2}{\rm d}x.

(199)

Taking the time derivative of (199) along with (4.8)-(196) yields

$\displaystyle\dot{W}(t)=$	$\displaystyle-D_{\mathrm{s},+}\int_{0}^{s(t)}\left(\frac{\partial\widetilde{w}}{\partial x}(x,t)\right)^{2}{\rm d}x-\lambda\int_{0}^{s(t)}\widetilde{w}(x,t)^{2}{\rm d}x$
	$\displaystyle+\dot{s}(t)\int_{0}^{s(t)}\widetilde{w}(x,t)\left[\int_{0}^{x}q^{\prime}(\overline{x},\overline{y})\right.$
	$\displaystyle\left.\left(\widetilde{w}(y,t)+\int_{0}^{y}P(\overline{y},\overline{z})\widetilde{w}(z,t){\rm d}z\right)dy\right]{\rm d}x.$	(200)

Applying Young’s, Cauchy Schwartz, and Poincare’s inequalities with the help of the properties (182) and (183), one can show that there exists a constant $a>0$ such that the following inequality holds:

\displaystyle\dot{W}(t)\leq

\displaystyle-bW(t)+a\dot{s}(t)W(t),

(201)

where $b=\frac{D_{\mathrm{s},+}}{4R_{\rm p,+}^{2}}+\lambda$ . With the help of (182) and (183), it yields the exponential decay of $W(t)$ as

\displaystyle W(t)\leq e^{aR_{\rm p,+}}W(0)e^{-bt}.

(202)

Hence, the origin of $\widetilde{w}$ -system is shown to be exponentially stable, from which we conclude Theorem 3.

Observer for Negative Electrode

The observer design for lithium ion concentration in the negative electrode is constructed by the copy of the dynamics (143)-(145) together with the output injection of the positive electrode

\displaystyle\frac{\partial\widehat{c}_{\mathrm{s},-}}{\partial t}(r,t)=

\displaystyle\frac{D_{\rm s,-}}{r^{2}}\frac{\partial}{\partial r}\left[r^{2}\frac{\partial\widehat{c}_{\mathrm{s},-}}{\partial r}(r,t)\right]+P_{-}(r_{\rm{p}}(t))\widetilde{c}_{\mathrm{s,+}}(R_{\rm p,+},t),

(203)

for $r\in(0,R_{\rm p,-})$ , $t>0$ with boundary conditions

	$\displaystyle\frac{\partial\widehat{c}_{\mathrm{s},-}}{\partial r}(0,t)=$	$\displaystyle 0,$		(204)
	$\displaystyle D_{\mathrm{s},-}\frac{\partial\widehat{c}_{\mathrm{s},-}}{\partial r}(R_{\mathrm{p},-},t)=$	$\displaystyle-j_{\mathrm{n},-}(t)+Q_{-}\left(r_{\rm p}(t)\right)\widetilde{c}_{\mathrm{s,+}}(R_{\rm p,+},t).$		(205)

Observer gains in the negative electrode are computed to conserve the total amount of lithium ions in the state observer defined as

\displaystyle\widehat{n}_{{\rm Li}}(t)=

\displaystyle\frac{3\epsilon_{\rm s,+}L_{+}}{R_{\rm{p},+}^{3}}\int_{0}^{R_{\rm{p},+}}\widehat{c}_{\rm{s},+}(r,t)r^{2}{\rm d}r+\frac{3\epsilon_{\rm s,-}L_{-}}{R_{\rm{p},-}^{3}}\int_{0}^{R_{\rm{p},-}}\widehat{c}_{\rm{s},-}(r,t)r^{2}{\rm d}r.

(206)

Taking the time derivative of (206) along with the dynamics (4.8)–(169) and (203)–(205) leads to

\displaystyle\frac{d\widehat{n}_{{\rm Li},+}}{dt}=

\displaystyle-a_{\rm s,+}L_{+}j_{\mathrm{n},+}(t)-a_{\rm s,-}L_{-}j_{\mathrm{n},-}(t)+F\widetilde{c}_{\mathrm{s,+}}(R_{\rm p,+},t),

(207)

where $F$ is defined by

	$\displaystyle F=$	$\displaystyle\leavevmode\nobreak\ a_{\rm s,+}L_{+}\left(\frac{\kappa}{R_{\rm{p},+}^{2}}\widehat{r}_{\rm p}(t)^{2}+Q(\widehat{r}_{\rm p}(t))\right)+a_{\rm s,-}L_{-}Q_{-}(\widehat{r}_{\rm p}(t))$
		$\displaystyle+\frac{3\epsilon_{\rm s,+}L_{+}}{R_{\rm{p},+}^{3}}\int_{\widehat{r}_{\rm p}(t)}^{R_{\rm{p},+}}r^{2}P(\widehat{r}_{\rm p}(t),r)dr+\epsilon_{\rm{s},-}L_{-}P_{-}(\widehat{r}_{\rm p}(t)).$		(208)

By the balance of the ionic molar fluxes given in (142), the first line in the right hand side of (207) is canceled. Therefore, by designing the observer gains as

	$\displaystyle Q_{-}(r_{\rm{p}}(t))$	$\displaystyle=-\frac{a_{s,+}L_{+}}{a_{s,-}L_{-}}\left(Q(r_{\rm{p}}(t))+\frac{\kappa}{R_{\rm{p},+}^{2}}\widehat{r}_{\rm p}(t)^{2}\right),$		(209)
	$\displaystyle P_{-}(r_{\rm p}(t))$	$\displaystyle=-\frac{\epsilon_{\rm{s},+}L_{+}}{\epsilon_{\rm{s},-}L_{-}}\frac{3}{R^{3}_{\rm{p},+}}\left[\int_{\widehat{r}_{\rm{p}}(t)}^{R_{\rm{p,+}}}P(r_{\rm{p}}(t))r^{2}dr\right],$		(210)

one can show that $\frac{d{\widehat{n}_{{\rm Li},+}}}{dt}=0$ from (207). Hence, the observer error in the negative electrode approaches to zero uniformly in space with the help of Theorem 3.

4.9 Numerical Simulation

Test with Constant Discharge Input

To test the observer we run a numerical example with a constant discharge current of $5$ [C-rate]. We are assuming $c_{\mathrm{ss},+}$ is available directly from measurements to be used as output error injection in the observer. In practice, this quantity could be estimated from measurements. Figure 15 shows the estimated concentration of lithium ions in $\beta$ -phase in the positive particle; one can compare this to the true concentration in Figure 14. Figure 16 shows the averaged concentration in the positive particle, both true value (black) and estimated value (blue). Convergence of the estimate to the true value is achieved within 0.8 [min], a relatively short time. Furthermore, Fig. 17 shows the time evolution of the moving interface of the both true value (black) and estimated value (blue), which also illustrates the convergence of the estimate to the true value. Note that SoC is directly proportional to the averaged concentrations; then the importance to evaluate the estimation of this quantity.

Comparison with the Extended Kalman Filter

Since the spatial discretization of the diffusion equations is performed for computation of the electrochemical model, we can apply the Extended Kalman Filter (EKF) to the reduced-order model as another approach for SoC estimation. Here, we compare the performance of SoC estimation between the proposed backstepping (BKS) observer and the EKF with incorporating a measurement noise.

Fig. 18 shows a simulation result of SoC estimation via the BKS observer (blue) and the EKF (red). The initial SoC in the model is around 66 %, while the initial SoC in the estimator is around 46 % for both BKS and EKF. Fig. 18 (a) depicts the result under the noise-free measurement, which illustrates that the BKS estimate converges and almost stays at the true value within 5 minutes, while the estimate by the EKF converges quickly first but does not stay at the true value even after 10 minues. Next, we incorporate the Gaussian noise in the measured value of the surface concentration, and compare the performance in Fig. 18 (b). The result illustrates that the BKS estimate still converges to the true value but it accompanies a noisy signal, while the EKF estimate has less noisy signal. Hence, in this one sample simulation, we observe that the BKS estimator is superior in convergence speed, while the EKF estimator is superior in noise attenuation. However, by lowing the gain parameters $(\lambda,\kappa)$ in the BKS estimator, we observe that the amplitude of the noisy estimate can be reduced in exchange for the convergence speed. Moreover, the EKF estimate also can be improved to accelerate the convergence speed by appropriately tuning the free parameters. Hence, there is an essential tradeoff between the convergence speed and the noise attenuation for both methods, and it is not appropriate to address which method is superior in general. Nevertheless, one of the advantages of the proposed BKS estimator is to have only two free parameters to tune for any given discretization number, while the EKF algorithm requires a tuning of covariance matrices in which the number of free parameters is increased as the discretization number is increased.

Table 3: Parameters of LFP used in the simulation.

Parameters
	Negative	Separator	Positive
$L[\rm m]^{a}$	$50\times 10^{-6}$	$25\times 10^{-6}$	$74\times 10^{-6}$
$c^{\rm max}_{\rm s}[\rm{mol}/m^{3}]^{a}$	$27760$		$20950$
$c_{\mathrm{s},\alpha}[\rm{mol}/m^{3}]^{b}$			0.0480 $\times c^{\rm max}_{\rm s,+}$
$c_{\mathrm{s},\beta}[\rm{mol}/m^{3}]^{b}$			0.8920 $\times c^{\rm max}_{\rm s,+}$
$R_{\rm p}\leavevmode\nobreak\ [\rm m]^{a}$	$11\times 10^{-6}$		$52\times 10^{-9}$
$D_{\rm s}\leavevmode\nobreak\ [\rm m^{2}/s]^{a}$	$9\times 10^{-14}$		$8\times 10^{-18}$
$\epsilon_{\rm s}\leavevmode\nobreak\ [-]^{\rm{a}}$	$0.33$		$0.27$
$R_{\rm f}\leavevmode\nobreak\ [\rm\Omega m^{2}]^{b}$	$1\times 10^{-5}$		$0$
$R_{\rm c}\leavevmode\nobreak\ [\rm\Omega m^{2}]^{b}$	$0$		$6.5\times 10^{-3}$
$k\leavevmode\nobreak\ [\rm m^{2.5}/mol^{0.5}s]^{a}$	$3\times 10^{-5}$		$3\times 10^{-17}$

Other Parameters and Physical Constants
$A\leavevmode\nobreak\ [\rm m]^{b}$	1
$F\leavevmode\nobreak\ [\rm As/mol]$	$96487$
$R\leavevmode\nobreak\ [\rm J/Kmol]$	8.314472
$T\leavevmode\nobreak\ [\rm K]^{b}$	298
$c_{\textrm{e}}\leavevmode\nobreak\ [\rm mol/m^{3}]^{a}$	$1\times 10^{3}$
$\alpha_{\rm a},\alpha_{\rm c}\leavevmode\nobreak\ [\rm-]^{a}$	0.5

^a borrowed from Srinivasan and Newman (2004a)
^b assumed

5 Conclusion and Open Problems

In this tutorial article, we have presented the state estimation of the Stefan PDE system modeling the dynamics of phase change phenomena of a material, and its applications to polar ice and lithium-ion batteries. The estimator is designed by an backstepping observer, which is given by a copy of the plant plus output error injection, where the observer gain is derived explicitly by solving a gain kernel equation of the state transformation. The convergence of the designed observer to the plant state from the Stefan system is ensured by Lyapunov analysis. The simulation results conducted for the thermodynamic model of Arctic sea ice illustrate the robust performance of the designed observer with respect to the neglected salinity effect and parameter uncertainty, where the convergence of the estimated temperature distribution to the true temperature is achieved within three days. The simulation performed for the electrochemical model of the lithium-ion batteries with phase transition materials has shown that the reduction of the error of more than 15 % in the SoC estimate is achieved within five minutes even in the presence of sensor noise.

There are various exciting open problems for the state estimation of the Stefan system, both from control-theoretic and application-driven perspectives. We summarize them in the following list for the control-theoretic problems:

1.

sensor-delay compensation in the observer for the Stefan system (see Ahmed-Ali et al. (2019)),
2.

adaptive observer design for the Stefan system (see Ahmed-Ali et al. (2015, 2017a)),
3.

observer for two-phase Stefan system under a single-boundary measurement (see Liu et al. (2016)),
4.

sampled-data observer for the Stefan system (see Ahmed-Ali et al. (2017b); Karafyllis et al. (2019)),
5.

prescribed-time observer design for the Stefan system (see Steeves et al. (2020, 2019)), and for the application-driven problems:
6.

cancer treatment via cryosurgery (see Rabin and Shitzer (1998); Rabin and Stahovich (2003); Kumar et al. (2017)),
7.

spreading of invasive species in ecology (see Du and Lin (2010)),
8.

information diffusion on online social networks (see Wang et al. (2020)),
9.

domain walls in ferroelectric thin films (see McGilly et al. (2015)),
10.

Black-Scholes model of American option pricing (see Chen et al. (2008)).

While the list of control-theoretic problems have all been considered for PDE and PDE-ODE systems on fixed domains, virtually all of them are unexplored for PDEs with moving boundaries. Moreover, the list of application-driven topics owns a significant impact in the real-world problems, by the utilization of the PDE-based estimation method. This tutorial review provides supporting technical materials to tackle those challenging topics.

References

Agarwala et al. (1995) Agarwala, M., Bourell, D., Beaman, J., Marcus, H., Barlow, J., 1995. Direct selective laser sintering of metals. Rapid Prototyping Journal 1, 26–36.
Ahmed-Ali et al. (2019) Ahmed-Ali, T., Fridman, E., Giri, F., Kahelras, M., Lamnabhi-Lagarrigue, F., Burlion, L., 2019. Observer design for a class of parabolic systems with large delays and sampled measurements. IEEE Transactions on Automatic Control 65, 2200–2206.
Ahmed-Ali et al. (2017a) Ahmed-Ali, T., Giri, F., Krstic, M., Burlion, L., Lamnabhi-Lagarrigue, F., 2017a. Adaptive observer design with heat pde sensor. Automatica 82, 93–100.
Ahmed-Ali et al. (2015) Ahmed-Ali, T., Giri, F., Krstic, M., Lamnabhi-Lagarrigue, F., Burlion, L., 2015. Adaptive observer for a class of parabolic pdes. IEEE Transactions on Automatic Control 61, 3083–3090.
Ahmed-Ali et al. (2017b) Ahmed-Ali, T., Karafyllis, I., Giri, F., Krstic, M., Lamnabhi-Lagarrigue, F., 2017b. Exponential stability analysis of sampled-data ODE–PDE systems and application to observer design. IEEE Transactions on Automatic Control 62, 3091–3098.
Alessandri et al. (2018) Alessandri, A., Bagnerini, P., Gaggero, M., 2018. Optimal control of propagating fronts by using level set methods and neural approximations. IEEE transactions on neural networks and learning systems 30, 902–912.
Alexiades and Solomon (2018) Alexiades, V., Solomon, A.D., 2018. Mathematical modeling of melting and freezing processes. Routledge.
Barbu et al. (1996) Barbu, V., Kunisch, K., Ring, W., 1996. Control and estimation of the boundary heat transfer function in stefan problems. ESAIM: Mathematical Modelling and Numerical Analysis 30, 671–710.
Bernauer and Herzog (2011) Bernauer, M.K., Herzog, R., 2011. Optimal control of the classical two-phase stefan problem in level set formulation. SIAM Journal on Scientific Computing 33, 342–363.
Bitz et al. (2001) Bitz, C., Holland, M., Weaver, A., Eby, M., 2001. Simulating the ice-thickness distribution in a coupled climate model. Journal of Geophysical Research: Oceans 106, 2441–2463.
Bitz and Lipscomb (1999) Bitz, C.M., Lipscomb, W.H., 1999. An energy-conserving thermodynamic model of sea ice. Journal of Geophysical Research: Oceans 104, 15669–15677.
Bonnerot and Jamet (1977) Bonnerot, R., Jamet, P., 1977. Numerical computation of the free boundary for the two-dimensional stefan problem by space-time finite elements. Journal of Computational Physics 25, 163–181.
Boyd et al. (1994) Boyd, S., El Ghaoui, L., Feron, E., Balakrishnan, V., 1994. Linear matrix inequalities in system and control theory. SIAM.
Braatz (2002) Braatz, R.D., 2002. Advanced control of crystallization processes. annual reviews in control 26, 87–99.
Cannon (1984) Cannon, J.R., 1984. The one-dimensional heat equation. Cambridge University Press.
Cannon and Primicerio (1971a) Cannon, J.R., Primicerio, M., 1971a. A two phase stefan problem with flux boundary conditions. Annali di Matematica Pura ed Applicata 88, 193–205.
Cannon and Primicerio (1971b) Cannon, J.R., Primicerio, M., 1971b. A two phase stefan problem with temperature boundary conditions. Annali di Matematica Pura ed Applicata 88, 177–191.
Chaturvedi et al. (2010) Chaturvedi, N.A., Klein, R., Christensen, J., Ahmed, J., Kojic, A., 2010. Algorithms for advanced battery-management systems. IEEE Control systems magazine 30, 49–68.
Chen et al. (2008) Chen, X., Chadam, J., Jiang, L., Zheng, W., 2008. Convexity of the exercise boundary of the american put option on a zero dividend asset. Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics 18, 185–197.
Chen et al. (2018) Chen, Z., Bentsman, J., Thomas, B.G., 2018. Bang-bang free boundary control of a stefan problem for metallurgical length maintenance, in: 2018 Annual American Control Conference (ACC), IEEE. pp. 116–121.
Chen et al. (2019) Chen, Z., Bentsman, J., Thomas, B.G., 2019. Optimal control of free boundary of a stefan problem for metallurgical length maintenance in continuous steel casting, in: 2019 American Control Conference (ACC), IEEE. pp. 3206–3211.
Chen et al. (2020) Chen, Z., Bentsman, J., Thomas, B.G., 2020. Enthalpy-based output feedback control of the stefan problem with hysteresis, in: 2020 American Control Conference (ACC), IEEE. pp. 2661–2666.
Colli et al. (1999) Colli, P., Grasselli, M., Sprekels, J., 1999. Automatic control via thermostats of a hyperbolic stefan problem with memory. Applied Mathematics and Optimization 39, 229–255.
Crepeau (2007) Crepeau, J., 2007. Josef stefan: His life and legacy in the thermal sciences. Experimental thermal and fluid science 31, 795–803.
Demir et al. (2021) Demir, C., Koga, S., Krstic, M., 2021. Neuron growth control by pde backstepping: Axon length regulation by tubulin flux actuation in soma, in: 2021 IEEE Conference on Decision and Control (CDC), IEEE.
Dey et al. (2015) Dey, S., Ayalew, B., Pisu, P., 2015. Nonlinear robust observers for state-of-charge estimation of lithium-ion cells based on a reduced electrochemical model. IEEE Transactions on Control Systems Technology 23, 1935–1942.
Di Domenico et al. (2010) Di Domenico, D., Stefanopoulou, A., Fiengo, G., 2010. Lithium-ion battery state of charge and critical surface charge estimation using an electrochemical model-based extended kalman filter. Journal of dynamic systems, measurement, and control 132.
Diehl et al. (2014) Diehl, S., Henningsson, E., Heyden, A., Perna, S., 2014. A one-dimensional moving-boundary model for tubulin-driven axonal growth. Journal of theoretical biology 358, 194–207.
Du and Lin (2010) Du, Y., Lin, Z., 2010. Spreading-vanishing dichotomy in the diffusive logistic model with a free boundary. SIAM Journal on Mathematical Analysis 42, 377–405.
Dunbar et al. (2003a) Dunbar, W.B., Petit, N., Rouchon, P., Martin, P., 2003a. Boundary control of a nonlinear stefan problem, in: 42nd IEEE International Conference on Decision and Control (IEEE Cat. No. 03CH37475), IEEE. pp. 1309–1314.
Dunbar et al. (2003b) Dunbar, W.B., Petit, N., Rouchon, P., Martin, P., 2003b. Motion planning for a nonlinear stefan problem. ESAIM: Control, Optimisation and Calculus of Variations 9, 275–296.
Ecklebe et al. (2020) Ecklebe, S., Woittennek, F., Winkler, J., 2020. Control of the vertical gradient freeze crystal growth process via backstepping. arXiv preprint arXiv:2002.11447 .
Ecklebe et al. (2019) Ecklebe, S., Woittennek, F., Winkler, J., Frank-Rotsch, C., Dropka, N., 2019. Towards model based control of the vertical gradient freeze crystal growth process. arXiv preprint arXiv:1908.02519 .
Fasano and Primicerio (1977) Fasano, A., Primicerio, M., 1977. General free-boundary problems for the heat equation, i. Journal of Mathematical analysis and applications 57, 694–723.
Fenty and Heimbach (2013) Fenty, I., Heimbach, P., 2013. Coupled sea ice–ocean-state estimation in the labrador sea and baffin bay. Journal of Physical Oceanography 43, 884–904.
Fenty et al. (2017) Fenty, I., Menemenlis, D., Zhang, H., 2017. Global coupled sea ice-ocean state estimation. Climate Dynamics 49, 931–956.
Friedman (1959) Friedman, A., 1959. Free boundary problems for parabolic equations i. melting of solids. Journal of Mathematics and Mechanics , 499–517.
Fujiwara et al. (2005) Fujiwara, M., Nagy, Z.K., Chew, J.W., Braatz, R.D., 2005. First-principles and direct design approaches for the control of pharmaceutical crystallization. Journal of Process Control 15, 493–504.
Gupta (2017) Gupta, S.C., 2017. The classical Stefan problem: basic concepts, modelling and analysis with quasi-analytical solutions and methods. volume 45. Elsevier.
Hall et al. (2004) Hall, D.K., Key, J.R., Casey, K.A., Riggs, G.A., Cavalieri, D.J., 2004. Sea ice surface temperature product from modis. IEEE transactions on geoscience and remote sensing 42, 1076–1087.
Hinze and Ziegenbalg (2007a) Hinze, M., Ziegenbalg, S., 2007a. Optimal control of the free boundary in a two-phase stefan problem. Journal of Computational Physics 223, 657–684.
Hinze and Ziegenbalg (2007b) Hinze, M., Ziegenbalg, S., 2007b. Optimal control of the free boundary in a two-phase stefan problem with flow driven by convection. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics 87, 430–448.
Hoffmann and Sprekels (1982) Hoffmann, K.H., Sprekels, J., 1982. Real-time control of the free boundary in a two-phase stefan problem. Numerical Functional Analysis and Optimization 5, 47–76.
Javierre et al. (2006) Javierre, E., Vuik, C., Vermolen, F., Van der Zwaag, S., 2006. A comparison of numerical models for one-dimensional stefan problems. Journal of Computational and Applied Mathematics 192, 445–459.
Kalman (1960) Kalman, R.E., 1960. A new approach to linear filtering and prediction problems. ASME Journal of Fluids Engineering 82, 35–45.
Kang and Zabaras (1995) Kang, S., Zabaras, N., 1995. Control of the freezing interface motion in two-dimensional solidification processes using the adjoint method. International Journal for Numerical Methods in Engineering 38, 63–80.
Karafyllis et al. (2019) Karafyllis, I., Ahmed-Ali, T., Giri, F., 2019. Sampled-data observers for 1-d parabolic PDEs with non-local outputs. Systems & Control Letters 133, 104553.
Khandelwal et al. (2015) Khandelwal, A., Hariharan, K.S., Gambhire, P., Kolake, S.M., Yeo, T., Doo, S., 2015. Thermally coupled moving boundary model for charge–discharge of lifepo4/c cells. Journal of Power Sources 279, 180–196.
Khandelwal et al. (2014) Khandelwal, A., Hariharan, K.S., Kumar, V.S., Gambhire, P., Kolake, S.M., Oh, D., Doo, S., 2014. Generalized moving boundary model for charge–discharge of lifepo4/c cells. Journal of Power Sources 248, 101–114.
Klein et al. (2012) Klein, R., Chaturvedi, N.A., Christensen, J., Ahmed, J., Findeisen, R., Kojic, A., 2012. Electrochemical model based observer design for a lithium-ion battery. IEEE Transactions on Control Systems Technology 21, 289–301.
Koga et al. (2020) Koga, S., Bresch-Pietri, D., Krstic, M., 2020. Delay compensated control of the stefan problem and robustness to delay mismatch. International Journal of Robust and Nonlinear Control 30, 2304–2334.
Koga et al. (2021a) Koga, S., Camacho-Solorio, L., Krstic, M., 2021a. State estimation for lithium-ion batteries with phase transition materials via boundary observers. Journal of Dynamic Systems, Measurement, and Control 143.
Koga et al. (2018) Koga, S., Diagne, M., Krstic, M., 2018. Control and state estimation of the one-phase stefan problem via backstepping design. IEEE Transactions on Automatic Control 64, 510–525.
Koga et al. (2021b) Koga, S., Karafyllis, I., Krstic, M., 2021b. Towards implementation of PDE control for Stefan system: Input-to-state stability and sampled-data design. Automatica 127, 109538.
Koga and Krstic (2020a) Koga, S., Krstic, M., 2020a. Arctic sea ice state estimation from thermodynamic PDE model. Automatica 112, 108713.
Koga and Krstic (2020b) Koga, S., Krstic, M., 2020b. Single-boundary control of the two-phase stefan system. Systems & Control Letters 135, 104573.
Koga and Krstic (2021) Koga, S., Krstic, M., 2021. Control of the Stefan system and applications: A tutorial. Annual Review of Control, Robotics, and Autonomous Systems 5.
Koga et al. (2020) Koga, S., Krstic, M., Beaman, J., 2020. Laser sintering control for metal additive manufacturing by PDE backstepping. IEEE Transactions on Control Systems Technology 28, 1928–1939.
Koga et al. (2020a) Koga, S., Makihata, M., Chen, R., Krstic, M., Pisano, A.P., 2020a. Energy storage in paraffin: A PDE backstepping experiment. IEEE Transactions on Control Systems Technology 29, 1490–1502.
Koga et al. (2020b) Koga, S., Straub, D., Diagne, M., Krstic, M., 2020b. Stabilization of filament production rate for screw extrusion-based polymer three-dimensional-printing. Journal of Dynamic Systems, Measurement, and Control 142.
Kolodner (1956) Kolodner, I., 1956. Free boundary problem for the heat equation with applications to problems of change of phase. Comm. Pure Appl. Math 9, 1–31.
Krstic and Koga (2020) Krstic, M., Koga, S., 2020. Materials Phase Change PDE Control and Estimation: From Additive Manufacturing to Polar Ice. Springer Nature.
Kumar et al. (2017) Kumar, A., Kumar, S., Katiyar, V., Telles, S., 2017. Phase change heat transfer during cryosurgery of lung cancer using hyperbolic heat conduction model. Computers in biology and medicine 84, 20–29.
Kutluay et al. (1997) Kutluay, S., Bahadir, A., Özdeş, A., 1997. The numerical solution of one-phase classical stefan problem. Journal of computational and applied mathematics 81, 135–144.
Kwok (2015) Kwok, R., 2015. Sea ice convergence along the arctic coasts of greenland and the canadian arctic archipelago: Variability and extremes (1992–2014). Geophysical Research Letters 42, 7598–7605.
Kwok and Rothrock (2009) Kwok, R., Rothrock, D., 2009. Decline in arctic sea ice thickness from submarine and icesat records: 1958–2008. Geophysical Research Letters 36.
Li et al. (2014) Li, J., Barillas, J.K., Guenther, C., Danzer, M.A., 2014. Sequential monte carlo filter for state estimation of lifepo4 batteries based on an online updated model. Journal of Power Sources 247, 156–162.
Liu et al. (2016) Liu, B.N., Boutat, D., Liu, D.Y., 2016. Backstepping observer-based output feedback control for a class of coupled parabolic pdes with different diffusions. Systems & Control Letters 97, 61–69.
Luenberger (1971) Luenberger, D., 1971. An introduction to observers. IEEE Transactions on automatic control 16, 596–602.
Maidi and Corriou (2014) Maidi, A., Corriou, J.P., 2014. Boundary geometric control of a linear stefan problem. Journal of Process Control 24, 939–946.
Marshall et al. (1997) Marshall, J., Adcroft, A., Hill, C., Perelman, L., Heisey, C., 1997. A finite-volume, incompressible navier stokes model for studies of the ocean on parallel computers. Journal of Geophysical Research: Oceans 102, 5753–5766.
Maykut and Untersteiner (1971) Maykut, G.A., Untersteiner, N., 1971. Some results from a time-dependent thermodynamic model of sea ice. Journal of Geophysical Research 76, 1550–1575.
McGilly et al. (2015) McGilly, L., Yudin, P., Feigl, L., Tagantsev, A., Setter, N., 2015. Controlling domain wall motion in ferroelectric thin films. Nature nanotechnology 10, 145–150.
Meng and Thomas (2003) Meng, Y., Thomas, B.G., 2003. Heat-transfer and solidification model of continuous slab casting: Con1d. Metallurgical and materials transactions B 34, 685–705.
Moura et al. (2016) Moura, S.J., Argomedo, F.B., Klein, R., Mirtabatabaei, A., Krstic, M., 2016. Battery state estimation for a single particle model with electrolyte dynamics. IEEE Transactions on Control Systems Technology 25, 453–468.
Moura et al. (2014) Moura, S.J., Chaturvedi, N.A., Krstic, M., 2014. Adaptive partial differential equation observer for battery state-of-charge/state-of-health estimation via an electrochemical model. Journal of Dynamic Systems, Measurement, and Control 136.
Nagy and Braatz (2012) Nagy, Z.K., Braatz, R.D., 2012. Advances and new directions in crystallization control. Annual review of chemical and biomolecular engineering 3, 55–75.
Padhi et al. (1997) Padhi, A.K., Nanjundaswamy, K.S., Goodenough, J.B., 1997. Phospho-olivines as positive-electrode materials for rechargeable lithium batteries. Journal of the electrochemical society 144, 1188.
Pawlow (1987) Pawlow, I., 1987. Optimal control of two-phase stefan problems?numerical solutions, in: Optimal Control of Partial Differential Equations II: Theory and Applications. Springer, pp. 179–206.
Pawłow (1990) Pawłow, I., 1990. Optimal control of dynamical processes in two-phase systems of solid-liquid type. Banach Center Publications 24, 293–319.
Perez et al. (2015) Perez, H., Shahmohammadhamedani, N., Moura, S., 2015. Enhanced performance of li-ion batteries via modified reference governors and electrochemical models. IEEE/ASME Transactions on Mechatronics 20, 1511–1520.
Perez and Moura (2015) Perez, H.E., Moura, S.J., 2015. Sensitivity-based interval pde observer for battery soc estimation, in: 2015 American Control Conference (ACC), IEEE. pp. 323–328.
Pernsteiner et al. (2021) Pernsteiner, D., Schirrer, A., Kasper, L., Hofmann, R., Jakubek, S., 2021. State estimation concept for a nonlinear melting/solidification problem of a latent heat thermal energy storage. Computers & Chemical Engineering 153, 107444.
Petit (2010) Petit, N., 2010. Control problems for one-dimensional fluids and reactive fluids with moving interfaces, in: Advances in the theory of control, signals and systems with physical modeling. Springer, pp. 323–337.
Petrus et al. (2010) Petrus, B., Bentsman, J., Thomas, B.G., 2010. Feedback control of the two-phase stefan problem, with an application to the continuous casting of steel, in: 49th IEEE Conference on Decision and Control (CDC), IEEE. pp. 1731–1736.
Petrus et al. (2012) Petrus, B., Bentsman, J., Thomas, B.G., 2012. Enthalpy-based feedback control algorithms for the stefan problem, in: 2012 IEEE 51st IEEE Conference on Decision and Control (CDC), IEEE. pp. 7037–7042.
Petrus et al. (2014) Petrus, B., Bentsman, J., Thomas, B.G., 2014. Application of enthalpy-based feedback control methodology to the two-sided stefan problem, in: 2014 American Control Conference, IEEE. pp. 1015–1020.
Petrus et al. (2017) Petrus, B., Chen, Z., Bentsman, J., Thomas, B.G., 2017. Online recalibration of the state estimators for a system with moving boundaries using sparse discrete-in-time temperature measurements. IEEE Transactions on Automatic Control 63, 1090–1096.
Rabin and Shitzer (1998) Rabin, Y., Shitzer, A., 1998. Numerical solution of the multidimensional freezing problem during cryosurgery. Journal of biomechanical engineering 120, 32–37.
Rabin and Stahovich (2003) Rabin, Y., Stahovich, T.F., 2003. Cryoheater as a means of cryosurgery control. Physics in Medicine & Biology 48, 619.
Rothrock et al. (2008) Rothrock, D., Percival, D., Wensnahan, M., 2008. The decline in arctic sea-ice thickness: Separating the spatial, annual, and interannual variability in a quarter century of submarine data. Journal of Geophysical Research: Oceans 113.
Rubinšteĭn (2000) Rubinšteĭn, L., 2000. The stefan problem. volume 8. American Mathematical Soc.
Schwunk et al. (2013) Schwunk, S., Armbruster, N., Straub, S., Kehl, J., Vetter, M., 2013. Particle filter for state of charge and state of health estimation for lithium–iron phosphate batteries. Journal of Power Sources 239, 705–710.
Semtner Jr (1976) Semtner Jr, A.J., 1976. A model for the thermodynamic growth of sea ice in numerical investigations of climate. Journal of Physical Oceanography 6, 379–389.
Sharma et al. (2009) Sharma, A., Tyagi, V.V., Chen, C., Buddhi, D., 2009. Review on thermal energy storage with phase change materials and applications. Renewable and Sustainable energy reviews 13, 318–345.
Spangler et al. (2016) Spangler, B., Fontaine, S.D., Shi, Y., Sambucetti, L., Mattis, A.N., Hann, B., Wells, J.A., Renslo, A.R., 2016. A novel tumor-activated prodrug strategy targeting ferrous iron is effective in multiple preclinical cancer models. Journal of medicinal chemistry 59, 11161–11170.
Srinivasan and Newman (2004a) Srinivasan, V., Newman, J., 2004a. Design and optimization of a natural graphite/iron phosphate lithium-ion cell. Journal of the Electrochemical Society 151, A1530.
Srinivasan and Newman (2004b) Srinivasan, V., Newman, J., 2004b. Discharge model for the lithium iron-phosphate electrode. Journal of the Electrochemical Society 151, A1517.
Steeves et al. (2019) Steeves, D., Krstic, M., Vazquez, R., 2019. Prescribed-time h 1-stabilization of reaction-diffusion equations by means of output feedback, in: 2019 18th European Control Conference (ECC), IEEE. pp. 1932–1937.
Steeves et al. (2020) Steeves, D., Krstic, M., Vazquez, R., 2020. Prescribed–time estimation and output regulation of the linearized schrödinger equation by backstepping. European Journal of Control .
STEFAN (1891) STEFAN, J., 1891. Uber die theorie der eisbilfung, insbesondere uber die eisbildung im polarmeere. Annal. Phy. Chem. 42, 269–286.
Tang et al. (2017) Tang, S.X., Camacho-Solorio, L., Wang, Y., Krstic, M., 2017. State-of-charge estimation from a thermal–electrochemical model of lithium-ion batteries. Automatica 83, 206–219.
Thomas et al. (2002) Thomas, K.E., Newman, J., Darling, R.M., 2002. Mathematical modeling of lithium batteries, in: Advances in lithium-ion batteries. Springer, pp. 345–392.
Untersteiner (1961) Untersteiner, N., 1961. On the mass and heat budget of arctic sea ice. Archiv für Meteorologie, Geophysik und Bioklimatologie, Serie A 12, 151–182.
Valkenaers et al. (2013) Valkenaers, H., Vogeler, F., Ferraris, E., Voet, A., Kruth, J., 2013. A novel approach to additive manufacturing: Screw extrusion 3d-printing, in: Proceedings of the 10th International Conference on Multi-Material Micro Manufacture, Research Publishing. pp. 235–238.
Wang et al. (2020) Wang, H., Wang, F., Xu, K., 2020. Modeling information diffusion in online social networks with partial differential equations. volume 7. Springer Nature.
Wang et al. (2014) Wang, Y., Fang, H., Sahinoglu, Z., Wada, T., Hara, S., 2014. Adaptive estimation of the state of charge for lithium-ion batteries: Nonlinear geometric observer approach. IEEE Transactions on Control Systems Technology 23, 948–962.
Winton (2000) Winton, M., 2000. A reformulated three-layer sea ice model. Journal of atmospheric and oceanic technology 17, 525–531.
Wunsch and Heimbach (2007) Wunsch, C., Heimbach, P., 2007. Practical global oceanic state estimation. Physica D: Nonlinear Phenomena 230, 197–208.
Zabaras et al. (1988) Zabaras, N., Mukherjee, S., Richmond, O., 1988. An analysis of inverse heat transfer problems with phase changes using an integral method. ASME Journal of Heat Transfer 110, 554–561.
Zabaras et al. (1992) Zabaras, N., Ruan, Y., Richmond, O., 1992. Design of two-dimensional stefan processes with desired freezing front motions. Numerical Heat Transfer, Part B Fundamentals 21, 307–325.
Zhang and White (2007) Zhang, Q., White, R.E., 2007. Moving boundary model for the discharge of a LiCoO2 electrode. Journal of The Electrochemical Society 154, A587.

State Estimation of the Stefan PDE: A Tutorial on Design and Applications to Polar Ice and Batteries

Abstract

keywords:

1 Introduction

1.1 Phase transitions

1.2 Stefan model

1.3 Control of Stefan system

1.4 Backstepping designs for Stefan system

1.5 State estimators for Stefan model

2 Stefan System

Assumption 1

Lemma 1

3 State Estimator Design by PDE Backstepping

3.1 Estimation of temperature profile by measuring boundary temperature and interface position

3.1.1 Estimator design and main theorem

Theorem 1

3.1.2 Fixed domain design

3.1.3 Analogous observer design on moving boundary domain

3.2 Estimation of Both Temperature Profile and Moving Interface by Measuring Only a Boundary Temperature

4 Sea Ice

4.1 Importance of the Arctic Sea Ice for Global Climate Modeling

4.2 Thermodynamic Model of Arctic Sea Ice

4.3 Annual Cycle Simulation of Sea Ice Thickness

Input Parameters

Simulation Test of MU71

4.4 Temperature Profile Estimation

Simplification of the Model

Assumption 2

Assumption 3

Observer Structure

Theorem 2

Remark 1

Remark 2

Gain Derivation via State Transformation

Stability Analysis

4.5 Numerical Tests of the Sea Ice Estimation

Initial conditions

Tuning method for gain parameters

Numerical simulation of state estimation

4.6 Battery Management Systems

4.7 Electrochemical Model with Phase Change Electrode

Single Particle Model

Remark 3

Mass Conservation

Lemma 2

4.8 State-of-Charge Estimation

Observer for Phase Transition Positive Electrode

Stability Analysis of the Estimation Error System with Known Interface Position

Theorem 3

Observer for Negative Electrode

4.9 Numerical Simulation

Test with Constant Discharge Input

Comparison with the Extended Kalman Filter

5 Conclusion and Open Problems

References

State Estimation of the Stefan PDE:
A Tutorial on Design and Applications to Polar Ice and Batteries