Anastasia Bolotnikova

Sébastien Courtois

Abderrahmane Kheddar

Compliant Robot Motion Regulated via Proprioceptive Sensor Based Contact Observer

In this paper, we present developments of a realtime compliant motion control for a personal humanoid robot. Our approach allows to interpret and react to human guidance through touch using only joint encoders measurements to monitor contact direction and intensity on both static and moving links. This novel method is developed with consideration of minimal sensor requirement of the hardware platform to meet high affordability criteria. We demonstrate performances in experiments with a humanoid robot Pepper.

I. INTRODUCTION

In various fields, such as manufacturing or health care, an efficient automation may require work tasks to be performed by a robot in close-contact with a human. To achieve safe and efficient collaboration, it is important that agents understand each other's intentions and react to them appropriately. In this work, we focus on physical human robot interaction (pHRI) and, foremost, on the interpretation and correct online reaction of a humanoid robot to human touch.

The human-robot collaboration must ideally be safe, efficient and optimal in terms of cost. Thus, a minimal sensor set-up is preferable in order to reduce robot manufacturing and maintenance costs. As a result, we focus on more costefficient, but in practice more challenging way of contact sensing by using minimum of proprioceptive sensors. That is to say, we use only robot's joint encoders.

Recent developments in the area of proprioceptive sensor based contact monitoring have shown promising results for serial manipulators [START_REF] Haddadin | Robot collisions: A survey on detection, isolation, and identification[END_REF], humanoids [START_REF] Flacco | Residual-based contacts estimation for humanoid robots[END_REF] and also a special case of affordable personal robots, such as Pepper, used in this work, with high motor-joint backlash and friction [START_REF] Flacco | Contact detection and physical interaction for low cost personal robots[END_REF]. In our previous work [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF], we proposed a contact observer method that uses data-driven machine learning techniques to detect a colliding link and monitor contact direction and intensity. This method was developed to fit well into minimal sensor set-up, as only joint encoder measurements are exploited for contact detection for both moving and static links.

In this paper, we describe the modifications made to our initial method [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF] in order to increase its accuracy and generalization properties in the case of multi-joint motions; we demonstrate its experimental performances in real-time feedback compliant motions regulation on the arm of a humanoid robot Pepper.

II. BACKGROUND

In order to achieve active compliant motion of a robot, it is necessary to sense or estimate the external forces applied (anywhere) on the robot's links. A sophisticated yet expensive solution would require to cover the entire robot structure with a haptic skin. Such an approach, however, is not appropriate due to the cost of adaptation of additional sensors (even with cheap components) and their maintenance; the introduction of additional weight to the robot structure; additional power consumption... see eg. [START_REF] Kheddar | A tactile matrix for whole-body humanoid haptic sensing and safe interaction[END_REF]. Thus, researches have directed their efforts towards achieving contact sensing by using must-be proprioceptive sensors only (i.e. those usually available on a robot platform), such as joint encoders, motors' electric current or torque sensors, etc.

A thorough overview of such contact sensing methods is presented in [START_REF] Haddadin | Robot collisions: A survey on detection, isolation, and identification[END_REF]. Methods, which were initially developed for serial manipulators, have been further developed to be applied on humanoid robots, e.g. [START_REF] Flacco | Residual-based contacts estimation for humanoid robots[END_REF], [START_REF] Mattioli | Interaction force reconstruction for humanoid robots[END_REF], [START_REF] Vorndamme | Collision detection, isolation and identification for humanoids[END_REF]. In [START_REF] Mattioli | Interaction force reconstruction for humanoid robots[END_REF] a proprioceptive sensor based external force reconstruction was used as a feedback in the control loop to regulate robot interaction forces with the environment -under assumption of static conditions for external force reconstruction, whereas we aim at developing a method that works equally well for static and dynamic motions. Compliant reaction to human touch using proprioceptive sensors has been demonstrated in [START_REF] Flacco | Contact detection and physical interaction for low cost personal robots[END_REF] in the experiment with a Romeo humanoid robot arm. However, this method could not be adapted for the current version of the Pepper platform. This is due to the inability to measure motor-link backlash and to measure or accurately estimate motor torque, which in turn prohibits to estimate external torque by means of the momentum observer method [START_REF] Luca | Sensorless robot collision detection and hybrid force/motion control[END_REF]. The latter method was used in robot control to achieve safe reaction to external collisions in pHRI context, see examples in [START_REF] Luca | Sensorless robot collision detection and hybrid force/motion control[END_REF], [START_REF] Haddadin | Collision detection and reaction: A contribution to safe physical humanrobot interaction[END_REF], [START_REF] Luca | Exploiting robot redundancy in collision detection and reaction[END_REF], [START_REF] Parusel | Modular state-based behavior control for safe human-robot interaction: A lightweight control architecture for a lightweight robot[END_REF].

In order to address a more challenging task of contact sensing without using cover haptics, joint torque or even motors' electric current sensors; we developed a position tracking error based contact observer in [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF]. This method is based on monitoring the discrepancy between the measured position tracking error and the expected one for a given desired free motion trajectory. First, a machine learning model is trained on a sample dataset with appropriate set of desired trajectory related features to predict the expected position tracking error value for collision-free joint motion. Then, for every iteration (of the control loop) the trained model is used to predict the expected position tracking error, assuming no contact has occurred. The contact is detected whenever there is significant discrepancies between measured and predicted (and therefore expected) error values. In our previous work, this method has shown promising results in contact detection experiments, and in monitoring both the contact intensity and direction. Now, we exploit our contact observer in realtime control for regulating compliant motion of Pepper. The motivation behind this goal is to enable efficient physical interaction between human and a robot. In this work, we demonstrate that enhanced version of our contact observer can be exploited for the online interpretation of the human touch and regulation of compliant motion.

The use of disturbance observers for the regulation of robot compliant motion has been explored is several works. In [START_REF] Hyon | Full-body compliant humanhumanoid interaction: balancing in the presence of unknown external forces[END_REF] whole-body compliant motion in human-humanoid interaction settings was proposed for maintaining balance of the torque-controlled bipedal humanoid platform in the presence of unknown external forces. A momentum based disturbance observer was applied to the floating-base model of a humanoid robot in [START_REF] Ott | Kinesthetic teaching of humanoid motion based on whole-body compliance control with interactionaware balancing[END_REF] to detect external forces and integrate them into a whole-body control for kinesthetic teaching and simultaneous compliant balancing. These methods, however, require a measurement of the joint torques, which may not be available and may also be challenging to estimate in case of highly geared and flexible joints encountered on many humanoid platforms.

In [START_REF] Kim | Disturbance observer based linear feedback controller for compliant motion of humanoid robot[END_REF] disturbance observer based compliant motion control scheme, that can compensate for high joint elasticity, was proposed including the modeling of flexible joints that was performed using the measurement of motorjoint backlash angle. On the current position-controlled Pepper platform, access to the motor side encoder measurements is not available through robot's centralized memory, prohibiting to handle motor-joint backlash and joint flexibility. There are also no motor torque sensors on the platform. Furthermore, estimation of the motor torque from electric current measurements is not feasible, due to the measurements being passed to the centralized memory in absolute and downsampled form. Nevertheless, in our work we try to overcome those constraints by exploiting the machine learning datadriven approach for disturbance observer.

We describe modifications of the contact observer w.r.t our previous work, with details about machine learning prediction model training process (Sec. III), and present the compliance control scheme and controller implementation details (Sec. IV). We demonstrate compliant motion results achieved with Pepper in pHRI experiments (Sec. V).

III. CONTACT OBSERVER FOR MULTI-JOINT MOTIONS

First we review the idea behind our methodology in [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF]; then we present the details of the updated framework we propose in order to enhance both robustness and performance of our initial method for the case of multi-joint motions.

A. Review of the Contact Observer

In our settings, a robot is controlled in a task-space closedloop by an acceleration resolved quadratic programming controller (QP) [START_REF] Bouyarmane | Multirobot and task-space force control with quadratic programming[END_REF]. As for now, state feedback is not used to compute desired trajectories -assuming that no contact has occurred, because when contact occurs, resultant force estimation is not yet integrated in the QP model. Moreover, at the test time, when external collisions are introduced, the QP control without state feedback computes the desired trajectories that in its nature match those trajectories that were present in the collision-free training set. At each control loop iteration desired joint accelerations qd are computed as an output of QP. We obtain values of desired velocity qd and position q d through numerical integration of qd , and compute desired torques τ d from inverse dynamics. Note, that τ d can also be obtained as an output of QP controller, by keeping torques as decision variables (not necessary in this context).

We define an expected position tracking error ǫ exp to be the kind of tracking error that occurs when joint moves freely and no external forces are applied to any of the robot links (i.e. external joint torques are zero, τ exp = 0). Our goal is to predict the value of ǫ exp given the information about the desired trajectory of the robot joints.

We have demonstrated in [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF] that -under some conditions, ǫ exp can be computed for a joint, actuated by a DC motor with PD control, as a linear combination of qd , qd , τ d and ǫexp (assuming no motor-link backlash). However, due to modeling errors and the lack of possibility to model some critical non-linear effects, such as a large motor-link backlash -which in practice cannot be ignored, we use a non-linear model, such as binary tree, to learn the relationship between the appropriate set of available features (q d , qd , τ d ) and the value of position tracking error, ǫ, from a sample free-motion data. As a result, we obtain: ǫ i exp = binary tree(q i d , qi d , τ i d ) for every joint i1 . The contact observer signal is defined as r = ǫǫ exp , where ǫ = q dq, with q being a vector of measured joint positions (joint encoder readings).

When a contact occurs at a link in the chain after joint i, the value of r i exceeds a predefined threshold ±δ. The magnitude and the sign of r i signal encode the intensity and the direction of the external force respectively.

Interested reader may refer to original work in [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF] for a more detailed explanation of the contact observer method.

B. Robustification of the Contact Observer

The use of r signal for contact monitoring has been demonstrated in several experiments with Pepper left and right arm joints, with one joint moving at a time, i.e. 1 degree of freedom (DoF) case [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF]. Occasionally false positive (FP) and false negative (FN) contact detections occurred during the experiments. Furthermore, when the system was tested on more complex motions, e.g. with all arm joints moving simultaneously with arbitrary speed to random set-points, it turned out to be not general enough to adapt to such motions. The reason for that is mainly the fact that sudden motions of one joint can have significant influence on the tracking error value of some other joints, especially in the presence of a significant motor-link backlash.

