Sampled-Data Control of the Stefan System

# Sampled-Data Control of the Stefan System

[    [    [
###### Abstract

This paper presents results for the sampled-data boundary feedback control to the Stefan problem. The Stefan problem represents a liquid-solid phase change phenomenon which describes the time evolution of a material’s temperature profile and the interface position. First, we consider the sampled-data control for the one-phase Stefan problem by assuming that the solid phase temperature is maintained at the equilibrium melting temperature. We apply Zero-Order-Hold (ZOH) to the nominal continuous-time control law developed in [23] which is designed to drive the liquid-solid interface position to a desired setpoint. Provided that the control gain is bounded by the inverse of the upper diameter of the sampling schedule, we prove that the closed-loop system under the sampled-data control law satisfies some conditions required to validate the physical model, and the system’s origin is globally exponentially stable in the spatial norm. Analogous results for the two-phase Stefan problem which incorporates the dynamics of both liquid and solid phases with moving interface position are obtained by applying the proposed procedure to the nominal control law for the two-phase problem developed in [30]. Numerical simulation illustrates the desired performance of the control law implemented to vary at each sampling time and keep constant during the period.

UCSD]Shumon Koga, iasson]Iasson Karafyllis, UCSD]Miroslav Krstic

Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, CA 92093-0411 USA

Department of Mathematics, National Technical University of Athens, 15780 Athens, Greece

Key words:  Sampled-data system, Stefan problem, moving boundaries, distributed parameter systems, nonlinear stabilization

## 1 Introduction

### 1.1 Background

Liquid-solid phase transitions are physical phenomena which appear in various kinds of science and engineering processes. Representative applications include sea-ice melting and freezing [28], continuous casting of steel [37], cancer treatment by cryosurgeries [38], additive manufacturing for materials of both polymer [26] and metal [6], crystal growth [7], lithium-ion batteries [25], and thermal energy storage systems [43]. Physically, these processes are described by a temperature profile along a liquid-solid material, where the dynamics of the liquid-solid interface is influenced by the heat flux induced by melting or solidification. A mathematical model of such a physical process is called the Stefan problem[14], which is formulated by a diffusion PDE defined on a time-varying spatial domain. The domain’s length dynamics is described by an ODE dependent on the gradient of the PDE state. Apart from the thermodynamical model, the Stefan problem has been employed to model several chemical, electrical, social, and financial dynamics such as tumor growth process [11], domain walls in ferroelectric thin films [36], spreading of invasive species in ecology [9], information diffusion on social networks [34], and optimal exercise boundary of the American put option on a zero dividend asset [5].

While the numerical analysis of the one-phase Stefan problem is broadly covered in the literature, their control related problems have been rarely addressed. In addition to it, most of the proposed control approaches are based on finite dimensional approximations with the assumption of an explicitly given moving boundary dynamics [8, 2]. For control objectives, infinite-dimensional approaches have been used for stabilization of the temperature profile and the moving interface of a 1D Stefan problem, such as enthalpy-based feedback [37] and geometric control [35]. These works designed control laws ensuring the asymptotical stability of the closed-loop system in the norm. However, the results in [35] are established based on the assumptions on the liquid temperature being greater than the melting temperature, which must be ensured by showing the positivity of the boundary heat input.

Recently, boundary feedback controllers for the Stefan problem have been designed via a “backstepping transformation” [31, 41] which has been used for many other classes of infinite-dimensional systems. For instance, [21] designed a state feedback control law by introducing a nonlinear backstepping transformation for moving boundary PDE, which achieved the exponentially stabilization of the closed-loop system in the norm without imposing any a priori assumption. Based on the technique, [22] designed an observer-based output feedback control law for the Stefan problem, [23] extended the results in [21, 22] by studying the robustness with respect to the physical parameters and developed an analogous design with Dirichlet boundary actuation, [24] designed a state feedback control for the Stefan problem under the material’s convection, [27] developed a control design with time-delay in the actuator and proved a delay-robustness, [29] investigated an input-to-state stability of the control of Stefan problem with respect to an unknown heat loss at the interface, and [30] developed a control design for the two-phase Stefan problem.

