I D Gil
email: idgilc@unal.edu.co

J C Vargas

J P Corriou

Optimal nonlinear control of an industrial emulsion polymerization reactor

Keywords: Dynamic optimization, Nonlinear geometric control, Optimal temperature profile, Optimal feed policies, Polymerization, Vinyl acetate

. The case study corresponds to an industrial reactor operated at a chemical company in Bogotá (Colombia). An industrial scale reactor (11 m 3 of capacity) is simulated. Three different dynamic optimization problems are solved from the more simplistic (only one control variable: reactor temperature) to the more complex (three control variables: reactor temperature, initiator flow rate and monomer flow rate) in order to minimize the reaction time. The results show that it is possible to minimize the reaction time while some polymer desired qualities (conversion, molecular weight and solids content) satisfy defined constraints. The optimal temperature profile and optimal feed policies of the monomer and initiator, obtained in a dynamic optimization step, are used as optimal set points for reactor control. A nonlinear geometric controller based on input/output linearization is implemented for temperature control.

Introduction

The purpose of dynamic optimization studies is to determine a set of variables of a dynamic system, such as flow rates, temperatures, pressures, heat duties, . . . , that optimize a given cost function or criterion (costs, productivity, time, energy, selectivity) subject to specific constraints (dynamic model, operating conditions, safety and environmental restrictions). Some of the common problems of chemical engineering addressed by means of dynamic simulation and optimization include startup, upset, shutdown and transient analysis, safety studies, control and scheduling of batch and semi-batch processes, and the validation of control schemes [START_REF] Biegler | An overview of simultaneous strategies for dynamic optimization[END_REF][START_REF] Cervantes | Optimization strategies for dynamic systems[END_REF]. In all cases, it is important to possess a dynamic model sufficiently representative of the real process by means of mass and energy balances, and algebraic equations for physical and thermodynamic relations, but with a moderate complexity in order to get a mathematical and numerical solution without difficulty [START_REF] Corriou | Process Control -Theory and Applications[END_REF][START_REF] Corriou | Commande des Procédés, 3rd Edition[END_REF]. Some of the most important objectives in resins and polymer production plants are related to the improvement of safety, quality and productivity, minimum operating costs and respect of environmental constraints [START_REF] Gentric | Optimization and nonlinear control of a batch emulsion polymerization reactor[END_REF]. These make the optimization and control of polymerization reactors of great interest. In most cases, an optimization problem for a polymerization system requires the definition of an objective function and constraints which are defined by the reaction time and/or polymer molecular characteristics, together with operating conditions. In terms of polymerization reactors, the main contributions concern homogeneous reactions and some multiphase considerations trying to minimize the batch period, improve quality control and minimize the molecular weight distribution. In these cases, nonlinear models are essential to accurately describe the dynamics of the process. The solution of this kind of optimal control problems can be obtained by means of various optimization methods such as variational calculus, Hamilton-Jacobi equations, Pontryagin's maximum principle for continuous time systems and Bellman dynamic programming for discrete time systems, among others [START_REF] Corriou | Process Control -Theory and Applications[END_REF][START_REF] Corriou | Commande des Procédés, 3rd Edition[END_REF][START_REF] Kameswaran | Simultaneous dynamic optimization strategies: recent advances and challenges[END_REF][START_REF] Biegler | An overview of simultaneous strategies for dynamic optimization[END_REF].

In the case of emulsion polymerization, several studies deal with dynamic optimization. For example, [START_REF] Jang | Dynamic optimization of batch emulsion polymerization of vinyl acetate -an orthogonal polynomial initiator policy[END_REF] report the dynamic minimization of the final time of a batch emulsion polymerization of vinyl acetate using initiator flow rate as control variable, and the maximum allowable reac-tion rate together with the total amount of initiator as constraints. [START_REF] Gentric | Optimization and nonlinear control of a batch emulsion polymerization reactor[END_REF] calculate the optimal temperature profile that minimizes the batch time of a copolymerization reactor of styrene and α-methylstyrene using orthogonal collocation coupled with a sequential quadratic programming method. As constraints, they used the final conversion and the final number average molecular weight. [START_REF] Sayer | Dynamic optimization of semicontinuous emulsion copolymerization reactions: composition and molecular weight distribution[END_REF] and [START_REF] Vicente | Dynamic optimization of non-linear emulsion copolymerization systems. Open-loop control of composition and molecular weight distribution[END_REF] calculated the optimal monomer and chain-transfer agent feed profiles for the semi-batch methylmethacrylate (MMA)/n-butylacrylate (n-BA) emulsion copolymerization, using iterative dynamic programming with an objective function that included a term for the copolymer composition and also a term for the molecular weight distribution, in a way close to multiobjective optimization. [START_REF] Araújo | Optimization of semicontinuous emulsion polymerization reactions by IDP procedure with variable time intervals[END_REF] used variable time intervals with an iterative dynamic programming procedure to minimize the reaction time while composition and molecular weight are controlled at specific values. [START_REF] Paulen | Dynamic optimization of a hybrid system: emulsion polymerization reaction[END_REF] worked on the dynamic optimization of the emulsion copolymerization of styrene and α-methylstyrene applying control vector parameterization (CVP) method in order to minimize the total reaction time.

Recently, batch and semibatch operation of copolymerization of styrene and MMA [START_REF] Ibrahim | Optimisation of emulsion copolymerization of styrene and MMA in batch and semi-batch reactors[END_REF] were studied in order to maximize the monomer conversion in one case and the average molecular weight in a second case by means of CVP techniques solved by successive quadratic programming. Multiobjective optimization refers to simultaneous optimization of more than one objective function, which is typical in most real-life optimization problems encountered in industry [START_REF] Benyahia | Multicriteria dynamic optimization of an emulsion copolymerization reactor[END_REF]. Multiobjective dynamic optimization has been also studied for a semibatch styrene polymerization process in order to establish optimal operating temperature and feeding policies, which maximize monomer conversion and minimize the residual initiator in the final product [START_REF] Silva | Multi-objective dynamic optimization of semi-batch polymerization processes[END_REF].

In the present work, the dynamic optimization of the industrial emulsion polymerization of vinyl acetate is performed with respect to three different optimization scenarios. In the three cases, the objective is to minimize the reaction time by varying separately or simultaneously the reactor temperature, initiator flow rate and monomer flow rate. After optimization, some open loop optimal results obtained are used for comparison with closed loop simulations. The nonlinear geometric controller coupled with state estimation is used for tracking the optimal reaction temperature profile found by dynamic optimization. The results show the potential of dynamic optimization in finding optimal feed policies and operating temperature to improve the productivity. Also, it is demonstrated that the optimal temperature trajectory is well followed by means of a nonlinear controller.

Dynamic optimization of vinyl acetate emulsion polymerization

Experimental validation of the model

The dynamic state space model of the process has been described in a previous work [START_REF] Gil | Nonlinear geometric temperature control of a vinyl acetate emulsion polymerization reactor[END_REF]. Nine states are used [I, M t , M M , V pol , µ 0 , µ 1 , µ 2 , T, T j ] standing respectively for the initiator concentration I, total added monomer M t , remaining monomer M M , total volume of polymer V pol , three first moments of polymerization µ i , reactor temperature T and jacket temperature T j . The first seven equations correspond to the reaction kinetics and the two following equations to the energy balances for the reactor and jacket. Kinetic model:

dI dt = q I -k I I (1) dM t dt = q M (2) dM M dt = q M -R pol (3) dV p pol dt = R pol MW M ρ pol (4) dµ 0 dt = k f m [M] p + k f p µ 0 + k t λ 0 αλ 0 - k f p λ 0 µ 0 -(1 -α) 2 αλ 0 + 0.5k t λ 2 0 (5) dµ 1 dt = λ 0 1 -α ( k f m [M] p + k f p µ 0 + k t λ 0 α (2 -α) + k t λ 0 )- k f p λ 2 0 1 -α (1 -α) 2 (6) dµ 2 dt = λ 0 (1 -α) 2 (2α k f m [M] p + k f p µ 0 + k t µ 0 + k t λ 0 (2α + 1)) -2k f p λ 2 0 1 -α (1 -α) 3 1 -α + dµ 1 dt (7)
Reactor dynamics:

dT dt = q i C p,i (T i -T ) -∆H r R pol + UA(T j -T ) -Q cond m i C p,i (8)
dT j dt = F j T j in -T j m w - UA m w C p,w T j -T (9)
The major part of the model has been validated previously [START_REF] Araújo | Optimization of semicontinuous emulsion polymerization reactions by IDP procedure with variable time intervals[END_REF][START_REF] Arora | Model based operation of emulsion polymerization reactors with evaporative cooling: Application to vinyl acetate homopolymerization[END_REF] based on the experimental data obtained by [START_REF] Penlidis | Dynamic modeling of emulsion polymerization reactors[END_REF][START_REF] Penlidis | Polymer reactor design, optimization and control in latex production technology[END_REF]. The main modifications to the model are related to the energy balances for the reactor and the jacket. In this section, a simulation of a pilot plant reactor from a Colombian chemical company is reported, in order to reproduce its industrial operation, to compare the results with some measurements of the solids content and to have an idea of the reaction conversion. Some tests in a pilot scale reactor were carried out and the results were compared with a simulation of the same system using Matlab. A semi-batch emulsion polymerization reaction of vinyl acetate was performed. The used recipe is given in Table 1. The reactor temperature was fixed taking into account the preheating step. The monomer and initiator flow rates were adjusted manually by the operator during the whole operation. The goal is to maintain a nearly constant value of the temperature in the reactor. In the current operation of the industrial reactor, a jacket is used but only for the initial preheating step. During the reaction, temperature is maintained almost constant by manipulating only initiator and monomer flow rates. Consequently, the temperature cannot accurately follow a given set point and this limits the feed flow rate of the reactants in view of a more efficient process. For that reason, a different control is desirable with efficient use of the jacket. To that intent, energy balances of the reactor will be used to take into account the behavior of the jacket and reactor contents. This will be useful in the next stages of the work to propose a control strategy associated also to a dynamic optimization of the system. The solids content φ S is calculated theoretically by summing the weight of polymer formed and the weight of the polyvinyl alcohol and dividing it by the total weight of the latex