Such false detection cases have to be eliminated or minimized in order to use r as a feedback signal in the active compliance control. We identified the main reasons for undesired false detection cases. First, due to the insufficient prediction accuracy of the binary tree prediction model for some new data (not used for the model training), r value can exceed threshold ±δ even though external collision is not present, meaning that prediction model fails to generalize to "unseen" data and predict value of ǫ exp correctly in those cases. Secondly, in some circumstances (e.g. singularity or near joint limits configurations), external contact force does not cause r to exceed ±δ, thus the contact remains unnoticed. Similarly, in case of light contacts, r does not exceed ±δ, which can be addressed by lowering the value of δ (which was set to 2.5 • in previous work). However, to do that the overall prediction accuracy has to be improved, for both seen and unseen data, which may not be achievable due to the accuracy-generalization trade-off. And, finally, as mentioned previously, the prediction of expected tracking error of a joint depends on the motion of other joints. Therefore, relevant variables related to other joints' desired trajectories must be identified and included in the prediction feature vector.

The aim is to ensure that the trained models generalize well to a wide range of complex motions, while encountering fewest possible false positive contact detections. On the other hand, it is equally important that prediction is computed within the control time constraint for all joints. Therefore, prohibiting the use of very complex models and potentially sacrificing the maximum achievable prediction accuracy. The trade-off between model complexity and feature vector size and model accuracy must be properly handled.

In several application areas of machine learning, so called, ensemble based techniques have often shown to perform better in practice compared to single classification or regression models [START_REF] Dietterich | Ensemble methods in machine learning[END_REF]. In particular, boosting is a machine learning technique for ensemble model training that is suitable for our goal to battle false positive contact detection cases [START_REF] Drucker | Boosting decision trees[END_REF].

According to the boosting methodology, a cascade of several models is trained successively; every new model k = 2, . . . , N (N being the total number of models in the ensemble) is trained taking into consideration the performance of the previously trained models 1, . . . , k -1. At every new iteration k, previously known inaccurate predictions of ǫ exp are given more weight before training of k th model takes place. As a result, every new model in the ensemble is trained to make a better prediction, especially in those cases where previous models performed poorly. This results in minimization of the maximum prediction error, and thus less false positives in the r signal. We chose a cascade of boosted decision trees as our prediction model; we use XGBoost library implementation of this technique [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF].

To achieve the best possible prediction model accuracy and generalization, while minimizing the required computation time, we make several important choices about model hyperparameters. The following developments are presented for LShoulderRoll joint, also referred to as joint nr. 6 or LSRoll for short. Analysis for other joints is performed similarly.

1) Feature Vector Components: First of all, on a sample recording of 45 minutes of collision-free data (with time step of 12 ms, i.e 225000 samples in total) we analyze the importance of different variables for ǫ 6 exp prediction to select the most relevant variables for the feature vector. We train a boosted ensemble of decision trees with a feature vector that includes the desired motion related variables for all the main left arm joints. The results of relative variable importance computed from the trained model are summarized in Fig. 1.

⋅⋅

q ⋅⋅ q ⋅⋅ q ⋅⋅

Variable importance (%)

Fig. 1: Variable importance for ǫ 6 exp prediction.

As can be seen from the plot, the most important variables for ǫ 6 exp prediction are desired speed and torque of the target joint, which confirms our previous findings [START_REF] Bolotnikova | Contact observer for humanoid robot pepper based on tracking joint position discrepancies[END_REF]. However, desired speed, acceleration and torque of the nearby joint, LShoulderPitch (joint nr. 5), also have significant weight for prediction of ǫ 6 exp , which even overweight importance of q6 d . The variables after the orange dashed line on the plot are considered of medium-importance and variables after red dashed line are considered as non-important.

To make a more informed decision on the final size of the feature vector, we perform k-fold cross-validation with k = 4 (75% samples for training and 25% for validation) and analyze how the root-mean-square error (RMSE) and maximum absolute prediction error (MAX) on validation sets are affected as we incrementally include more of most important variables into prediction feature vector. The mean RMSE evaluated over 4 cross-validation folds only improves significantly if 3 to 4 most important variables are used (see Fig. 2). The mean MAX improves significantly until 6 most important variables are used and stops improving and even deteriorates after this point. Thus, we chose to include 6 most important variables in the final feature vector. Feature vector size Prediction error (deg) Fig. 2: Prediction error with various feature vector sizes.

2) Ensemble Size and Learning rate: Now that the best components of the feature vector have been set, we identify an optimal value for the ensemble size N , and a learning rate η of the model. Since we intend to use the trained model in the real-time application, we must select the model with the smallest possible ensemble size. Selecting a higher learning rate for the training allows to achieve good performance with fewer number of sub-models in the ensemble. However, higher learning rate may result in lower model accuracy, as in this case the training process is only allowed to make "big steps" towards an optimal solution. We analyze the mean cross-validation RMSE for various sizes of ensemble and for different values of the learning rate. The Fig. 3 illustrates the performance of models of various sizes trained with three different learning rates. We see that lower η values do not result in significant improvement of performance in terms of the mean RMSE over k cross-validation folds as the size of the ensemble increases. The difference in performance between the best and worst performing models is < 0.02 • . Therefore, we choose to use η = 0.3 with model size N = 50 as our optimal choice, since it provides a good performance while keeping ensemble size N relatively low.

3) Individual Decision Tree Maximum Depth: Now that other model hyperparameters are selected, we identify an optimal size of the individual decision trees. We train several models with N = 50 and η = 0.3 and analyze mean RMSE and MAX computed for training and validation folds. We also compute the rate of false positive detections as a total amount of absolute validation prediction errors exceeding the threshold ±δ divided by the total number of the validation samples. The results are shown in Fig. 4.

The result is such that with maximum tree depth of 7 the false positive rate decreases to < 0.1%. The values of RMSE for both the validation and the training data are acceptably low, always below soft and hard thresholds indicated as orange and red dashed lines on the plot respectively. The maximum absolute prediction error decreases for both seen and unseen data until maximum tree depth reaches 8, after this point, trained model becomes too complex and overfits the training data (training error decreases, while validation error remains the same or increases). Adapting regularization strategies did not help to improve the performance of the model on validation data. The final value for the maximum tree depth is, therefore, set to 8.

With this final decision on hyperparameters, we are able to achieve the performance of the contact observer that is acceptable for the use of r as a feedback in a real-time control for active compliance. The next section presents the compliant motion control implementation details.

IV. COMPLIANT MOTION REGULATED VIA MONITORING

OF JOINT POSITION DISCREPANCIES Here we describe the QP controller design that enables compliant motion using the contact observer signal as a feedback on collisions. The overview of the entire process is presented by pseudocode in Algorithm 1. For simplicity, consider a QP controller with one posture task in the objective function, and typical kinematics constraints, such as self-collision avoidance and joint limits. The posture task specifies a target posture, q target2 . The output, qd , of the QP controller generates a whole-body motion that brings the robot joints closer to q target while satisfying the kinematics constraints. The rate of convergence of the joints to q target is regulated via the QP posture task stiffness gain.

Before the QP controller starts, we read the initial robot sensor values, q init . For the first iteration of the controller execution, timestamp t = 0, we assume that no collision with the environment occurred, i.e. r(0) = 0. At this point, the posture task target position is equal to the values of the initial sensor reading (q target = q init ). The QP controller computes qd , that is integrated twice to have q d that is sent to the low-level PD servo and brings progressively the robot joints closer to q target . In the first iteration, since q target = q init the robot does not move, i.e. qd ≈ 0.

Using known robot model and QP output, the desired torque τ d is computed. Now the ǫ exp can be computed using the trained prediction models described in the previous section. At the same time, the measured position tracking error is computed as ǫ = q dq init . Finally, the contact observer signal, r = ǫǫ exp , is computed and passed to the next QP iteration along with the new sensor readings.

Starting from the second iteration, if |r| < δ the posture task remains unchanged. Otherwise, if a discrepancy between ǫ and ǫ exp has occurred (i.e. |r| ≥ δ), before the QP controller is inferred, the posture task has to be readjusted to comply in the direction of estimated disturbance. The compliant motion of every joint is expressed as follows:

q target (t) = q d (t -1) -K c r + K v qd (t -1) (1)
where K c and K v are compliance and velocity gains respectively, which are tuned for every joint. The posture task is updated in the QP objective function with new q target , and the entire process repeats until the controller is stopped. The compliance gain, K c , regulates how much influence contact observer signal r has on updated target posture. Lower values of K c allow the target posture to be updated only slightly, which means that QP will not compute high qd and robot will comply stiffly (little amount with lower speed). On the contrary, a high K c causes the q target of the posture task to change significantly resulting in high compliance and faster motion, as the qd is proportional to the QP task error. The velocity gain, K v , regulates the proportional contribution of joint desired speed to facilitate compliance of moving joints (those that started to comply in the previous iterations).

Note that this basic strategy can be augmented with any other QP tasks (e.g. set point, visual servoing, force tasks, etc.), by adding them to the QP objective function via higher level planning algorithm either along with compliant posture task or within periods when |r| < δ. In such case, a proper tuning of QP tasks' weight and stiffness values can influence greatly the overall performance of the complex QP controller.

V. EXPERIMENTAL RESULTS

We use our contact observer methodology in two pHRI experimental setups. The threshold for contact detection δ is set to 2.5 • for all joints in all experiments.

In the first scenario, a robot is making a motion that initiates the process of assisting a human to stand up from a sitting position. The robot moves its arm towards the human back and stops moving further if it detects a contact (i.e. K c and K v gains of Eq. 1 are set to zero). If no contact is detected the robot continues to move its arm until joint limits. The motion is repeated several times with varying values of the target LShoulderPitch (LSPitch for short) joint position, allowing to take contact lower or higher on the back of a human. The values of the desired joint positions Algorithm 1: The pseudocode of algorithm for the QP controller regulating joint compliance via r signal