The aforementioned results assumed the control input to be varying continuously in time; however, in practical implementation of the control systems it is impossible to dynamically change the control input continuously in time due to limitations of the sensors, actuators, and software. Instead, the control input can be adjusted at each sampling time at which the measured states are obtained or the actuator is manipulated. One of the most fundamental and well known method to design such a “sampled-data” control is the so-called “emulation design” that applies “Zero-Order-Hold” (ZOH) to the nominal “continuous-time” control law. A general result for nonlinear ODEs to guarantee the global stability of the closed-loop system under such a ZOH-based sampled-data control was studied in [16], and the sampled-data observer design under discrete-time measurement is developed in [17] by introducing inter-sampled output predictor. As further extensions, the stability of the sampled-data control for general nonlinear ODEs under actuator delay is shown in [18, 19] by applying predictor-based feedback developed in [32], and results for a linear parabolic PDE are given in [20] by employing Sturm-Liouville operator theory. The sampled-data control for parabolic PDEs has been intensively developed by Fridman and coworkers by utilizing linear matrix inequalities [1, 12, 13, 39]. However, none of the existing work on the sampled-data control has studied the class of the Stefan problem described by a parabolic PDE with state-dependent moving boundaries “(a nonlinear system)”.

### 1.2 Contributions and results

This paper presents the first theoretical result for the sampled-data boundary feedback control for the Stefan problem. The approach employed in this paper is distinct from the methodology developed in literature. Namely, we solve the growth of the system’s energy analytically in time under the proposed sampled-data feedback control that is in the form of an energy-shaping design. Then, a perturbation that is incorporated in the closed-loop system due to the error between the continuous-time design and the sampled-data design can be represented analytically, and the closed-loop stability is proven by using Lyapunov method.

First, we consider the one-phase Stefan problem by assuming that the solid phase temperature is maintained at the melting temperature and focusing on the single melting process. We employ ZOH to the nominal continuous-time feedback controller for the one-phase Stefan problem developed in [23], and prove the required conditions for the model validity and the global exponential stability of the closed-loop system under explicit conditions for the setpoint position and the control gain with respect to the sampling scheduling. Next, we consider the two-phase Stefan problem by incorporating the dynamics of the solid phase temperature and prove the analogous results for the sampled-data control for the two-phase Stefan problem. The results established in this paper hold for arbitrary sampling schedules, and not necessarily uniform sampling schedules.

### 1.3 Organization

The mathematical model the one-phase Stefan problem for a single phase change is presented in Section 2 with stating some important properties. The sampled-data control law and the stability proof of the closed-loop system is given in Section 3. The extension of the presented procedure to the two-phase Stefan problem is described in Section 4. The numerical simulation of the proposed control law is provided in Section 5. The paper ends with the concluding remarks in Section 6.

## 2 Description of the One-Phase Stefan Problem

Consider a physical model which describes the melting or solidification mechanism in a pure one-component material of length in one dimension. In order to mathematically describe the position at which phase transition occurs, we divide the domain into two time-varying sub-domains, namely, the interval which contains the liquid phase, and the interval that contains the solid phase. A heat flux enters the material through the boundary at (the fixed boundary of the liquid phase) which affects the liquid-solid interface dynamics through heat propagation in liquid phase. As a consequence, the heat equation alone does not provide a complete description of the phase transition and must be coupled with the dynamics that describes the moving boundary. This configuration is shown in Fig. 1.

Assuming that the temperature in the liquid phase is not lower than the melting temperature of the material , the energy conservation and heat conduction laws yield the heat equation of the liquid phase as follows

 Tt(x,t) =αTxx(x,t),α:=kρCp,0≤x≤s(t), (1)

with the boundary conditions

 −kTx(0,t) =qc(t), (2) T(s(t),t) =Tm, (3)

and the initial values

 T(x,0)=T0(x),s(0)=s0 (4)

where , , , , and are the distributed temperature of the liquid phase, the manipulated heat flux, the liquid density, the liquid heat capacity, and the liquid heat conductivity, respectively. Moreover, the local energy balance at the liquid-solid interface yields

 ˙s(t)=−βTx(s(t),t),β:=kρΔH∗ (5)