φ S = (M t -M M )MW M + M P V OH M t MW M + M P V OH + ρ w V w (10)
The viscosity of the reactor contents η is calculated from the expression (11) proposed by [START_REF] Chylla | Temperature control of semibatch polymerization reactors[END_REF]) that depends on the solids content and the reactor temperature. Some of the parameters of (11) were proposed by other authors [START_REF] Graichen | Feedforward control with online parameter estimation applied to the Chylla-Haase reactor benchmark[END_REF][START_REF] Hvala | Modelling, simulation and control of an industrial, semi-batch, emulsion-polymerization reactor[END_REF] or fitted using experimental information from this test

η = c 0 e (c 1 φ S ) 10 c 2 a 0 T -c 3 (11)
where a 0 , c 0 , c 1 , c 2 and c 3 are model parameters, and T is the reactor temperature. Figure 1 shows monomer and initiator flow rates used during all the reaction. As it was mentioned previously, these flow rates are imposed manually by the operator to guarantee a nearly constant temperature in the reactor.

In the same way, Figure 2a and Figure 2b show the corresponding quantities of initiator and monomer respectively remaining in the reactor. It can be seen that the quantity of remaining monomer in the reactor changes accordingly to variations of the monomer flow rate, and also it is influenced by additional injections of initiator that increase the polymerization rate. The initiator reacts in the first part of the batch rapidly as a result of the high increase of the temperature and then, when the initiator feed starts, there is an accumulation of initiator (Figure 2a). It is also evident that almost all the monomer is consumed in the polymerization reaction (Figure 2b) verifying the large conversions typical of emulsion polymerization processes.

Reactor temperature during the total batch time is shown in Figure 3.

The average temperature in the reactor is close to 343 K. Fron Figures 2 and3, it can be noted that, when reactor temperature increases, reaction rate also increases thus decreasing the number of moles of monomer remaining in the reactor. Average molecular weight, dispersity and conversion are presented in Figure 4. The average molecular weight and, in consequence, the dispersity are varying in function of the initiator and monomer injections during all the batch while the conversion increases up to 350 min (approx.), then decreases slowly, and finally increases until the end of the run where the conversion is higher than 98%. This is due to the effect of temperature on the reaction rate and also it can be explained by the constant initiator flow rate at the end of the batch. and viscosity were plotted with those data of experimental pilot test (Figure 5). As can be verified, there is a good representation of the solids content along the run. This is an indication of the adequacy of the model and its parameters, taking into account that the solids content is directly related to the reaction and the formation of polymer particles. On the other way, it should be noted also that, apart from a single point which is erroneous, viscosity is well represented by the model. The comparison made with the experimental data obtained shows that a good approximation with respect to typical values of viscosity of this kind of emulsions is achieved. In addition, viscosity measurement is spoilt by an experimental error that was not quantified here, but depends on the measurement technique and the operator. Of course, the experimental error has an effect on the fitting of the viscosity model to the experimental values determined in the pilot reactor.

Industrial case study

In the following, the dynamic optimization of an industrial emulsion polymerization reactor to produce poly-vinyl acetate will be presented. The case study corresponds to the industrial reactor operated in a chemical company in Colombia. An industrial scale reactor (11 m 3 of capacity) is simulated in the case where a semi-batch emulsion polymerization reaction of vinyl acetate is performed. A schematic of the reactor is shown in Figure 6. The used industrial recipe is shown in Table 2.

Three different dynamic optimization problems are solved with a piecewise constant control using different discretization scenarios. In the three problems, three different variables were considered as control variables u(t): solves constrained NLP problems. The mathematical model for the emulsion polymerization was reported in a previous work [START_REF] Gil | Nonlinear geometric temperature control of a vinyl acetate emulsion polymerization reactor[END_REF] and was summarized in section 2.1. In the dynamic optimization problem, the energy balances representing the reactor temperature dynamics are not considered because the reactor temperature T will be assumed as a control vari- for dynamic optimization then corresponds to the following seven states

[I, M t , M M , V pol , µ 0 , µ 1 , µ 2 ].
Here, the goal is to calculate an optimal reactor temperature profile based on the knowledge of the polymerization kinetics.

In the second part of this study, devoted to closed loop control, the dynamics of the reactor as energy balances will be considered.

Process operation

An emulsion polymerization process displays different behaviors according to the relative rates of initiation, propagation and termination, which at the same time depend on the monomer flow rate, initiator flow rate and reaction conditions. Typically, semi-batch emulsion polymerizations are divided in two steps: batch and fed-batch (Figure 7). At initial time t = 0, specific quantities of monomer, initiator, water and protective colloid, representing a fraction of the recipe, are charged to the reactor. In the process studied here, according to the procedure followed by the chemical company, monomer, ini- In this study, the pre-heating step is not taken into account for the dynamic optimization calculations. At the end of the pre-heating stage, when the reactor temperature reaches 340 K, the reaction is assumed to start and it corresponds to the initial reaction time t = 0 which is thus difficult to determine exactly. Later, the reactor temperature will take a value between 348 and 355 K, as it will be explained in next section.

Three way valve

Minimization of batch time with T as control variable

In the first case, in order to maximize the productivity of the industrial polymerization reactor, i.e. to minimize the final batch time, the optimal temperature profile is calculated. Temperature is chosen because of its large influence on the polymerization reaction and polymer properties. The opti-mization problem can be formulated as min T (t) The monomer and initiator flow rates are constant according to the recipe of Table 2. The value of monomer flow rate corresponds approximately to the ratio of the monomer quantity of the industrial recipe over the duration of the fed-batch time, similarly for the initiator flow rate. Indeed, as the final time is not known before dynamic optimization, the monomer is fed from t = 70 min with a flow rate of 2.1 mol/s, whereas the maximum authorized flow rate in the industrial operation is 2.3 mol/s. The initiator flow rate is 1×10 -3 mol/s also from t = 70 min (Figures 9a and9b). Then, the operation is performed in two stages, the first one in batch mode (without addition of reactants) and the rest of the operation in fed-batch mode where monomer and initiator are fed to the reactor. Consequently, two different models are used, one for the batch operation, followed by another one for the fed-batch.

t f t 0 dt = t f -t 0 s.t. ẋi = f i (x(t), T (t), t) , i = 1, ...,
Dynamic optimization is applied to the entire operation, i.e. including the batch mode and fed-batch mode stages (Figure 8). The influence of the number N u of discrete time segments used during the total reaction time has been studied. Thus, four different discretization scenarios were calculated considering N u equal to 3, 5, 10 and 20, the temperature taken as the control variable being piecewise constant on each time interval. Table 3 shows the influence of the number of discrete time increments on the most important quality indicators of the final polymer. As it can be observed, the total batch time decreases little with the increase of the number of piecewise controls. At the same time, the overall quality of the polymer such as the conversion and the average molecular weight is very little influenced by the increase of the number of time increments. The chain length decreases when temperature increases due to the transfer reactions, which induces the entry of radicals to the particles. This results finally in instantaneous termination reaction of these radicals inside the particles. It must be noted that the same strategy of temperature increase, at the end of the operation, is currently used by the operators of the reactor at the company. In the same way, Figures 11a, 11c and11e show that the rate of change of conversion in time is directly proportional to the temperature. Figures 10a and10c show that the temperature often switches from the lower bound to the upper bound and viceversa, in particular at the start and close to the end of the total reaction time. This behaviour is typical from minimum time problems and is known as bang-bang control [START_REF] Corriou | Process Control -Theory and Applications[END_REF][START_REF] Corriou | Commande des Procédés, 3rd Edition[END_REF][START_REF] Chachuat | Nonlinear and Dynamic Optimization: From Theory to Practice[END_REF]. Many systems in chemical engineering and in other domains are controlled in on-off mode, in a way similar to bang-bang control.

Minimization of batch time with reactor T and q I as control variables

In the second optimization problem, one additional optimization variable is considered. Here, the initiator flow rate is also used as a control variable due to its large effect on the monomer conversion and molecular properties of the final product. Again, the objective is to minimize the final reaction time, and in consequence the optimal temperature and initiator flow rate profiles are calculated. In this case, the optimization problem is formulated as min T (t),q I (t)