q ← ROBOT.GET SENSORS() r ← 0 CNTRL ← INITIALIZE CONTROLLER() POSTURE TASK ← INITIALIZE POSTURE(q) CNTRL ← INSERT TASK(POSTURE TASK) CNTRL.RUN() while CNTRL.RUNNING() do if |r| ≥ δ then qtarget(t) ← COMPLY JOINTS(r, q d (t -1), qd (t -1)) CNTRL.POSTURE TASK ← UPDATE TASK(qtarget(t)) end if qd ← SOLVE QP(POSTURE TASK) qd ← INTEGRATE(q d ) q d ← INTEGRATE( qd ) ROBOT.SETJOINTANGLES(q d ) τ d ← INV. DYNAMICS(q d , qd , q d , ROBOT MODEL) ǫexp ← PREDICT(q d , qd , q d , τ d ) ǫ = q d -q r = ǫ -ǫexp q ← ROBOT.GET SENSORS() end while
and the contact observer of the LSRoll are shown in Fig. 5. When the contact observer signal exceeds the threshold, the commanded position for the LSRoll joint is altered such that robot stops instead of continuing to move forward.

In the second experimental scenario, the robot moves 4 joints of the left arm randomly. When the contact is detected, Eq. 1 with positive K c and K v gains, is used to guide the robot motions through human touch and comply the joints of the arm. The values for K c and K v gains are set manually for every joint; exact values used in the experiments are presented in Tab. I. The plot segment of desired position of the shoulder roll joint and corresponding contact observer signals are shown in Fig. 6. Due to the random set-points being commanded to the robot during the free motion -to show that the method can be applied for movement of arbitrary position and speed, the compliance of the joints may not be obvious from the plot. However, looking closely one sees that at the moments of contacts (indicated with dashed blue lines) the desired joint position is being altered taking into account the direction in which the contact is applied.

No special handling of the contact/no-contact transition was done in aforementioned experiments, although the control method would benefit from meticulous handling of such transition, which is left for the future work. Extended presentation of the results of both experiments is included in the video accompanying this paper.

VI. CONCLUSION

We have presented the integration of a joint position tracking discrepancy based contact observer as a feedback signal into the compliant motion controller. The theoretical and implementation details of contact observer robustification and compliance controller design have been described. Finally, we demonstrated the performance of the proposed method in pHRI experiments with humanoid robot Pepper.

Future work will focus on enhancing the overall contact observer and compliance behaviors to the entire robot body including also an omnidirectional mobile base or a floating base of a humanoid and provide it as an additional functionality to the commercial package of Pepper.

3 Fig. 3 :

 33 Fig. 3: Mean validation RMSE for various ensemble sizes and learning rates.

Fig. 4 :

 4 Fig. 4: Optimal max individual tree depth selection from false positive rate (top) analysis and training/test error (bottom).

Fig. 5 :

 5 Fig. 5: Robot stops moving when in contact. The dotted blue lines show start of the contacts, red dashed lines show ±δ.

Fig. 6 :

 6 Fig.6: Robot moves joints to random set-points. When contact is detected, joint position is set according to Eq. 1.

TABLE I :

 I Gain values for main left arm joints.

SoftBank Robotics Europe, Paris, France

University of Montpellier-CNRS LIRMM, Interactive Digital Humans, Montpellier, France

x i being the i th element of vector x.

Note that qtarget is the final target robot posture, whereas q d is intermediate desired joints position command which iteratively and incrementally brings a robot from its initial state to qtarget.


P Jaque Gonzalez

D Dubuc #

K Grenier #
email: 3grenier@laas.fr

T Veronese

Millifluidic Sensor Dedicated to the Microwave Dielectric Spectroscopy of Liquids

Keywords: microwave spectroscopy, coplanar waveguides, biosensors, millifluidic, complex permittivity, liquids

This paper presents a microwave sensor, which has been designed to dielectrically characterize liquids in the frequency range of 0.1 to 4GHz for a fluid volume of the order of 1 ml. Accuracy of the technique has been evaluated for two case studies. One is based on BSA molecules in aqueous solution at different concentrations, whereas the second case study is related to the dielectric measurements of eggs constituents. Even if these different liquids may present intrinsic variability due to the living, a weak standard deviation of 0.01 for BSA solutions and 0.05 for extra-embryonic liquids have been reached respectively. This accuracy enables the future use of the sensor and measurement technique in various applications and notably for the agriculture of precision, which mostly requires ubiquitous and affordable analytic equipment.

I. INTRODUCTION

The microwave spectroscopy is a non-destructive technique, which enables to obtain the dielectric properties of a solid, liquid or gaseous material. This technique has been developed of many applicative contexts such as soil moisture measurement or pipeline detection because of its high sensitivity to water, or to verify the quality of materials in agriculture, food and others [START_REF] Kupfer | Elektromagnetische Aquametrie (Electromagnetic Aquametry)[END_REF]. More recently, the technique has also made possible to study molecules and biological cells due to their dielectric properties in the microwave range, associated with the use of micrometric and microfluidic technologies [START_REF] Grenier | Integrated broadband microwave and microfluidic sensor dedicated to bioengineering[END_REF]- [START_REF] Zedek | Microwave permittivity extraction of individual biological cells submitted to different stimuli[END_REF]. This scale is indeed compatible with fluids available in small quantities for performing tests and analysis. Nevertheless, applications exist for which the quantity of liquid to be analysed does not constitute a limiting factor with up to millilitres available for analysis. It is consequently interesting to develop new analytical tools, which enables the study of the dielectric properties of such liquids, while relaxing the technological constraints, the cost and finally ubiquitous.

This paper presents a microwave device, which enables to perform the dielectric spectroscopy of millilitres' liquids. This requires the combination of both microwave sensing capability with millifluidic technologies. This makes possible to broaden the sensing capabilities of the microwave spectroscopy and its applications, while using a low cost and less complicated technology for microdevices as well as a transportable Vector Network Analyzer (VNA). The developed technique makes it suitable for various exploitation areas such as health, biology, chemistry or agriculture.

As far as chemistry is concerned, the quantification of molecules in aqueous solution has already been demonstrated at the microscale. Regarding the agriculture of precision, there are strong interests in the development of analytic techniques, which could provide quantitative parameters on the quality of biological entities. Freshness of fishes, fruits or eggs [START_REF] Trabelsi | Sensing Egg Quality during Storage by Radiofrequency Complex Permittivity Measurement[END_REF]- [START_REF] Ates | The freshness analysis of an apple and a potato using dielectric properties at the microwave frequency region[END_REF] have already been investigated with various microwave-based sensing techniques, as well as the moisture content in cereal grains or meat [START_REF] Nelson | Using cereal grain permittivity for quality sensing by moisture determination[END_REF], [START_REF] Ji | Study of water-injected pork and pork freshness using microwave techniques[END_REF].

Pertinent examples of chemical and biological samples will consequently be employed for the evaluation of the developed microwave and millifluidic sensor.

This paper first presents the architecture and the fabrication of the microwave and millifluidic sensor. The next section is then dedicated to the modelling of the structure and the extraction of the complex permittivity of the liquid under test. Finally, 2 case studies will be investigated in the fourth section, in order to evaluate the sensitivity of the technique depending on the type of liquids and the associated applicative domain.

II. ARCHITECTURE AND FABRICATION OF THE MICROWAVE

AND MILLIFLUIDIC SENSOR The millifluidic and microwave sensor has been designed to extract the complex permittivity of a liquid between 0.1 to 4 GHz. As shown in Fig. 1, it is composed of two parallel coplanar waveguides with identical electrical and dimensional characteristics, except length. The characteristic impedance is Z0 = 50Ω, the widths of the central conductor and the coplanar slot correspond to w=0.4 mm and s= 0.2 mm respectively. The coplanar waveguides only differ from their length, which are L1 = 43 mm and L2 = 56 mm respectively.

The microwave structure is fabricated on the Rogers RO4003 substrate with a thickness h of 0.8 mm and dielectric constant εr=3.38 and dissipation factor tanδ=0.0027. Above these two waveguides are placed two fluidic channels also exhibiting different lengths of 13 and 26 mm respectively. These channels are made of PDMS, a traditional material for fluidic experiments. These channels enables to precisely localize the fluids to analyze on the coplanar lines (see Fig. 1). A block of PEEK material is then placed on the PDMS layer and allows the liquid injection within the fluidics channels. The system is screwed to an acrylic support placed under the device in order to maintain by pressure the adhesion of the fluidic channels on the coplanar lines (see Figure 2). The required volume of liquids corresponds to 1.3 ml. Next section is related to the modelling of this device to extract the complex permittivity of the liquid under test.

III. MODELLING

The modelling method is based on the one developed by Bianco and Parodi [START_REF] Bianco | Determination of the propagation constant of uniform microstrip lines[END_REF], and next applied by Reynosos-Hernandez et al. [START_REF] Reynoso-Hernández | An improved method for the wave propagation constant γ estimation in broadband uniform millimeterwave transmission line[END_REF] using two micromachined lines transmission with millimetric lengths. This L-L method (for Line-Line) allows extracting the propagation constant from the transition matrix (T) of two nonreflecting lines exhibiting the same characteristic impedance, but different lengths, as shown in Figure 3. The two-scattering matrix of the input and ouput ports Ti and To correspond to the connections between the lines and the cables connected to the VNA. Based on the schematic of Fig. 3, the matrix of the complete structure can be written as:

𝑇 1 = 𝑇 𝑖 𝑇 𝐿 1 𝑇 𝑜 (1) And 𝑇 2 = 𝑇 𝑖 𝑇 𝐿 2 𝑇 𝑜 , (2)
Where 𝑇 𝐿 1 and 𝑇 𝐿 2 is the transition matrix of the two lines. After algebraic manipulation of equations ( 1) and ( 2), the products 𝑇 1 𝑇 2 -1 and 𝑇 2 𝑇 1 -1 are given respectively by:

𝑇 1 𝑇 2 -1 = 𝑇 𝑜 -1 𝑇 𝐿 2 -1 𝑇 𝐿 1 𝑇 𝑜 = 𝑇 𝑖 𝑇 𝐿 1 𝑇 𝐿 2 -1 𝑇 𝑖 -1
(3)