where represents the latent heat of fusion.

###### Remark 1

As the moving interface depends on the temperature, the problem defined in (2)–(2) is nonlinear.

There are two underlying assumptions to validate the model (2)-(2). First, the liquid phase is not frozen to the solid phase from the boundary . This condition is ensured if the liquid temperature is greater than the melting temperature. Second, the material is not completely melt or frozen to single phase through the disappearance of the other phase. This condition is guaranteed if the interface position remains inside the material’s domain. In addition, these conditions are also required for the well-posedness (existence and uniqueness) of the solution in this model. Taking into account of these model validity conditions, we emphasize the following remark.

###### Remark 2

To maintain the model (2)-(2) to be physically validated, the following conditions must hold:

 T(x,t)≥ Tm,∀x∈(0,s(t)),∀t>0, (6) 0< s(t)0. (7)

Based on the above conditions, we impose the following assumption on the initial data.

###### Assumption 1

, for all , and is continuously differentiable in .

The existence and uniqueness of the solution of the one-phase Stefan problem (2)–(2) is presented in [3] as follows.

###### Lemma 1

Under Assumption 1, if is a bounded piecewise continuous function with generating nonnegative heat for a time interval, i.e., for all , then there exists a unique solution for the Stefan problem (2)-(2) with satisfying the condition (2) for all . Moreover, it holds

 ˙s(t)>0,∀t∈[0,¯t]. (8)

## 3 Sampled-Data Control for the One-Phase Stefan Problem

### 3.1 Problem statement and main result

The steady-state solution of the system (2)–(2) with zero manipulating heat flux yields a uniform melting temperature and a constant interface position given by the initial data. In [21], the authors developed the exponential stabilization of the interface position at a desired reference setpoint through the design of as

 qc(t)=−c⎛⎜⎝kαs(t)∫0(T(x,t)−Tm)dx+kβ(s(t)−sr)⎞⎟⎠, (9)

where is the controller gain. However, in practical implementation, the actuation value cannot be changed continuously in time. Instead, by obtaining the measured value as signals discretely in time, the control value needs to be implemented at each sampling time. One of the most typical design for such a sampled-data control is the application of ”Zero-Order-Hold”(ZOH) to the nominal continuous time control law. Through ZOH, during the time intervals between each sampling, the control maintains the value at the previous sampling time. Let be the -th sampling time for , and be defined by

 τj=tj+1−tj. (10)

The application of ZOH to the nominal control law (3.1) leads to the following design for the sampled-data control

 qc(t)= −c⎛⎜ ⎜⎝kαs(tj)∫0(T(x,tj)−Tm)dx+kβ(s(tj)−sr)⎞⎟ ⎟⎠, ∀t∈[tj,tj+1), (11)

of which the right hand side is constant during the time interval . Let us denote for . Hereafter, all the variables with subscript denote the variables at . First, we introduce the following assumptions on the setpoint and the sampling scheduling.

###### Assumption 2

The setpoint is chosen to verify

 (12)
###### Assumption 3

The sampling schedule has a finite upper diameter and a positive lower diameter, i.e., there exist constants such that

 supj∈Z+{τj}≤ R, (13) infj∈Z+{τj}≥ r. (14)

Our main theorem is given next.

###### Theorem 1

Consider the closed-loop system (2)–(2), (2), (3.1) under Assumptions 1, 2. Then for every , there exists a constant for which the following property holds: for every sequence with for which Assumption 3 holds, the initial-boundary value problem (2)–(2) with (3.1) has a unique solution satisfying (2), (2) as well as the following estimate:

 Ψ(t)≤MΨ(0)exp(−bt), (15)

where , for all , in the norm .

The proof of Theorem 1 is established through several steps in the next sections. The positive constant in (15) has a dependency on as

 M(r)=M1+M21−(1−cr)2ecr8, (16)

for some positive constants and that are not dependent on .

### 3.2 Some key properties of the closed-loop system

We first provide the following lemma.

###### Lemma 2