t f t 0 dt = t f -t 0 s.t. ẋi = f i (x(t), T (t), q I (t), t) , i = 1, ..., 7 and ∀t ∈ [t 0 , t f ] state model x 1 (t 0 ) = 5 , initiator moles
x 2 (t 0 ) = 4000 , total monomer moles x 3 (t 0 ) = 4000 , residual monomer moles x i (t 0 ) = 0 , i = 4, ..., 7 initial conditions x f ≥ 0.99 , final conversion Mn f ≥ 2.1 × 10 5 , final number average molecular weight φ S ≥ 49% , final solids content 348K ≤ T (t) ≤ 355K , temperature interval q I (t) ≤ 0.8 × 10 -3 mol/s , initiator flow rate interval piecewise discretization using 3, 5, 10 and 20 control segments was studied.

274

Table 4: Dynamic optimization case using T and q I as control variables: Optimization results Table 4 shows the results for the four optimization runs. For this second dynamic optimization problem, the constraint of the minimum conversion was increased from 95% to 99%. For this reason, the final times are slightly larger than those obtained previously. However, the conversion has been increased importantly which is beneficial for the process efficiency. The reaction times calculated here are approximately 9% lower than the typical reaction time in the plant, thus improving the productivity of the process.

N u t f (s) x f Mn,f × 10 -5 Mw,f × 10 -5 D
✲ ✻ ✻ ✻ t = 0 t = 70 t = t f inal ❄ Reaction starts ✲ ✛ Monomer flowrate q M = 0 ✲ ✛ Fed-batch mode Constant monomer flow rate ✲ ✛
Calculation of optimal T and q I profiles Figure 12: Dynamic optimization case using T and q I as control variables: Scheduling of operations Again, it can be noted that the total reaction time decreases with the increase of N u , but the final values of the average molecular weight and the dispersity of the polymer do not vary significantly.

Figures 13 and14 show the optimal temperature profile and optimal initiator flow rate profile, respectively. Again, the bang-bang control tendency is observed for the temperature and a similar behaviour is observed for the initiator when N u = 20 (Figure 14d). Also, the multivariable nature of the system is observed. The initiator flow rate and the reaction temperature appear to be correlated. In the first half of the operation, the tendency is reactor increase the polymerization rate and the conversion, and reduce the final molecular weight (Figure 15).

Minimization of batch time with T , q I and q M as control variables

The third optimization problem involves three control variables: temperature, initiator flow rate and monomer flow rate. These three variables are easily controlled in the industrial reactor and therefore are susceptible to be changed at the same time to achieve a desired performance of the reactor. Taking into account that vinyl acetate has a high rate of radical transfer to polymer, the monomer feed flow rate has an important effect on the molecular weight. For that reason, it is also considered in this optimization case as a control variable. The optimization problem is formulated as min T (t),q I (t),q M (t)

t f t 0 dt = t f -t 0 s.t.
ẋi = f i (x(t), T (t), q I (t), q M (t), t) i = 1, ..., 7 and ∀t ∈ [t 0 , t f ] st. model x 1 (t 0 ) = 5 , initiator moles x 2 (t 0 ) = 4000 , total monomer moles x 3 (t 0 ) = 4000 , residual monomer moles x i (t 0 ) = 0 , i = 4, ..., 7 initial conditions x f ≥ 0.992 , final conversion Mn f ≥ 2.2 × 10 5 , final number average molecular weight

φ S ≥ 50% , final solids content 348K ≤ T (t) ≤ 355K , temperature interval q I (t) ≤ 0.8 × 10 -3 mol/s , initiator flow rate interval
The monomer flow rate, initiator flow rate and temperature profiles were determined according to the dynamic optimization problem formulated in equation ( 14). In this case, fed-batch operation mode is used as the unique operation mode for the dynamic optimization calculations (Figure 16). As in the two last cases, a piecewise discretization using 3, 5, 10 and 20 control segments was studied. Table 5 shows the results for the four optimization runs. The final conversion constraint here is slightly higher than the conversion constraint used in the dynamic optimization based on temperature and initiator flow rate as control variables. This is mentioned because the final time obtained here is lower than the minimum time calculated in the two previous cases where the conversion constraint was lower as well as the constraints of molecular weight and solids content. The most interesting result of this optimization case is the final time obtained when using 10 and 20 piecewise controls.

✲ ✻ ✻ t = 0 t = t f inal ❄ Reaction starts

Monomer and initiator flowrates are variable

Fed-batch mode ✲ ✛ Calculation of optimal T , q I and q M profiles Figure 16: Dynamic optimization case using T , q I and q M as control variables: Scheduling of operations

In this case, the influence of the number of time intervals or piecewise controls used is more evident because more degrees of freedom exist that promote interactions during the reaction. This multivariable problem with three control variables allows the optimization solver to more easily find optimal operating values of the process variables which minimize the total reaction time. Specifically, in the case of 20 piecewise controls, a total time of 23762 seconds is obtained (Figure 17). This optimal final time is at least 20% lower than the current batch time used in the plant to perform this polymerization.

Optimal profiles are shown in Figures 17, 18 and 19 for temperature, initiator flow rate and monomer flow rate, respectively. The influence of during the major part of the reaction, and finally is decreased when the polymerization is ending (Figure 19). In the real plant operation, the initiator is fed almost continuously and this differs largely from the optimal policy obtained by dynamic optimization where some initiator is fed at the end of the reaction. Also, in the real plant, reactor temperature is maintained Figure 19: Dynamic optimization case using T , q I and q M as control variables: Optimal monomer flow rate profile approximately constant by manual operation during the reaction without using the jacket, but only by means of the sensible heat of the monomer and initiator fed to the reactor. In the present study, the proposition will be to use the jacket for heat exchange, but to feed the reactor in monomer and initiator according to the optimal policy and control the temperature by means of a nonlinear controller. In the same way as the initiator takes a maximum value just before the final time, it can also be noticed that the temperature increases to its upper bound during the final control segment also trying to reduce the residual monomer of the polymer (Figure 17). In this section, a nonlinear geometric temperature controller is applied to 354 the industrial emulsion polymerization reactor. Temperature is controlled by 355 a nonlinear controller designed in our previous work [START_REF] Gil | Nonlinear geometric temperature control of a vinyl acetate emulsion polymerization reactor[END_REF]. The optimal temperature profile calculated in the last section is used here as the set point of the temperature control loop. The optimal feed policies of initiator and monomer are also used here as set points supposing an automatic regulatory control for these two flow rates. Figure 21 shows a schematic representation of the control loops and the use of dynamic optimization results.

Conversion and polymer quality results obtained are compared with those corresponding to the current operation data. The controlled output is the temperature of the reactor contents T . The position u of a three-way valve is used as a manipulated input that imposes the respective flow rates through the cold and hot heat exchangers so that the inlet coolant temperature T j in could also be considered as the manipulated variable. Thus, the controlled system is reduced to Single Input Single Output.

Dynamic optimization (open loop) ✲ ✲ Automatic flow rate controller ✲ ✲ q ref M q ref I q M q I ✲ ✍✌ ✎☞ ✻ T ref ✲ Nonlinear geometric controller ✲ Emulsion Polymerization Reactor ✲ T u (control valve)
Two models are used, on one side a complete detailed model considered as the plant describing the dynamic reactor behavior and the polymerization reaction, including the moments of the polymer chains described in eqs (5-7). On the other side, in order to simplify the nonlinear geometric control law and the state estimation, a reduced model is built with the following reduced state vector for control and estimation

x = (I, M t , M M , V pol , T, T j ) ( 15
)
where the three moments of the polymer chains are not taken into account [START_REF] Gil | Nonlinear geometric temperature control of a vinyl acetate emulsion polymerization reactor[END_REF]. Among these states, I, M t and V pol are not observable and they are only predicted, i.e. they are obtained by simple integration of the differential equations without correction. The continuous-discrete extended Kalman filter is implemented to estimate the three states M M , T, T j which are used in the nonlinear control law. The controller and the observer share the same input and output. It must be noted that the state estimations or predictions which are provided are useful also for monitoring of the reactor.

The details of the extended Kalman filter as well as the nonlinear control law are presented in [START_REF] Gil | Nonlinear geometric temperature control of a vinyl acetate emulsion polymerization reactor[END_REF].

Control under optimal conditions

Now, it is important to show that a nonlinear geometric controller is capable of tracking the optimal temperature trajectories with the combined use of the optimal feed policies calculated for the initiator and the monomer. It 3.2. Optimal temperature control with optimal feed policies (q I and q M )

The most interesting result of the dynamic optimization study is the minimization of the reaction time (maximization of productivity). Now, it will be supposed that the reactor is operated under the operating conditions found in section 2.6 and following the recipe defined for the industrial operating conditions in Table 2. The case of the dynamic optimization using also the three control variables and N u equal to 20 is discussed here. The nonlinear geometric controller presents a good performance for tracking the temperature trajectory calculated by means of dynamic optimization (Figure 22a). At the end of the operation, the temperature increases up to the upper limit (355 K) and also the initiator flow rate is increased in order to end the reaction and satisfy the final constraints (Figure 22a).