𝑇 2 𝑇 1 -1 = 𝑇 𝑜 -1 𝑇 𝐿 1 -1 𝑇 𝐿 2 𝑇 𝑜 = 𝑇 𝑖 𝑇 𝐿 2 𝑇 𝐿 1 -1 𝑇 𝑖 -1 . ( 4
)
The addition of equations ( 3) and ( 4) is expressed as

[𝑇 1 𝑇 2 -1 + 𝑇 2 𝑇 1 -1 ] ⏞ 𝐴 = 𝑇 𝑜 -1 [𝑇 𝐿 2 -1 𝑇 𝐿 1 + 𝑇 𝐿 1 -1 𝑇 𝐿 2 ] ⏞ 𝐵 𝑇 𝑜 = 𝑇 𝑖 -1 [𝑇 𝐿 1 𝑇 𝐿 2 -1 + 𝑇 𝐿 2 𝑇 𝐿 1 -1 ] ⏟ 𝐶 𝑇 𝑖 (5)
Due to equation ( 5), the terms A, B and C present similar matrices, with the same trace.

𝑡𝑟(𝐴) = 𝑡𝑟(𝐵) = 𝑡𝑟(𝐶) (6)
If we consider the adapted case with no transmission loss for both L1 and L2 lines, their transition matrix of the Li can write as:

𝑇 𝐿 𝑖=1,2 = ( 𝑒 -𝛾𝑙 𝑖 0 0 𝑒 𝛾𝑙 𝑖 ) , ( 7
)
Where γ is the propagation constant of the lines and li their length. After manipulating equations ( 6) and ( 7), we can express the propagation constant as a function of transition matrices of the system as shown in equation [START_REF] Redzwan | Study of Maturity Fruit Assessment using Permittivity and Microwave Reflectivity Measurements for Quality Classification[END_REF]. Note the dependence of γ as a function of 𝛥𝑙 = 𝑙 2 -𝑙 1 .

𝛾 = 1 ∆𝑙 cosh -1 ( 𝑡𝑟(𝑇 1 𝑇 2 -1 +𝑇 2 𝑇 1 -1 ) 4 ) (8)
On the other hand, the propagation constant may also be written as a function of the frequency f, the light velocity c, and the relative permittivity εr as shown in equation ( 9)

𝛾 = 𝑗𝛽 = 𝑗 2𝜋𝑓√𝜀 𝑟 𝐶 . ( 9
)
Finally, due to equality of equation ( 8) and ( 9), the relative permittivity of the propagation constant may be obtained and is described in equation [START_REF] Nelson | Using cereal grain permittivity for quality sensing by moisture determination[END_REF]:

𝜀 𝑟 = -( 𝛾𝑐 2𝜋𝑓 ) 2 . ( 10
)
This method has then been used to obtain the permittivity of different liquids placed on top of the two coplanar waveguides for two applicative case studies for sensitivity evaluation.

IV. EXPERIMENTAL EVALUATION OF THE SENSING CAPABILITIES FOR TWO CASE STUDIES

A. Molecule in aqueous solution

The first case study corresponds to the detection of bovine serum albumin (BSA) dissolved in water and measured at different concentrations: 10, 50 and 100 mg/ml. Each solution has been measured several times, starting from the low BSA concentration to the highest one. The measurements enable a clear discrimination of the liquids with different concentrations of BSA.

The method described in the precedent section has been applied to the measurements. The obtained results, shown in Fig. 4 and Fig. 5, correspond to the real and imaginary parts respectively of the effective permittivity obtained for the different concentrations of BSA as a function of frequency.

A quasi linear relationship between the effective permittivity parameters may be observed in the Fig. 6, which provides the effective permittivity (real and imaginary parts in blue and red respectively) versus the BSA concentration at a fixed frequency of 3 GHz.

The standard deviation of the effective permittivity values (both for real and imaginary parts) is presented in Table 1 for the different concentrations of BSA. The dispersion was calculated based on 5 measurements for each solution. One may observe a correlation between the permittivity values and the concentration of the measured liquids. A low standard deviation is obtained for the concentrations of 10 and 50 mg/ml, whereas it slightly increases for the 100 mg/ml concentration. This increase can be attributed to a pollution of the fluidic channel due to the viscosity of the liquid at high concentrations. This standard deviation is however still considered reasonable for rapid, ubiquitous and low-cost experimentations.

B. Microwave measurements of egg constituents

The second case study is related to the agriculture domain. The chosen material corresponds the extra embryonic liquid of eggs. A previous study indeed targeted to link the moisture content of the yolk or albumen of stored eggs in order to determine the egg freshness as an indicator of quality [START_REF] Trabelsi | Sensing Egg Quality during Storage by Radiofrequency Complex Permittivity Measurement[END_REF].

Before performing any quality assessment on eggs through microwave spectroscopy, it is mandatory to evaluate the sensitivity of the technique with the targeted samples. Inversely to BSA molecules in aqueous solution, yolk and albumen materials are extracted from various eggs, which are naturally presenting intrinsic dispersions. To do so, different liquids of fertilized duck eggs at 12-day of embryonic development have therefore been measured. The liquid corresponds to the mixing of the extra-embryonic constituents of eggs (removal of the duck embryo and shell egg), while including the membranes, the allantoic, amniotic, and vitellin fluids, and the albumin.

In order to evaluate the repeatability of our measurements of extra-embryonic fluids, 10 samples of different eggs with the same origins and incubation conditions have been analyzed. Due to the high viscosity and strong attachment of egg yolk and albumen to any other materials, a specific cleaning protocol has been developed and applied between each measurement. The extracted effective permittivity was referenced to water as described in equation [START_REF] Ji | Study of water-injected pork and pork freshness using microwave techniques[END_REF]. The standard deviation σ of extra-embryonic fluids measured at 2 GHz is 0.05 for both real and imaginary parts. The differences of Δε values is partially because of the different biochemical characteristics of different eggs, however the standard deviation remains low. These results allow to validate the fluidic device for the evaluation of extra-embryonic fluids, and constitute the first step toward the analysis of the quality of the eggs fluids.

∆𝜀 𝑒𝑓𝑓 = 𝜀 𝑒𝑓𝑓 𝐿𝑖𝑞𝑢𝑖𝑑 -𝜀 𝑒𝑓𝑓 𝑤𝑎𝑡𝑒𝑟 (11)

V. CONCLUSIONS

The microwave and millifluidic sensor is presented. A modeling technique based on the L-L method permits the dielectric parameters of the structure loaded with different liquids on the frequency range of 0.1 to 4 GHz.

Its sensitivity has been first evaluated with BSA molecules in aqueous solutions with different concentrations. A weak standard deviation of 0.01 on the real part of the permittivity for low concentration has been reached, which validates the developed device and modeling.

These results have been completed with the measurements of extra-embryonic fluids of 10 fertilized duck eggs at 12-day of embryonic development. Even if such a living material presents an intrinsic variability with high viscosity and attachment to any surface, the millifluidic and microwave sensor only exhibits a standard deviation of 0.05 on the real part of permittivity contrast. This result allows us to validate the device and the methodology with the measurement of the dielectric properties of milliliter liquids using a low cost and affordable technology, compatible with many application areas. It notably opens the door towards further investigations for agricultural quality assessment.

Fig. 1 .

 1 Fig. 1. Schematic of the millifluidic and microwave sensor, which includes 2 coplanar waveguides of different lengths with 2 PDMS fluidic channels on top.

Fig. 2 .

 2 Fig. 2. Photograph of the millifluidic and microwave device, seen from above, with coplanar waveguides, SMA connectors, PEEK block and inlet and outlet tubes.

Fig. 3 .

 3 Fig. 3. Schematic of the two transmission lines with the input and output access transmission matrix for the L-L method.

Fig. 4 .

 4 Fig. 4. Real part of the effective permittivity of the structure when loaded with liquids presenting different concentrations of BSA.

Fig. 5 .

 5 Fig. 5. Imaginary part of the effective permittivity of the structure when loaded with liquids presenting for different concentrations of BSA.

Fig. 6 .

 6 Fig. 6. Real and imaginary parts of the effective permittivity of the structure versus BSA concentration at 3 GHz.

  The real and imaginary parts of the effective permittivity contrast (Δε) of 10 extra-embryonic fluids samples are shown in Fig.7and Fig.8respectively.

Fig. 7 .

 7 Fig. 7. Contrast of the real part of the effective permittivity of 10 extraembryonic samples fluids referenced to water at 2 GHz.

Fig. 8 .

 8 Fig. 8. Contrast of the imaginary part of the effective permittivity of 10 extraembryonic samples fluids referenced to water at 2 GHz. Curves are quite close to each other. For further visibility, the Cole-Cole plot is given in Fig. 9 at 2 GHz.

Fig. 9 .

 9 Fig. 9. Cole-Cole plot for 10 extra-embryonic fluids referenced to water at 2 GHz.

Table 1 .

 1 Standard deviation of the real and imaginary parts of effective permittivity for different BSA concentrations at 3 GHz.

	BSA	σ ε'	σ ε''
	Concentrations
	10 mg/ml	0.02	0.003
	50 mg/ml	0.01	0.003
	100 mg/ml	0.057


A Le Coënt

L Fribourg
email: fribourg@lsv.fr

Guaranteed optimal reachability control of reaction-diffusion equations using one-sided Lipschitz constants and model reduction

We show that, for any spatially discretized system of reactiondiffusion, the approximate solution given by the explicit Euler timediscretization scheme converges to the exact time-continuous solution, provided that diffusion coefficient be sufficiently large. By "sufficiently large", we mean that the diffusion coefficient value makes the one-sided Lipschitz constant of the reaction-diffusion system negative. We apply this result to solve a finite horizon control problem for a 1D reactiondiffusion example. We also explain how to perform model reduction in order to improve the efficiency of the method.

Introduction

Guaranteed reachability analysis

