Thérèse Azar

Laetitia Perez

Christophe Prieur

Emmanuel Moulay

Laurent Autrique

Quasi-Online Disturbance Rejection for Nonlinear Parabolic PDE using a Receding Time Horizon Control

Null controllability of nonlinear partial differential equation is a very complex challenge. The context underlying this study is to improve the behavior of the plasma in a tokamak reactor in order to lengthen the duration of the nuclear fusion process. Considering the class of a specific parabolic PDE, the well known heat equation is nonlinear if thermal properties are temperature dependent. In such a context a numerical method based on the resolution of inverse heat conduction problem is proposed. It aims to provide identified control laws quasi-online in order to guarantee that the thermal state is kept close to its equilibrium state at zero. The iterative conjugate gradient method is implemented in order to control the temperature in the one-dimensional spatial domain despite several disturbances (time-dependent or thermo-dependent). The proposed strategy is based on successive numerical resolutions of ill-posed inverse problem on receding time horizons which are adapted considering the previous evolution of the system. Numerical results in the investigated configuration highlight that identified control laws are able to reject disturbances and to ensure null controllability.

I. INTRODUCTION

Partial differential equations (PDE) are widely investigated in physics and the development of efficient control strategies is a key-requirement for numerous processes. A classic objective is to seek to maintain the state of the system in the neighbourhood of a desired target regardless of the evolution of the uncontrolled inputs. In the following a parabolic nonlinear PDE is studied. A one dimensional domain is considered and the study is motivated by the control of nuclear fusion for which both magnetic flux density and thermal state is described by such a PDE system [START_REF] Attar | Plasma thermal characterization in nuclear fusion context[END_REF], [START_REF] Gaye | Robust stabilization of the current profile in tokamak plasmas using sliding mode approach in infinite dimension[END_REF], [START_REF] Bribiesca Argomedo | Safety Factor Profile Control in a Tokamak[END_REF] . Numerous approaches have been developed in the past decades: sliding mode approach in infinite dimension, feedback control [START_REF] Wijnands | Feedback control of the current profile on Tore Supra[END_REF] , H ∞ control [START_REF] Gaye | H∞ stabilization of the current profile in tokamak plasmas via LMI approach[END_REF], model predictive control [START_REF] Ou | Model predictive control of parabolic PDE systems with dirichlet boundary conditions via galerkin model reduction[END_REF], control-oriented model [START_REF] Witrant | A control-oriented model of the current in tokamak plasma[END_REF], . . . However, disturbances rejection is complex to ensure from a theoretical point of view if mathematical model nonlinearities are not neglected. Moreover, if the location of the actuators is different from the spatial area where disturbances occur, the problem of optimal control is quite difficult to deal with.

In the following the method presented in [START_REF] Azar | Stabilization using in-domain actuator: A Numerical Method for a Non Linear Parabolic Partial Differential Equation[END_REF] is modified in order to provide a quasi-online numerical strategy for the identification of control laws. The goal is to identify the control that keeps the system state close to zero (null controllability) using the conjugate gradient problem. This method is relevant for ill-posed problem such as Inverse Heat Conduction Problem IHCP [START_REF] Alifanov | Extreme Methods for Solving Ill-Posed Problems with Applications to Inverse Heat Transfer Problems[END_REF], [START_REF] Morozov | Methods for solving incorrectly posed problems[END_REF]. In thermal context, adaptive selection of relevant sensors in a network for unknown mobile heating flux estimation has been proposed in [START_REF] Vergnaud | Adaptive selection of relevant sensors in a network for unknown mobile heating flux estimation[END_REF]. The approach presented in the following is different from model predictive control but corresponds to a delayed command [START_REF] Yue | Delayed feedback control of uncertain systems with time-varying input delay[END_REF], [START_REF] Roy | Adaptive-robust Control With Limited Knowledge on Systems Dynamics: An Artificial Input Delay Approach and Beyond[END_REF]. From the theoretical point of view, the null controllability problem for the heat equation has been investigated for example in [START_REF] Bárcena-Petisco | Null controllability of the heat equation in pseudo-cylinders by an internal control[END_REF].