The controller works well to remove the reaction heat (Figure 22e) and maintain the temperature close to the optimal set point previously established (Figure 22a). After preheating stage (the first 80 min) where the seed is formed, the monomer conversion increases up to 99.2% as in the dynamic optimization problem. It is interesting to note that, according to the dynamic optimization results, the temperature is maintained at the minimum tor flow rate also increases at the end of the reaction to satisfy the final constraints related with the final molecular weight. This demands an additional effort from the controller which is compensated, in principle, by the control valve (Figure 22b) but also it compensates the release of an important additional quantity of heat (Figure 22e). However, it is clear that the controller is capable of following the temperature setpoint. Consequently, the desired conversion and quality results are obtained. Another interesting observation from Figure 22f is that the average number of radicals per particle n increases importantly at the end of the operation, due to the final large injection of initiator calculated by the optimization. The low radical desorption rate with respect to radical entry rate produces a large increase of the heat released by the reaction and, simultaneously, the viscosity at the end of the operation increases up to 1200 P, limiting the mobility of radicals and favoring the increase of the temperature. All these observations could indicate a gel effect phenomenon which also accompanies final conversions as high as 99.2%. However, again, this is satisfactorily managed by the controller as it can be observed in the last 20 minutes in Figure 22a. Table 6 summarizes the constraints established in two of the dynamic optimization cases studied and the results obtained with the application of the nonlinear control to these optimal scenarios. It can be noted that, in all the cases, the constraints previously defined in dynamic optimization calculations are also satisfied when applying the nonlinear control to the simulated plant.

Robustness tests

The robustness of the nonlinear controller coupled with state estimation is studied with respect to typical modeling errors or real situations found in the normal plant operation. For that purpose, the rigorous model of the plant, representing the process, remains the same as it was defined in section 2.1 while the simplified model used in control calculations is modified by introducing modeling errors. Four different scenarios were supposed: errors about the propagation reaction constant, heat of reaction and heat transfer coefficient values, and a change in the temperature of the coolant used in the reactor cooling. In the first robustness test, it was assumed that the propagation constant is reduced by 50% with respect to its real value. Thus, the reaction rate calculated in the reduced model of the state estimator is lower than the real rate of the process, and therefore also the estimation of the reaction heat produced. The second test consisted in increasing by 50% the value of the heat of the reaction reducing the exothermicity of the reaction.

In the third robustness test, it was assumed that the heat transfer coefficient U is decreased by 40% with respect to its real value. This is representative of changes due to equipment fouling and variations in the viscosity and the solids content of the reacting mixture that also affect the value of U. Finally, in the last test, the coolant temperature was increased from 293 to 300 K simulating a failure in the cooling tower of the plant. Thus, the value of the inlet jacket temperature considered in the reduced model is erroneous and this influences the energy balances.

The main effect of these robustness tests is observed about the estimated monomer conversion (Figure 23). In the first 300 minutes, the difference between the estimated monomer conversion and the real value is slight for the cases where errors about heat transfer coefficient (Figure 23a) and heat of reaction (Figure 23b) were introduced. In spite of important differences with respect to the real values, the controller works well and the monomer conversion is well estimated. On the other way, the errors due to the change in the coolant temperature and propagation rate constant are more important with differences around 10% at maximum compared with the actual monomer conversion. Globally, the same tendency is maintained up to 300 minutes, but, in the last part of the operation (between 300 and 500 minutes), the deviation between the estimated and calculated values is large in three of the four cases. This deviation is caused by the addition of an important quantity of initiator at time t = 300 min. The initiator accelerates the reaction, an additional quantity of heat is released and the average number of particles also increases promoving a gel effect in the system. Also, it should be mentioned that, in the last part of the operation, the solids content and viscosity increase importantly modifying heat and mass transfer coefficients.

However, a good temperature control is achieved in all cases (Figure 24)

demonstrating the efficiency of nonlinear control and its robustness.

Conclusion

The dynamic optimization of the emulsion polymerization of vinyl acetate was studied. Three different optimization scenarios were established from the more simplistic (only one control variable) to the more complex (three control variables) in order to minimize the reaction time. Constraints are imposed with respect to some polymer desired qualities (conversion, molecular weight and solids content) as well as allowed flow rates and reactor temperature. Piecewise constant profiles were assumed and the influence of time discretization was studied. The influences of initiator, temperature and monomer were identified. In all the cases, the control variables often change during the batch according to the well-known bang-bang effect, typical of A nonlinear controller was used to track the temperature in the polymer-492 ization reactor in spite of typical disturbances such as initiator and monomer 493 injections. The optimal temperature profile, obtained by a dynamic opti-494 mization study, was used as the set point for the nonlinear control. In the 495 same time, the optimal feed policies of monomer and initiator were followed 496 by means of a regulatory control of their flow rates. The results show that 497 the nonlinear controller is appropriate to track the optimal temperature trajectories calculated previously. Also, the final temperature increase due to the initiator injection is rapidly corrected by the controller action making the operation of the reactor safer while, at the same time, the productivity is improved satisfactorily.

Acknowledgement

The financial support of the Embassy of France, Colfuturo and Universidad Nacional de Colombia is greatly acknowledged.

Nomenclature

C pj

Specific heat of component j [J.K -1 .kg Concentration of first moment for dead chains [-]

Figure 1 :

 1 Figure 1: Flow rates for the pilot reactor test. (a) Initiator, (b) Monomer

Figure 2 :

 2 Figure 2: Number of moles in the pilot reactor test. a) Initiator, b) Monomer

FinallyFigure 3 :

 3 Figure 3: Temperature profile for the pilot reactor test

Figure 4 :

 4 Figure 4: Quality results for the pilot reactor test. a) Number average molecular weight, b) Weight average molecular weight, c) Dispersity, d) Conversion

Figure 5 :

 5 Figure 5: Quality results of the experimental pilot reactor test. a) Solids content, b) Viscosity

Figure 6 :Figure 7 :

 67 Figure 6: Schematic of the industrial emulsion polymerization reactor

Figure 9 :

 9 Figure 8: Dynamic optimization case using T as control variable: Scheduling of operations

FiguresFigure 10 :

 10 Figures 10 and 11 show the influence of the number of time increments

Figure 11 :

 11 scheme of the operation is shown in Figure12. As in the first problem, a

Figure 13 :Figure 14 :

 1314 Figure13: Dynamic optimization case using T and q I as control variables: Optimal temperature profile

Figure 15 :

 15 Figure 15: Dynamic optimization case using T and q I as control variables: Quality results. Left column: monomer conversion. Right column: average molecular weight.

Figure 17 :Figure 18 :

 1718 Figure17: Dynamic optimization case using T , q I and q M as control variables: Optimal temperature profile

Figure 20 :

 20 Figure20: Dynamic optimization case using T , q I and q M as control variables: Quality results. Left column: monomer conversion. Right column: average molecular weight

Figure 21 :

 21 Figure 21: Schematic representation of the control of the emulsion polymerization reactor

  must be noted that the dynamic optimization was performed in open loop using only the kinetic model without the energy balances of the reactor. The nonlinear controller of course makes use of the same kinetic model but with the energy balances. Thus, apart from the study of temperature tracking, it will be interesting to compare the various characteristics of the polymer obtained in open loop and in closed loop, such as viscosity, conversion, number average molecular weight.

Figure 22 :

 22 Figure 22: Nonlinear geometric control of the emulsion polymerization reactor with time minimization using q I , q M and T as optimization variables and N u = 20. a) Optimal temperature profile b) Valve position; c) Dispersity d) Monomer conversion; e) Reaction power f) Average number of radicals per particle

Figure 23 :

 23 Figure 23: Influence of the robustness tests on the monomer conversion. a) Heat transfer coefficient error; b) Reaction heat error; c) Coolant temperature error; d) Propagation constant error

Figure 24 :

 24 Figure 24: Influence of the robustness tests on the reactor temperature control. a) Heat transfer coefficient error; b) Reaction heat error; c) Coolant temperature error; d) Propagation constant error

  491

Table 1 :

 1 Recipe used in the pilot emulsion polymerization reactor

	Component	Load (kg)
	Water	36
	Vinyl acetate	35
	Potassium persulfate	0.12
	Polyvinyl alcohol	3.5
	In order to follow the reaction, solids content and viscosity were measured
	by withdrawing samples at specific reaction times. The procedure established

Table 2 :

 2 Recipe used for the simulation of the industrial reactor

	Component	Load (kg)
	Water	5400
	Vinyl acetate	4651
	Potassium persulfate	12.8
	Polyvinyl alcohol	701

T T

  7 and ∀t ∈ [t 0 , t f ] , state model x 1 (t 0 ) = 5 , initiator moles x 2 (t 0 ) = 4000 , total monomer moles x 3 (t 0 ) = 4000 , residual monomer moles x i (t 0 ) = 0, i = 4, ..., 7 , initial conditions x f ≥ 0.95 , final conversion Mn f ≥ 1.8 × 10 5 , final number average molecular weight φ

S ≥ 46% , final solids content 348K ≤ T (t) ≤ 355K , temperature interval (12)

Table 3 :

 3

	x f	Mn,f × 10 -5	Mw,f × 10 -5 D
	3 26463 0.9625	2.6217	5.7279	2.18
	5 26458 0.9655	2.6341	5.7907	2.19
	10 26347 0.9678	2.6161	5.7509	2.19
	20 26214 0,9706	2.6586	5.8062	2.18

Dynamic optimization case using T as control variable: Influence of the number of time increments on the optimization results N u t f (s)