Given a system of Ordinary Differential equations (ODEs) of dimension n satisfying standard conditions of existence and uniqueness of the solution, the area of Numerical Analysis makes use of numerical tools in order to compute the approximate value of the solution, starting at an initial point of R n , with high accuracy: 1st order methods (explicit/implicit Euler method, trapezoid rule), higher-order Runge-Kutta methods, etc. In contrast, the area of Guaranteed (or Symbolic) Analysis is devoted to the construction of an overapproximation of the set of solutions that start, not at a single point of R n , but from a dense compact set of initial points. Guaranteed analysis, in its modern form, has been initiated in the 60's by R.E. Moore and his creation of Interval Arithmetic [START_REF] Moore | Interval Analysis[END_REF]: the set of solutions (or trajectories) are overapproximated by a sequence of "rectangular sets", i.e., cross-product of intervals of R. A set of arithmetic and differential calculus has been created for manipulating such sets. An overapproximation of the set of trajectories is computed using a Taylor development up to some order and an overestimation of the "Lagrange remainder". The method has been considerably refined in the 90's [START_REF] Berz | Computation and application of taylor polynomials with interval remainder bounds[END_REF][START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF][START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF][START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF][START_REF] Nedialko | Interval arithmetic, affine arithmetic, taylor series methods: Why, what next?[END_REF]. These recent techniques make use of different convex data structures such as parallelepipeds [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF] or zonotopes [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF][START_REF] Kühn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF] instead of rectangular sets in order to enclose the flow of ODEs.

Such methods are typically applied to the formal proof of correctness of ODE integration, and more generally, to guarantee that the solutions of the ODEs satisfy some desired properties. Guaranteed reachability analysis generally treats linear systems. Extensions to nonlinear systems have been proposed, e.g., in [4], using local linearizations (see also [START_REF] Mitchell | Validating a hamiltonjacobi approximation to hybrid system reachable sets[END_REF][START_REF] Mitchell | Overapproximating reachable sets by hamiltonjacobi projections[END_REF]).

Guaranteed optimal control

In presence of inputs, we can use guaranteed analysis to describe a law that allows the system to satisfy a desired property. This corresponds to the topic of guaranteed (or correct-by-design) control synthesis. Several works have reecently applied guaranteed analysis to optimal control synthesis. Thus, in [START_REF] Schürmann | Optimal control of sets of solutions to formally guarantee constraints of disturbed linear systems[END_REF][START_REF] Schürmann | Reachset model predictive control for disturbed nonlinear systems[END_REF], the authors focus on a (finite time-horizon) optimal control procedure with a formal guarantee of safety constraint satisfaction, using zonotopes as state set representations. In [START_REF] Silva | Synthesis of safe controllers for nonlinear systems using dynamic programming techniques[END_REF], the authors focus on (periodically) sampled systems, and perform reachability analysis using convex polytopes as state set representations. In [START_REF] Kapela | A lohner-type algorithm for control systems and ordinary differential inclusions[END_REF][START_REF] Maidens | Reachability analysis of nonlinear systems using matrix measures[END_REF][START_REF] Fan | Simulation-driven reachability using matrix measures[END_REF][START_REF] Reissig | Symbolic optimal control[END_REF][START_REF] Rungger | Arbitrarily precise abstractions for optimal controller synthesis[END_REF], the authors construct an over-approximation of the set of trajectories using a growth bound (bounding the distance of neighboring trajectories) exploiting the notion of one-sided Lipschitz constant (also called "logarithmic norm" or "matrix norm"). The notion of "one-sided Lipschitz (OSL) constant" has been introduced independently by Dahlquist [START_REF] Dahlquist | Stability and error bounds in the numerical integration of ordinary differential equations[END_REF] and Lozinskii [START_REF] Mikhailovich | Error estimate for numerical integration of ordinary differential equations. i[END_REF] in order to derive error bounds in initial value problems (see survey in [START_REF] Söderlind | The logarithmic norm. history and modern theory[END_REF]). We used ourselves OSL constants in the context of symbolic optimal control in [START_REF] Le | Guaranteed control of sampled switched systems using semi-Lagrangian schemes and one-sided Lipschitz constants[END_REF]. The main difference with previous work [START_REF] Kapela | A lohner-type algorithm for control systems and ordinary differential inclusions[END_REF][START_REF] Maidens | Reachability analysis of nonlinear systems using matrix measures[END_REF][START_REF] Fan | Simulation-driven reachability using matrix measures[END_REF][START_REF] Reissig | Symbolic optimal control[END_REF][START_REF] Rungger | Arbitrarily precise abstractions for optimal controller synthesis[END_REF] is that our method makes use of explicit Euler's algorithm for ODE integration (cf. [START_REF] Le Coënt | Distributed control synthesis using Euler's method[END_REF][START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]) instead of sophisticated algorithms such as Lohner's algorithm [START_REF] Kapela | A lohner-type algorithm for control systems and ordinary differential inclusions[END_REF] or interval Taylor series methods [START_REF] Nedialko | Interval arithmetic, affine arithmetic, taylor series methods: Why, what next?[END_REF]. This leads us to a simple implementation of just a few hundred lines of Octave (see [START_REF] Le | OSLator 1.0[END_REF]).

As explained in [START_REF] Saluzzi | Error estimates for a tree structure algorithm solving finite horizon control problems[END_REF], using the Dynamic Programming (DP) [START_REF] Bellman | Dynamic Programming[END_REF] one can approximate the "value" of the solution of Hamilton-Jacobi-Bellman (HJB) equations. In [START_REF] Falcone | An approximation scheme for evolutive hamilton-jacobi equations[END_REF][START_REF] Saluzzi | Error estimates for a tree structure algorithm solving finite horizon control problems[END_REF], the authors thus show how to use finite difference schemes, Euler time integration and DP for solving finite horizon control problems. Furthermore, they give a priori errors estimates which are first-order in the size ∆t of the time discretization step; however, the error involves a constant CpT q which depends exponentially on the length T of the finite horizon3 . We solve here finite horizon control problems along the same lines (using finite difference, explicit Euler and DP) but, under the hypothesis of OSL negativity (see section 1.3), we obtain an error upper bound that is linear in T (see Section 2.4, Theorem 2).

Reaction-diffusion equations

It is natural to adapt the optimal control methods of ODEs to the control of Partial Differential Equations (PDEs). This can be done by transforming the PDE into (a vast system of) ODEs, using space discretization techniques such as finite difference or finite element methods. In the present work, we focus on a particular class of non-linear PDEs called "reaction-diffusion" equations. Reactiondiffusion equations cover a variety of particular cases with important applications in mathematical physics, and in biological models such as the Schlögl model or the FitzHugh-Nagumo system [START_REF] Casas | Optimal control of a class of reaction-diffusion systems[END_REF]. The problem of optimal control of reaction-diffusion equations has been recently the topic of many works of (classical) numerical analysis: see, e.g., [START_REF] Barthel | Optimal boundary control of a system of reaction diffusion equations[END_REF][START_REF] Court | Hybrid optimal control problems for a class of semilinear parabolic equations[END_REF][START_REF] Finotti | Optimal control of advective direction in reaction-diffusion population models[END_REF][START_REF] Griesse | A primal-dual active set strategy for optimal boundary control of a nonlinear reaction-diffusion system[END_REF][START_REF] Moura | Optimal boundary control & estimation of diffusion-reaction pdes[END_REF][START_REF] Moura | Optimal boundary control of reactiondiffusion partial differential equations via weak variations[END_REF].

The notion OSL constant can be naturally extended to PDEs and reactiondiffusion equations in particular, as shown in [START_REF] Arcak | Certifying spatially uniform behavior in reaction-diffusion pde and compartmental ode systems[END_REF][START_REF] Aminzare | Logarithmic lipschitz norms and diffusioninduced instability[END_REF][START_REF] Aminzare | Guaranteeing spatial uniformity in reaction-diffusion systems using weighted l 2 norm contractions[END_REF][START_REF] Aminzare | Some remarks on spatial uniformity of solutions of reaction-diffusion pdes[END_REF]. In these works, the authors focus on the case where the OSL constant associated with the reactiondiffusion equation is negative. In this case, the system has a contractivity (or "incremental stability") property which expresses the fact that all solutions converge exponentially to each other (see [START_REF] Sontag | Contractive systems with inputs[END_REF]).

In this work, we also study reaction-diffusion equations with negative OSL constants, but the equations are equipped with control inputs, and the problem of controlling these inputs in an optimal way is here considered.

Model reduction