The paper is organized as follows. In addition to this introductory section, Section II is devoted to the mathematical description of the physical system. Direct problem is numerically solved in order to describe temperature evolution of the uncontrolled nonlinear system from a given initial state. Effect of several disturbances is highlighted and analyzed considering a relevant cost-function. In Section III, an iterative regularization method is described in order to solve the IHCP for which the unknown parameter which has to be identified is the control (depending on time and space). Adaptation of such method to an online implementation is presented in Section IV and numerical results are shown in order to discuss the method performances. Comparison between several online strategies are presented in Section V. In the last section concluding remarks and several outlooks are proposed.

II. PROBLEM STATEMENT A. Mathematical Description

Let us consider the direct problem modeled by the following parabolic PDE:

       ρC ∂θ ∂t - ∂ ∂x λ(θ) ∂θ ∂x = u [a,b] + f [c,d] θ(0, t) = θ(L, t) = 0 θ(x, 0) = θ 0 (x) (1)
The space variable is x ∈ Ω = [0, L], t ∈ [0, t f ] is the time variable where t f is the final time, ρC is the volumic heat in J.m -3 .K -1 , θ(x, t) in K is the temperature at point x at the instant t. It is important to specify that θ is the temperature variation compared to a reference temperature (in what follows, negative temperatures therefore mean that the temperature is lower than the reference temperature). λ(θ) is the thermal conductivity in W.m -1 .K -1 .

B. Direct Problem Resolution and Disturbance Impact

In this work, input parameters of (1) are taken into account: L = 0.1m, t f = 900s and ρC = 10 6 J.m -3 .K -1 . Thermal conductivity is assumed to depend on temperature:

λ(θ) = 4 exp -θ 2 10 4 + 1.
In order to describe the impact of the disturbance, let us consider the following disturbance function localized at x ∈ [0.04, 0.08] such as:

f (x, t) =        10 6 (t -240) 60 if t ∈ [240, 300] 1.4 × 10 4 × θ(x, t) + 10 4 if t > 600 0 if not
Such a disturbance which depends on the temperature introduces an interesting non-linearity. It could be also suitable to consider random disturbances. Initial temperature is:

θ 0 (x) = 25 1.5exp -(x -0.03) 2 5 × 10 -5 + exp -(x -0.05) 2 10 -4 -0.5exp -(x-0.08) 2 5×10 -5
.

This configuration aims to show how the control law acts to drive the state of the system more quickly from its initial state to the equilibrium state (zero without disturbances) and how the two types of disturbances (time-dependent for t ∈ [240, 300] and thermo-dependent for t ∈ [600, 900]) are rejected.

It is important to notice that the orders of magnitude of all the previous parameters are quite realistic. With the thermophysical parameters previously defined, ( 1) is solved numerically using finite element method with Comsol-Multiphysics solver interfaced with Matlab program. Without control i.e. u [a,b] (x, t) = 0, temperature evolution is shown in Fig. 1. for x ∈ [0.04, 0.08] and t > 600 then

Based on Fig. 1., it is obvious that if f (x, t) = 1.4 × 10 4 × θ(x, t) + 10 4
lim t→∞ θ(x, t) 2 L 2 (Ω) = lim t→∞ L 0 [θ(x, t)] 2 dx = +∞.
In such configuration and due to the disturbance, temperature θ(x, t) governed by ( 1) is not stable. In order to illustrate how the temperature is affected by the disturbances f , it is relevant to compute the following criterion (2) which is shown in Fig. 2.: Considering previous figures from 0 to 240 seconds, the system evolves from its initial temperature and naturally tends towards its equilibrium state (θ(x, t) = 0, ∀x ∈ Ω) without disturbances. Then between 240 and 300 seconds, the system is affected by a local (x ∈ [0.04, 0.08]) disturbance whose amplitude increases linearly as a function of time. State system is thus quickly moved away from zero. When this first disturbance vanishes at t = 300s, temperature naturally converges towards zero. Finally, from 600 seconds, the system is once more disturbed by a function which depends on the temperature. As it has been previously mentioned, such disturbance leads to non stability. This numerical study illustrates the disturbance impact and will be used as a try-model for our control problem.