Table 5 :

 5 Dynamic optimization case using T , q I and q M as control variables: Optimization results

			N u t f (s)	x f		Mn,f × 10 -5	Mw,f × 10 -5 D
			3 29338 0.9920	2.4041	5.6376	2.34
			5 28426 0.9920	2.6186	6.0685	2.32
			10 26678 0.9920	2.1999	5.1265	2.33
			20 23762 0.9920	2.2000	5.1171	2.32
		360
	Temperature (K)	350 355
		0 345	100	200	300	400	500
				Time (min)

Table 6 :

 6 Results for the constraints established in the dynamic optimization. CDO: Constraint in Dynamic Optimization, CS: Control simulation

		Number average molecular weight Final conversion Solids content (%)
	N u	CDO	CS	CDO	CS	CDO	CS
	5 2.2 × 10 5	2.259 × 10 5	0.992	0.9921	50	53.1
	20 2.2 × 10 5	2.269 × 10 5	0.992	0.992	50	50

  -1 ] I Moles of initiator in the reactor [mol] k f m Rate constant for chain transfer to monomer [m 3 .mol -1 .s -1 ] k f p Rate constant for chain transfer to polymer [m 3 .mol -1 .s -1 ] k t Termination rate constant [m 3 .mol -1 .s -1 ] k I Overall initiation rate constant [s -1 ] M M Moles of monomer in the reactor [mol] Mn Number average molecular weight [g.mol -1 ] Mw Weight average molecular weight [g.mol -1 ] M t Total moles of monomer fed to the reactor [mol] MW M Monomer molecular weight [kg.mol -1 ] m w Mass of water in the reactor jacket [kg] n Average number of radicals per particle [-] q IFlow rate of initiator fed to the reactor [mol.s -1 ] q M Flow rate of monomer fed to the reactor [mol.s-1 ] R pol Overall reaction rate [mol.s -1 ] T Reactor temperature [K] T j Jacket temperature [K] U Overall heat transfer coefficient [W.m -2 .K -1 ] V pol Total volume of polymer generated in the reaction [m 3 ] V w Total volume of aqueous phase [m 3 ] α Probability of propagation [-] ∆H r Heat of reaction [J.kg -1 ] λ 0Total concentration of zeroth moment for growing chains [-] µ 0Concentration of zeroth moment for dead chains [-] µ 1


Ahmed Maidi

Jean-Pierre Corriou
email: jean-pierre.corriou@univ-lorraine.fr

Boundary Geometric Control of a Nonlinear Diffusion System with Time-Dependent Spatial Domain

Keywords: distributed parameter system, Cole-Hopf tangent transformation, Stefan problem, boundary control

Boundary Geometric Control of a Nonlinear Diffusion

System with Time-Dependent Spatial Domain

Introduction

Physical systems involving phase change phenomenon are characterized by a moving interface that separates the solid and liquid phases [START_REF] Crank | Free and moving boundary problems[END_REF]. To cite a few, melting, solidification, casting, chemical reaction, heat and mass transfer belong to this class of systems, which is commonly termed Stefan problems in the literature [START_REF] Crank | Free and moving boundary problems[END_REF][START_REF] Gupta | The Classical Stefan Problem[END_REF]. The characteristic variables of this kind of systems depend both on space and time variables. These systems are classified as distributed parameter systems (DPSs), which are of infinite dimension and characterized by a time-dependent spatial domain [START_REF] Dunbar | Motion planning for a nonlinear Stefan problem[END_REF][START_REF] Maidi | Boundary control of linear Stefan problem[END_REF].

The dynamic behavior modeling of Stefan problems leads to a hybrid model [START_REF] Crank | Free and moving boundary problems[END_REF][START_REF] Mattheij | Partial Differential Equations[END_REF] involving coupled partial and ordinary differential equations together with initial and boundary conditions. The partial differential equations (PDEs) describe the system dynamic behavior while the movement of the solid-liquid interface, that forms actually a moving boundary of the system, is modeled by an ordinary differential equation (ODE).