In order to reduce the large dimension of ODE systems originating from the PDE space discretization, Model Order Reduction (MOR) techniques are often used in conjunction with the analysis of ODE systems. The idea is to first infer the optimal control at a reduced level, then apply it at the original level. In the field of guaranteed analysis, the MOR technique of "balanced truncation" was used to treat linear systems (e.g., [3,[START_REF] Han | Reachability analysis of hybrid control systems using reduced-order models[END_REF][START_REF] Han | Reachability analysis of large-scale affine systems using low-dimensional polytopes[END_REF][START_REF] Le Coënt | Guaranteed control synthesis of switched control systems using model order reduction and state-space bisection[END_REF]). In [START_REF] Kalise | Reduced-order minimum time control of advectionreaction-diffusion systems via dynamic programming[END_REF], a MOR technique based on spectral element method was coupled to an HJB approach for application to advection-reaction-diffusion systems (cf. [START_REF] Kalise | Polynomial approximation of high-dimensional hamilton-jacobi-bellman equations and applications to feedback control of semilinear parabolic pdes[END_REF] for application to semilinear parabolic PDEs). The MOR technique of "Proper Orthogonal Decomposition (POD)" was coupled to an HJB approach in [START_REF] Alla | Error analysis for POD approximations of infinite horizon problems via the dynamic programming approach[END_REF][START_REF] Alla | A HJB-POD approach for the control of nonlinear PDEs on a tree structure[END_REF][START_REF] Kunisch | Hjb-pod-based feedback design for the optimal control of evolution problems[END_REF]. Here, we couple our HJBbased method to a simple ad hoc reduction method (see Section 2.5).

The plan of the paper is as follows: We explain how to convert the reactiondiffusion equation into a system of ODEs by domain discretization in Section 2.1, and how to approximate the solution of the latter system using the explicit Euler scheme of time integration in Section 2.2. Our procedure for solving finite horizon control problems is explained in Section 2.3. In Section 2.4, we give an upper bound to the error between the approximate value thus computed and the exact optimal value. In Section 2.5, we explain how to perform MOR in order to treat systems of larger dimension. We conclude in Section 3.

Optimal Reachability Control of Reaction-Diffusion Equations

Let us consider the special class of PDEs called "reaction-diffusion" equations.

For the sake of notation simplicity, we focus on 1D reaction-diffusion equations with Dirichlet boundary conditions (the domain Ω is of the form r0, Ls Ă R), but the method applies to 2D or 3D reaction-diffusion equations with other boundary conditions. A 1D reaction-diffusion system with Dirichlet boundary conditions is of the form:

Bypt, xq Bt " σ B 2 ypt, xq Bx 2 `f pypt, xqq, t P r0, T s, x P Ω " r0, Ls. ypt, 0q " u 0 ptq, ypt, Lq " u L ptq, t P r0, T s, yp0, xq " y 0 pxq, x P Ω " r0, Ls.
Here, y " ypt, xq is an R-valued unknown function, Ω is a bounded domain in R with boundary BΩ :" t0, Lu, and f is a function from r0, T s ˆΩ to r0, 1s. Also y 0 pxq is a given function called "initial condition", and σ a positive constant, called "diffusion constant".

The boundary control up¨q :" pu 0 p¨q, u L p¨qq that we consider here, is a piecewise constant (or "staircase") function from r0, T s to a finite set U Ă r0, 1s r0, 1s.The control uptq changes its value periodically at t " τ, 2τ, . . . . We assume that T " kτ for some positive integer k. The constant τ is called the "switching (or sampling) period".

Given an initial condition y 0 p¨q such that y 0 pxq P r0, 1s for all x P r0, Ls, we assume that, for any boundary control up¨q, the solution yp¨, ¨q of the system exists, is unique, and ypt, xq P r0, 1s for all pt, xq P r0, T s ˆr0, Ls.

Domain discretization

A well-known approach in numerical analysis of PDEs (see, e.g., [START_REF] Koto | IMEX Runge-Kutta schemes for reaction-diffusion equations[END_REF]) is to discretize in space by finite difference or finite element methods in order to transform the PDE into a system of ODEs.

Let M be a positive integer, h " L{pM `1q, and let Ω h be a uniform grid with nodes x j " jh, j " 1, . . . , M . By replacing the 2nd order spatial derivative with the second order centered difference, we obtain a space-discrete approximation:

dy dt " σL h y `σϕ h pt, uq `f pt, yq,
with yptq " ry 1 ptq, . . . , y M ptqs T , y j ptq « ypt, x j q, and

L h " 1 h 2 » - - - - - ´2 1 0 ¨¨¨0 1 ´2 1 ¨¨¨0 0 1 ´2 ¨¨¨0 ¨¨0 0 ¨¨¨1 ´2fi ffi ffi ffi ffi fl ϕ h pt, uq " 1 h 2 ru 0 ptq, 0, . . . , 0, u L ptqs J .
The point yptq, often abbreviated as y, is thus an element of S " r0, 1s M .

Explicit Euler time integration

Let us abbreviate the equation dy dt " σL h y `σϕ h pt, uq `f pt, yq by: dy dt " f u pt, yq.

We denote by Y u t,y0 , the solution y of the system at time t P r0, τ q controlled by mode u P U , for initial condition y 0 . Given a sequence of modes (or "pattern") π :" u k ¨¨¨u 1 P U k , we denote by Y π t,y0 the solution of the system for mode u k on t P r0, τ q with initial condition y 0 , extended continuously with the solution of the system for mode u k´1 on t P rτ, 2τ q, and so on iteratively until mode u 1 on t P rpk ´1qτ, kτ s.

Let us now approximate the solution of the system by performing time integration with the explicit Euler scheme. This yields:

y n`1 " y n `τ f u pt n , y n q,
Here y n is an approximate value of ypt n q. Given a starting point z P X and a mode u P U , we denote by Ỹ u t,z the Euler-based image of z at time t via u for t P r0, τ q. We have: Ỹ u t,z :" z `t f u pzq. We denote similarly by Ỹ π t,z the Euler-based image of z via pattern π P U k at time t P r0, kτ s.

Finite horizon control problems

Let us now explain the principle of the method of optimal control of ODEs used in [START_REF] Le | Guaranteed control of sampled switched systems using semi-Lagrangian schemes and one-sided Lipschitz constants[END_REF], in the present context. We consider the cost function: J k : r0, 1s M ˆU k Ñ R ě0 defined by: J k py, πq " }Y π kτ,y ´yf }, where }¨} denotes the Euclidean norm in R M , and y f P r0, 1s M is a given "target" state.

We consider the value function v k : r0, 1s M Ñ R ě0 defined by:

v k pyq :" min πPU k tJ k py, πqu " min πPU k t}Y π kτ,y ´yf }u.
Given k P N and τ P R ą0 , we consider the following finite time horizon optimal control problem: Find for each y P r0, 1s M the value v k pyq, i.e. min πPU k t}Y π kτ,y ´yf }u,

and an optimal pattern:

π k pyq :" arg min πPU k t}Y π kτ,y ´yf }u.
In order to solve such optimal control problems, a classical "direct" method consists in spatially discretizing the state space S " r0, 1s M (i.e., the space of values of y). We consider here a uniform partition of S into a finite number N of cells of equal size: in our case , this means that interval r0, 1s is divided into K subintervals of equal size, and N " K M . A cell thus corresponds to a M -tuple of subintervals. The center of a cell coresponds to the M -tuple of the subinterval midpoints. The associated grid X is the set of centers of the cells of S. The center z P X of a cell C is considered as the ε-representative of all the points of C. We suppose that the cell size is such that }y ´z} ď ε, for all y P C (i.e. K ě ? M {2ε). In this context, the direct method proceeds as follows (cf. [START_REF] Le | Guaranteed control of sampled switched systems using semi-Lagrangian schemes and one-sided Lipschitz constants[END_REF]): we consider the points of X as the vertices of a finite oriented graph; there is a connection from z P X to z 1 P X if z 1 is the ε-representative of the Euler-based image pz `τ f u pzqq of z, for some u P U . We then compute using dynamic programming the "path of length k with minimal cost" starting at z: such a path is a sequence of k `1 connected points z z k z k´1 ¨¨¨z 1 of X which minimizes the distance }z 1 ´yf }. This procedure allows us to compute a pattern π ε k pzq of length k, which approximates the optimal pattern π k pyq.

Definition 1. The function next u : X Ñ X is defined by:

-next u pzq " z 1 , where z 1 is the ε-representative of Ỹ u τ,z .
Definition 2. For all point x P X , the spatially discrete value function v ε k : X Ñ R ě0 is defined by: for k " 0, v ε k pzq " }z ´yf }, for k ě 1, v ε k pzq " min uPU tv ε k´1 pnext u pzqqu.

Definition 3. The approximate optimal pattern of length k associated to z P X , denoted by π ε k pzq P U k , is defined by:

-if k " 0, π ε k pzq " nil, -if k ě 1, π ε k pzq " u k pzq ¨π1 where u k pzq " arg min uPU tv ε k´1 pnext u pzqqu and π 1 " π ε k´1 pz 1 q with z 1 " next u k pzq pzq.
It is easy to construct a procedure P ROC ε k which takes a point z P X as input, and returns an approximate optimal pattern π ε k P U k .

Remark 1. The complexity of P ROC ε k is Opm ˆk ˆN q where m is the number of modes (|U | " m), k the time-horizon length (T " kτ ) and N the number of cells of X (N " K M with K " ? M {2ε).

Error upper bound

Given a point y P S of ε-representative z P X , and a pattern π ε k returned by P ROC ε k pzq, we are now going to show that the distance } Ỹ π ε k kτ,z ´yf } converges to v k pyq as ε Ñ 0. We first consider the ODE: dy dt " f u pyq, and give an upper bound to the error between the exact solution of the ODE and its Euler approximation (see [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]). Definition 4. Let µ be a given positive constant. Let us define, for all u P U and t P r0, τ s, δ u t,µ as follows:

if λ u ă 0 : δ u t,µ " ˆµ2 e λut `C2 u λ 2 u ˆt2 `2t λ u `2 λ 2 u `1 ´eλut ˘˙˙1 2 if λ u " 0 : δ u t,µ " `µ2 e t `C2 u p´t 2 ´2t `2pe t ´1qq ˘1 2 if λ u ą 0 : δ u t,µ " ˆµ2 e 3λut `C2 u 3λ 2 u ˆ´t 2 ´2t 3λ u `2 9λ 2 u `e3λut ´1˘˙˙1 2
where C u and λ u are real constants specific to function f u , defined as follows:

C u " sup yPS L u }f u pyq},
where L u denotes the Lipschitz constant for f u , and λ u is the OSL constant associated to f u , i.e., the minimal constant such that, for all y 1 , y 2 P S:

xf u py 1 q ´fu py 2 q, y 1 ´y2 y ď λ u }y 1 ´y2 } 2 , where x¨, ¨y denotes the scalar product of two vectors of S.

Proposition 1. [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF] Consider the solution Y u t,y0 of dy dt " f u pyq with initial condition y 0 of ε-representative z 0 (hence such that }y 0 ´z0 } ď ε), and the approximate solution Ỹ u t,z0 given by the explicit Euler scheme. For all t P r0, τ s, we have:

}Y u t,y0 ´Ỹ u t,z0 } ď δ u t,ε .
Proposition 2. Consider the system dy dt " f u pyq with f u pyq :" σL h y`σϕ h pt, uqf pyq. For a diffusion coefficient σ ą 0 sufficiently large, the OSL constant λ u associated to f u is such that: λ u ă 0.

Proof. Consider the ODE: dy dt " f u pyq " σL h y `σϕ h pt, uq `f pyq. For all y 1 , y 2 P S, we have: xf py 2 q´f py 1 q, y 2 ´y1 y ď λ f }y 2 ´y1 } 2 , where λ f is the OSL constant of f . Hence: xf u py 2 q ´fu py 1 q, y 2 ´y1 y " xσL h py 2 ´y1 q `f py 2 q ´f py 1 q, y 2 ´y1 y ď py 2 ´y1 q J pσL h `λf qpy 2 ´y1 q.