The closed-loop system consisting of the plant (2)–(2) under the sampled-data control law (3.1) has a unique classical solution which is equivalent to the open-loop solution of (2)–(2) with the control law of

 qc(t)=qj= q0j−1∏i=0(1−cτi),∀t∈[tj,tj+1),∀j∈Z+ (17)

where

 q0=−c⎛⎜⎝kαs0∫0(T0(x)−Tm)dx+kβ(s0−sr)⎞⎟⎠. (18)

PROOF. We introduce the following reference error states:

 u(x,t)=T(x,t)−Tm,X(t)=s(t)−sr. (19)

The governing equations (2)–(2) are rewritten as the following reference error system

 ut(x,t)= αuxx(x,t), (20) ux(0,t)= −qc(t)/k, (21) u(s(t),t)= 0, (22) ˙X(t)= −βux(s(t),t). (23)

Define the internal energy of the reference error system as follows:

 ˜E(t)=kαs(t)∫0u(x,t)dx+kβX(t). (24)

Taking the time derivative of (24) along the solution of (20)–(23) leads to the following energy conservation law

 ddt˜E(t)=qc(t). (25)

Noting that is constant for as under ZOH-based sampled-data control, taking the integration of (25) from to yields

 ˜Ej+1−˜Ej=τjqj, (26)

where and . The sampled-data control (3.1) and the internal energy (24) at each sampling time satisfy the following relation:

 qj=−c˜Ej. (27)

Substituting (27) into (26), we obtain

 ˜Ej+1=(1−cτj)˜Ej, (28)

which leads to the explicit solution as follows:

 ˜Ej= ˜E0j−1∏i=0(1−cτi). (29)

Substituting (29) into (27) yields as (17). Therefore, the closed-loop system under the sampled-data feedback control (3.1) is equivalent to the open-loop solution with the control input (17). Moreover, under Assumptions 2, 3, and the fact that , the input (17) is shown to be a bounded piecewise continuous function and for all . Thus, the existence and uniqueness of the solution is ensured by Lemma 1, from which we conclude Lemma 2.

###### Lemma 3

The closed-loop system satisfies the following properties:

 ˙s(t)>0,∀t≥0, (30) s0

PROOF. Combining Lemma 1 with Lemma 2, one can deduce (30), and for all . We show for all . Integrating (25) from to leads to

 ˜E(t)−˜Ej=(t−tj)qj,∀t∈[tj,tj+1). (32)

With the help of (27) and (29), equation (32) yields

 ˜E(t)=(1−c(t−tj))˜Ej,∀t∈[tj,tj+1). (33)

By Assumption 3 and since , we have for all . In addition, for all and for all , it holds . Hence, we have , for all and for all . Applying this to (33) and noting that

 ˜Ej<0,∀j∈Z+, (34)

deduced from (29) and Assumption 2, one can obtain

 ˜E(t)<0,∀t≥0. (35)

Substituting (35) into (24) and applying for all and , we have

 X(t)<0,∀t≥0, (36)

which leads to for all .

### 3.3 Stability analysis

To conclude Theorem 1, this section is devoted to the stability proof of the closed-loop system under the designed sampled-data control law. First, we introduce the backstepping transformation developed in [23] for the continuous-time design, and apply the transformation to the closed-loop system under the sampled-data control in this paper.

#### 3.3.1 State transformation

Introduce the following backstepping transformation

 w(x,t)= u(x,t)−βαs(t)∫xϕ(x−y)u(y,t)dy −ϕ(x−s(t))X(t), (37)

which maps into

 wt(x,t)= αwxx(x,t)+˙s(t)ϕ′(x−s(t))X(t), (38) wx(0,t)= βαϕ(0)u(0), (39) w(s(t),t)= εX(t), (40) ˙X(t)= −cX(t)−βwx(s(t),t). (41)

The objective of the transformation (3.3.1) is to add a stabilizing term in (41) of the target -system which is easier to prove the stability than -system. By taking the derivative of (3.3.1) with respect to and respectively, to satisfy (38), (40), (41), we derive the conditions on the gain kernel solution, and they leads to the following solution:

 ϕ(x)= cβx−ε. (42)

