Flatness Based Trajectory Planning for Two Heat Conduction Problems Occurring in Crystal Growth Technology

Publication e-STA e-STA 2004-1
OAI : oai:www.see.asso.fr:545:2004-1:20068


Flatness Based Trajectory Planning for Two Heat Conduction Problems Occurring in Crystal Growth Technology


278.27 Ko


Creative Commons Aucune (Tous droits réservés)
<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/545:2004-1/20068</identifier><creators><creator><creatorName>Joachim Rudolph</creatorName></creator><creator><creatorName>Jan Winkler</creatorName></creator><creator><creatorName>Frank Woittennek</creatorName></creator></creators><titles>
            <title>Flatness Based Trajectory Planning for Two Heat Conduction Problems Occurring in Crystal Growth Technology</title></titles>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><dates>
	    <date dateType="Created">Sun 1 Oct 2017</date>
	    <date dateType="Updated">Sun 1 Oct 2017</date>
            <date dateType="Submitted">Sun 24 Mar 2019</date>
	    <alternateIdentifier alternateIdentifierType="bitstream">7c18c0caaa7799d3982bc28b9f43146be1f8aa5d</alternateIdentifier>
            <description descriptionType="Abstract"></description>

1 Flatness Based Trajectory Planning for Two Heat Conduction Problems Occurring in Crystal Growth Technology Joachim Rudolph Jan Winkler, and Frank Woittennek Abstract— Two heat conduction problems from crystal growth technology are considered. The growth process studied as the first problem comprises two phases (solid crystal and liquid melt) which are separated by a moving interface; affine temperature dependencies of the parameters are taken into account. Planning trajectories for finite time transitions between two regimes is considered. This is possible using recently developed “flatness based” methods leading to series representations. The series are directly parameterized by the flat output trajectories. A proof of series convergence is given. The interest of the second problem, called tempering process, resides in the fact that here a fixed temperature at a moving point is prescribed. An “indirect” solution is proposed where a flat output trajectory is designed such as to approximately achieve that goal. Index Terms— Boundary control, crystal growth, distributed parameter systems, flatness, free boundary, heat equation, tra- jectory planning. I. INTRODUCTION Two key steps during the production of semi-conductors are the growth of mono-crystals from melt and the tempering of the grown crystals for improving their quality. Especially III- V-semi-conductors like Gallium-Arsenide (GaAs) are grown using the so called Vertical-Gradient-Freeze (VGF) technique, where the melt is cooled down in a crucible in such a way that the crystallization starts at the bottom and moves upwards until the whole melt is solidified. Both, crystal growth and tempering, are typical representa- tives of processes where the heat conduction problems play a crucial role. However, while the first one involves a two phase system (melt/crystal) with a free (moving) boundary and two control inputs at the top and bottom of the crucible, the latter one consists of one phase and one control input only. In both cases a fixed temperature should be assigned at a moving point, the trajectory of which is prescribed. In order to achieve this the trajectories of the control inputs can be calculated using “flatness based” methods, partially combined with optimization techniques. Differential flatness is a concept which is very useful in the trajectory planning and feedback design for nonlinear finite dimensional systems, i.e., systems described by ordinary differential equations. The flatness based control methods place an emphasis on trajectory design and open-loop control The authors are with the Institut für Regelungs- und Steuerungstheorie, TU Dresden, 01069 Dresden, Germany; Jan Winkler also with the Institut für Kristallzüchtung Berlin. [1]–[4]. This aspect gains even more importance in infinite dimension, namely for distributed parameter systems with boundary control action, the mathematical models of which comprise partial differential equations, including the subclass of (linear and nonlinear) time delay systems. As a consequence, flatness based methods have been devel- oped for distributed parameter systems during the last years; for a survey see [5], [6], e.g. In the particular case of heat conduction problems the system trajectories are described by series which are parameterized by trajectories of a finite number of lumped variables, called a flat output. The flat output trajectory occurs together with all its derivatives, and it must be a Gevrey function of an order less than 2 in order to get series convergence. Such representations of solutions of partial differential equations have been known for a long time [7], [8]. Trajectory planning problems modeled with the heat equa- tion and other second order linear or quasi-linear parabolic equations have been treated using this flatness based approach in [9]–[16] and others. Similar methods have been applied on fourth order beam equations, see e.g. [17]–[20]. II. VERTICAL-GRADIENT-FREEZE CRYSTAL GROWTH A. Process overview Gallium-Arsenide (GaAs) crystals are mainly produced us- ing the Vertical-Gradient-Freeze (VGF) method. A crucible (cf. Fig. 1) is filled with liquid GaAs. The temperatures at the top and the bottom, possibly also the jacket, of the crucible are manipulated in such a way that the crystallization starts at the bottom of the crucible, slowly moving upwards until the whole liquid is solidified. The plant sketched in the left scheme in Fig. 1 has only two heaters at the top and the bottom of the crucible while the middle one is equipped with four heaters surrounding the crucible in addition. As an alternative the whole crucible may be moved in a fixed temperature field as shown in the right sketch in Fig. 1 (known as the Vertical- Bridgman plant). The crystal growth process consists of two stages. First we reach a temperature profile in the melt from which crystalliza- tion can start immediately. Crystallization then starts at the beginning of the second phase, which ends once the whole melt has become solid. B. Mathematical Model We assume that, due to the small Peclet-number of GaAs and its opaque nature, the convection and the radiation within 2                                                                                                                                                             Crystal Crystal Crystal Melt Heater Melt T op − Heater Shaft Bottom Heater Melt Interface P hase Interface P hase Interface P hase Fig. 1. Different types of crystal growth plants.           H TM TI TC Crystal Front Crystallization Crucible Melt R vI r ζ = −zI (t) ζ = H − zI (t) r Bottom Heater Top Heater ζ z Fig. 2. Sketch of a Vertical-Gradient-Freeze plant. the crystal are negligible. Under the additional assumptions that the crucible is a cylinder of height H and inner radius R (cf. Fig. 2), and its jacket is ideally isolated we may neglect angular and radial gradients and can, therefore, model the process using the one-dimensional heat equation ∂ ∂t  ρ T(z, t)  T(z, t)  = ∂ ∂z  λ(T(z, t)) ∂T ∂z (z, t)  . (2.1) Here ρ is the volume-specific heat capacity, λ the thermal conductivity, T is the temperature, z the position, and t is the time. At the bottom and the top of the crucible, where the heaters reside, we have radiation. The boundary conditions resulting from the heat flow, thus, read λ(T) ∂T ∂z = αb(T − ub) + εbσS T4 − u4 b  , z = 0 (2.2a) −λ(T) ∂T ∂z = αt(T − ut) + εtσS T4 −u4 t  , z = H, (2.2b) where t 7→ ub(t) and t 7→ ut(t) are the temperatures of the top ρC temperature dependent volume-specific heat capacity of solid GaAs: ρC (TC ) = ρC,0 + ρC,1TC , (ρC,0, ρC,1 ∈ R+) λC temperature dependent heat conduction coefficient of solid GaAs: λC (TC ) = λC,0 +λC,1TC , (λC,0, λC,1 ∈ R+) ρM,0 volume-specific heat capacity of liquid GaAs λM,0 heat conduction coefficient of liquid GaAs αb, αt heat transfer coefficient (bottom, top) b, t emissivity (bottom, top) σS Stephan-Boltzmann number h volume-specific crystallization enthalpy of GaAs (h < 0) TABLE I PARAMETERS OF THE CRYSTAL GROWTH PROCESS MODEL. and the bottom heater, respectively. The constant parameters given in Tab. I are used1 . The temperature at the interface is equal to the melting (or crystallization) temperature TI: T(zI(t), t) = TM (zI(t), t) = TC(zI(t), t) = TI, (2.3a) where zI denotes the actual location of the interface. Due to the energy release resulting from the crystallization, at the interface the temperature gradient is discontinuous λM (TI) ∂TM ∂z (zI(t), t) − λC(TI) ∂TC ∂z (zI(t), t) = vI(t)h. (2.3b) Here the indices M and C are used in order to distinguish the variables and parameters of the melt from those of the crystal2 . The initial conditions for the process are given by ∂TM ∂t (z, 0) = 0, TM (z, 0) = TM,0 ∈ R, z ∈ [0, H]. (Remember that in the beginning there is no crystal.) 1Note, that the above equations are valid for both the melt and the crystal. They will be adapted in the following sections depending on the stage of the process considered. 2These compatibility conditions at the crystallization front will be used only during the second stage, when the crystal is present. 3 C. Adjusting a temperature profile in the melt As mentioned above, in the first stage we want to adjust a temperature profile from which crystallization can start, i.e., a profile increasing with z such that at the bottom of the crucible the temperature equals TI. 1) Adapted model for the melt: As the temperature in the melt is almost constant during the process, we neglect the temperature dependency of the parameters λ and ρ. Therefore, the differential equation (2.1) may be simplified to ρM,0 ∂TM ∂t (z, t) = λM,0 ∂2 TM ∂z2 (z, t), (z, t) ∈ [0, H] × R+ . (2.4a) From (2.2) we obtain adapted boundary conditions λM,0 ∂TM ∂z = αb(TM − ub) + εbσS T4 M −u4 b  , z = 0 (2.4b) −λM,0 ∂TM ∂z = αt(TM − ut) + εtσS T4 M −u4 t  , z = H (2.4c) with constant parameters ρM,0 and λM,0. 2) Formal power series solution: The temperature profile in the melt can be expressed using the power series TM (z, t) = ∞ X k=0 ak(t) zk k! , ak+2(t) = ρM,0 λM,0 ȧk(t), (2.5) where the recursion for the coefficients is obtained by sub- stituting the series into the differential equation (2.4a) and comparing the coefficients of like powers of z. 3) Parameterization of the process: Up to now we have not yet fully specified the coefficients of the series; nothing has been said about the initialization of the recursion for the series coefficients (i.e., about a0 and a1). This freedom is now used to parameterize a transition by introducing a flat output. We start the process at a time t = 0 where the whole GaAs is molten, assuming a uniform constant initial temperature TM,0 > TI. At t = t1 we want to reach another stationary profile, at which the temperature at the bottom is equal to the crystallization temperature TI, and the temperature gradient has a desired value ∆TM,1. Hence, the stationary profiles to be connected are given by TM (z, 0) = TM,0 = const., z ∈ [0, H], TM,0 > TI (2.6a) TM (z, t1) = TI + ∆TM,1z, z ∈ [0, H], ∆TM,1 > 0. (2.6b) Since we have no constraints resulting from boundary conditions, the first two coefficients of the power series (2.5) can be freely chosen and, therefore, form a flat output y(t) =  y1(t) y2(t)  =  a0(t) a1(t)  =  TM (0, t) ∂TM ∂z (0, t)  . (2.7) From given trajectories for the flat output we can calculate the temperatures and the temperature gradients using the power series (2.5). This yields TM = ∞ X n=0  ρM,0 λM,0 n  z2n (2n)! y (n) 1 + z2n+1 (2n + 1)! y (n) 2  . (2.8a) Now the trajectories for the controls ub and ut can be calculated using (the algebraic) equations (2.4b) and (2.4c), which are solved numerically. Since the system trajectories depend on derivatives of the flat output of arbitrary order, those trajectories have to be chosen C∞ . As an analytic function would be completely defined by these initial conditions, we must use non-analytic functions to parameterize the solution. The convergence of the series can be ensured using functions of Gevrey order α ≤ 2 [7], [12], [21]. We may use the Gevrey function Φσ,T (t) = 1 2  1 + tanh  2t/T − 1 (4t/T (1 − t/T)) σ  , (2.9) which realizes a smooth transition from 0 to 1 on the interval [0, T]. The parameter σ allows us to vary the slope of the transition. We have α = 1 + 1/σ for the Gevrey order. D. The crystallization We will now indeed study the second stage of the process where the crystal is growing and we have a two-phase system. 1) Adapted mathematical model: Starting from equations (2.1) – (2.3), we have to adapt the model equations for the melt and the crystal, being valid during the growth. As in the previous section, we assume the heat capacity ρM,0 and the heat conduction coefficient λM,0 in the melt do not depend on the temperature. For the solid phase we use an affine approximation for these parameters: ρC(TC) = ρC,0 + ρC,1TC, λC(TC) = λC,0 + λC,1TC. Using a moving frame the origin of which is located at the interface between solid and liquid GaAs, we apply the transformation ζ = z − zI(t), ϑ(ζ, t) = T(z, t) which, with żI = vI, leads to  ∂ ∂t − vI(t) ∂ ∂ζ  ρ ϑ(ζ, t)  ϑ(ζ, t)  = ∂ ∂ζ  λ ϑ(ζ, t)  ∂ ∂ζ ϑ(ζ, t)  . Using the constant parameter approximations for the melt and the affine ones for the crystal we obtain ρM,0  ∂ϑM ∂t − vI ∂ϑM ∂ζ  = λM,0 ∂2 ϑM ∂ζ2 (2.10a) (ρC,0 +2ρC,1ϑC)  ∂ϑC ∂t − vI ∂ϑC ∂ζ  = (λC,0 +λC,1ϑC) ∂2 ϑC ∂ζ2 + λC,1  ∂ϑC ∂ζ 2 . (2.10b) 4 The boundary conditions for the second stage of the process are given by λM,0 ∂ϑ ∂ζ = αb(ϑM − ub) + εbσS ϑ4 M − u4 b  (2.11a) at ζ = −zI, and −(λC,0 + λC,1ϑC) ∂ϑ ∂ζ = αt(ϑC − ut) + εtσS ϑ4 C − u4 t  (2.11b) at ζ = H − zI. Finally, the adapted compatibility conditions read ϑC = ϑM = TI, ζ = 0 (2.12a) vIh = λM,0 ∂ϑM ∂ζ − (λC,0 + λC,1TI) ∂ϑC ∂ζ , ζ = 0. (2.12b) 2) Formal power series solution: Again we use power series for planning trajectories of the growth process: ϑM (ζ, t) = ∞ X k=0 bk(t) ζk k! (2.13a) ϑC(ζ, t) = ∞ X k=0 ck(t) ζk k! . (2.13b) Substituting (2.13a) into the differential equation (2.10a) we obtain the recursion formula bk+2 = %M,0 λM,0 ḃk − vIbk+1  , k ≥ 0. (2.14a) Furthermore, using (2.10b) the recursion for the coefficients ck in (2.13b) reads3 ck+2 = 1 λC,0 + λC,1c0 " ρC,0 (ċk − vIck+1) + k X l=0  k l  2ρC,1 (clċk−l − vIclck+1−l) − λC,1 k + 1 l + 1 cl+1ck+1−l # , k ≥ 0. (2.14b) 3) Parameterization of the process: During the second stage of the process the temperature at ζ = 0 has to be equal to the crystallization temperature of GaAs (recall that our frame is moving with the crystallization front) and, therefore, b0 = c0 ≡ TI. Moreover, the coefficients b1, c1 and the velocity żI = vI have to be chosen such that the second compatibility condition (2.12b) is satisfied. Thus, the process can be parameterized, for instance, by ỹ(t) =  ỹ1(t) ỹ2(t)  =  zI(t) b1(t)  =  zI(t) ∂ϑM ∂ζ (0, t)  . (2.15) The trajectories for the other variables can be calculated from the flat output trajectory t 7→ ỹ(t) using the compatibility con- dition (2.12b), the series solutions (2.13) with the recursions (2.14a) and (2.14b), and the boundary conditions (2.11). 3We assume, that the heat conduction coefficient is non-zero (λC,0 + λC,1c0 > 0). We assume that when crystallization starts the crystallization front is stationarily located at the bottom of the crucible. The initial conditions are now those given in (2.6b) completed by those for the crystal, which due to the compatibility condition and the zero initial velocity are equal to those for the melt. We want to realize a transition to another “stationary” solution4 of (2.10), to be reached at t = t2, which satisfies the conditions żI(t2) = v2, v2 > 0 (2.16a) ∂ϑM ∂ζ (0, t2) = ∆TM,2, ∆TM,2 > 0 (2.16b) ∂ϑC ∂ζ (0, t2) = λM,0∆TM,2 − v2h λC,0 + λC,1TI . (2.16c) This can be achieved using the trajectories ỹ(t) = Φσ,t2 (t)  v2t ∆TM,2  for the flat output where Φσ,t is the function defined in (2.9). Alternatively, it is possible to start the crystallization at t = t1 with a constant initial velocity v2. To this end, we have to change the parameterization of the first stage. The ordinary differential equation for the melt describing the desired “stationary” regime is obtained from (2.10a) by setting ∂ϑM ∂t to zero: −v2ρM,0 ∂ϑ̃M ∂ζ (ζ) = λM,0 ∂2 ϑ̃M ∂ζ2 (ζ). The “stationary” solution of the second stage is given by ϑ̃M (ζ) = TI + λM,0 ρM,0v2 ∆TM,2  1 − exp  −ζv2 ρM,0 λM,0  . (2.17) In order to start crystallization with a prescribed constant velocity, we have to meet this solution at the end of the first stage of the process. Using fixed coordinates (2.17) reads ϑ̃M (ζ) = TM (z, t) = TI + λM,0 ρM,0v2 ∆TM,2  1 − exp  −(z − zI(t))v2 ρM,0 λM,0  . (2.18) Since żI = v2 is constant and zI(t1) = 0, at z = 0 we have TM = TI + λM,0 ρM,0v2 ∆TM,2  1−exp  (t − t1)v2 2 ρM,0 λM,0  | {z } =:b y1(t) ∂TM ∂z = ∆TM,2 exp  (t − t1)v2 2 ρM,0 λM,0  | {z } =:b y2(t) , which are the ficticious analytic trajectories to be reached by y1 and y2 at time t = t1. In order to connect the profiles (2.6a) and (2.18) we have to choose trajectories for y1 and y2 all the derivatives of which 4Here “stationary” means that the temperature profile is constant for an observer moving with the crystallization front. 5 at t = t1 are equal to those of b y1 and b y2, respectively. For instance, we may choose y1(t) = (b y1(t) − TM,0) Φσ,t1 (t) + TM,0 y2(t) = b y2(t)Φσ,t1 (t). Again Φσ,T is the function defined in (2.9). 4) Convergence: It can be shown, as in the previous sec- tion, that the series giving the solution for the melt converges everywhere provided the trajectories are chosen Gevrey of order less than 2. For the series solution of the solid phase, which results from a nonlinear heat equation, convergence can be ensured by using Gevrey functions of order α < 2, but contrary to the linear case we can prove a finite radius of convergence only. Assume that the trajectories for vI and c1 are Gevrey of order α < 2, i.e., |v (l) I | ≤ mv (l!)α γl , |c (l) 1 | ≤ m1 (l!)α γl , (2.19) where l ∈ N, 1 ≤ α ≤ 2, m1, mv, γ ∈ R+ . Then the radius of convergence is larger than the unique positive root of the polynomial (cf. the Appendix) p(R̄) = R̄3 4m1|ρ1| γ + R̄2  4 3 |b ρ0| γ + 4m1|ρ1|mv  + R̄  3 2 mv|b ρ0| + 4m1|λ1|  − |b λ0| (2.20) with b ρ0 := ρ0 + 2ρ1TI, b λ0 := λ0 + 2λ1TI. This additional condition bounds the maximal values admis- sible for v2 and ∆TC,2 (and therefore for ∆TM,2). Moreover, the transition time t2 must be chosen large enough. E. Simulation result 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 1.2 1.4 Melt Crystal z t ub ut T Fig. 3. Temperature profile during a crystallization process. Fig. 1 shows the temperature profile which is obtained using the flatness based motion planning strategies for the parameter- ization of the crystallization process. The computations have been done using normalized coordinates. The bold lines on the left and the right of the plot show the associated trajectories of the temperatures of the bottom heater (ub ) and the top heater (ut), respectively, while the third bold line marks the location of the crystallization front zI(t). III. TEMPERING OF CRYSTALS A. Process overview and control problem                             u(t) z = −L z = 0 z∗(t) Crystal T(z∗(t), t) T(z, t) v∗(t) Fig. 4. Crystal with ideal isolation and one boundary control u. Grown single-crystals have to be reworked to reduce stress in the crystalline structure. Therefore, the whole crystal is re-heated up to a temperature slightly lower than its melting temperature. Then the crystal is cooled down from one side to the other. During the process it is required that the temperature at a moving point is prescribed. The heat flow enters the crystal only from one of its front ends. As in the previous section we assume that the crystal is ideally isolated (except at the front where the heater resides), and therefore, the radial temperature gradients are negligible. Under these assumptions the model reads (cf. (2.1)) λ0 ∂2 T ∂z2 (z, t) = ρ0 ∂T ∂t (z, t), t ≥ 0, z ∈ [−L, 0] (3.1) with heat conduction coefficient λ0 and volume-specific heat capacity ρ0. The boundary conditions are given by ∂T ∂z (0, t) = 0 (3.2a) λ0 ∂T ∂z (−L, t) = αL (T(−L, t) − u(t)) , (3.2b) where the first one is due to the ideal isolation at z = 0, and the second one results from the control input u(t) being the temperature of the heater. Here αL denotes the heat transfer coefficient. Starting from stationary initial conditions T(z, 0) = T0, z ∈ [−L, 0], (3.3) we want to heat up the crystal to a constant predefined temperature T1: T(z, t) = T1, for t = 0 (3.4) 6 Having reached this stationary profile the crystal has to be cooled down to T2 while prescribing the temperature T2 at a moving point z = z∗(t): T(z, t) = T2, for t ≥ t∗ , (3.5) with T(z∗(t), t) = T2 = const. (3.6) B. Series solution and parameterization From the differential equation (3.1) and the boundary con- dition (3.2a) we obtain the power series solution T(z, t) = ∞ X k=0  ρ0 λ0 k y(k) (t) z2k (2k)! , (3.7) where y(t) = T(0, t) is the flat output. The re-heating of the crystal can be parameterized by choosing the trajectory [−t̄, 0] 3 t 7→ y(t) Gevrey of order α < 2 with y(−t̄) = T0, y(0) = T1, y(k) (−t̄) = y(k) (0) = 0 (k > 0). The second part of the control problem (on [0, t∗ ]) is more difficult to handle since we cannot parameterize the system by the trajectory t 7→ T(z∗(t), t), but by the temperature T(0, t) rather. To overcome this problem we use an optimization strategy, detailed in the next section. C. Solution via optimization Obviously, it is not possible to achieve T(−L+v∗t, t) = T2 for v∗ = const. and t ∈ [0, t∗ ]. The initial/boundary condition T(−L, 0) = T1 forms an obstruction if T1 6= T2. We may, however, assign T(−L + v∗t, t) = T2 with v∗ = const. on the interval [ε, t∗ ]. Since the system can be parameterized by the flat output and its derivatives, the control problem which has to be solved consists in choosing a trajectory for the flat output y such that a point having a fixed temperature T2 moves through the crystal with prescribed velocity v∗. Having found such a trajectory we are in a position to calculate the associated control input5 . A possible choice for the flat output trajectory is y(t) = T1 + (T2 − T1)Φ(0) (t) | {z } transient + P(t; γ)Φ(1) (t) | {z } for parameterizing , (3.8) where Φ is a Gevrey function realizing a smooth transi- tion from 0 to 1 in time t∗ , i.e., with all the initial and final derivatives being equal to zero. Moreover, P(t; γ) is a function depending linearly on a parameter vector γ = (γ0, γ1, . . . , γR−1), R ∈ N+ , i.e., P(t; γ) = R−1 X r=0 γrBr(t), where C∞ 3 Br : [0, t∗ ] → R, (r = 0, . . . , R − 1) is the rth function of an appropriate function base. Note that, for any choice of the parameter vector γ we obtain a solution 5Note, that this problem is easier to handle than the “direct” approach which consists in choosing a trajectory for the input u(t) such that the additional constraints are satisfied. Here we still avoid integration. of the differential equation (3.1) satisfying the initial an final conditions as well as the boundary condition (3.2a). For simplicity we set T1 = 1 and T2 = 0 and obtain y(t) = Φ(t) |{z} transient + Φ(1) (t) R−1 X r=0 γrBr(t) | {z } for parameterizing (3.9) the kth derivative of which can be written as y(k) (t) = Φ(k) (t) + k X i=0  k i  Φ(i+1) (t) R−1 X r=0 γrB(k−i) r (t). (3.10) Using (3.10) in (3.7), after some elementary calculations, we obtain T(z, t; γ) = ∞ X k=0  ρ0 λ0 k Φ(k) (t) z2k (2k)! + R−1 X r=0 γr ∞ X k=0 k X i=0  k i   ρ0 λ0 k z2k (2k)! Φ(i+1) (t)B(k−i) r (t). (3.11) Using a least squares method to perform the optimization, we have to minimize the integral Q̄(γ) = Z t∗ ε T(z∗(t), t; γ) − T2 2 dt with respect to the parameter vector γ. For numerical eval- uation we discretize this integral using N ≥ R points and obtain Q(γ) = N X n=1 T(z∗(tn), tn; γ) − T2 2 . D. Simulation results Fig. 5 and Fig. 6 illustrate our optimization strategy. As function base for optimization we have used Legendre poly- nomials up to order 30. All the calculations have been done using normalized coordinates (including the fact that the final value of T is +1). The diagonal line in Fig. 5 marks the point z∗(t) where the temperature is prescribed. An alternative, maybe even better, illustration is given in Fig. 6 which shows several snapshots of the temperature profile. Here the desired temperatures at the point z∗(t) are marked by circles. After the beginning of the process, which has been excluded from optimization, the temperature at z∗(t) is indeed equal to the desired one. IV. CONCLUSION Two heat equation problems have been solved using the flatness based approach leading to series representations. One is a free boundary problem for a system with two phases (liquid and solid) — see [8], [16] for a similar approach to a related problem. The other one requires an approximation of a trajectory of a variable which does not play the role of a flat output. 7 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 z t Fig. 5. Temperature profile during the second stage of a tempering process. −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 z T t = 0.1 t = 0.2 t = 0.3 t = 0.4 t = 0.5 t = 0.6 Fig. 6. Snapshots of the temperature profile during the second stage of the tempering process. An interesting extension is the consideration of a non-planar phase boundary. Again a transformation in a frame moving with the boundary is useful; this leads to curved coordinates. A flatness based parameterization seems being possible if distributed control is admitted, see [22] where a formal series solution is derived. However, while series convergence has been proven here for the planar case, this proof is outstanding for the case of a bent interface. APPENDIX The following lemmas are required in the proof of conver- gence of the series (2.13) in Section II. (Most of the proofs of these Lemmas are trivial, references are given for the others.) Lemma 1.1: For real numbers α > 1, Lk ≥ 0, l ∈ N (cf. [7]) l X k=0 Lα k ≤ l X k=0 Lk !α . Lemma 1.2: For i, j, l ∈ N (cf. [12]) i!j!(i + j + l + 1)! (i + j + 1)! = l X r=0  l r  (j + r)!(i + l − r)!. Lemma 1.3: For 0 < a ≤ b, k > 0 a b ≤ a + k b + k . Lemma 1.4: For real numbers a > b ≥ 0, the function f : R+ → R f(x) = x + a x + b is monotonically decreasing, i.e., for c ∈ R+ f(x) > f(x+c). Lemma 1.5: For real α ≥ 0 and non-negative integers l, k, c  (l + k)! k! α ≤  (l + k + c)! (k + c)! α . Lemma 1.6: For integers n1 ≤ n2, i > −n1 n2 X k=n1 1 (k + i)(k + i + 1) = n2 X k=n1  1 k + i − 1 k + i + 1  = n2 X k=n1 1 k + i − n2+1 X k=n1+1 1 k + i = 1 n1 + i − 1 n2 + i + 1 . We are now in a position to state the convergence result, the proof of which is inspired by those in [12]. To this end, it is advantageous to use the relative temperature w.r.t. TI. This yields exactly the same recursion formula as given in (2.14b) but with the first coefficient in the series being equal to zero and the adapted constant parameters b ρ0 := ρ0 + 2ρ1TI, b λ0 := λ0 + 2λ1TI. Theorem 1.1: Consider ρ0, ρ1, λ0, λ1 ∈ R, γ, m, mv ∈ R+ , 1 ≤ α ≤ 2, and a trajectory t 7→ vb(t) satisfying |v (l) I | ≤ mv (l!)α γl , 1 ≤ α ≤ 2, mv, γ ∈ R+ . Define the formal power series6 X n≥0 anxn (1.1) via the recursion an+2 = 1 λ0 " ρ0 ȧn (n + 1)(n + 2) − ρ0vIan+1 (n + 2) + n X k=0  2ρ1  akȧn−k − vIakan+1−k(n + 1 − k) (n + 1)(n + 2)  − λ1 ak+1an+1−k(n − k + 1) (n + 2) # , (1.2) initialized by a0 and a1 which satisfy |a (l) 1 | ≤ m 2 (l!)α γl , l ∈ N, 1 ≤ α ≤ 2, γ > R+ . 6The recursion formula used here is obtained from the one given in (2.14b) in Section II by setting ak = ck/k!. 8 Then, the unique positive root of the polynomial7 R3 2m|ρ1| γ + R2  4 3 |ρ0| γ + 2m|ρ1|mv  + R  3 2 mv|ρ0| + 2m|λ1|  − |λ0| (1.3) determines a lower bound of the radius of convergence. Proof: Assuming that the coefficients an satisfy a0 = 0, |an| ≤ m Rn−1 1 n(n + 1) , n ≥ 1 (1.4) we calculate the radius of convergence R̄ of (1.1) using the Cauchy-Hadamard formula: 1 R̄ = lim n→∞ |an| 1 n ≤ lim n→∞  m Rn−1 1 n(n + 1)  1 n = 1 R . Hence, if (1.4) holds (which will be shown below), the series (1.1) converges absolutely, for |x| < R. Next, we prove by induction that, provided R is the positive root of (1.3), the inequality (1.4) holds. More precisely, we show that the l-th derivative (l ∈ N) of the coefficient an satisfies a(l) n ≤ An,l with An,l = m Rn−1γl ((l + n − 1)!)α ((n − 1)!)α 1 n(n + 1) . (1.5) By the assumption, the coefficients a0 and a1 satisfy the above inequality. This initializes the proof. Assuming that (1.5) holds for 0 ≤ k ≤ n + 1, we show that the assertion is true for a (l) n+2 as well. To this end, we re-write the recurrence relation (1.2) as an+2 = 1 λ0 5 X j=1 Bj,n+2, where B1,n+2 = ρ0ȧn (n + 1)(n + 2) B4,n+2 = −2ρ1 n X k=1 (n − k + 1)vIakan+1−k (n + 1)(n + 2) B2,n+2 = − vIan+1ρ0 n + 2 B5,n+2 = −λ1 n X k=0 n + 1 − k (n + 2) ak+1an+1−k B3,n+2 = 2ρ1 n−1 X k=1 akȧn−k (n + 1)(n + 2) . We use the triangle inequality |a (l) n+2| ≤ 1 |λ0| 5 X j=1 |B (l) j,n+2|, and majorize the terms |B (l) 1,n+2|, . . . , |B (l) 5,n+2|. 7We use R here, there should be no confusion with the notation used in Section II where it denotes the radius R of the cylinder. For n ≥ 1 the term B1,n+2 satisfies |B (l) 1,n+2| ≤ |ρ0|An,l+1 (n + 1)(n + 2) ≤ m|ρ0| Rn−1γl+1 ((l + n)!)α ((n − 1)!)α 1 n(n + 1)2(n + 2) = m|ρ0| Rn−1γl+1 ((l + n)!)α (n!)α nα n(n + 1)2(n + 2) . We proceed using α ≤ 2, and the Lemmas 1.5, 1.3, and 1.4: |B (l) 1,n+2| ≤ An+2,l |ρ0|R2 γ  n(n + 3) (n + 1)2 ff ≤ An+2,l |ρ0|R2 γ  (n + 3) (n + 2) ff ≤ 4 3 |ρ0|R2 γ An+2,l. Note that B1,2 = 0 and, therefore, it satisfies the above inequality trivially. To bound the term |B2,n+2| we have to expand the deriva- tives of a product using the Leibnitz formula. In order to calculate a bound for the resulting sum, we make use of the Lemmas 1.1 and 1.2. We have |B (l) 2,n+2| = ˛ ˛ ˛ ˛ ˛ ρ0 l X r=0 l r ! v (l−r) I a (r) n+1 n + 2 ˛ ˛ ˛ ˛ ˛ ≤ ˛ ˛ ˛ ˛ ˛ ρ0 l X r=0 l r ! |v (l−r) I |An+1,r n + 2 ˛ ˛ ˛ ˛ ˛ ≤ m|ρ0|mv γlRn 1 (n + 2)2(n + 1) × l X r=0 l r ! ((r + n)!)α ((l − r)!)α (n!)α ≤ m|ρ0|mv γlRn 1 (n + 2)2(n + 1) × ( l X r=0 l r ! (r + n)!(l − r)! n! )α = m|ρ0|mv γlRn ((l + n + 1)!)α ((n + 1)!)α 1 (n + 2)2(n + 1) ≤ 3 2 R|ρ0|mvAn+2,l. For n = 0, 1, the terms B3,2 and B3,3 are easily verified to be zero. Using the Lemmas 1.1, 1.2, and 1.5 successively, for B3,n+2 (n ≥ 2) we obtain: |B (l) 3,n+2| ≤ 2|ρ1| n−1 X k=1 1 (n+1)(n+2) l X r=0 l r ! Ak,rAn−k,l−r+1 ≤ 2m2 |ρ1| Rn−2γl+1 1 (n + 1)(n + 2) × n−1 X k=1 1 (k + 1)k(n − k)(n − k + 1) × ( l X r=0 l r ! (k + r − 1)!(n−k+l−r)! (k − 1)!(n − k − 1)! )α = 2m2 |ρ1| Rn−2γl+1 1 (n + 1)(n + 2) × n−1 X k=1 (n − k)α (k + 1)k(n − k)(n − k + 1) ( (l + n) n! )α ≤ 2m|ρ1|R3 γ n+3 n+1 n−1 X k=1 (n − k) k(k+1)(n−k+1) ! An+2,l. 9 Applying the Lemmas 1.3 (two times) and 1.6 |B (l) 3,n+2| ≤ 2m|ρ1|R3 γ ( (n + 3) (n + 2) n−1 X k=1 1 k(k + 1) ) An+2,l, = 2m|ρ1|R3 γ  (n + 3) (n + 2) „ 1 − 1 n «ff An+2,l ≤ 2 m|ρ1|R3 γ An+2,l. In B4,n+2 we have a product of three time dependent variables. Expanding this product, we obtain a triple sum, which can be handled by applying the Lemmas 1.1 and 1.2 twice. Note, that for n = 0, B4,2 is zero as well. For n ≥ 1 we have |B (l) 4,n+2| ≤ 2|ρ1| n X k=1 (n − k + 1) (n + 2)(n + 1) × l X r=0 l r ! An−k+1,l−r r X s=0 r s ! Ak,s ˛ ˛ ˛v (r−s) I ˛ ˛ ˛ ≤ 2 |ρ1|m2 mv γlRn−1 n X k=1 1 (n + 2)(n + 1)k(k + 1)(n−k+2) × ( l X r=0 l r ! (n−k+l−r)! (n − k)! r X s=0 r s ! (s+k−1)!(r−s)! (k − 1)! ! )α = 2 m2 mv|ρ1| Rn−1γl n X k=1 1 (n+2)(n+1)k(k + 1)(n − k + 2) ×  (n + l + 1)! (n + 1)! ffα = 2R2 mmv|ρ1| (n + 1) ( n X k=1 (n + 3) (k+1)k(n−k+2) ) An+2,l = 2R2 mmv|ρ1| (n + 1) ( n X k=1 1 k „ 1 (k+1) + 1 (n−k+2) « ) An+2,l = 2R2 mmv|ρ1| (n + 1)  n (n + 1) + n X k=1 1 (n + 2) „ 1 k + 1 k + 1 «ff An+2,l ≤ 2R2 mmv|ρ1|  n (n + 1)2 + 3 2 n (n + 1)(n + 2) ff An+2,l ≤ R2 mmv|ρ1|An+2,l. Finally, we calculate a bound for B5,n+2: |B (l) 5,n+2| ≤ |λ1| n X k=0 (n − k + 1) (n + 2) l X r=0 l r ! Ak+1,rAn−k+1,l−r ≤ m2 |λ1| Rn+2γl 1 (n + 2) × n X k=0 1 (k + 2)(k + 1)(n − k + 2) ( (n + l + 1)! (n + 1)! )α = Rm|λ1| ( n X k=0 (n + 3) (k + 1)(k + 2)(n − k + 2) ) An+2,l = Rm|λ1| ( n X k=0 1 (k + 1)(k + 2) + 1 (k + 2)(n − k + 2) ) An+2,l = Rm|λ1| ( n + 1 (n + 2) + 2 n + 4 n X k=0 1 k + 2 ) An+2,l ≤ 2Rm|λ1|An+2,l. Collecting |B1|, . . . , |B5| yields a (l) n+2 ≤ An+2,l |λ0|  4 3 R2 γ |ρ0| + 3 2 Rmv|ρ0| + 2 R3 m|ρ1| γ + 2R2 m|ρ1|mv + 2Rm|λ1|  . Hence, the bounds (1.5) hold, if R3 2m|ρ1| γ + R2  4 3 |ρ0| γ + 2m|ρ1|mv  + R  3 2 mv|ρ0| + 2m|λ1|  − |λ0| ≤ 0. The polynomial defined by the left hand side of the above inequality is negative for R = 0, strictly increasing for positive values of R, and positive for R large enough. Hence, it possesses exactly one positive root, which determines a lower bound on the radius of convergence. ACKNOWLEDGMENT The authors would like to thank Daniel Groß for early contributions to the results presented in this paper, and in particular for writing the simulation programs. REFERENCES [1] M. Fliess, J. Lévine, P. Martin, and P. Rouchon, “Flatness and defect of non-linear systems: introductory theory and examples,” Internat. J. Control, vol. 61, pp. 1327–1361, 1995. [2] ——, “A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems,” IEEE Trans. Automat. Control, vol. AC–44, pp. 922–937, 1999. [3] P. Martin, R. M. Murray, and P. Rouchon, “Flat systems, equivalence and feedback,” in Advances in the Control of Nonlinear Systems, ser. Lecture Notes in Control and Inform. Sci., A. Baños, F. Lamnabhi-Lagarrigue, and F. J. Montoya, Eds. Springer-Verlag, 2001, vol. 264, pp. 3–32. [4] ——, “Flat systems: open problems, infinite dimensional extension, symmetries and catalog,” in Advances in the Control of Nonlinear Systems, ser. Lecture Notes in Control and Inform. Sci., A. Baños, F. Lamnabhi-Lagarrigue, and F. J. Montoya, Eds. Springer-Verlag, 2001, vol. 264, pp. 33–57. [5] P. Rouchon, “Motion planning, equivalence, infinite dimensional sys- tems,” Int. J. Appl. Math. Comput. Sci., vol. 11, pp. 165–188, 2001. [6] J. Rudolph, Flatness based control of distributed parameter systems, ser. Berichte aus der Steuerungs- und Regelungstechnik. Aachen: Shaker Verlag, 2003. [7] M. Gevrey, “La nature analytique des solutions des équations aux dérivées partielles,” Ann. Sci. Ecole Norm. Sup., vol. 25, pp. 129–190, 1918. [8] C. D. Hill, “Parabolic equations in one space variable and the non- characteristic Cauchy problem,” Communications on Pure and Applied Mathematics, vol. XX, pp. 619–633, 1967. [9] M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Controlling the transient of a chemical reactor: a distributed parameter approach,” in Proc. Computational Engineering in Systems Application IMACS Multiconference, (CESA’98), Hammamet, Tunisia, 1998. [10] ——, “A distributed parameter approach to the control of a tubular reactor: a multi-variable case,” in Proc. 37th IEEE Conference on Decision and Control, Tampa, FL, 1998, pp. 439–442. [11] P. Rouchon and J. Rudolph, “Réacteurs chimiques différentiellement plats : planification et suivi de trajectoires,” in Automatique et procédés chimiques — Réacteurs et colonnes de distillation, J. P. Corriou, Ed. Hermès Science Publications, 2001, ch. 5, pp. 163–200. 10 [12] A. F. Lynch and J. Rudolph, “Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems,” Internat. J. Control, 2002. [13] B. Laroche, P. Martin, and P. Rouchon, “Motion planning for the heat equation,” Int. J. Robust and Nonlinear Control, vol. 10, pp. 629–643, 2000. [14] B. Laroche, Extension de la notion de platitude à des systèmes décrits par des équations aux dérivées partielles linéaires. Thèse de Doctorat, École Nationale Supérieure des Mines de Paris, 2000. [15] R. Rothfuß, U. Becker, and J. Rudolph, “Controlling a solenoid valve — a distributed parameter approach,” in Proc. 14th Int. Symp. Mathematical Theory of Networks and Systems — mtns 2000, Perpignan, France, 2000. [16] W. B. Dunbar, N. Petit, P. Rouchon, and P. Martin, “Motion planning for a nonlinear Stefan problem,” École Nationale Supérieure des Mines de Paris, Tech. Rep. [17] M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Systèmes linéaires sur les opérateurs de Mikusiński et commande d’une poutre flexible,” in ESAIM Proc., vol. 2, 1997, pp. 183–193, (http://www.emath.fr/proc). [18] Y. Aoustin, M. Fliess, H. Mounier, P. Rouchon, and J. Rudolph, “Theory and practice in the motion planning and control of a flexible robot arm using Mikusiński operators,” in Proc. 5th Symposium on Robot Control, Nantes, France, 1997, pp. 287–293. [19] W. Haas and J. Rudolph, “Steering the deflection of a piezoelectric ben- der,” in Proc. 5th European Control Conference, Karlsruhe, Germany, 1999. [20] J. Rudolph and F. Woittennek, “Flachheitsbasierte Randsteuerung von elastischen Balken mit Piezoaktuatoren,” at — Automatisierungstechnik, vol. 50, pp. 412–421, 2002. [21] H. Chen and L. Rodino, “General theory of pde and gevrey classes,” in General Theory of Partial Differential Equations and Microlocal Analysis, Q. Min-you and L. Rodino, Eds. Harlow, Essex: Addison Wesley Longman, 1996, pp. 6–81. [22] J. Rudolph, J. Winkler, and F. Woittennek, Flatness based control of distributed parameter systems: Examples and computer exercises from various technological domains, ser. Berichte aus der Steuerungs- und Regelungstechnik. Aachen: Shaker Verlag, 2003. Joachim Rudolph Biography text here. Jan Winkler Biography text here. Frank Woittennek Biography text here.