Since y J L h y ă 0 for all y P S (negativity of the quadratic form associated to L h ), we have: λ u }y 1 ´y2 } 2 ď py 2 ´y1 q T pσL h `λf qpy 2 ´y1 q ă 0, for σ ą 0 sufficiently large. Hence λ u ă 0. l Lemma 1. Consider the system dy dt " f u pyq where the OSL constant λ u associated to f u is negative, and initial error e 0 :" }y 0 ´z0 } ą 0. Let G u :"

? 3e0|λu| Cu
. Consider the (smallest) positive root

α u :" 1 `|λ u |G u {4 ´a1 `pλ u G u {4q 2 of equation: ´1 2 |λ u |G u `p2 `1 2 |λ u |G u qα ´α2 " 0. Suppose: |λu|Gu 4 ă 1.
Then we have 0 ă α u ă 1, and, for all t P r0, τ s with τ ď G u p1 ´αu q: δ u e0 ptq ď e 0 . Proof. See Appendix 1.

Remark 2. In practical case studies |λ u | is often small, and the term pλ u G u {4q 2 can be neglected, leading to α u « |λ u |G u {4 and G u p1 ´αu q « G u p1 ´|λu|Gu

4 q « G u .
Remark 3. It follows that, for τ ď G u p1´α u q, the Euler explicit scheme is stable, in the sense that initial errors are damped out. Remark 4. If τ ą G u p1´α u q, we can make use of subsampling, i.e., decompose τ into a sequence of elementary time steps ∆t with ∆t ď G u p1 ´αu q in order to be still able to apply Lemma 1 (see Example 1). Let us point out that Lemma 1 (and the use of subsampling) allows to ensure set-based reachability with the use of procedure P ROC ε k . Indeed, in this setting, the explicit Euler scheme leads to decreasing errors, and thus, point based computations performed with the center of a cell can be applied to the entire cell.

We suppose henceforth that the system dy dt " f u pyq satisfies:

pHq : λ u ă 0, |λ u |G u 4 ă 1
and τ ď G u p1 ´αu q, for all u P U.

From Proposition 1 and Lemma 1, it easily follows:

Theorem 1. Consider a system dy dt " f u pyq satisfying pHq, and a point y P S of ε-representative z P X . We have:

}Y π t,y ´Ỹ π t,z } ď ε,
for all π P U k and t P r0, kτ s.

Proposition 3. Let z P X and π ε k be the pattern of U k returned by P ROC ε k pzq. For all π P U k , we have:

} Ỹ π ε k kτ,z ´yf } ď } Ỹ π kτ,z ´yf } `2kε.
Proof. W.l.o.g., let us suppose that y f is the origin O. Let us prove by induction on k:

} Ỹ π ε k kτ,z } ď } Ỹ π kτ,z } `2kε.
Let π ε k :" u k ¨¨¨u 1 . The base case k " 1 is easy. For k ě 2, we have:

} Ỹ π ε k kτ,z } " } Ỹ u k´1 ¨¨¨u1 pk´1qτ,z k } with z k " Ỹ u k τ,z with u k " argmin uPU tv ε k´1 pnext u pzqqu ď } Ỹ u k´1 ¨¨¨u1 pk´1qτ,next u k pz k q } `ε ď } Ỹ π 1 pk´1qτ,next u k pz k q } `p2k ´1qε for all π 1 P U k´1 by induction hypothesis, ď } Ỹ π 1
pk´1qτ,z 1 } `2kε for all π 1 P U k´1 and all z 1 P tnext u pzq |u P U u ď } Ỹ π τ,z } `2kε for all π P U k . l Theorem 2. Let y P S be a point of ε-representative z P X . Let π ε k be the pattern returned by P ROC ε k pzq, and π ˚:" argmin πPU k }Y π kτ,y ´yf }. The discretization error E ε pT q :" |} Ỹ π ε k kτ,z ´yf } ´vk pyq|, with v k pyq :" }Y π kτ,y ´yf } and T " kτ , satisfies:

E ε pT q ď p2k `1qε. It follows that } Ỹ π ε k kτ,z ´yf } converges to v k pyq as ε Ñ 0.
Proof. W.l.o.g., let us suppose that y f is the origin O. For all π P U k , we have by Proposition 3 and Theorem 1:

} Ỹ π ε k kτ,z } ď } Ỹ π kτ,z } `2kε ď }Y π kτ,y } `p2k `1qε. Hence } Ỹ π ε k kτ,z } ď min πPU k }Y π kτ,y } `p2k `1qε " }Y π kτ,y } `p2k `1qε.
On the other hand, for all π P U k , it follows from Theorem 1:

}Y π kτ,y } ď }Y π kτ,y } ď } Ỹ π kτ,z } `ε.
Hence:

}Y π kτ,y } ď } Ỹ π ε k kτ,z } `ε.
Therefore we have:

|} Ỹ π ε k kτ,z } ´}Y π kτ,y }| ď p2k `1qε. l
Remark 5. The error bound E ε pT q is thus linear in k " T {τ . In order to decrease k, one can apply consecutively p ě 2 modes in a row (without intermediate ε-approximation); this is equivalent to divide k by p, at the price of considering m p "extended" modes instead of just m modes. (see Example 1, Figure 2). An alternative for decreasing k is to increase τ (which may require in turn to decrease ∆t for preserving assumption ∆t ď G u p1 ´αu q, see Remark 4).

Example 1. Consider the 1D reaction-diffusion system with Dirichlet boundary condition (see [START_REF] Pouchol | Phase portrait control for 1D monostable and bistable reaction-diffusion equations[END_REF], bistable case):

Bypt, xq Bt " σ B 2 ypt, xq
Bx 2 `f pypt, xqq, t P r0, T s, x P r0, Ls ypt, 0q " u 0 , ypt, Lq " u L , yp0, xq " y 0 pxq, x P r0, Ls with σ " 1, L " 4 and f pyq " yp1 ´yqpy ´θq with θ " 0.3. The control switching period is τ " 0.1. The values of the boundary control u " pu 0 , u L q are in U " tp0, 0q, p0.2, 0.2q, p0.4, 0.4q, p0.6, 0.6q, p0.8.0.8q, p1, 1qu. 4We discretize the domain Ω " r0, Ls of the system with M 1 " 5 discrete points, using a finite difference scheme. Our program returns an OSL constant λ u " ´0.322 for all u P U . Constant C u varies between 10.33 and 11.85 depending on the values of u.

We then discretize each interval component of the space S " r0, 1s M1 of values of y into 15 points with spacing η " 1{15 « 0.066. The grid X is of the form t0, η, 2η, . . . , 15ηu M1 , and the initial error e 0 equal to ε " ? M 1 η{2. This leads to G u varying between 0.00155 and 0.00178 depending on the value of u P U . One checks: |λu|Gu 4 ă 1 for all u P U . The time step upper bound required by Theorem 1 for ensuring numeric stability is 0.00155. Since the switching period is τ " 0.1, we perform subsampling (see, e.g., [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]) by decomposing every time step riτ, pi `1qτ q (1 ď i ď k ´1) into a sequence of elementary Euler steps of length ∆t " τ {100 ă 0.00155. This ensures that the system satisfies pHq, hence, by Theorem 1, the explicit Euler scheme is stable and error }Y π t,y0 ´Ỹ π t,z0 } never exceeds ε.

For objective with y f " p0.3, 0.3, 0.3, 0.3, 0.3q and horizon time T " kτ " 2 (i.e., k " 20), our program5 returns an approximate optimal controller in 2 minutes. Let z 0 be the ε-representative of y 0 " 0.8x{L `0.1p1 ´x{Lq. Let π ε k be the pattern output by P ROC ε k pz 0 q. A simulation of zptq :" Ỹ π ε k t,z0 is given in Figure 1 with T " 2, τ " 0.1 (k " 20), ∆t " τ 100 . We have }zpT q ´yf } « 0.276. The simulation presents some similarity with simulations displayed in [START_REF] Pouchol | Phase portrait control for 1D monostable and bistable reaction-diffusion equations[END_REF] (see, e.g., lower part of Figure 6), with a phase control u 0 " u L ą θ (here, u 0 " u L " 0.4) alternating with a phase control u 0 " u L ă θ (here, u 0 " u L " 0.2). The discretization error E ε pT q is smaller than p2k `1qε " 41 ? 5{30 ă 3.1. Let us now proceed with extended modes of length p " 2 and p " 4, as explained in Remark 5. For p " 2 (i.e., k " 10), the control is synthesized in 7mn of CPU time. The controller simulation is given in the left part of Figure 2; we have: }zpT q ´yf } « 0.445 with E ε pT q ă 1.57. For p " 4 (i.e., k " 5), the computation of the control requires 8h of CPU time. The corresponding simulation is given in the right part of Figure 2; we now have: }zpT q´y f } « 0.164 with E ε pT q ă 0.82.

Model reduction

Let us consider the system S 2 on space S h2 " r0, 1s M2 (with M 2 even). The differential equation can be written under the form:

dy 2 dt
" σL h2 y 2 `ϕh2 puq `f py 2 q.

where L h2 corresponds to the pM 2 ˆM2 q Laplacian matrix, and h 2 " L M2`1 . Let us consider the "reduced" system S 1 defined on S h1 " r0, 1s M1 with M 1 " M 2 {2, defined by:

dy 1 dt " σL h1 y 1 `ϕh1 puq `f py 1 q,
where L h1 is the pM 1 ˆM1 q Laplacian matrix and h 1 " L M1`1 .

With M 1 " M 2 {2, we have h 2 " L 2M1`1 (" h1pM1`1q 2M1`1 ). Let us consider the pM 1 ˆM2 q reduction matrix:

Π :" 1 ? 2 » - - - 1 1 0 ¨¨¨0 0 0 0 1 1 ¨¨¨0 ¨¨0 0 ¨¨¨0 1 1 fi ffi ffi fl
Note that ΠΠ J " I M . Let us consider a point w 0 P S h2 , and let z 0 " Πw 0 P S h1 . Theorem 3. Consider the system S 2 and a point w 0 P S h2 , and let z 0 " Πw 0 P S h1 . Let Y h2 w0 and Y h1 z0 be the solutions of S 2 and S 1 with initial conditions w 0 P S h2 and z 0 P S h1 respectively. We have:

@t ě 0 }ΠY h2 w0 ptq ´Y h1 z0 ptq} ď K 2 σ |λ h1 | , where K 2 :" sup wPS h 2 }pΠL h2 ´Lh1 Πqw},
and L h2 (resp. L h1 ) is the Laplacian matrix of size M 2 ˆM2 (resp. M 1 ˆM1 ).

Proof. Let us consider the system S 2 : dy 2 dt " σL h2 y 2 `ϕh2 puq `f py 2 q.

By application of the projection matrix Π, we get: dΠy 2 dt " σΠL h2 y 2 `ϕh1 puq `f pΠy 2 q.

By substracting pairwise with the sides of S 1 , we have:

dΠy 2 dt ´dy 1 dt
" σpΠL h2 y 2 ´Lh1 y 1 q `f pΠy 2 q ´f py 1 q " F h1 pΠy 2 q ´Fh1 py 1 q `σpΠL h2 ´Lh1 Πqy 2 ,

where F h1 pyq " σL h1 pyq `f pyq for y P S h1 . On the other hand, we have:

1 2 d dt p}Πy 2 ´y1 } 2 q " x d dt pΠy 2 ´y1 q, Πy 2 ´y1 y " xF h1 pΠy 2 q ´Fh1 py 1 q `σpΠL h2 ´Lh1 Πqy 2 , Πy 2 ´y1 y " xF h1 pΠy 2 q ´Fh1 py 1 q, Πy 2 ´y1 y `σxpΠL h2 ´Lh1 Πqy 2 , Πy 2 ´y1 y ď λ h1 }Πy 2 ´y1 } 2 `σxpΠL h2 ´Lh1 Πqy 2 , Πy 2 ´y1 y ď λ h1 }Πy 2 ´y1 } 2 `K2 σ}Πy 2 ´y1 } with K 2 :" sup wPS h 2 }pΠL h2 ´Lh1 Πqw} ď λ h1 }Πy 2 ´y1 } 2 `K2 σ 1 2 pα}Πy 2 ´y1 } 2 `1 α q, for all α ą 0. Choosing α ą 0 such that K 2 σα " ´λh1 , i.e.: α " ´λh 1 K2σ , we have:

1 2 d dt p}Πy 2 ´y1 } 2 q ď λ h1 2 }Πy 2 ´y1 } 2 ´pK 2 σq 2 2λ h1 .
Since y 2 p0q " w 0 and y 1 p0q " z 0 , we get by integration: }Πy 2 ptq ´y1 ptq} when the same control modes are applied to both systems. 6Let y 0 2 P S 2 and y f 2 P S 2 be an initial and objective point respectively. Let y 0 1 :" Πy 0 2 P S 1 and y f 1 :" Πy f 2 P S 1 denote their projections. Suppose that π ε is the pattern returned by P ROC ε k py 0 1 q for the reduced system S 1 . Then, from Theorem 3, it follows that, when the same control π ε is applied to the original system S 2 with y 2 p0q " y 0 2 P S 2 , it makes the projection Πy π ε 2 ptq P S 1 reach a neighborhood of y f 1 at time t " T . Formally, we have:

}P y π ε 2 pT q ´yf 1 } ď }y π ε 1 pT q ´yf 1 } `K2 σ |λ h1 | .
Example 2. Let us take the system defined in Example 1 as reduced system S 1 (M 1 " 5), and let us take as "full-size" system S 2 the system corresponding to M 2 " 10. Since the size of the grid X 2 associated to S 2 is exponential in M 2 , the size X 2 is multiplied by p1{ηq M2´M1 " 15 5 « 7.6 ¨10 5 w.r.t. the size of the grid X 1 associated to S 1 . The complexity for synthesizing directly the optimal control of S 2 thus becomes intractable. On the other hand, if we apply to S 2 the optimal strategy π ε P U k found for S 1 in Example 1, we obtain a simulation depicted in Figure 3 for extended mode of length 1, which is the counterpart of Figure 1 with M 2 " 10 (instead of M 1 " 5), and has a very similar form. Likewise, if we apply to S 2 the optimal strategy π ε P U k found for S 1 in Example 1, we obtain a simulation depicted in Figure 4 for extended modes of length 2 and 4, which is the counterpart of Figure 2, and very similar to it. As seen above, we have:

}Πy π ε 2 pT q ´yf 1 } ď }y π ε 1 pT q ´yf 1 } `K2 σ |λ h1 | ,
where y f 1 " p0.3, 0.3, 0.3, 0.3, 0.3q, and the reduction error is bounded by K2σ |λ h 1 | " 17.9 σ. The subexpression }y π ε 1 pT q´y f 1 } can be computed a posteriori by simulation: see Table 1 of Appendix 2, with σ " 1, σ " 0.5. The value of }y π ε 2 pT q´y f 2 } for S 2 is also given in Table 1 for comparison.

The upper bound }y π ε 1 pT q ´yf 1 } `K2σ |λ h 1 | of the distance }P y π ε 2 pT q ´yf 1 } is very conservative, due to a priori error bound K2σ |λ h 1 | . On can obtain a posteriori a much sharper estimate of }P y π ε 2 pT q´y f 1 } by simulation: see Table 2, Appendix 2.

Final Remarks

Using the notion of OSL constant, we have shown how to use the finite difference and explicit Euler methods in order to solve finite horizon control problems for reaction-diffusion equations. Furthermore, we have quantified the deviation of this control with the optimal strategy, and proved that the error upper bound is linear in the horizon length. We have applied the method to a 1D bi-stable reaction-diffusion equation, and have found experimental results similar to those of [START_REF] Pouchol | Phase portrait control for 1D monostable and bistable reaction-diffusion equations[END_REF]. We have also given a simple and specific model reduction method which allows to apply the method to equations of larger size. In future work, we plan to apply the method to 2D reaction-diffusion equations (e.g., Test 1 of [START_REF] Alla | A HJB-POD approach for the control of nonlinear PDEs on a tree structure[END_REF]).

Fig. 1 .

 1 Fig. 1. Simulation of the system of Example 1 discretized with M1 " 5 points, for initial condition y0 " 0.8x{L `0.1p1 ´x{Lq, objective y f " 0.3 and horizon time T " 2 (τ " 0.1, ∆t " τ 100 ).

Fig. 2 .

 2 Fig. 2. Simulation of the system of Example 1 discretized with M1 " 5 points, with extended modes of length 2 (left) and extended modes of length 4 (right).

CpT q " Ope L f T q where L f is the Lipschitz constant associated with vector field f .

Note that, in[START_REF] Pouchol | Phase portrait control for 1D monostable and bistable reaction-diffusion equations[END_REF], the values of the boundary control are in the full interval r0, 1s, not in a finite set U as here. In[START_REF] Pouchol | Phase portrait control for 1D monostable and bistable reaction-diffusion equations[END_REF], they focus, not on the bounding of computation errors during integration as here, but on a formal proof that the objective state y f " θ (0 ă θ ă 1) is reachable in finite time iff L ă L ˚for some threshold value L ˚.

The program, called "OSLator"[START_REF] Le | OSLator 1.0[END_REF], is implemented in Octave. It is composed of 10 functions and a main script totalling

lines of code. The computations are realised in a virtual machine running Ubuntu 18.06 LTS, having access to one core of a 2.3GHz Intel Core i5, associated to 3.5 GB of RAM memory.

By comparison, in[START_REF] Alla | A HJB-POD approach for the control of nonlinear PDEs on a tree structure[END_REF], the error term originating from the POD model reduction is exponential in T (see C1pT, |x|q in the proof of Theorem 5.1).

Appendix 1: Proof of Lemma 1

Proof. It is easy to check that 0 ă α u ă 1 when |λu|Gu 4 ă 1. Let t ˚:" G u p1 ´αu q. Let us first prove δ e0 ptq ď e 0 for t " t ˚. We have:

Hence:

˚λu pt ˚q2 `2α u ´α2 u " 0.

We have: ´1 4G 2 u t ˚λu pt ˚q4 e λut ˚ě 0. It follows:

˚λu pt ˚q4 e λut ˚ě 0.

Hence:

By multiplying by t

Since G " ? 3|λ u |e 0 {C u :

By multiplying by λ u :

Note that, in the above formula, the subexpression λ u t ˚`1 2 λ 2 u pt ˚q2 is such that:

On the other hand, the subexpression ´1 3 λ u pt ˚q3 ´1 12 λ 2 u pt ˚q4 e λut ˚is such that:

u pt ˚q3 ´1 24 λ 4 u pt ˚q4 e λuθ q for some 0 ď θ ď t

It follows:

i.e.

pδ u e0 pt ˚qq 2 ď e 2 0 . Hence: δ u e0 pt ˚q ď e 0 . It remains to show: δ u e0 ptq ď e 0 for t P r0, t ˚s. Consider the 1rst and 2nd derivative δ 1 p¨q and δ 2 p¨q of δp¨q. We have: δ 1 ptq " λ u e 2 0 e λut `C2 u λ 2 u p2t `2 λu ´2 λu e λut q δ 2 ptq " λ 2 u e 2 0 e λut `C2 u λ 2 u p2 ´2e λut q. Hence δ 2 ptq ą 0 for all t ě 0. On the other hand, for t " 0, δ 1 ptq " λ u e 2 0 ă 0, and for t sufficiently large, δ 1 ptq ą 0. Hence, δ 1 p¨q is strictly increasing and has a unique root. It follows that the equation δptq " e 0 has a unique solution t ˚˚for t ą 0. Besides, δptq ď e 0 for t P r0, t ˚˚s, and δptq ě e 0 for t P rt ˚˚, `8q. Since we have shown: δpt ˚q ď e 0 , it follows t ˚ď t ˚˚and δptq ď e 0 for t P r0, t ˚s. l

Appendix 2: Numerical results

Dimension

Extended mode length }y π ε i pT q ´yf i } for σ " 1 }y π ε i pT q ´yf i } for σ " 0. i pT q ´yf i } for σ " 1 and σ " 0.5 (T " 2, i " 1, 2).

Extended mode length }P y π ε 2 pT q ´yf 1 } for σ " 1 }P y π ε 2 pT q ´yf 1 } for σ " 0. 2 pT q ´yf 1 } for σ " 1, σ " 0.5 (T " 2).