The control of DPSs represents a very challenging field and occupies an important place in control theory [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF]. Note that most contributions are based on the early lumping approach [START_REF] Ray | Advanced Process Control[END_REF][START_REF] Christofides | Nonlinear and robust control of PDE systems: methods and applications to transport-reaction processes[END_REF], i.e. the controller design is achieved using an equivalent ODEs model, denoted as a reduced model, that approximates the original PDEs one. This reduced model is derived using some powerful reduction techniques that achieve an accurate approximation characterizing the dominant dynamic behavior of the distributed parameter system (DPS) described by the PDEs model [START_REF] Li | Modeling of distributed parameter systems for applications-A synthesized review from time-space separation[END_REF]. This reduced model is obtained by approximating either the PDEs (for instance by discretization using either finite differences or finite elements methods) or their solutions (using for example Galerkin's method). This design procedure, denoted as early lumping approach, leads to a high dimension controller with a limited performance due to the approximation process that masks the spatially distributed nature of the DPS [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Christofides | Control of nonlinear distributed process systems: Recent developments and challenges[END_REF]. Nevertheless, during these last years, significant efforts have been made to develop sophisticated control strategies that enhance significantly the closed-loop performance following the early lumping approach [START_REF] Pourkargar | Modification to adaptive model reduction for regulation of distributed parameter systems with fast transients[END_REF][START_REF] Pourkargar | Geometric output tracking of nonlinear distributed parameter systems via adaptive model reduction[END_REF].

An interesting alternative to the early lumping approach is the late lumping approach that consists in designing the controller using directly the PDEs model without any approximation [START_REF] Ray | Advanced Process Control[END_REF] The late lumping approach preserves the spatially distributed nature and the fundamental control theoretical properties of the DPS, which allows to increase the controller performance [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Christofides | Nonlinear and robust control of PDE systems: methods and applications to transport-reaction processes[END_REF]. Several successful applications of the late lumping are reported in the literature [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Maidi | Distributed control of nonlinear diffusion systems by input-output linearization[END_REF][START_REF] Shang | Feedback control of hyperbolic distributed parameter systems[END_REF][START_REF] Wu | Output regulation of nonisothermal plug-flow reactors with inlet perturbations[END_REF][START_REF] Maidi | Boundary control of linear Stefan problem[END_REF]. Nevertheless, most developed control strategies and applications assume a fixed spatial domain, hence control of DPSs with a time-dependent domain constitutes a very challenging field. In this area, few contributions are reported in the literature [START_REF] Armaou | Nonlinear feedback control of parabolic partial differential equation systems with time-dependent spatial domains[END_REF][START_REF] Armaou | Robsut control of parabolic PDE systems with time-dependent spatial domains[END_REF][START_REF] Dunbar | Motion planning for a nonlinear Stefan problem[END_REF][START_REF] Ng | Optimal boundary control of a diffusion-convectionreaction PDE model with time-dependent spatial domain: Czochralski crystal growth process[END_REF][START_REF] Wang | Stabilization and control of distributed parameter system with time-dependent spatial domains[END_REF][START_REF] Wang | Feedback control of a heat diffusion system with timedependent spatial domains[END_REF] and most of them are based on the early lumping approach [START_REF] Armaou | Nonlinear feedback control of parabolic partial differential equation systems with time-dependent spatial domains[END_REF][START_REF] Armaou | Robsut control of parabolic PDE systems with time-dependent spatial domains[END_REF][START_REF] Ng | Optimal boundary control of a diffusion-convectionreaction PDE model with time-dependent spatial domain: Czochralski crystal growth process[END_REF].

The Stefan problem has been investigated in literature and, in most contributions, it is formulated as an inverse geometry problem [START_REF] Dunbar | Motion planning for a nonlinear Stefan problem[END_REF][START_REF] Khosravifard | Simultaneous control of solidus and liquidus lines in alloy solidification[END_REF][START_REF] Nowak | Application of bezier surfaces to the 3-d inverse geometry problem in continuous casting[END_REF] that is solved in open loop. The approach developed here is based on closed-loop control and thus differs by a large extent. To the best of our knowledge, the first application of the late lumping approach to the DPSs with timedependent spatial domain has been done by the authors of the present paper in the case of a linear Stefan problem [START_REF] Maidi | Boundary control of linear Stefan problem[END_REF]. In their contribution, the authors tackle the closed-loop control of a linear Stefan problem, that is, a linear diffusion process with a time-dependent spatial domain. Thus, a boundary geometric control law that enforces the solid-liquid interface position tracking was developed using the concepts of characteristic index [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Christofides | Nonlinear and robust control of PDE systems: methods and applications to transport-reaction processes[END_REF] and the closed-loop stability was demonstrated using the Lyapunov stability test. The performance of the developed control law was evaluated by simulation in the case of a zinc melting process.

In recent years, geometric control emerged as an interesting and suitable approach for designing controllers for DPSs using the late lumping approach [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Maidi | Boundary geometric control of a countercurrent heat exchanger[END_REF][START_REF] Maidi | Distributed feedback design for systems governed by the wave equation[END_REF][START_REF] Maidi | Distributed control of nonlinear diffusion systems by input-output linearization[END_REF][START_REF] Maidi | Boundary control of linear Stefan problem[END_REF]. Geometric control presents the following advantages:

• the controller can be easily designed using the PDE model by performing only successive time derivatives of the controlled variable [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Maidi | Distributed control of nonlinear diffusion systems by input-output linearization[END_REF][START_REF] Maidi | Boundary control of linear Stefan problem[END_REF],

• it allows the use of the full potential of an existing control theory for lumped linear systems that offers powerful controller design techniques to enhance the performance of the obtained linearized system [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Maidi | Boundary geometric control of a countercurrent heat exchanger[END_REF][START_REF] Maidi | Distributed feedback design for systems governed by the wave equation[END_REF],

• in the case of linear DPSs, it is easy to show, using some concepts from semigroup and Lyapunov theories, the stability of closed-loop system [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF][START_REF] Maidi | Distributed control of nonlinear diffusion systems by input-output linearization[END_REF][START_REF] Maidi | Boundary control of linear Stefan problem[END_REF].

The advantages cited above and the results developed by the authors of the present paper in the case of linear Stefan problem [START_REF] Maidi | Boundary control of linear Stefan problem[END_REF] motivate the extension of the geometric control to a nonlinear diffusion process with timedependent domain. Thus, based on their original results, the authors address the control of the solid-liquid interface position in the case of a nonlinear Stefan problem. It is shown that by using the Cole-Hopf tangent transformation [START_REF] Meleshko | Methods for Constructing Exact Solutions of Partial Differential Equations[END_REF][START_REF] Vadasz | Analytical solution to nonlinear thermal diffusion: Kirchhoff versus Cole-Hopf transformations[END_REF], an equivalent PDE linear model of the nonlinear Stefan problem can be easily derived. The resulting control law, that enforces position tracking of the liquid-solid interface, of the assumed nonlinear Stefan problem can be derived directly from the control law of the linear Stefan problem. Then, based on the fact that the Cole-Hopf tangent transformation is bijective, it is demonstrated that the obtained controller yields a stable closed-loop system. The performance of the developed controller is shown by numerical simulation in the case of melting of stainless steel characterized by a temperature-dependent thermal conductivity, which makes it a nonlinear diffusion process. The objective is to control the position of the liquid-solid interface by manipulating a heat flux applied at the boundary.

The paper is structured as follows: Section 2 is devoted both to the general Stefan problem and to the control problem formulation. The mathematical model of the Stefan problem is presented in Section 3. Section 4 gives the geometric control law in the case of a linear Stefan problem. The boundary geometric of the nonlinear Stefan problem is investigated in Section 5. The performance of the controller is evaluated, by numerical simulation, in Section 6 in the case of stainless steel melting. Section 7 concludes the paper.

Note that the present study addresses the control of a nonlinear diffusion with a time-dependent space domain, thus without loss of generality, the terminology from heat conduction transfer will be used in this paper to simplify the presentation.

Stefan problem and control problem formulation

A schematic representation of a melting process of a finite strip of a solid material of length L is given in Figure 1. The heat flux q(t) applied at the left boundary z = 0 causes the melting process that moves in the positive z direction. The system is thus composed of two phases, liquid and solid, separated by a moving liquid-solid interface. At the position z = l(t) corresponding to the liquid-solid interface, the temperature of the liquid is equal to the melting temperature T f of the considered material.

In the Stefan problem, the aim is to move the liquid-solid interface l(t) to a desired position l d (t) by manipulating the heat flux q(t) at the boundary z = 0. Thus, the control problem consists in designing a control profile q(t) that ensures this objective in a finite time. The design of this control law needs to set the formulated control problem under a mathematical model form which is given in the following Section.

Liquid phase

T (z, t) Solid phase Liquid-Solid interface q(t) z 0 z = l(t) L

Mathematical model of Stefan problem

Assuming that the heat conductivity of the liquid phase is temperaturedependent, the evolution of the liquid phase temperature is modeled by the following nonlinear diffusion equation [START_REF] Mattheij | Partial Differential Equations[END_REF]

ρ c p ∂T (z, t) ∂t = ∂ ∂z k T (z, t) ∂T (z, t) ∂z (1)
with the following fixed boundary conditions

-k T (0, t) ∂T (z, t) ∂z z=0 = q(t) (2)
T l(t), t = T f (3)
and the initial condition

T (z, 0) = T 0 (z) (4)
with 0 ≤ z ≤ l(t) ≤ L. The PDE ( 1) is quasilinear, that is a kind of nonlinear PDEs [START_REF] Salsa | A Primer on PDEs. Models, Methods, Simulations[END_REF].

The motion of the liquid-solid interface is described by the following ODE [START_REF] Crank | Free and moving boundary problems[END_REF][START_REF] Mattheij | Partial Differential Equations[END_REF]

ρ ∆H f dl(t) dt = -k T l(t), t ∂T (z, t) ∂z z=l(t) (5)
In this model, T (z, t) is the temperature of liquid phase, z the space variable, t ∈ [0, ∞[ the time variable, l(t) the position of the moving boundary (liquid-solid interface). The spatial domain at time t is denoted by Ω t = [0, l(t)] ⊂ Ω with Ω = [0, L] the largest spatial domain, i.e., 0 ≤ l(t) ≤ L. q(t) denotes the manipulated heat flux at boundary z = 0. The parameters ρ, c p , k, ∆H f and T f are the liquid density, liquid specific heat, liquid heat conductivity, latent heat of fusion and melting temperature of the material, respectively.

Note that, in the case where both c p and k are constant, i.e. independent of the temperature T (z, t), the nonlinear PDE (1) is reduced to a linear one. In this case, the control problem formulated in Section 2 has been solved by Maidi [4] in the framework of geometric control. The main results of their work used in the present study are summarized in the following section.

Geometric control of linear Stefan problem

Assuming that k is temperature-independent, i.e. constant, the nonlinear model ( 1)-( 5) is reduced to the following linear model

∂T (z, t) ∂t = α ∂ 2 T (z, t) ∂z 2 , α = k ρ c p (6) -k ∂T (z, t) ∂z z=0 = q(t) (7)
T (l(t), t) = T f (8) T (z, 0) = T 0 (z) (9) ρ ∆H f dl(t) dt = -k ∂T (z, t) ∂z z=l(t) (10)
which corresponds to the linear Stefan problem. In this case, the control of the position of liquid-interface by manipulating the heat flux q(t) has been solved by the authors of the present paper [START_REF] Maidi | Boundary control of linear Stefan problem[END_REF]. To give more generality to their results, a dimensionless model of ( 6)-( 10) has been used in both control design and closed-loop stability analysis. The dimensionless variables used are

w(ξ, τ ) = T (z, t) -T min T max -T min , τ = α t L 2 , ξ = z L (11)
which yields the following dimensionless model

∂w(ξ, τ ) ∂τ = ∂ 2 w(ξ, τ ) ∂ξ 2 , 0 ≤ ξ ≤ l a (τ ) ≤ 1 (12) - ∂w(ξ, τ ) ∂ξ ξ=0 = qa (τ ), qa (τ ) = L q(t) k (T max -T min ) (13) w(l a (τ ), τ ) = w f (14) w(ξ, 0) = w 0 (ξ) (15) ∆H f c p dl a (τ ) dτ = -(T max -T min ) ∂w(ξ, τ ) dξ ξ=la(τ ) (16)
where the dimensionless variables w f , w 0 , qa and l a are the solid melting temperature, initial solid strip temperature, heat flux and liquid-solid interface position, respectively. T min and T max are arbitrary temperatures adequately chosen for dimensionless reasons.

The control law developed by Maidi [START_REF] Maidi | Boundary control of linear Stefan problem[END_REF], in the framework of geometric control, using the concept of the characteristic index introduced by Christofides [START_REF] Christofides | Feedback control of hyperbolic PDE systems[END_REF], is given as follows

qa (τ ) = 1 γ w d m (τ ) -w m (τ ) -1 - c p w f (T max -T min ) ∆H f ∂w(ξ, τ ) ∂ξ ξ=la(τ ) (17)
where w m (τ ) is an auxiliary output defined as the dimensionless spatial average temperature of the liquid phase, along the spatial domain Ω τ , defined as follows

w m (τ ) = Ωτ w(ξ, τ ) dξ (18)
which is introduced to ensure the controllability condition, i.e. to have a finite characteristic index σ of the input-output system with input q(τ ) and the output w m (τ ) since, in the case of the output l a (τ ), the characteristic index σ does not exist (σ → ∞). γ and w d m (τ ) are the dimensionless closedloop time constant and desired set point of the controlled output w m (τ ), respectively.

The developed control (17) yields the following closed-loop system

γ dw m (τ ) dτ + w m (τ ) = w d m (τ ) ( 19
)
and it is shown that, using the Lyapunov stability test, the internal dynamics is stable. Note that the developed control law [START_REF] Armaou | Robsut control of parabolic PDE systems with time-dependent spatial domains[END_REF] ensures the output tracking of the auxiliary variable w m (τ ) but our objective is to control the position of the liquid-solid interface l a (τ ). Hence, to solve the original control problem, a simple control strategy has been proposed based on the steady state relation between the spatial average temperature of the liquid phase w m (τ ) and the position of the liquid-solid interface l a (τ ).

In this section, the authors summarize only the necessary results that will be used later in the present study. Geometric control of a linear Stefan problem has been studied extensively by Maidi [START_REF] Maidi | Boundary control of linear Stefan problem[END_REF], thus more details about the derivation of the control law [START_REF] Armaou | Robsut control of parabolic PDE systems with time-dependent spatial domains[END_REF], the stability analysis of the resulting closed-loop system and the general control strategy proposed can be found in this original work.

By assuming the dimensionless variables [START_REF] Pourkargar | Modification to adaptive model reduction for regulation of distributed parameter systems with fast transients[END_REF], the dimensional control law q(t) derived from the dimensionless control law [START_REF] Armaou | Robsut control of parabolic PDE systems with time-dependent spatial domains[END_REF] is

q(t) = ρ c p µ T d m (t) -T m (t) - µ α 1 - c p T f ∆H f ∂T (z, t) ∂z z=l(t) ( 20
)
where µ is the dimensional closed-loop time constant and T m (t) is the dimensional spatial average temperature of the liquid phase T (z, t) defined as

T m (t) = CT (z, t) = Ωt T (z, t) dz (21)
and T d m (t) is its desired set point. C is a bounded linear operator. The results presented in this section are exploited to derive the geometric control law of the nonlinear Stefan problem (1)-( 5), which is discussed in the following section.

Geometric control of nonlinear Stefan problem

If the control theory of the linear DPSs has reached a certain level of maturity thanks to the well-established and powerful semi-group theory [START_REF] Pazy | Semigroup of Linear Operators and Applications to Partial Differential Equations[END_REF], the control theory of nonlinear DPSs is a challenging field with many interesting open questions [START_REF] Christofides | Control of nonlinear distributed process systems: Recent developments and challenges[END_REF]. Concerning the control problem of nonlinear DPSs, most contributions come from the early lumping approach [START_REF] Baker | Finite-dimensional approximation and control of non-linear parabolic PDE systems[END_REF][START_REF] Dubljevic | Distributed nonlinear control of diffusion-reaction processes[END_REF]. Recently, the late lumping approach has been applied, in the framework of geometric control, to control nonlinear DPSs and particularly nonlinear diffusion systems [START_REF] Maidi | Distributed control of nonlinear diffusion systems by input-output linearization[END_REF].

Tangent transformations play a key role in studying and solving nonlinear PDEs. They allow to convert nonlinear PDEs to linear ones [START_REF] Meleshko | Methods for Constructing Exact Solutions of Partial Differential Equations[END_REF]. For instance, under a certain assumption, the nonlinear diffusion equation ( 1) can be easily converted to a linear diffusion equation using the Cole-Hopf transformation [START_REF] Vadasz | Analytical solution to nonlinear thermal diffusion: Kirchhoff versus Cole-Hopf transformations[END_REF]. This transformation has been exploited in the case of a nonlinear diffusion process, with a distributed control, to design a controller following the late lumping approach [START_REF] Maidi | Distributed control of nonlinear diffusion systems by input-output linearization[END_REF] and to study the controllability propriety of a nonlinear diffusion system [START_REF] Maidi | Controllability of a nonlinear diffusion system[END_REF]. In the following, the Cole-Hopf transformation is used to deduce the control of the nonlinear Stefan problem (1)-( 5) from the control law [START_REF] Wang | Feedback control of a heat diffusion system with timedependent spatial domains[END_REF] obtained in the case of the linear Stefan problem.

Cole-Hopf Transformation

The Cole-Hopf transformation consists in converting the nonlinear diffusion equation ( 1) into a linear one if the thermal diffusivity defined by the following ratio

α = k (T (z, t)) ρ c p (22)
is approximately constant. This hypothesis is accepted in many cases since the variation of α with T (z, t) is much less important than that of k (T (z, t)), so that this approximation is reasonable [START_REF] Carslaw | Conduction of Heat in Solids[END_REF].

To linearize the nonlinear diffusion equation ( 1), using the Cole-Hopf technique, one seeks a transformation of the form [START_REF] Vadasz | Analytical solution to nonlinear thermal diffusion: Kirchhoff versus Cole-Hopf transformations[END_REF] T (z, t) = ϕ(θ(z, t)) [START_REF] Maidi | Boundary geometric control of a countercurrent heat exchanger[END_REF] where ϕ( . ) is a nonlinear function.

From a physical point of view, the liquid phase temperature T (z, t) cannot undergo discontinuities in a single-phase medium and, to each temperature T (z, t), corresponds one temperature θ(z, t). Consequently, the nonlinear function ϕ( . ) is a bijective mapping and both temperatures T (z, t) and θ(z, t) are continuous functions with respect to z, that is, they are spatially differentiable.

By using the transformation [START_REF] Maidi | Boundary geometric control of a countercurrent heat exchanger[END_REF], the evaluation of the derivatives of the left-hand side of (1) gives

∂T (z, t) ∂t = dϕ(θ(z, t)) dθ(z, t) ∂θ(z, t) ∂t (24)
and for the right-hand side of (1), one obtains

∂ ∂z k (T (z, t)) ∂T (z, t) ∂z = M ∂θ(z, t) ∂z 2 + k(ϕ(θ(z, t))) dϕ(θ(z, t)) dθ(z, t) ∂ 2 θ(z, t) ∂z 2 (25)
with

M =k(ϕ(θ(z, t))) d 2 ϕ(θ(z, t)) dθ 2 (z, t) + dk(ϕ(θ(z, t))) dϕ(θ(z, t)) dϕ(θ(z, t)) dθ(z, t) 2 (26)
To make the right-hand side of (25) linear, the term M is set equal to zero. This term is a differential equation, which can be expressed in the following integrable form

d dθ(z, t) k(ϕ(θ(z, t))) dϕ(θ(z, t)) dθ(z, t) = 0 (27)
Integrating [START_REF] Salsa | A Primer on PDEs. Models, Methods, Simulations[END_REF] gives

k(ϕ(θ(z, t))) dϕ(θ(z, t)) dθ(z, t) = c 1 (28) hence k(ϕ(θ(z, t))) dϕ(θ(z, t)) = c 1 dθ(z, t) = c 1 θ(z, t) + c 2 (29)
or equivalently θ(z, t) =ϕ -1 (T (z, t))

= 1 c 1 k(T (z, t)) dT (z, t) - c 2 c 1 (30)
where ϕ -1 ( . ) is the inverse function of ϕ( . ) and c 1 = 0. According to the discussion in Subsection 5.1, as the mapping ϕ( . ) is a continuous bijective function, it follows that its inverse ϕ -1 ( . ) exists and is also a continuous function. The two arbitrary constants c 1 and c 2 can be chosen as follows [START_REF] Vadasz | Analytical solution to nonlinear thermal diffusion: Kirchhoff versus Cole-Hopf transformations[END_REF]

c 1 = k(0) ; c 2 = - k(T ) dT T =0 (31)
Thus c 2 is the primitive integral of the heat conductivity k(T ) taken at T = 0. This leads, as demonstrated by Vadasz [START_REF] Vadasz | Analytical solution to nonlinear thermal diffusion: Kirchhoff versus Cole-Hopf transformations[END_REF], after some elementary mathematical manipulations, to the following well-known Kirchhoff's transformation

θ(z, t) = ϕ -1 (T (z, t)) = 1 k(0) T (z, t) 0 k(ζ) dζ (32)
which is actually a particular case of Cole-Hopf transformation. In summary, using the transformation [START_REF] Maidi | Boundary geometric control of a countercurrent heat exchanger[END_REF], with ϕ(θ(z, t)) satisfying ( 27) and considering the expression [START_REF] Nowak | Application of bezier surfaces to the 3-d inverse geometry problem in continuous casting[END_REF] of diffusivity α, the nonlinear state equation ( 1) will be converted to the following linear one

∂θ(z, t) ∂t = α ∂ 2 θ(z, t) ∂z 2 (33)
and by taking into account both equations ( 28) and ( 31), the following new boundary and initial conditions can be easily deduced

-k(0) ∂θ(z, t) ∂z z=0 = q(t) (34) θ(l(t), t) = ϕ -1 (T f ) = θ f (35) θ(z, 0) = ϕ -1 (T 0 (z)) = θ 0 (z) (36)
and the equation of the interface motion ( 5) takes the following form

ρ ∆H f dl(t) dt = -k(0) ∂θ(z, t) ∂z z=l(t) (37)
Equations ( 33)-(37) represent a linear Stefan problem obtained by linearizing the nonlinear Stefan model ( 1)-( 5) by means of Cole-Hopf transformation. Therefore, the control law of the nonlinear Stefan problem (1)-( 5) can be deduced based on its equivalent linear model ( 33)-(37). This is explained in the following subsection.

Control law and closed-loop stability analysis

Based on the assumption made with respect to the diffusivity α, given by [START_REF] Nowak | Application of bezier surfaces to the 3-d inverse geometry problem in continuous casting[END_REF], being nearly constant compared to the variation of heat conductivity k, the equivalent linear model ( 33)-(37) provided by the Cole-Hopf transformation can be used to design the control law of the nonlinear Stefan problem (1)- [START_REF] Mattheij | Partial Differential Equations[END_REF]. Hence, according to the development given in Section 4, the control law that enforces the output tracking for the linear equivalent model ( 33)-( 37) is

q(t) = ρ c p µ θ d m (t) -θ m (t) -µ α 1 - c p θ f ∆H f ∂θ(z, t) ∂z z=l(t) (38)
where the measured output is

θ m (t) = Ωt θ(z, t) dz ( 39
)
and θ d m (t) is its corresponding set point. Thereafter, considering relations ( 30), ( 35) and (39), the control law (38) can be expressed according to the temperature T (z, t) as follows

q(t) = ρ c p µ T d m (t) - l(t) 0 ϕ -1 T (z, t) dz -µ α 1 - c p T f ∆H f ∂ ϕ -1 T (z, t) ∂z z=l(t) (40)
which leads to the following closed-loop system

µ dT m (t) dt + T m (t) = T d m (t) (41)
which is externally or BIBO-stable. Nevertheless, without the stability of the unobservable part, termed internal dynamics, the closed-loop stability of (41) is not guaranteed. In geometric control, the behavior of the internal dynamics, i.e. of the state of the system, plays a central role in proving the stability.

In the case of linear Stefan problem [START_REF] Maidi | Boundary control of linear Stefan problem[END_REF], based on the fact that the latent heat overcomes the sensible heat by a large extent for most physical systems, the authors demonstrated the stability of the internal dynamics of the closed-loop system using Lyapunov stability test. They concluded that the closed-loop system is exponential stable and θ(z, t) → θ f when t → ∞. Consequently, the exponential stability of the closed-loop system (41), that is the stability of its internal dynamics, can be deduced from the stability of its equivalent linear model ( 33)-(37) using the fact that Cole-Hopf transformation ( 23) is bijective. Lyapunov stability test for the equivalent linear model ( 33)-(37) yields that θ(z, t) → θ f when t → ∞, therefore using the Cole-Hopf transformation [START_REF] Maidi | Boundary geometric control of a countercurrent heat exchanger[END_REF], it can be concluded that T (z, T ) → ϕ(θ f ) when t → ∞, i.e. according to (35), it follows that lim t→∞ T (z, t) = T f . This means that the closed loop system (41) is internally stable.

l d (t) External Controller T d m (t) State feedback (41) q(t) Nonlinear Stefan process T (z, t) C T m (t) f( . ) l(t) + - + -

Global control strategy

Recall that our objective consists in solving the initial boundary control problem, i.e., controlling the position of the liquid-solid interface l(t) by manipulating the heat flux q(t) at the boundary z = 0. Thus, it is proposed to keep the control law (40) derived for the auxiliary output [START_REF] Khosravifard | Simultaneous control of solidus and liquidus lines in alloy solidification[END_REF] and to define its set point T d m (t) by mean of an external controller taking as input the error e(t) = l d (t)l(t), where l d (t) is the corresponding set point of the position l(t) of the liquid-solid interface. The proposed global control strategy is shown in Fig 2.

The external controller in Fig. 2 is designed for the following nonlinear system

Ṫm (t) = - 1 µ T m (t) + 1 µ T d m (t) (42)
l(t) = f (T m (t)) (43)
Note that the bijective function f ( . ) can be identified by assuming several values of the desired set point T d m (t) and their corresponding set points l d (t) (see the application example given is Section 6).

Melting of the stainless steel

The tracking performance of the proposed control strategy is evaluated in the case of melting of a strip of stainless steel of length l max = 0.5 m the physical properties of which are given in Table 1. The thermal conductivity k is a nonlinear function of the temperature given as [START_REF] Taler | Solving Direct and Inverse Heat Conduction Problems[END_REF] k(T (z, t)) = 14.6 + 1.27 10 -2 T (z, t)

where T is expressed in • C and k in W.m -1 .K -1 . Thereafter, eq. ( 32) gives

ϕ -1 (T (z, t)) = T (z, t) + 4.3493 10 -4 T 2 (z, t) (45)
Note that, for some complex thermal conductivity k, the analytical expression of ϕ -1 ( . ) cannot be calculated so that, in that case, numerical integration methods should be used.

For the closed-loop system simulation, the boundary immobilization method [START_REF] Caldwell | Mathematical Modelling. Case Studies and Projects[END_REF][START_REF] Maidi | Boundary control of linear Stefan problem[END_REF] is used to derive an equivalent nonlinear model with an immobilized right-hand boundary condition for the original nonlinear system (1)- [START_REF] Mattheij | Partial Differential Equations[END_REF]. Then, the equivalent nonlinear model obtained is simulated using the method of lines,based on finite differences, by taking 200 discretization points. The integral term of the control law (40) is evaluated using the trapezoidal method.

The relation, at steady-state, between the position of the liquid-solid interface position and the spatial average temperature T m (t) is given by Fig. 3. From this figure, it follows that the relation between l(t) and T m (t) is approximately linear, that is,

T m (t) = a l(t) (46) Designation Symbol Value Density ρ 6820 kg • m -3 Melting temperature T f 1648 K Heat of fusion ∆H f 260000 J • kg -1 Heat capacity c p 490 J • kg -1 • K -1
T d m (t) = K c l d (ζ) -l(ζ) + 1 τ i t 0 l d (ζ) -l(ζ) dt (48)
where K c , τ i and l d (t) are respectively the proportional gain, integral time constant and desired output of the liquid-solid interface position l(t).

The controller tuning parameter µ is taken equal to 5 min. The PI parameters are tuned so that the closed-loop system l d (t)l(t) is characterized by a damping factor η = 1 and a natural pulsation ω n = 0.18 rad • min -1 . The obtained parameters are K c = 1318.8 and τ i = 4.9383 min. In order to have a smooth control signal, the desired set point l d (t) has been filtered by a second-order filter with time constants 1 min and 2 min.

In the following, the tracking capabilities of the proposed control strategy and its robustness against the parametric variation are investigated through numerical simulations.

Tracking problem

The objective is to melt a portion of strip of length 0.3 m, initially solid. To show the effectiveness of the developed controller, the melting process is achieved in two stages. First, a portion of 0.15 m is melted, once this is achieved, a new melting is started to melt the rest of the strip. Thus, two step set points l d (t) = 0.15 m and l d (t) = 0.30 m are imposed at t = 0 and t = 80 min, respectively. Fig. 4 shows that the imposed set points are perfectly tracked, i.e. the liquid-solid interface is moved to the specified positions. From Fig. 5, it follows that the heat flux evolution is physically acceptable and it is observed that the heat flux, initially zero, increases rapidly just after the variation of set point, then returns to its steady value as the liquid-solid interface reaches the specified positions.

Robustness

To evaluate the robustness of the controller, a perturbed model is simulated for a variation of the liquid specific heat c p assumed as parametric uncertainty. The robustness test performed consists in assuming that the control law (40) is calculated by assuming a constant liquid specific heat c p = 454 J • kg -1 • K -1 whereas for the process is assumed to be temperaturedependent, that is, time-varying parameter. In this case, the thermal con- ductivity c p is a nonlinear function of the temperature T (z, t) given as [START_REF] Taler | Solving Direct and Inverse Heat Conduction Problems[END_REF] c(T (z, t)) =454 + 0.388 T (z, t) + 3.22 10 -4 T 2 (z, t) + 1.1 10 -7 T 3 (z, t)

(T in • C) (49)
For this test, the objective is to melt a portion of strip, 0.25 m in length. Figures 67show that in spite of the variation of c p , the control objective is correctly achieved with admissible moves of the manipulated heat flux.

Conclusion

In this paper, a control strategy is proposed to control the position of the liquid-solid interface in a nonlinear Stefan problem. This is modeled by coupled partial and ordinary differential equations. The PDE describes the temperature evolution of the phase liquid, which is a nonlinear diffusion equation as the thermal conductivity and the specific heat are temperaturedependent. The motion of the liquid-solid interface is modeled by an ODE. The objective is to control the interface position by manipulating the heat flow at the opposite boundary.

Based on the tangent Cole-Hopf transformation, the control law of the nonlinear Stefan problem is deduced directly from the control law developed by the authors, in their previous work, for the control of Stefan linear problem. Thereafter, based on the fact that the Cole-Hopf transformation is a bijective mapping, it is shown that the developed control law yields a stable closed-loop system.

Finally, to achieve the initial objective control, i.e., to move the liquidsolid interface to a desired position, a control strategy is proposed where an external controller is introduced to define the set point of the spatial average temperature of the liquid phase. The error between the liquid-solid interface position and its desired set point is defined as the input of the external controller.

This control strategy is successfully tested by numerical simulation of the melting of a stainless steel strip. The obtained results show the ability of the controller to track the desired position of the interface with physically acceptable moves of the applied heat flux.

Figure 1 :

 1 Figure 1: Schematics of material melting.

Figure 2 :

 2 Figure 2: Boundary control strategy of the liquid-solid interface.

Figure 4 :

 4 Figure 4: Tracking: Evolution of the position of the liquid-solid interface l(t).

2 ]Figure 5 :

 25 Figure 5: Tracking: Evolution of the manipulated heat flux q(t).

Figure 6 :

 6 Figure 6: Robustness: Evolution of the position of the liquid-solid interface l(t).

2 ]Figure 7 :

 27 Figure 7: Robustness : Evolution of the manipulated heat flux q(t).

Table 1 :

 1 Stainless steel properties.Figure 3: Steady-state relation between T m (t) and l(t).

		400	Interpolation points
		350	Function f ( . )
		300
	T ss m (t) [ K • m ]	150 200 250
		100
		50
		0 0	0.05	0.1 l ss (t) [ m ] 0.15	0.2
	Hence, combining (41) and (46), the following transfer function between
	T d m (t) and l(t) results	L(s) T d m (s)	=	a µ s + 1	(47)
	In this application, as the relation (46) is linear, hence the set point T d m (t)
	is defined by means of a PI controller (external controller in Fig 2) as follows