By taking the derivative of the transformation (3.3.1) in and substituting , we have

 wx(0,t)= −qc(t)k−βαεu(0,t)−cαs(t)∫0u(y,t)dy−cβX(t). (43)

Substituting the design of the sampled-data control for all and for all , and recalling the definition of in (24), the boundary condition (43) can be written as

 wx(0,t)= −ck(˜E(t)−˜Ej)−βαεu(0,t). (44)

Moreover, substituting (33), we can describe (44) as

 wx(0,t)= f(t)−βαεu(0,t), (45)

where is an explicit function in time defined by

 f(t)=c2k˜Ej⋅(t−tj),∀t∈[tj,tj+1),j∈Z+. (46)

The closed form representation of (45) using variables is given after the inverse transformation is obtained in the next section.

#### 3.3.2 Inverse transformation

Consider the following inverse transformation

 u(x,t)= w(x,t)−βαs(t)∫xψ(x−y)w(y,t)dy −ψ(x−s(t))X(t). (47)

Taking the derivatives of (3.3.2) in and along (38)-(41), we obtain the gain kernel solution as

 ψ(x)= eλx(psin(ωx)+εcos(ωx)), (48)

where , , , and is to be chosen later. Finally, using the inverse transformation, the boundary condition (39) is rewritten as

 wx(0,t)= f(t)−βαε[w(0,t) −βαs(t)∫0ψ(−y)w(y,t)dy−ψ(−s(t))X(t)⎤⎥ ⎥⎦. (49)

Therefore, the closed form of the target -system is described by (38), (40), (41), and (3.3.2).

#### 3.3.3 Lyapunov method

To show the stability of the original system, first we show the stability of the target system (38), (40), (41), and (3.3.2). For a given , we define the most recent sampling number as

 n:={n∈Z+|tn≤t

and we firstly apply Lyapunov method for the time interval for all , and next for the interval from to . For both cases, we consider the following functional

 V=12α||w||2+ε2βX(t)2, (51)

where denotes norm defined by . Note that Poincare’s and Agmon’s inequalities for the system (38)–(40) with lead to

 ||w||2≤2srε2X(t)2+4s2r||wx||2, (52) w(0,t)2≤2ε2X(t)2+4sr||wx||2. (53)

Taking the time derivative of (51) along with the solution of (38)–(41), (3.3.2), we have

 ˙V= −||wx||2−εβcX(t)2−w(0,t)f(t)+βαεw(0,t)2 −βαεw(0,t)⎡⎢ ⎢⎣βαs(t)∫0ψ(−y)w(y,t)dy+ψ(−s(t))X(t)⎤⎥ ⎥⎦ +˙s(t)α⎛⎜⎝ε22X(t)2+cβs(t)∫0w(x,t)dxX(t)⎞⎟⎠. (54)

Applying Young’s inequality to the second line of (3.3.3) twice, we get

 −w(0,t)f(t)≤γ12w(0,t)2+12γ1f(t)2, (55) −w(0,t)⎡⎢ ⎢⎣βαs(t)∫0ψ(−y)w(y,t)dy+ψ(−s(t))X(t)⎤⎥ ⎥⎦ ≤ +γ2ψ(−s(t))2X(t)2, (56)

where and are parameters to be determined. Applying (55), (56), (52), (53), and Cauchy Schwarz inequalities to (3.3.3) with choosing and , we have

 ˙V≤ −(12−2βsrα(64cs2rα+3)ε)||wx||2 −ε(c8β+g(ε))X(t)2+2srf(t)2 +˙s(t)2α⎛⎜⎝ε2X(t)2+2cβ∣∣ ∣ ∣∣s(t)∫0w(x,t)dxX(t)∣∣ ∣ ∣∣⎞⎟⎠, (57)

where . Since and for all , there exists such that for and . Thus, setting , the inequality (3.3.3) leads to

 ˙V≤ −bV+2srf(t)2+a˙s(t)V, (58)

where

 b=18min{αs2r,c},a=2βεαmax{1,αc2sr2β3ε3}. (59)

Consider the following functional

 W=Ve−as(t). (60)

Taking the time derivative of (60) with the help of (58), we deduce

 ˙W ≤−bW+2srf(t)2e−as(t) ≤−bW+2srf(t)2. (61)

(i) For , for all
Applying comparison principle to (61) for leads to

 W(t)≤W(tj)e−b(t−tj)+2sre−btt∫tjebτf(τ)2dτ. (62)

Setting and recalling , we get

 Wj+1≤ Wje−bτj+2c4srk2e−bτj˜E2jIj, (63)

where , and is defined by

 Ij:=tj+1∫tjeb(τ−tj)(τ−tj)2dτ. (64)

Then, by introducing the variable and integration by substitution, with the help of for all derived by (59), Assumption 3 and the fact that , one can derive the following inequality:

 Ij= 1b3bτj∫0ess2ds≤Jb3, (65)

where is defined by . Applying (65) to (63) yields

 Wj+1≤ Wje−bτj+Bj, (66)

where is defined by

 Bj=2Jc4srk2b3e−bτj˜E2j. (67)

Applying (66) from to inductively, we get

 Wn≤W0e−b∑n−1i=0τi+Bn−1+n−2∑i=0Bie−b∑n−1j=i+1τj. (68)

By (67) and the solution of given in (29), we have

 n−2∑i=0Bie−b∑n−1j=i+1τj ≤ 2Jc4sr˜E20e−b∑n−1j=0τjk2b3(1+n−2∑i=1(i−1∏k=0(1−cτk)2)eb∑i−1j=0τj) ≤ 2Jc4sr˜E20e−b∑n−1j=0τjk2b3(1+n−2∑i=1(i−1∏k=0(1−cτk)2ebτk)). (69)

Since , by using given in Assumption 3, the following inequality holds

 (1−cτi)2ebτi≤(1−cr)2ecr8:=δ<1,∀j∈Z+. (70)

Thus, the inequality (69) leads to

 n−2∑i=0Bie−b∑n−1j=i+1τj≤ 2Jc4sr˜E20e−b∑n−1j=0τjk2b3(1+n−2∑i=1δi) ≤ 2Jc4sr˜E20k2b3(1−δ)e−b∑n−1j=0τj. (71)

In the similar way, we get

 Bn−1≤2Jc4sr˜E20k2b3(1−δ)e−b∑n−1j=0τj. (72)

Recalling that and , we get . Applying (71) and (72) to (68), we arrive at

 Wn≤(W0+A˜E20)e−btn. (73)

where .

(ii) For ,
Applying comparison principle to (61) from to , we get

 W(t)≤ Wne−b(t−tn)+Bne−b(t−tn+1) ≤ Wne−b(t−tn)+A˜E20e−bt. (74)

Finally, combining (73) and (74), the following bound is obtained

 W(t)≤(W0+2A˜E20)e−bt. (75)

Recalling the relation defined in (60), and applying , the norm estimate for in (75) leads to the following estimate for :

 V(t)≤easr(V0+2A˜E20)e−bt. (76)

We consider the -norm of -system defined by

 Ψ(t)=s(t)∫0u(x,t)2dx+X(t)2. (77)

Due to the invertibility of the transformation from to together with the boundedness of the domain , there exist positive constants and such that the following inequalities hold:

 M–––Ψ(t)≤V(t)≤¯¯¯¯¯¯MΨ(t). (78)

Moreover, due to the definition of the reference energy given in (24), using Young’s and Cauchy Schwarz inequalities one can show that

 ˜E20≤KΨ0, (79)

where . Applying (78) and (79) to (76), we deduce that there exists positive constant such that the following inequality holds

 Ψ(t)≤MΨ0e−bt, (80)

which completes the proof of Theorem 1.

## 4 Sampled-Data Design for Two-Phase Stefan Problem

In this section, we extend the results we have established in the previous section to the ”two-phase” Stefan problem, where the temperature dynamics in the solid phase is governed by the heat equation with different physical parameters from the liquid phase, following the work in [30]. This configuration is depicted in Fig. 2.

### 4.1 Problem statement

The governing equations are descried by the following coupled PDE-ODE-PDE system:

 ∂Tl∂t(x,t)= αl∂2Tl∂x2(x,t),0