J(θ) = 1 2 θ(x, t) 2 L 2 (Ω) . (2)
In order to reduce the effect of disturbances, an inverse problem is formulated in the next section.

III. IHCP RESOLUTION

In this section, an optimization method is proposed in order to minimize the cost-function defined in (2). This method is based on the Conjugate Gradient Method (CGM) which acts as an iterative regularization method (convergence of this well-known minimization method is discussed for linear systems in [START_REF] Strakos | On the real convergence rate of the conjugate gradient method[END_REF], [START_REF] Shi | The convergence of conjugate gradient method with nonmonotone line search[END_REF] and for specific non-linear system in [START_REF] Hager | A survey of nonlinear conjugate gradient methods[END_REF], [START_REF] Sabi'u | Two optimal Hager-Zhang conjugate gradient methods for solving monotone nonlinear equations[END_REF]). Examples of numerical implementation for identification purposes in a thermal context are given in [START_REF] Vergnaud | Quasi-online parametric identification of moving heating devices in a 2D geometry[END_REF], [START_REF] Beddiaf | Heating source localization in a reduced time[END_REF] or [START_REF] Vergnaud | Quasi on-line identification of temperature-dependent thermal conductivity[END_REF]. However, in the works cited above, the authors never broached the question of control.

Minimization algorithm implemented for IHCP resolution is presented hereafter according to the following notations:

time interval T = [t * o , t * f ] is considered, control law u(x, t
) is discretized as a piecewise continuous linear function defined by u(x i , t j ) = u ij . Temperature is measured by sensors equally spaced in domain Ω = [0, L].

A. Conjugate gradient algorithm

At each iteration of the algorithm, three well-posed problems have to be solved:

• The direct problem in order to determine the temperature distribution θ k (x, t * f ), and then to estimate the criterion J(θ k );

• The adjoint problem to determine the gradient of the cost function J(θ k (x, t * f )) and thus to define the next descent direction d k ;

• The sensitivity problem to estimate the descent depth γ k (in the descent direction).

Optimization algorithm is briefly described as follows in order to minimize criterion (2) (in [START_REF] Azar | Stabilization using in-domain actuator: A Numerical Method for a Non Linear Parabolic Partial Differential Equation[END_REF] a similar algorithm is presented in offline context):

1) Initialization of the unknown parameter (control flux) at the first iteration k = 0: u k ij = 0; 2) Resolution of the direct problem (1) in order to determine θ k (x, t), 3) Determination of the criterion J(θ k (x, t * f )) according to (2); 4) Solve the adjoint problem (5) in order to determine the Lagrange multiplier ψ k (x, t) and the cost function gradient according to [START_REF] Ou | Model predictive control of parabolic PDE systems with dirichlet boundary conditions via galerkin model reduction[END_REF]. Determination of the descent direction d k according to (7); 5) Resolution of the sensitivity problem (3) in the descent direction d k to calculate the sensitivity function δθ k (x, t) and determination of the descent depth γ k according to (4); 6) Update new estimations of the control according to

[u ij ] k+1 = [u ij ] k -γ k [d ij ] k 7)
Increment of the iteration k = k + 1 and back to step 2. In the following section, sensitivity problem is briefly presented.

B. Sensitivity Problem

In order to calculate at iteration k the descent depth γ k in the descent direction d k , the sensitivity problem has to be solved. Let us consider temperature variation:

θ(x, t) + ε 0 δθ(x, t)
induced by a variation of the control (thermal flux):

u(x, t) + ε 0 δu(x, t).
Sensitivity function δθ is solution of this so-called sensitivity problem:

       ρC ∂δθ k ∂t - ∂ 2 λ(θ)δθ k ∂x 2 = ξδθ k + δu k δθ k (0, t) = δθ k (L, t) = 0 δθ k (x, t * 0 ) = 0 (3)
where

ξ(x, t) = 1.4 × 10 4 if t > 600 and x ∈ [0.04, 0.08] 0 if not
It should be noticed that in the previous equation, the coefficient ξ(x, t) is due to the disturbance f (x, t).

The optimal descent depth γ k is determined as follow:

γ k = arg min γ∈R J(u k a -γd k ),
and estimated as in the pioneer work [START_REF] Jarny | A general optimization method using adjoint equation for solving multidimensional inverse heat conduction[END_REF]:

γ k = - L 0 θ k (x, t * f )δθ k (x, t * f )dx L 0 δθ k (x, t * f ) 2 dx (4)
where θ k (x, t * f ) is the solution of the direct problem and δθ k (x, t * f ) is the solution of the sensitivity problem (3) (solved in descent direction d k ).

In the following section, adjoint problem is developed in order to determine the descent direction d k .

C. Adjoint Problem

In order to calculate at each iteration k, the gradient ∂J(θ k ) ∂u k ij and the descent direction d k , an adjoint problem is formulated. Let us denote by a Lagrangian formulation which is a function of u k ij (.), θ k (.) and ψ k (.) where ψ k (.) is the adjoint function:

(u k , θ k , ψ k ) = J(θ k (x, t * f )) + t * f t * 0 L 0 F (.)dxdt
where F is defined as:

F (.) = ρC ∂θ k ∂t - ∂ ∂x λ(θ k ) ∂θ k ∂x -u k -f ψ k
The adjoint function ψ k (.) is fixed such as ∂ ∂θ k δθ k = 0. Thus ψ k (x, t) is solution of the following adjoint problem:

         -ρC ∂ψ k (.) ∂t -λ(θ) ∂ 2 ψ k (.) ∂x 2 = ξψ k (.) ψ k (0, t) = ψ k (L, t) = 0 ψ k (x, t * f ) = - 1 ρC θ k (x, t * f ) (5)
Cost function gradient is then obtained as follows:

∂J ∂u k ij = - t * f t * 0 b a ψ k (x, t)s i (x)s j (t)dxdt (6)
where s i (x) is the basis function in space and s j (t) is the basis function in time for piecewise continuous linear functions (in space and time). Descent direction can be estimated at each new iteration k from the previous gradient (at iteration k -1), as follows:

d k = - ∂J ∂u k ij + β k d k-1 (7)
with

β k = ∂J ∂u k ij 2 ∂J ∂u k-1 ij 2 and . is the Euclidean norm.
In [START_REF] Witrant | A control-oriented model of the current in tokamak plasma[END_REF] the previous optimization method has been successfully implemented in order to determine a control strategy which ensures the convergence to zero. The problem was solved offline since T = [0, t f ] i.e. t * 0 = 0 and t * f = t f . In the following the online adaptation is presented.

IV. QUASI-ONLINE IMPLEMENTATION

Three quasi-online strategies based on an adaptation of the iterative regularization method of the CGM are implemented in order to determine a relevant control u which steers the temperature to zero. The proposed approaches are based on the choice of receding intervals T . Let us introduce general comments:

• For the three strategies, it is important to have enough observations contained in T in order to understand the evolution of the system and to be able acting on it; • Temperature is measured every centimeter of the domain. Such observation strategy could be easily modified considering experimental constraints; • The control law is discretized with a spatial step of 3 millimeters inside [0.03; 0.06]. This arbitrary choice leas to determine at each time step nine unknown parameters. If this spatial step is too large, then unknown required control distribution might be not possible to describe. If this spatial step is too small, a large amount of unknown parameters is not useful and overparametrization problems could appear; • The choice of the time step has obviously a significant effect on the performance of the method. It must be realistic and consistent with the technologies and based on a knowledge of the dynamics of heat transfers. A time step equal to 1 second for both the control and the observation has been chosen in this paper. • Last but not the least, spatial interval, time discretization and space discretization for disturbances, measurements and control may be different. The formulation of the methodology is not affected by such choices (but it is obvious that results are). In the next sections, the three strategies are briefly presented.

A. Strategy #1 -constant offset

For the first strategy, τ in seconds is a constant offset. Strategy #1 acts sequentially in order to provide updated control laws at constant intervals. For example, let us denote by T m the interval [mτ, (m + 1)τ ].

• for m = 0: the process starts at t 0 = 0, the system evolves "naturally" from its initial state and under the effect of any external disturbances. After τ seconds, the temperature is measured throughout the plate and allows us to know θ(x, τ ); • for m = 1: the temperature continues to evolve freely without controller and simultaneously the first control law is estimated considering the measurements collected during the time interval T 0 ; • for m = 2: the first control validated for the time interval T 0 is applied and simultaneously the second control law is estimated considering the measurements collected during the time interval T 1 ; • for m = 3, the second control validated for the time interval T 1 is applied and simultaneously the third control law is estimated considering the measurements collected during the time interval T 2 ; • etc. IHCP numerical resolution requires a computational time (lower than τ ) to identify the command which acts with a delay of 2τ . Results obtained for the strategy #1 are shown for τ = 6 seconds in Fig. 3., Fig. 4. and Fig. 5.

With the identified control law, disturbances effect is reduced as can be seen in Fig. 3. and Fig. 5. The second disturbances (after 600 s) is rejected more efficiently because its effect is proportional to the temperature. In Fig. 4., control evolution is shown at location x = 0.045.

B. Strategy #2 -adaptive offset

For the second strategy new control is updated only when the final temperature distribution deviates too far from a previously fixed threshold related to a given J. Temperatures are measured every second. If J(θ) > J then, the identification procedure is launched taking into account the previous measurements on interval T m whom size is equal to τ seconds. This strategy is called adaptive offset. Results obtained with strategy #2, are shown for τ = 6 and J = 0.5 in the next Fig. 6., Fig. 7. and Fig. 8. It is shown in Fig. 6., Fig. 7. and Fig. 8 that disturbance is attenuated even if control laws are calculated less often than according to strategy #1.

C. Strategy #3 -adapted duration

As for the previous strategy, the temperatures are measured every second and the criterion is calculated every second. Strategy #3 differs from the other previous strategies by the modification of the duration of receding time interval T m = [t * 0m , t * f m ]. In fact, control laws are estimated only if temperature observations leads to a criterion J(θ(x, t)) > J max . In such a case, t * f m = t and t * 0m is the previous time where J(θ(x, .)) > J min . Thus J max is considered as a prohibitive threshold (for which temperature distribution is too far from zero) and J min is considered as significant threshold (for which temperature distribution begins to move slightly away from its equilibrium state). Strategy #3 authorizes the recovery of the control laws which can reduce delays. Since the control laws have variable duration, this strategy is called adapted duration strategy.

The request to determine a new control law is triggered as soon as the criterion exceeds prohibitive threshold J max and the new control is implemented as soon as it is determined. Numerical result obtained with J min = 0.1 and J max = 0.5 are shown in Fig. 9., Fig. 10. and Fig.

V. COMPARISON OF DIFFERENT STRATEGIES

In this section, different tables are proposed in order to compare the previous strategies. Obviously, the closer the criterion is to zero, the better the equilibrium state is obtained (Identified control laws are relevant).

A. Numerical Results

Table I corresponds

Without control: M ≈ 3.

Table I shows that the strategy #3 with a low trigging threshold J max = 0.05 ensures the best null-controllability. It has the closest temperature to zero. Strategies #1 and #2 with a small constant offset, provides good results and keeps the domain temperature close to zero. These values are calculated from the distributions θ(x, t) which are calculated for each of the three strategies. Without control T = 900 seconds (obtained from the distributions plotted in Fig. 1). Indeed, without control,

u

  [a,b] (x, t) is the control applied to the interval [a, b] ⊂ ]0, L[ which aims to keep the system state θ close to zero despite disturbances f [c,d] (x, t) acting on [c, d] ⊂ ]0, L[. Boundaries condition for x ∈ {0, L} are Dirichlet conditions and θ 0 (x) is the initial temperature at t = 0.

Fig. 1 .

 1 Fig. 1. Temperature evolution without control.

Fig. 2 .

 2 Fig. 2. Criterion evolution without control.

Fig. 3 .

 3 Fig. 3. Temperature evolution with control (strategy #1; τ = 6).

Fig. 4 .

 4 Fig. 4. Control law evolution at x = 0.045 (strategy #1; τ = 6).

Fig. 5 .

 5 Fig. 5. Criterion evolution with and without control (strategy #1; τ = 6).

Fig. 6 .

 6 Fig.6. Temperature evolution with control (strategy #2; τ = 6; J = 0.5).

Fig. 7 .

 7 Fig. 7. Control law evolution at x = 0.045 (strategy #2; τ = 6; J = 0.5).

Fig. 8 .

 8 Fig.8. Criterion evolution with and without control (strategy #2; τ = 6; J = 0.5).

  .

Fig. 9 .

 9 Fig. 9. Temperature evolution with control (strategy #3; J min = 0.1; Jmax = 0.5).

Fig. 10 .

 10 Fig. 10. Control law evolution at x = 0.036 (strategy #3; J min = 0.1; Jmax = 0.5).

Fig. 11 .

 11 Fig. 11.Criterion evolution with and without control (strategy #3; J min = 0.1; Jmax = 0.5).

  max x∈Ω | θ(x, t) | > 1 • C.

  max x∈Ω | θ(x, t) | > 4.39 • C.

Table

  II provides the duration T during which

This work was supported by the regional programme Atlanstic 2020, funded by the French Region Pays de la Loire and the European Regional Development Fund.

Table II shows that in the studied configuration and according to the parameters chosen, the strategy #1 is the one which makes it possible to remain most often under the threshold temperature fixed here at 1 • C. Strategy #3 could give better results by lowering the trigger threshold J max .

Table III provides the cost control U for each strategy such as:

Table III shows that strategy #3 requires less energy to maintain the temperature near zero (system state equilibrium θ = 0).

Table IV provides the time t act in which the control is not null.

Table IV shows that the strategy #1 requires a lot of calculations to identify the control laws. This could be obviously explained by the fact that the controls are calculated even when it is not useful. We can also notice that strategy #2 requires few calculations but the three previous tables have shown that overall it was less effective than strategies #1 and #3.

VI. CONCLUDING REMARKS

In this paper, the identification of control law for nonlinear parabolic PDE with homogeneous Dirichlet boundaries conditions is formulated as an inverse problem. It has been shown how disturbances considerably move the temperature away from its equilibrium state. In order to ensure the null-controllability for such nonlinear PDE where internal controls and disturbances are not collocated, three quasi-online strategies have been proposed. The third strategy based on a receding time interval with an adapted duration (related to prohibitive and significant threshold) is the most relevant. Disturbances effect is divided by 10 and this strategy requires less energy comparing to the others one.

In the next future this methodology will be implemented and tested for the control of the nuclear fusion in order to improve the duration of the plasma. In such a way a coupled system of two parabolic nonlinear PDE will be investigated as investigated in [START_REF] Bribiesca Argomedo | Safety Factor Profile Control in a Tokamak[END_REF].