Charles-Henri Clerget

Thèse De Doctorat

Président M Yann Le Gorrec Ensmm

M Lorentz

T Biegler

M Gabriele Pannocchia

Jean-Philippe Grimaldi

Examinateur Total Rc

Mme Mériam Chebre

M Nicolas Petit

Charles-Henri Clerget Le

Spécialité Mathematiques

Automatique Dirigée

Nicolas Petit

Contributions to the control and dynamic optimization of processes with varying delays

Keywords: Optimisation dynamique, Retards variables, Contrôle de procédés Dynamic optimization, Variable delays, Process control Dynamic optimization, Variable delays, Process control

Dans cette thèse, nous avons étudié le contrôle et l'optimisation de systèmes dynamiques sujets à des retards variables.

L'existence de retards, de commande ou d'état, est un problème classique en automatique, susceptible de réduire les performances du système en régime transitoire, voire de remettre en cause la stabilité de contrôleurs en boucle fermée. De tels phénomènes de retards variables jouent un rôle important dans de nombreuses applications en génie des procédés.

Dans une première partie, nous avons étudié la régulation en boucle fermée d'un système soumis à des retards de métrologie variables et incertains. Nous avons établi de nouveaux résultats garantissant la stabilité robuste sous certaines conditions explicites sur le gain du contrôleur. Dans une seconde partie, nous avons abordé le problème de l'optimisation dynamique de systèmes présentant des retards variables dépendant de la commande liés à des phénomènes de transport dans des réseaux hydrauliques. Nous avons proposé un algorithme itératif d'optimisation et garanti sa convergence grâce à une analyse détaillée.

Introduction to varying delays and related control problems

In this thesis, we study some aspects of the general problem of controlling systems subject to delays. As is well-known, delays are very frequent in process industries, our primary focus for applications. Interestingly, this is not the only area of engineering where delays are ubiquitous. Delays are a central topic in many control applications, for instance in networked systems (see [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF]), aerospace (see [START_REF] Asher | Optimal control of systems with state-dependent time delay[END_REF], [START_REF] Hirzinger | Telerobotics with large time delays-the rotex experience. In Intelligent Robots and Systems '94[END_REF], [START_REF] Bresch-Pietri | Delay-adaptive predictor feedback for systems with unknown long actuator delay[END_REF] or [START_REF] Efremov | The aerospace vehicle control characterizing by considerable time delay[END_REF]), automotive engine (see [START_REF] Depcik | One-dimensional automotive catalyst modeling[END_REF], [START_REF] Lepreux | Warm-up strategy for a diesel oxidation catalyst[END_REF], [START_REF] Lepreux | Model-based Temperature Control of a Diesel Oxidation Catalyst[END_REF], [START_REF] Lepreux | Model-based temperature control of a Diesel Oxidation Catalyst[END_REF] or [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF], [GPP + 17], [START_REF] Sabatini | Design and experimental validation of a physics-based oxygen storage -thermal model for three way catalyst including aging[END_REF]), manufacturing industries (see [START_REF] Grimble | LQG controllers for state-space systems with pure transport delays: Application to hot strip mills[END_REF]), population dynamics (see [START_REF] Hrinca | An optimal control problem for the lotka-volterra system with delay[END_REF] or [START_REF] Lopes | Optimal release strategies for the biological control of aphids in melon greenhouses[END_REF]) or biomedical systems (see [START_REF] Rundeil | Systematic method for determining intravenous drug treatment strategies aiding the humoral immune response[END_REF], [START_REF] Ledzewicz | Antiangiogenic therapy in cancer treatment as an optimal control problem[END_REF] or [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]), to name but a few.

In process industries, delays (also referred to as dead-times or lags) are usually caused by the process itself, its instrumentation or its actuation and control technology. The delays can have several root causes. A first cause is data acquisition, the limitations of the Information Technology (IT) structure and the induced mis-synchronization between networks that lead to communication delays and errors in measurements time-stamping (see [START_REF] Jankovic | Developments in control of time-delay systems for automotive powertrain applications[END_REF], [START_REF] Noble | Loss, latency, and speed -TN0021[END_REF], [START_REF] Magnis | Impact of imprecise dating of measurements for bulk material flow network monitoring[END_REF] or [START_REF] Petit | Analysis of problems induced by imprecise dating of measurements in oil and gas production[END_REF]). Another common issue in process control is that the analysis of material samples may require a significant time to be performed. As evidenced for instance in [START_REF] Findeisen | Computational delay in nonlinear model predictive control[END_REF] or [START_REF] Zavala | The advanced-step NMPC controller: Optimality, stability and robustness[END_REF], an important source of lags in the control is the computation times of advanced process control algorithms such as Model Based Control (MPC) and its recent variants (adaptive MPC, economic MPC, distributed and cooperative MPC among others) which consider fairly advanced models of plants (e.g. [START_REF] Richalet | Model predictive heuristic control: Application to industrial processes[END_REF], [START_REF] Mayne | Constrained model predictive control: stability and optimality[END_REF], [DBS + 02], [START_REF] Qin | A survey of industrial model predictive control technology[END_REF], [START_REF] Pannocchia | Disturbance models for offset-free model-predictive control[END_REF], [START_REF] Pannocchia | Combined design of disturbance model and observer for offset-free model predictive control[END_REF], [START_REF] Rawlings | Model Predictive Control: Theory and Design[END_REF], [SVR + 10], [START_REF] Scheu | Sensitivity-based coordination in distributed model predictive control[END_REF], [START_REF] Maestre | Distributed model predictive control based on agent negotiation[END_REF] or [START_REF] Christofides | Distributed model predictive control: A tutorial review and future research directions[END_REF]). As for the intrinsic sources of delays in the processes themselves, complex chemical schemes involving activation times and the associated reaction lags are well-documented problems in chemical engineering (see [START_REF] Roussel | The use of delay differential equations in chemical kinetics[END_REF]). At a more macroscopic scale, material transport phenomena give birth to delays in plants where several parts of the process are not collocated. This is typically the case when some feeds are relatively far from the core of the reactor or when complex piping or recycling architectures are present (see [START_REF] Harmand | The optimal design of two interconnected (bio)chemical reactors revisited[END_REF], [START_REF] Chèbre | Feed control for an hydrodesulphurization unit using Anamel blend optimizer[END_REF], [START_REF] Cristea | Nonlinear predictive control of processes with variable time delay. a temperature control case study[END_REF], [START_REF] Roca | Solar field control for desalination plants[END_REF] or [START_REF] Zenger | Modelling and control of a class of timevarying continuous flow processes[END_REF] for instance). Such situations, with dead-times in pipes used to transfer material among process units, have been documented in many oil refining and petrochemical processes such as distillation columns, hydrodesulfurization units (see [START_REF] Barraud | Commande de procédés à paramètres variables[END_REF] which provides typical values of the delay and its variability) or fluid catalytic cracking (FCC e.g. [START_REF] Khandalekar | Control and optimization of fluid catalytic cracking process[END_REF] or [START_REF] Antoniades | Robust control of nonlinear time-delay systems[END_REF]). Finally, delays are also frequently used to build simple models of high-order, or even distributed, systems. For instance, a rich body of literature has studied the identification of First Order Plus Chapter 1. Introduction Dead Time (FOPDT) models (see [START_REF] Fedele | A new method to estimate a first-order plus time delay model from step response[END_REF], [START_REF] Niemann | A simple method for estimation of parameters in first order systems[END_REF]) and Second Order Plus Dead Time (SOPDT) models (see [START_REF] Neshasteriz | Generalized predictive control and tuning of industrial processes with second order plus dead time models[END_REF]). Such models can be used in combination with popular industrial PID auto-tune algorithms (e.g. [START_REF] Liu | A tutorial review on process identification from step or relay feedback test[END_REF]). Generalization to higher order models have also been studied (see [START_REF] Wang | Robust identification of continuous systems with dead-time from step responses[END_REF]).

In view of the general open-loop and closed-loop paradigms of automatic control, the presence of delays in the dynamics of a system raises a number of challenges. Closed-loop feedback-control schemes are usually employed for disturbance rejection and stabilization purposes. Having a delay in the loop causes several well-documented problems: degraded transient behaviours, oscillations around steady state and even loss of closed-loop stability (see [START_REF] Halanay | Differential equations. Stability, oscillations, time lags[END_REF], [START_REF] Rugh | Linear System Theory[END_REF], [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF] or [START_REF] Krstic | Delay Compensation for Nonlinear, Adaptive, and PDE Systems[END_REF], [START_REF] Aström | Feedback Systems: An Introduction for Scientists and Engineers[END_REF]), water-bed effects that limit both tracking and disturbance rejection performances and impose trade-offs (see [START_REF] Leong | Effects of delays, poles, and zeros on time domain waterbed tradeoffs and oscillations[END_REF]). Numerous strategies have been developed to address these issues, either extending the seminal idea of the Smith predictor through Internal Model Control (IMC) or using distributed control laws. The reader can refer to [START_REF] Morari | Robust Process Control[END_REF], [START_REF] Niculescu | Survey on recent results in the stability and control of time-delay systems[END_REF] or [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems: An Eigenvalue-Based Approach[END_REF] for further details. The classic techniques are limited to the case of constant delays. The case of time-varying delays have recently been studied following the early work of [START_REF] Krstic | Backstepping boundary control for firstorder hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF], treating the delay under the form of an hyperbolic equation stabilized using the tools of backstepping developed in the field of Partial Differential Equation (PDE) control. In the numerous extensions that have been considered following [START_REF] Krstic | Backstepping boundary control for firstorder hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF], the future evolution of the delay is either exactly known a priori or unknown but assumed to be smooth and bounded in a favourable sense (see [START_REF] Zhang | Delay-dependent stabilization of linear systems with time-varying state and input delays[END_REF], [START_REF] Bresch-Pietri | Delay-adaptive predictor feedback for systems with unknown long actuator delay[END_REF], [START_REF] Bekiaris-Liberis | Nonlinear Control Under Nonconstant Delays[END_REF], [START_REF] Bekiaris-Liberis | Robustness of nonlinear predictor feedback laws to time-and state-dependent delay perturbations[END_REF], [START_REF] Bekiaris-Liberis | Simultaneous compensation of input and state delays for nonlinear systems[END_REF] or [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF]).

In the context of open-loop motion planning and MPC control, the objective in process applications is chiefly to improve adverse transient behaviours leading to offspec production or degraded regimes. The key tools in this area are those of optimal control and its numerical implementations. In most cases, a trade-off must be made between accurate determination of optimal solutions and acceptable computation times. With constant delays, the mathematical formulations have long been investigated by the control community (see [START_REF] Kharatishvili | The maximum principle in the theory of optimal processes with delay[END_REF], [START_REF] Halanay | Optimal controls for systems with time lag[END_REF], [START_REF] Soliman | Optimal control of multivariable systems with pure time delays[END_REF], [START_REF] Malez-Zavarei | Suboptimal control of systems with multiple delays[END_REF], [START_REF] Basin | Optimal control for linear systems with multiple time delays in control input[END_REF] or [START_REF] Frederico | Noether's symmetry theorem for variational and optimal control problems with time delay[END_REF], [START_REF] Boccia | The maximum principle for optimal control problems with time delays[END_REF]). A detailed survey can be found in [START_REF] Göllmann | Optimal control problems with delays in state and control variables subject to mixed control-state constraints[END_REF] and [START_REF] Göllmann | Theory and applications of optimal control problems with multiple time-delays[END_REF]. These works cover cases of multiple input and state delays, with state constraints, in the framework of Pontryagin's maximum principle, [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. In their linear versions, industrial implementations in commercial MPC software are common place and routinely handle fixed delays. Following Richalet's early work (see [START_REF] Richalet | Pratique de la commande prédictive[END_REF]), AspenTech's DMCplus is extremely popular in the oil refining industry, along with the solutions of other vendors such as Honeywell or ABB. Interestingly enough, it appears that only little attention has been given to dynamic optimization problems under varying delays. Since the seminal work of [START_REF] Banks | Necessary conditions for control problems with variable time lags[END_REF], most research efforts have focused on closed-form solutions to LQR problems for dynamics impacted by time-varying delays, see [START_REF] Carravetta | Quadratic optimal control of linear systems with time-varying input delay[END_REF]. To the best of our knowledge, the only contributions regarding the case of a structured variability of the delay was laid out in the early work of [START_REF] Asher | Optimal control of systems with state-dependent time delay[END_REF]. This work treated the case of state dependant delays but does not seem to have received the appropriate attention, and is not implemented in any available software package.

In this thesis, we study cases of varying delays. We propose two contributions. Early in our work, we have realized that in many applications, the instrumentation delay is actually varying. We illustrate this phenomenon with an example of IMC where the delay is variable in an unstructured fashion and where instruments and IT limitations do not allow incoming measurements to be assigned an exact dating upon reception. This yields an unstructured varying delay in the underlying closed-loop discrete-time dynamics. After some steps of modelling, the problem boils down to studying the asymptotic stability of the following nonlinear dynamics

       u 0 = 0, δ 0 = 0 n 0, u n+1 = f -1 (c -δ n ) δ n+1 = δ n + α(f p (u n-Dn ) -f (u n-Dn-∆n ) -δ n )
where D n ≥ 0 is an uncertain delay and ∆ n is a mis-synchronization term while f and f p are close functions, c is a given parameter and α is a tuning parameter. We establish a proof of global convergence and robustness using an original convergence analysis. This contribution is presented in Chapter 2.

Then, we have spent most of our efforts on optimization problems. In a second part, we consider a class of structured delays denoted as hydraulic delays (see [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF] and [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF]). Such delays are defined by the following relation

t ru(t) φ(τ ) dτ = 1
where u is one of the system variable (typically some of the system inputs such as a flow rate), φ is a strictly positive, scalar-valued function and r u = t -D u (t) is the function of delayed time instants. This type of delay D u (t) is the exact solution of a plug-flow transport equation, hence its designation as hydraulic. Given a general objective function to minimize, we tackle the optimal control of systems whose dynamics is subject to such delays. Prior to our analysis, we conduct preliminary numerical investigations in Chapter 3. These outline both the surprisingly rich nature of the optimal control strategy of an elementary process along with the numerical challenges and limitations facing the state-of-the-art approaches and implementations.

Having met some difficulties in the numerics, we investigate the mathematical characterization of the problem and study its calculus of variations in a general context in Chapter 4. These investigations highlight the non-differentiability of the straightforward optimal control formulation. This serves as a basis to develop in Chapter 5 a complete numerical procedure to solve such problems. This procedure involves an explicit regularization and a specifically tailored fixed-point scheme. A complete convergence proof of this scheme is established. Numerical applications are treated in Chapter 6.

This work was carried out to support the development of the in-house tools used by TOTAL for process control and optimization in refining. It has been the subject of the following publications :

• C.-H. Clerget, J.-P. Grimaldi, M. Chèbre, and N. Petit

Introduction aux problèmes de retards variables en automatique

Dans cette thèse, nous étudions certains aspects du problème général de la commande des systèmes à retards. De tels phénomènes de retards sont bien connus dans les industries de procédés, qui constituent notre principal champs d'application. Ceci étant, il ne s'agit pas du seul domaine dans lequel la présence de retards est fréquente. Les retards constituent un problème central pour de nombreuses applications en automatique, par exemple pour le contrôle de systèmes en réseau (voir [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF]), en aéronautique (voir [START_REF] Asher | Optimal control of systems with state-dependent time delay[END_REF], [START_REF] Hirzinger | Telerobotics with large time delays-the rotex experience. In Intelligent Robots and Systems '94[END_REF], [START_REF] Bresch-Pietri | Delay-adaptive predictor feedback for systems with unknown long actuator delay[END_REF] ou [START_REF] Efremov | The aerospace vehicle control characterizing by considerable time delay[END_REF]), en contrôle moteur (voir [START_REF] Depcik | One-dimensional automotive catalyst modeling[END_REF], [START_REF] Lepreux | Warm-up strategy for a diesel oxidation catalyst[END_REF], [START_REF] Lepreux | Model-based Temperature Control of a Diesel Oxidation Catalyst[END_REF], [START_REF] Lepreux | Model-based temperature control of a Diesel Oxidation Catalyst[END_REF], [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF], [GPP + 17] ou [START_REF] Sabatini | Design and experimental validation of a physics-based oxygen storage -thermal model for three way catalyst including aging[END_REF]), dans les industries manufacturières (voir [START_REF] Grimble | LQG controllers for state-space systems with pure transport delays: Application to hot strip mills[END_REF]), en dynamique des populations (voir [START_REF] Hrinca | An optimal control problem for the lotka-volterra system with delay[END_REF] ou [START_REF] Lopes | Optimal release strategies for the biological control of aphids in melon greenhouses[END_REF]) ou encore en ingénierie biomédicale (voir [START_REF] Rundeil | Systematic method for determining intravenous drug treatment strategies aiding the humoral immune response[END_REF], [START_REF] Ledzewicz | Antiangiogenic therapy in cancer treatment as an optimal control problem[END_REF] ou [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]). Dans les industries chimiques, les retards peuvent avoir plusieurs causes. Ils sont en général causés par le procédé lui-même, son instrumentation ou les technologies de contrôle employées. Une première source de retards vient des processus d'acquisition de données, des défauts des architectures IT et de la désynchronisation induite entre réseaux qui conduisent à des retards de communication et des erreurs dans la datation des mesures (voir [START_REF] Jankovic | Developments in control of time-delay systems for automotive powertrain applications[END_REF], [START_REF] Noble | Loss, latency, and speed -TN0021[END_REF], [START_REF] Magnis | Impact of imprecise dating of measurements for bulk material flow network monitoring[END_REF] ou [START_REF] Petit | Analysis of problems induced by imprecise dating of measurements in oil and gas production[END_REF]). Un autre problème classique en contrôle de procédés est le fait que les mesures basées sur des analyses d'échantillons peuvent nécessiter un temps significatif pour être réalisées. Comme l'ont mis en évidence [START_REF] Findeisen | Computational delay in nonlinear model predictive control[END_REF] ou [START_REF] Zavala | The advanced-step NMPC controller: Optimality, stability and robustness[END_REF], les temps de calculs des algorithmes de contrôle avancé peuvent être importants. C'est en particulier le cas pour la commande prédictive et ses variantes récentes (commande prédictive adaptative, commande prédictive à critère économique ou commande prédictive distribuée par exemple) qui emploient des modèles relativement complexes des procédés (e.g. [START_REF] Richalet | Model predictive heuristic control: Application to industrial processes[END_REF] [START_REF] Maestre | Distributed model predictive control based on agent negotiation[END_REF] ou [START_REF] Christofides | Distributed model predictive control: A tutorial review and future research directions[END_REF]). En ce qui concerne les sources de retards propres aux procédés, les schémas réactionnels complexes présentant des temps d'activation et les retards associés sont des problèmes bien connus en génie chimique. A une échelle plus macroscopique, les phénomènes de transport de matière induisent également des retards dans les usines dont les différentes unités ne se trouvent pas au même endroit. C'est typiquement le cas lorsque certains flux d'alimentation sont loin du coeur du réacteur ou lorsque des architectures complexes de tuyauterie ou de recyclage sont présentes. De telles situations avec des retards dus à l'écoulement des produits entre unités est classique dans de nombreux procédés de raffinage et de pétrochimie comme la distillation, l'hydrodésulfuration (voir [START_REF] Barraud | Commande de procédés à paramètres variables[END_REF] qui fournit des ordres de grandeur du retard et de sa variabilité) ou le craquage catalytique (e.g. [START_REF] Khandalekar | Control and optimization of fluid catalytic cracking process[END_REF] ou [START_REF] Antoniades | Robust control of nonlinear time-delay systems[END_REF]). Enfin, des retards sont également souvent utilisés pour construire des modèles simples de systèmes d'ordres élevés (voire de dimension infinie). Ainsi, une vaste littérature a étudié l'identification de modèles du premier ordre avec retard (voire [START_REF] Fedele | A new method to estimate a first-order plus time delay model from step response[END_REF], [START_REF] Niemann | A simple method for estimation of parameters in first order systems[END_REF]) et du second ordre avec Chapter 1. Introduction retard (voir [START_REF] Neshasteriz | Generalized predictive control and tuning of industrial processes with second order plus dead time models[END_REF]). En particulier, de tels modèles peuvent être utilisés en combinaison avec des algorithmes de réglage automatique de PID (e.g. [START_REF] Liu | A tutorial review on process identification from step or relay feedback test[END_REF]). Des généralisations à des modèles d'ordres plus élevés ont également été étudiées (voir [START_REF] Wang | Robust identification of continuous systems with dead-time from step responses[END_REF]).

, [MRRS00], [DBS + 02], [QB03], [PR03], [PB07], [RM09], [SVR + 10], [SM11],
Du point de vue de l'automatique, la présence de retards dans la dynamique d'un système soulève de nombreuses difficultés, tant en boucle fermée qu'en boucle ouverte. Les algorithmes de contrôle en boucle fermée sont généralement utilisés pour rejeter des perturbations ou stabiliser un système. Dans ce cadre, la présence d'un retard peut causer plusieurs problèmes : un transitoire dégradé, des oscillations autour de l'état stable et même une perte de stabilité en boucle fermée (voir [START_REF] Halanay | Differential equations. Stability, oscillations, time lags[END_REF], [START_REF] Rugh | Linear System Theory[END_REF], [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF], [START_REF] Krstic | Delay Compensation for Nonlinear, Adaptive, and PDE Systems[END_REF] ou [START_REF] Aström | Feedback Systems: An Introduction for Scientists and Engineers[END_REF]), des effets de waterbed qui imposent des compromis entre les performances de suivi de trajectoire et de rejet de perturbations (voir [START_REF] Leong | Effects of delays, poles, and zeros on time domain waterbed tradeoffs and oscillations[END_REF]). De nombreuses stratégies existent pour traîter ces problèmes, soit en généralisant les idées du prédicteur de Smith par le contrôle à modèle interne, soit en utilisant des lois de contrôle intégrales. Le lecteur intéressé peut consulter [START_REF] Morari | Robust Process Control[END_REF], [START_REF] Niculescu | Survey on recent results in the stability and control of time-delay systems[END_REF] ou [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems: An Eigenvalue-Based Approach[END_REF] pour davantage de détails. Néanmoins, ces techniques classiques ne traitent que le cas de retards constants. Le cas de retards variables n'a que récemment été étudié, dans la foulée des travaux novateurs de [START_REF] Krstic | Backstepping boundary control for firstorder hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF], traitant le retard comme une équation de transport hyperbolique stabilisée en utilisant les outils du backstepping développés pour le contrôle d'équations aux dérivées partielles. Dans les nombreuses extensions qui ont été envisagées depuis, l'évolution future du retard est soit supposée être connue a priori, soit inconnue mais régulière et bornée dans un sens favorable (voir [START_REF] Zhang | Delay-dependent stabilization of linear systems with time-varying state and input delays[END_REF], [START_REF] Bresch-Pietri | Delay-adaptive predictor feedback for systems with unknown long actuator delay[END_REF], [START_REF] Bekiaris-Liberis | Nonlinear Control Under Nonconstant Delays[END_REF], [START_REF] Bekiaris-Liberis | Robustness of nonlinear predictor feedback laws to time-and state-dependent delay perturbations[END_REF], [START_REF] Bekiaris-Liberis | Simultaneous compensation of input and state delays for nonlinear systems[END_REF] ou [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF]).

Dans le cadre de la planification de trajectoire en boucle ouverte et de la commande par modèle prédictif, l'objectif est principalement de limiter les transitoires néfastes conduisant à la fabrication de produits hors spécifications ou à des régimes de fonctionnement dégradés. Les principaux outils utilisés dans ce domaine sont la théorie de la commande optimale ainsi que ses nombreuses implémentations numériques. Dans la plupart des cas, un compromis doit être trouvé entre la détermination exacte des solutions optimales et des temps de calcul acceptables en pratique. Le cas à retards constants a déjà fait l'objet d'une étude approfondie par les automaticiens (voir [START_REF] Kharatishvili | The maximum principle in the theory of optimal processes with delay[END_REF], [START_REF] Halanay | Optimal controls for systems with time lag[END_REF], [START_REF] Soliman | Optimal control of multivariable systems with pure time delays[END_REF], [START_REF] Malez-Zavarei | Suboptimal control of systems with multiple delays[END_REF], [START_REF] Basin | Optimal control for linear systems with multiple time delays in control input[END_REF], [START_REF] Frederico | Noether's symmetry theorem for variational and optimal control problems with time delay[END_REF] ou [START_REF] Boccia | The maximum principle for optimal control problems with time delays[END_REF]). Une revue détaillée peut être trouvée dans [START_REF] Göllmann | Optimal control problems with delays in state and control variables subject to mixed control-state constraints[END_REF] et [START_REF] Göllmann | Theory and applications of optimal control problems with multiple time-delays[END_REF]. Ces travaux couvrent le cas de retards multiples dans l'état et la commande, en présence de contraintes d'état, dans le cadre du principe du maximum de Pontryagin, [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. En outre, dans le cadre de systèmes linéaires, les logiciels commerciaux de commande à modèle prédictif utilisés dans l'industrie permettent généralement de prendre en compte des retards fixes. Dans la lignée des travaux pionniers de Richalet (voir [START_REF] Richalet | Pratique de la commande prédictive[END_REF]), l'application DMCplus d'AspenTech est ainsi devenu très populaire dans le secteur du raffinage, aux côtés de solutions d'autres fournisseurs comme Honeywell ou ABB. En revanche, il semble que peu d'attention ait été portée au cas de retards variables. Depuis les travaux fondateurs de [START_REF] Banks | Necessary conditions for control problems with variable time lags[END_REF], la plus grande partie des efforts de recherche ont porté sur la synthèse de solutions de problèmes LQR impactés par des retards variables, comme par exemple dans [START_REF] Carravetta | Quadratic optimal control of linear systems with time-varying input delay[END_REF]. A notre connaissance, la seule contribution à l'étude de retards variant de façon structurée (en fonction de l'état ou de la commande) a été établie dans les travaux précoces de [START_REF] Asher | Optimal control of systems with state-dependent time delay[END_REF]. Cette étude traitait le cas de retards dépendant de l'état mais ne semble pas avoir reçue l'attention qu'elle méritait, tant d'un point de vue théorique que pratique.

Dans cette thèse, nous étudions le cas de retards variables et proposons deux contributions. Dans de nombreuses applications de contrôle, les retards de métrologie liés à l'instrumentation du système sont en pratique variables. Dans la première partie de notre travail, nous étudions ce phénomène sur le cas d'un contrôleur à modèle interne soumis à un retard variable non structuré et dont le protocole d'acquisition des données ne permet pas de les dater de façon exacte. Cela conduit à une dynamique en temps échantillonné présentant un retard variable non structuré. Après quelques étapes de ré-écriture, le problème se ramène à l'étude de la stabilité asymptotique de la dynamique non linéaire suivante

       u 0 = 0, δ 0 = 0 n 0, u n+1 = f -1 (c -δ n ) δ n+1 = δ n + α(f p (u n-Dn ) -f (u n-Dn-∆n ) -δ n )
où D n ≥ 0 est un retard incertain, ∆ n est un terme de désynchronisation, f et f p sont des fonctions proches, c est un paramètre fixé et α est un paramètre de réglage. Nous établissons une preuve de convergence robuste en utilisant une technique de preuve originale. Cette contribution est présentée au Chapitre 2.

Par la suite, nous avons consacré l'essentiel de nos efforts à des problèmes d'optimisation dynamique. Dans la seconde partie de la thèse, nous considérons une classe de retards dits hydrauliques (voir [START_REF] Chèbre | Feedback control and optimization for the production of commercial fuels by blending[END_REF] et [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF]). De tels retards sont définis par la relation suivante

t ru(t) φ(τ ) dτ = 1
où u est un sous-ensemble de variables du système (typiquement une partie des variables de contrôle telles que des débits), φ une fonction à valeurs réelles strictement positives et r u = t -D u (t) la fonction des temps retardés. Ce type de retards est la solution exacte d'une équation d'écoulement piston, d'où sa désignation comme hydraulique. Étant donné un critère général à minimiser, nous étudions le problème du contrôle optimal de tels systèmes dont la dynamique présente ce type de retards. En amont de notre analyse, nous réalisons des investigations numériques préliminaires dans le Chapitre 3. Elles permettent de mettre en évidence la nature riche des solutions optimales au contrôle d'un procédé simple ainsi que les limitations de l'état de l'art sur ces problèmes.

Ayant rencontré des difficultés dans le traitement numérique de cet exemple, nous réalisons le calcul des variations de ce problème pour caractériser la structure de ses solutions optimales dans le Chapitre 4. Cette étude met en évidence qu'une formulation naïve du problème n'est pas différentiable. Ces résultats forment la base des développements du Chapitre 5 dans lesquelles une procédure complète est proposée pour résoudre ce type de problèmes. Cette procédure implique une régularisation du problème et sa résolution par un algorithme de point fixe. Une preuve de convergence de cet algorithme est présentée. Les applications numériques sont traitées dans le Chapitre 6.

Ce travail a été réalisé pour permettre l'améliorations d'outils propriétaires développés chez TOTAL pour le contrôle et l'optimisation de procédés de raffinage. Il a fait l'objet des publications suivantes : Tous les calculs réalisés dans cette thèse utilisent un processeur 2.60 GHz Intel(R) Core(TM) i7-4720HQ sur un système 64 bits avec 16.0 GB de RAM. Nous noterons CPU time le temps CPU nécessaire pour une accomplir une certaine tâche.

Introduction

In this chapter, we investigate the effects of delay variability and uncertainty on the classic Internal Model Controller (IMC, see e.g [START_REF] Morari | Robust Process Control[END_REF]) of a single-input single-output (SISO), static, nonlinear, sampled-data process with delayed measurements whose dating is uncertain. As is well-known, the uncertainty and the variability of delays lead to challenging control problems that may jeopardize closed-loop stability, see [START_REF] Krstic | Delay Compensation for Nonlinear, Adaptive, and PDE Systems[END_REF], [HOS + 16] and references therein. It is also known, see [WHQ + 05], that metrology delays coupled with inaccurate process models could lead to closed-loop instability. As discussed earlier in this thesis, the general treatment of these issues is still an open problem.

The process under consideration and its controller constitute a sampled-data system (following the terminology employed in e.g. [START_REF] Chen | Optimal sampled-data Control Systems[END_REF], [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF]) which can be reformulated using a classic discrete time representation. The specific case under consideration is actually also formally very similar to a scalar run-to-run controller, the robustness of which is not trivial. Run-to-run control is a popular and efficient class of techniques, originally proposed in [START_REF] Sachs | Process control system for VLSI fabrication[END_REF], specifically tailored for processes lacking in situ measurement for the quality of the production (see [START_REF] Wang | Survey on iterative learning control, repetitive control, and run-to-run control[END_REF]). Numerous examples of implementations have been reported in the semiconductor, and materials industry, in particular, see e.g. [START_REF] Wang | Survey on iterative learning control, repetitive control, and run-to-run control[END_REF][START_REF] Moyne | Run-to-Run Control in Semiconductor Manufacturing[END_REF] and references therein. Indeed, the field of run-to-run control encounters two of the practical problems addressed in this chapter: nonlinear model uncertainty and variable metrology delays. While these issues have often been reported (see, e.g. [WHQ + 05], [START_REF] Good | Stability analysis of double ewma run-to-run control with metrology delay[END_REF], [START_REF] François | Convergence analysis of run-to-run control for a class of nonlinear systems[END_REF], [START_REF] Wan | Double EWMA controller for semiconductor manufacturing processes with time-varying metrology delay[END_REF]), they have not received any definitive treatment from a theoretical viewpoint.

In the problem considered here, model uncertainty stems from the interactions between Chapter 2. A problem of robustness to uncertain metrology delay the input and the system states which can be rather complex, and, in turn, cause some non-negligible uncertainty on the quantitative effects of the input. On the other hand, the measurements are available after a long time lag covering the various tasks of sample collection, receipt, preparation, analysis and transfer of data through an information technology (IT) system to the control system. Measurements are thus impacted by large delays, which can be varying to a large extent, and in some applications be state-or input-dependant. As detailed in the introduction of this thesis, this variability of the delay builds up with the intrinsic IT dating uncertainty, because, in numerous implementations, no reliable timestamp can be associated to the measurements. The delay variability cannot be easily represented by Gaussian models (e.g. additive noise on the measurement), nor can it be fully described as deterministic input or state dependant delay, nor known varying delays that could be exactly compensated for by predictor techniques (as done in e.g. [BPCP12, BPCP14, BLK13b, BLK13c, BLK13a]).

In the absence of measurement dating uncertainty, robust stability in the presence of model mismatch can be readily established, using the monotonicity of the system and model which is formulated here as an assumption. The study of measurement dating uncertainty effects is more involved. Once expressed in the sampled time-scale, the control scheme exhibits a variable delay discrete-time dynamics. No straightforward eigenvalues or Nyquist criterion analysis (see [START_REF] Good | Stability analysis of double ewma run-to-run control with metrology delay[END_REF]) can be used to infer stability. Instead, a complete stability analysis in a space of sufficiently large dimension, with a well chosen norm, yields a proof of robust stability under a small gain condition. Interestingly, this small-gain bound is reasonably sharp, so that it can serve as guideline for practical implementation. The novelty of the approach presented in this chapter lies in the proof technique. It does not treat the uncertainty of the delay using the Padé approximation approach considered in [START_REF] Zhang | Minimum entropy based run-to-run control for semiconductor processes with uncertain metrology delay[END_REF], but directly uses an extended dimension of the discrete time dynamics. It is believed that these arguments of proof could be extended to address more general problems, in particular to higher dimensional forms (lifted forms) usually considered to recast general iterative learning control into run-to-run as is clearly explained in [START_REF] Wang | Survey on iterative learning control, repetitive control, and run-to-run control[END_REF].

The chapter has two objectives. After laying out a precise problem statement in Section 2.3, it establishes robust stability results with respect to model mismatch when measurements are delayed but exactly dated in Section 2.4.1. Then, it extends robust stability to small model errors when measurements are delayed and their dating is uncertain in Section 2.4.2. Those results are illustrated through simulations in Section 2.5.

Notations

Given I an interval of R, and f : I → R a smooth function, let us define

f ∞ = sup x∈I |f (x)|
For any vector X, note X 1 , X 2 and X ∞ its 1-norm, its Euclidean norm and its infinity norm, respectively. Note . * any of the vector norms above. For any square matrix A, note A * the norm of A, subordinate to . * . Classically (e.g. [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]), for all A, B, one has AB * ≤ A * B * We note x the floor value of x, mapping x to the largest previous integer.

For any matrix of dimension s, define E i the matrix of general term e k,l where δ is the Kronecker delta δ i,j = 1 if i = j and 0 otherwise.

∀(k, l), e k,l = δ k,s δ l,i (2.

Problem statement and proposed solution

Plant under consideration (delay-free)

We note y the controlled variable (output) of the considered plant and u the control variable (input). It is assumed that there exists f p a strictly monotonous smooth function such that

y = f p (u) (2.2)
Although f p is unknown, we can use a model of it, f , which is also smooth and monotonous1 , such that f p (0) = f (0). Usually, f is a rough estimate of f p . Typical models are represented in Figure 2.1. For the simulations considered in this chapter, the model error can be as large as 20-40%, which is representative of industrial applications requirements. The target value c for the controlled variable is assumed to be reachable by both the system and the model, i.e. there exists u c and ũc verifying

f p (u c ) = c, f (ũ c ) = c (2.3)

Metrology delay

A measurement system provides estimates of y with some time delay in a sampled manner. In many cases, this delay is time varying. Depending on the IT structure, measurements dating is usually done either using timestamping or an a priori estimation of the measurement delay. Either way, exact measurement dating is usually impractical, and some uncertainty on the measurement delay must be considered.

In the system considered in this chapter, the measurements available for feedback in a control loop thus have two specificities. They are delayed and the measurement delay 0 ≤ D itself is varying and uncertain. With 0 ≤ D the available estimation of D, we note ∆ D -D the mismatch. Next, we formulate two modelling assumptions Assumption 1. There exits D max such that D ≤ D max . Assumption 2. There exits ∆ max such that ∆ ≤ ∆ max . If Assumption 1 holds, it is clear from definition that -D max ≤ ∆.

Control problem

A closed-loop controller can be designed for system (2.2). Each time a measurement is received, the control is updated and the value of the control is kept constant until the next measurement is received, creating piece-wise constant control signals (with varying step-lengths). Repetitive application of this process generates a sequence of inputs and outputs. The delay results in shift of index in the measurement sequence.

Formally, the control design should aim at solving the following problem Problem 2.1

Create a sequence (u n ) using the approximate model f and the delayed measurements (f p (u n-Dn )) n∈N of y n such that lim

n→+∞ f p (u n ) = c
We propose a simple nonlinear IMC algorithm to address the problem. This algorithm adapts a bias term used in a model inversion. Assuming that one could estimate exactly the measurement delay D n , the implementation of such an algorithm would be

       u 0 = 0, δ 0 = 0, α ∈]0; 1] n 0, u n+1 = f -1 (c -δ n ) δ n+1 = δ n + α(y n-Dn -f (u n-Dn ) -δ n ) (2.4) with y n-Dn = f p (u n-Dn ).
(2.4) can be wrapped up in the following usual block diagram of Figure 2.2.

However, the uncertainties in measurements dating have an impact on the controller dynamics. Instead of (2.4), one is able to implement the following

       u 0 = 0, δ 0 = 0 n 0, u n+1 = f -1 (c -δ n ) δ n+1 = δ n + α(y n-Dn -f (u n-Dn ) -δ n ) When ∆ n = Dn -D n = 0, this becomes        u 0 = 0, δ 0 = 0 n 0, u n+1 = f -1 (c -δ n ) δ n+1 = δ n + α(y n-Dn -f (u n-Dn-∆n ) -δ n ) (2.5) f -1 D D f p f Low Pass Filter u c + y + - δ - Figure 2.2: Idealized closed-loop control scheme. f -1 D ∆ f p f D Low Pass Filter u c + y + - δ - Figure 2.3: Realistic closed-loop control scheme.
where ∆ n is a dating uncertainty term. The situation is pictured in Figure 2.3. To show that (2.5) constitutes a viable solution to our control problem 2.1, it is necessary to investigate the closed-loop stability of the controller in this case.

Convergence analysis

Convergence with model mismatch and delay but without measurement dating uncertainty

In the analysis, three problems must be treated: model mismatch, delayed measurements and measurement dating uncertainty. We first consider the system without the later. Used in closed loop, controller (2.4) gives

       u 0 = 0, δ 0 = 0, δ 1 = α(f p (0) -f (0)) n 0, u n+1 = f -1 (c -δ n ) δ n+2 = (1 -α)δ n+1 + α(δ n-D n+1 -c + f p • f -1 (c -δ n-D n+1 )) (2.6)
The asymptotic behaviour of (2.6) is determined by the extended dynamics of (δ n ) since convergence of (δ n ) clearly implies convergence of (u n ). If (u n ) and (δ n ) converge toward the limits u and δ respectively, then, necessarily,

u = u c and δ = c -f (u c )
We now define the sequence (d n δ n -δ, n ≥ 0). Equivalently, the error dynamics is represented by the equation

d n+2 = (1 -α)d n+1 + α(d n-D n+1 + f p • f -1 (f (u c ) -d n-D n+1 )) -αc
Chapter 2. A problem of robustness to uncertain metrology delay Applying the mean value theorem to the function x → x + f p • f -1 (f (u c ) -x), one easily deduces that there exists

a n ∈ [min(0, d n-D n+1 ); max(0, d n-D n+1 )] such that d n+2 = (1 -α)d n+1 + α 1 - f p • f -1 (f (u c ) -a n ) f • f -1 (f (u c ) -a n ) d n-D n+1
Gathering past values of d n over the range n -D max , ..., n + 1 into a single vector X n , the system can be written as a linear time varying system (LTV) of dimension p D max + 2

X n+1 = A n X n (2.7)
where

X n = d n-Dmax • • • d n+1 T with A n = C + αh(a n )F n (2.8)
where

C =                  0 1 0 • • • • • • • • • 0 . . . . .
• • • • • • • • • • • • 0 1 -α                 
using the notation (2.1),

F n = E Dmax+1-D n+1 and h(a n ) = 1 - f p • f -1 (f (u c ) -a n ) f • f -1 (f (u c ) -a n ) (2.9)
One shall note that, h can be interpreted as a metric of the model error: if f ≡ f p , we do indeed get h ≡ 0. Since (2.7) is a LTV system, establishing its convergence is non-trivial. Establishing the asymptotic (not to say exponential) convergence of a general LTV discrete time system is usually a difficult task. In particular, it is not sufficient to study its eigenvalues (see [START_REF] Rugh | Linear system Theory[END_REF]). Some results have long been available for slowly varying systems (see [START_REF] Rosenbrock | The stability of linear time-dependent control systems[END_REF] for instance) and have recently been refined in [START_REF] Hill | Exponential stability of time-varying linear systems[END_REF], in particular. However, in our present case, it is not necessary to use them. The particular structure of the varying term allows more straightforward investigations. Define the (infinite) set of possible transition matrices (2.8)

A = {C + αh(x)E i , x ∈ R, i ∈ 1; D max + 1 } (2.10)
Let us assume that h ∞ < 1. This assumption implies that the set A is bounded. Consider a sequence of n transition matrices (A k ) k∈ 0;n-1 ∈ A n and for i ∈ 0; n -1 note

A k = C + αh k E i k 2.4. Convergence analysis 25 Define ∀k ∈ 1; n , Π k = k i=1 A k-i =     L k 1 . . . L k p     (2.11)
where L k i designates the i th row of the product of the k matrices. Calculating the product Π k+1 , one obtains that

L k+1 Dmax+2 = (1 -α)L k Dmax+2 + αh k+1 L k i k+1 While for j ∈ 1; D max + 1 L j k+1 = L j+1 k It follows that L k+1 Dmax+2 1 ≤ ((1 -α) + α h ∞ <1 ) max( L k Dmax+2 1 , L k i k+1 1 )
Recursively, it is then straightforward to show that

Π Dmax+1 ∞ = max j∈ 1;Dmax+2 L Dmax+1 j 1 ≤ (1 -α) + α h ∞ (2.12)
Now, we establish a preliminary result. Consider a discrete linear time-varying system (2.13) of dimension s, and A a bounded set of possible transition matrices in M s (R) and initial condition X 0 ∀n ≥ 0, X n+1 = A n X n , A n ∈ A

(2.13)

For any vector norm . * and any N ∈ N * , we define

M N, * sup A N -i ∈A N i=1 A N -i * = sup A i ∈A N -1 i=0 A i * (2.14)
Proposition 2.1: Sufficient condition for exponential stability

Consider the system (2.13). If there exists N 0 ∈ N * such that M N 0 , * < 1, then the system (2.13) (globally) exponentially converges to 0. One has, for some K > 0,

∀n ∈ N * , X n * ≤ K X 0 * (M N 0 , * ) n N 0 (2.15)
Proof. The proof is relatively straightforward

∀n ∈ N, X n = n i=1 A n-i X 0
Hence, grouping terms in N 0 -size bundles starting from the right

X n * ≤ n-n N 0 N 0 i=1 A n-i * × n N 0 i=1 N 0 j=1 A n N 0 N 0 -(i-1)N 0 -j * X 0 *
Chapter 2. A problem of robustness to uncertain metrology delay and

X n * ≤ M n-n N 0 N 0 , * M n N 0 N 0 , * X 0 * Besides, ∀n ∈ N, 0 ≤ n - n N 0 < N 0
Hence, we get the desired result by defining

K max k∈ 0;N 0 -1 M k, *
As a consequence, using the notation (2.14)

M Dmax+1,∞ = sup (A i )∈A Dmax+1 Π Dmax+1 ∞ < 1
which, according to Proposition 2.1, means that the (X n ) sequence is exponentially convergent. This, in turn, allows us to formulate the following result

Theorem 2.1: Global convergence without measurement dating uncertainty

Let ∆ = 0. Consider any 0 < α ≤ 1.

If h ∞ < 1, then the closed loop error (2.6) converges exponentially and

lim n→+∞ f p (u n ) = c
Remark 2.1. In particular, one can notice that f and f p must have the same sign so that the condition h ∞ < 1 can be verified. In this case, if

0 < f p f ∞ < 2 (2.16)
then the sufficient condition is satisfied.

Remark 2.2. The result derived in Theorem 2.1 is a sufficient condition for the controller stability. We can still get some additional insight into the controller behaviour by studying the particular case in which D max = ∆ max = 0. Then, we derive a necessary stability condition from a straightforward eigenvalue analysis showing that the equilibrium point of the system is locally stable if and only if

0 ≤ f p (u c ) f (u c ) ≤ 1 + 1 α (2.

17)

If α = 1, this shows that the sufficient condition previously derived is also necessary. Otherwise, this shows that taking α < 1 small enough may allow one to stabilize systems where the ratio

f p
f is greater than 2 (actually, this is indeed observed in simulations).

Convergence with measurement dating error

We now consider the implementation of the controller (2.5) with measurement dating uncertainty causing the discussed mis-synchronization between measurement and prediction with ∆ = 0. Using the same transformation as in Section 2.4.1, we establish the closed-loop error

d n+2 = (1 -α)d n+1 +α f p (f -1 (f (u c ) -d n-D n+1 )) -f (u c ) + d n-D n+1 -∆ n+1 -α(c -f (u c ))
Applying the mean value theorem, we get

d n+2 = (1 -α)d n+1 -αρ(a n )d n-D n+1 + αd n-D n+1 -∆ n+1 where ρ = 1 -h and a n ∈ [min(0, d n-D n+1 ); max(0, d n-D n+1 )]
The system can be written as a LTV system of dimension p D max + ∆ max + 2

X n+1 = A n X n (2.18)
where

X n = d n-Dmax-∆max • • • d n+1 T with A n = C + αF n -αρ(a n )F n (2.19)
and, with the notation (2.1),

F n = E p-1-D n+1 -∆ n+1 and F n = E p-1-∆ n+1

Convergence analysis without model error

Let us first assume that there is no model error. Under this assumption ρ = 1 and the transition matrices A n of the dynamics (2.18) all belong to the finite set

A = {C + αE k -αE k , (k, k ) ∈ 1; p -1 × p -1 -D max ; p -1 } (2.20) Consider a sequence of n transition matrices (A i ) i∈ 0;n-1 ∈ A n . Similarly to Sec- tion 2.4.1, define ∀k ∈ 1; n , Π k = k i=1 A k-i =     L k 1 . . . L k p     (2.21)
Chapter 2. A problem of robustness to uncertain metrology delay where L k i designates the i th row of the product of the k matrices. The convergence analysis is built recursively upon the fact that there exists K > 0 such that for all n

∈ N * M n,∞ = n i=1 A n-i ∞ ≤ Kη n-1 p+Dmax (2.22) where η max r∈ 1;Dmax+1 1 - α (1 -α) r + 2α 1 -α (1 -(1 -α) r ) (2.23) For all n ≥ 2, it is clear that ∀j ∈ 1; p -1 , L n j = L n-1 j+1 (2.24)
and

∃(r n , m n ) ∈ 1; p -1 2 , L n p = (1 -α)L n-1 p -αL n-1 p-rn + αL n-1 p-mn (2.25)
We wish to prove that there exists K such that the following relation holds for all n ≥ 0

∀j ∈ 1; p , L n j 1 ≤ Kη n+j-2 p+Dmax (2.26) Let us define K max (A 0 ,...,A n-1 )∈A n max i∈ 1;Dmax+1 L i p 1
It is clear that (2.25) is true for all indexes from 2 to D max + 1. Given n ≥ D max + 1, let us assume that the property is true for this rank. One has

Π n+1 = n+1 i=1 A n+1-i =     L n+1 1 . . . L n+1 p     with ∀j ∈ 1; p -1 , L n+1 j = L n j+1 and L n+1 p = (1 -α)L n p -αL n p-r n+1 + αL n p-m n+1
Hence, according to (2.24)

L n+1 p = (1 -α)L n p -αL n-r n+1 p + αL n p-m n+1 (2.27)
To proceed to the induction, we choose to develop the second term of (2.27). According to (2.25) at rank n

L n-r n+1 +1 p = (1 -α)L n-r n+1 p -αL n-r n+1 p-r n-r n+1 +1 + αL n-r n+1 p-m n-r n+1 +1
It follows that

L n-r n+1 p = 1 1 -α (L n-r n+1 +1 p + αL n-r n+1 p-r n-r n+1 +1 -αL n-r n+1 p-m n-r n+1 +1 ) (2.28)
As a consequence, after substitution with (2.28), (2.27) gives

L n+1 p = (1 -α)L n p + αL n p-m n+1 2.4. Convergence analysis 29 - α 1 -α (L n-r n+1 +1 p + αL n-r n+1 p-r n-r n+1 +1 -αL n-r n+1 p-m n-r n+1 +1 )
Recursively, from rank n -r n+1 + 1 to n, we get

L n+1 p = [(1 -α) - α (1 -α) r n+1 ]L n p + αL n p-m n+1 -α 2 r n+1 -1 i=0 1 (1 -α) i+1 (L n-r n+1 +i p-r n-r n+1 +1 +i -L n-r n+1 +i p-m n-r n+1 +1 +i ) Then, if the following condition holds 0 ≤ (1 -α) - α (1 -α) r n+1
(2.29)

Using one cancellation and a careful succession of terms reorderings, one has

L n+1 p 1 ≤   1 - α (1 -α) r n+1 + 2α 2 r n+1 -1 i=0 1 (1 -α) i+1   max j∈ 1;p k∈ n-rmax;n L k j
And, finally, using the explicit summation of the geometric sequence

L n+1 p 1 ≤ 1 - α (1 -α) r n+1 + 2α 1 -α (1 -(1 -α) r n+1 ) max j∈ 1;p k∈ n-Dmax;n L k j
which leads, by induction with (2.26), to

L n+1 p 1 ≤ Kη n-Dmax+1-2 p+Dmax +1
and after a simplification

L n+1 p 1 ≤ Kη (n+1)+p-2 p+Dmax
This proves (2.26) at rank n + 1. As a consequence, (2.22) directly follows using the relation between the infinity norm of a matrix and the one norm of its rows

∀n ∈ N, Π n ∞ = max i∈ 1;p L n i 1

General case

Based on this first result, we introduce a small model error, and formulate an extension by continuity. This last result shows that the proposed controller solves the control problem at stake, in the presence of model mismatch, delayed measurements and dating error.

According to (2.22) there exists N 0 ∈ N such that if there is no model error

M N 0 ,∞ ≤ 1 2
With model error, any transition matrix of the dynamics A n can be written under the additive form

A n = A 0 n + P n (2.30)
where A 0 n is a matrix of the set (2.20)

A 0 n ∈ {C + αE kn -αE k n , (k n , k n ) ∈ 1; p -1 × p -1 -D max ; p -1 } (2.31)
Chapter 2. A problem of robustness to uncertain metrology delay and P n is a perturbation matrix

P n = αh(x n )E k n (2.32)
with x n a given real number. Consider any collection of N 0 such matrices (A i ) i∈ 0;N 0 -1 and assume that there exists > 0 such that h ∞ ≤ , then

N 0 -1 i=0 A i ∞ ≤ N 0 i=0 A 0 i ∞ + N 0 i=1 C N 0 -i N 0 (1 + α) N 0 -i α i i ≤ 1 2 + N 0 i=1 C N 0 -i N 0 (1 + α) N 0 -i α i i-1
By upper-bounding the (finite) sum appearing in the right-hand side, it follows that there exists a sufficiently small value of such that for any (A i ) i∈ 0;N 0 -1

N 0 -1 i=0 A i ∞ ≤ 3 4 < 1 (2.33)
Then, Proposition 2.1 yields the exponential convergence of X n and leads to the following (main) result Dmax+1 (2.34) and max

Theorem 2.2: Exponential convergence under measurement dating uncertainty and model mismatch

Let ∆ ≤ ∆ max . Consider any 0 < α ≤ 1 such that 0 ≤ (1 -α) - α (1 -α)
r∈ 1;Dmax+1 1 - α (1 -α) r + 2α 1 -α (1 -(1 -α) r ) < 1 (2.35)
There exists ∈ R + such that, if h ∞ ≤ , then the controller (2.5) is exponentially stabilizing and guarantees lim

n→+∞ f p (u n ) = c
Remark 2.3. In particular, one sees from (2.22) that the larger D max and ∆ max is, the slower the guaranteed convergence rate is.

Remark 2.4. One can easily check that there always exist a neighbourhood of α = 0 on which conditions (2.34)-(2.35) are verified. Indeed

(1 -α) - α (1 -α) Dmax+1 -→ α→0 + 1 and ∀r ≥ 1 1 - α (1 -α) r + 2α 1 -α (1 -(1 -α) r ) = 1 -α + O(α 2 )
While (2.34)-(2.35) are only sufficient conditions derived in a conservative fashion, we expect that as α increases, it may reach a bifurcation value beyond which the system is no longer stable. This behaviour is confirmed by numerical simulations in Section 2.5.1.

Numerical illustration

Remark 2.5. In terms of controller design specifications, we can draw the following statements from our analysis :

• the sign of the estimated gain must be correct and its value cannot be too small compared with the reality. On the other hand, taking an estimated gain larger than the true one will slow the controller down but cannot jeopardize its stability

• variable delays cause no specific problem of convergence if we assume that exact measurement dating is available

• if there is dating uncertainty of the measurement, stability can still be retained, provided that the measurements filtering is strong enough (α small enough)

Numerical illustration

In this section, we consider a setting where each sample is analysed during a certain lapse of time during which no other sample is taken. New control actions are only implemented when a new measurement result is received. In that way, the time-sampled measurement delay is always zero (D n = 0, D max = 0), i.e. the measurement we receive is always informative of the result of the last control action taken. This is a special case of Theorem 2.2 which is of practical importance in the implementation of many controllers.

In classic run-to-run cases, nonzero D max could be considered, without loss of generality. We will assume that the actual physical time required for the measurement to reach the controller, T p , depends on the measured value according to the following relation

T p (y) = 8.5 -0.75y mes + 10(ξ 1 + ξ 2 ) (2.36)
where ξ 1,2 are the realizations of two independent random variables taking the values 0 or 1 with respective probabilities of 1 2 . In some simulations, we will assume that timestamps are not available and that the times at which the measurements are taken have to be estimated using an approximate model T of T p . This estimation may be inexact, thus leading to dating uncertainty of the measurements, i.e. ∆ max ≥ 0. In all simulations, the target will be c = -10.

Influence of pure dating uncertainty

In this subsection, the response of the system, f p , is assumed to be a simple linear function of the control u. We further assume that it is perfectly known, so that

f p (u) = f (u) = -25u
(2.37) Furthermore, we consider a situation where no timestamp is available. As a consequence, T p must be estimated and we assume that the approximate model available to the controller is

T (y) = 17 -1.5y mes (2.38)
This mismatch results into nonzero values of ∆ which are plotted to provide a graphical view of its statistical distribution. We also corrupt the measurements with a zero-mean, uniform noise of small amplitude to excite the system. We simulate the response of the system for various values of the gain α. The results show that while the closed-loop system remains stable for α small enough, as α increases, Chapter 2. A problem of robustness to uncertain metrology delay destabilization can arise in absence of any model error on the function f p , simply because of measurement dating uncertainty. This illustrates the fact that Theorem 2.2 does not hold after a critical value of α, which is conservatively estimated by conditions (2.34) and (2.

Influence of the system response uncertainty

We further illustrate Theorem 2.2 by simulating the closed-loop response of the system with a model mismatch (but without any measurement noise). We assume that f p is given by f p (u) = -8.4339 -6 arctan(8u -6) (2.39) and f by

f (u) = -25u (2.40)
We first run the simulations without measurement dating error with three different values of α corresponding to a low (α = 0.38 discussed above), medium (α = 0.60) and high (α = 0.80) filtering. Figure 2.5 shows the results of the simulations. Despite model mismatch, all of them display asymptotically stable responses.

Figure 2.6 then shows the results of the simulations where measurement dating error is introduced (the previous expressions for T p and T are used). One sees that beyond a critical value for the gain α, the response becomes increasingly poor, and eventually becomes unstable for the high α scenario.

These simulations illustrate the merits of the theoretical results established in this chapter. A tuning of the controller gain following the (conservative) estimate provided by the small-gain condition gives satisfactory closed-loop responses even when the measurement dating error is not negligible. If the gain is chosen above the threshold, some divergence (or strong oscillations) can be observed. The situation is similar with reasonable levels of measurement noise.

Conclusions

As a static SISO control problem, the core problem tackled in this chapter appears, at first sight, as simple as it could be. However, the variability and uncertainty of the delay makes the problem particularly tricky. We have provided explicit robustness margins in regard of model error and asymptotic analysis on the consequences of uncertain measurement dating.

In the case where an underlying dynamical system should be considered to model the system, the preceding approach should be updated, significantly. Because the measurement will remain sampled by nature, the closed loop system will naturally become a sampleddata ordinary differential equation as considered in e.g. [START_REF] Fridman | Robust sampled-data stabilization of linear systems: an input delay approach[END_REF]. Also, it is known that the introduction of time-varying gains may improve the exponential convergence, when measurements are subjected to (known) delays. If estimates of the delay are available, such tuning rules could bring some performance improvement. While the problem becomes significantly harder due to the time-varying nature of the discretized system transition matrices, it would be interesting to investigate whether, in a more general context of multi-input multi-output (MIMO) dynamical systems, an event-triggered discretization approach such as the one developed here could be used to obtain results on the influence of measurement dating uncertainty.

Chapter 2. A problem of robustness to uncertain metrology delay

Introduction

In this chapter, we consider dynamic optimization of systems subject to hydraulic time delays. The field of dynamic optimization is primarily divided into two types of approaches: direct and indirect methods. In indirect methods, the problem is viewed as one of functional analysis (the unknown variables being functions) and its resolution is based upon the calculus of variations to determine stationarity conditions necessarily satisfied by locally optimal solutions. These conditions can be derived in a variety of settings using Pontryagin maximum principle, see [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. Classically, this involves the introduction of a dual (adjoint) variable. The stationarity conditions have the form of a two-point boundary value problem (TPBVP). The numerical resolution then boils down to finding solutions of this TPBVP. A detailed overview of methods can be found in [START_REF] Bryson | Applied Optimal Control[END_REF], [START_REF] Roberts | Two-point boundary value problems: shooting methods[END_REF], [START_REF] Ascher | Numerical solution of boundary value problems for ordinary differential equations[END_REF], [START_REF] Bonnans | The shooting algorithm for optimal control problems: A review of some theoretical and numerical aspects[END_REF] or [START_REF] Bonnans | The shooting approach to optimal control problems[END_REF]. Direct methods, on the other hand, first discretize the input sequence into a vector of finite dimension and approximate the continuous optimal control problem as a Non Linear Program (NLP) that can be handled using numerical solvers (see e.g. [START_REF] Betts | Practical Methods for optimal control using nonlinear programming[END_REF]). Direct simultaneous methods are a subclass of direct methods. They discretize both the control and state variables and rely on orthogonal collocations on finite time elements to achieve Chapter 3. Direct simultaneous strategies for hydraulic delays a high-order approximation of the dynamics of the system. The resulting set of discretized equations can be shown to be equivalent to a high-order implicit Runge-Kutta integration method (see [START_REF] Biegler | Nonlinear Programming, Concepts, Algorithms, and Applications to Chemical Processes[END_REF]). Solving the problem then requires the solution of large but very structured NLP (the Jacobian of the constraints being sparse). This can be efficiently done using dedicated software packages, such as IPOPT. For further details regarding the merit of the different types of direct methods and the current state of the art in terms of direct simultaneous strategies, the reader can refer to [START_REF] Biegler | An overview of simultaneous strategies for dynamic optimization[END_REF] for a comprehensive survey and [START_REF] Flores-Tlacuahuac | Global optimization of highly nonlinear dynamic systems[END_REF] or [START_REF] Shi | Optimization of grade transitions in polyethylene solution polymerization processes[END_REF] for detailed case-studies.

These approaches both have well-known advantages and shortcomings. The trade-off between them is hence very dependant on the applications, and in some cases they can even be advantageously combined (see [START_REF] Stryk | Direct and indirect methods for trajectory optimization[END_REF]). Indirect methods are recognized for their accuracy and their speed of convergence. However, they require good initial guesses to converge. More critically, the treatment of constraints is difficult in indirect methods and requires a good a priori knowledge of the solution structure. Specifically, the number of constrained arcs in the optimal solution needs to be estimated a priori in order to parametrize the resolution of the problem. They have found their main applications in the field of aerospace control, (see for instance [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF], [START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF], [START_REF] Aronna | A well-posed shooting algorithm for optimial control problems with singular arcs[END_REF] or [START_REF] Trélat | Optimal control and applications to aerospace: Some results and challenges[END_REF]). Direct methods, on the other hand, are much more robust to poor starting guesses and allow a straightforward treatment of inequality constraints. They can also be viewed as more straightforward to implement as no adjoint variable need to be explicitly defined prior to resolution. Among them, collocation-based simultaneous approaches have attractive numerical stability and approximation properties.

A problem worth discussing is constraints, both on input and state variables. In practical process applications, bounds are routinely imposed on variables. Thermodynamic equilibriums also frequently give rise to algebraic relations between variables that can not be easily written-off and can hamper the stability of the solutions.

Time constants in process industries are relatively larger than in most other control applications. Nowadays, the computational burden associated with the resolution of direct formulations of dynamic optimization is considered as tractable for many applications with standard (desktop) computers. Generally speaking, real-time implementations of MPC can take advantage of the results of the previous iterations to provide a good initialization point for the following run. A practical discussion regarding the implementation of such "warmstart" procedure for primal-dual interior point methods (such as the one implemented in IPOPT) can be found in [START_REF]Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT[END_REF]. Alternatively, strategies such as advanced-step MPC (see [START_REF] Zavala | The advanced-step NMPC controller: Optimality, stability and robustness[END_REF]) where an optimal solution is pre-computed between sampling instants based on a prediction of the future state and very quickly updated using a sensitivity analysis when a new measurement is received can be considered.

For these reasons, direct methods are usually favoured in this field and we have decided to treat the class of systems under consideration in this thesis using these state of the art solutions, specifically a collocation-based simultaneous strategy. This is the topic covered in this chapter. In Section 3.2, we briefly review some of the existing literature on the application of collocation-based simultaneous approaches for dynamic optimization of systems subject to time delay and emphasize the specific challenges raised by the case of input-dependant hydraulic delays. In Section 3.3, we introduce a tutorial example of a water heater subject to such delays. It will serve as illustrating example. In Section 3.4, we present a resolution strategy based on the use of a partial differential equation (PDE) to describe the delay equation and its approximation by a finite-order model. Several discretization strategies are presented, along with their numerical results. As will appear, the results stress the rich (complex) nature of the optimal solution, and the numerical difficulties. Their limits motivate the investigation of an alternative approach in the next chapters of the thesis.

Direct simultaneous methods for systems subject to delays

Delays cause specific concerns in the field of simultaneous methods, because they require specific tools in their numerical simulation. A rich body of literature has long studied the numerical simulation of delay-differential algebraic equations (DDAE). Useful references can be found in [START_REF] Banks | Spline approximations for functional differential equations[END_REF], [START_REF] Karoui | Computer solutions of state-dependent delay differential equations[END_REF] or [START_REF] Ascher | The numerical solution of delay-differentialalgebraic equations of retarded and neutral type[END_REF]. In particular, the implementation of collocation methods has been studied, for instance in [GH01]- [START_REF] Betts | Lobatto IIIA methods, direct transcription, and DAEs with delays[END_REF]. However, to the best of our knowledge, there exists no contributions on dynamic optimization using such schemes for input or state dependant delays. Indeed, the formulation of this case seems to raise specific difficulties. In order to illustrate this, let us consider a system defined by a state variable x, a control variable u and a hydraulic delay r u such that the dynamics be given by

∀t ∈ [0; T ] , ẋ(t) = f (x(t), x(r u (t)), u(t))) (3.1)
and

t ru(t) φ(u(τ ))dτ = 1
with f and φ > 0 smooth functions. Consider a collocation mesh defined by N c ≥ 1 finite elements of time of equal size h = T Nc , each containing K collocation points. We note t ij the time instant associated to the j th collocation point of the i th time element and t i the endpoint of the i th time element. In general, the time instants r u (t ij ) have no reason to match the discretization points of the collocation mesh and the transcription of (3.1) requires the local interpolation at each mesh point t ij of x(r u (t ij )). This might seem a minor difficulty as, typically, noting U = (u i ) 1;Nc the vector of discretized inputs and using a Lagrange interpolation formula, one has

ỹ(r u (t ij )) ≈ K k=1 L k r U (t ij ) -t i-d(U )-1 h y i-d(U ),k (3.2)
where d(U ) is defined as

r U (t ij ) ∈ t i-d(U )-1 ; t i-d(U )
However, (3.2) raises the problem of being non differentiable at every point where there exists integers (i, j, k) ≥ 0 such that r U (t ij ) = t k . This is made apparent by the equivalent representation using binary variables

ỹ(r u (t ij )) ≈ Nc r=r 0 δ r K k=1 L k r U (t ij ) -t r-1 h y r,k
where

δ r = 1 if r U (t ij ) ∈ [t r-1 ; t r [ 0 if r U (t ij ) ∈ [t r-1 ; t r [
The complexity of optimization problems involving binary variables (MINLP, e.g. [START_REF] Grossmann | Review of nonlinear mixed-integer and disjunctive programming techniques[END_REF]) is such that such a formulation is not suitable for a dynamic optimization applications.

Chapter 3. Direct simultaneous strategies for hydraulic delays q(t)

Example of application with hydraulic delay

The system under consideration here is a flow rate controlled water heater with downstream measurements (pictured on Fig. 3.1) first introduced in [MSP08] and [START_REF] Himpe | Design and advanced control of a process with variable time delay[END_REF] (where complete description of the test bench is given) to outline the challenges associated with closed-loop control of systems featuring hydraulic time delays. It is composed of a tank filled with a constant volume V of water and heated by a fixed thermal flux Q. A controlled flow rate of water q passes through the tank, coming in at a fixed inlet temperature T in . Since water gets heated as it flows through the tank, the outlet temperature of the tank is higher than T in . After having left the tank, water flows through a pipe of (constant) cross-section S over a length L. The variable one measures and seeks to control is the temperature T o at the outlet of this pipe. Neglecting heat losses, the average temperature in the tank T (t) satisfies the following balance equation

dT (t) dt = Q ρc p V + q(t) V (T in -T (t)) f (T, q) (3.3)
where ρ, c p and Q are the density of water, its specific heat and the power of supplied heat, respectively. Assuming instantaneous mixing in the tank, one has

T o (t) = T (t -D(t, q)) = T (r q (t)) with t t-D(t,q)
q(τ ) dτ = LS

In our numerical investigations, we will take Q = 10 7 J.s -1 , ρ = 10 3 kg.m -3 , c p = 4185 J.kg -1 .K -1 , V = 1 m 3 , L = 0.5 m, S = 1 m 2 , T in = 0°C . Given some desired 3.4. Numerical implementation and results 41 reference signals T r and q r (chosen to be consistent with the desired steady-state value T r ), the optimal control problem under consideration is

min q T 0 T o (t) -T r (t) 2 + w • q(t) -q r (t) 2 dt s.t. Ṫ = f (T, q) T o (t) = T (r q (t)) T (t ≤ 0) = T 0 q min ≤ q(t) ≤ q max
(3.4)

Numerical implementation and results

As discussed in Section 3.2, the transcription of the optimal control problem (3.4) into a smooth NLP using orthogonal collocations is directly not possible. As a consequence, instead of working directly with the delayed equation, a classic idea is to replace it with the underlying transport equation governing the system (e.g. [START_REF] Shi | Optimization of grade transitions in polyethylene solution polymerization processes[END_REF]). This leads to the following optimization problem

min q T 0 T o (t) -T r (t) 2 + w • q(t) -q r (t) 2 dt s.t. Ṫ = f (T, q) ∂Z ∂t (x, t) = - q(t) S ∂Z ∂x (x, t), x ∈ [0; L] T o (t) = Z(L, t) Z(0, t) = T (t) Z(x, 0) = T 0 , x ∈ [0; L] q min ≤ q(t) ≤ q max (3.5)
Formally, this change of representation does not generate any approximation (equation (3.3) is the exact solution of the PDE of (3.5), see [START_REF] Bresch-Pietri | Robust control of variable time-delay systems. Theoretical contributions and applications to engine control[END_REF]). The resolution of (3.5) requires the discretization of the transport PDE. It is well-known that good numerical schemes can be obtained for transport phenomena using finite volumes methods (see [START_REF] Leveque | Finite Volumes Methods for Hyperbolic Problems, Cambridge texts in applied mathematics[END_REF]). Classically, space is divided into a set of cells over which averaged properties are defined (with a running index j)

T j (t) = 1 ∆x x j+ 1 2 x j-1 2 Z(x, t) dx
We follow the approach described in [START_REF]Advanced Strategies for Optimal Design and Operation of Pressure Swing Adsorption Processes[END_REF] based on the use of the Method of Lines (MOL) by discretizing the PDE only with respect to space into a set of ordinary differential equations. Such a (second-order accurate) scheme defined over a mesh of N cells is given by

dT n dt = - q(t) ∆V T n (t) -T n-1 (t) + 1 2 (T n (t) -2T n-1 (t) + T n-2 (t)) T 0 (t) = T 1 (t) = T (t), T o (t) = T N (t), ∆V = V L N -1 (3.6)
Now, using (3.6) we have a standard optimal control problem, the resolution of which is straightforwardly implemented using a simultaneous transcription strategy (see [START_REF] Biegler | An overview of simultaneous strategies for dynamic optimization[END_REF]).

Chapter 3. Direct simultaneous strategies for hydraulic delays

The discretization of the ODEs is achieved using 3 Radau collocation points (see [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]) in each time finite element and the discretized optimisation problems are solved using IPOPT 3.11.8 through the algebraic modelling language AMPL (see [START_REF] Biegler | Nonlinear Programming, Concepts, Algorithms, and Applications to Chemical Processes[END_REF]). All simulations are initialized using a constant solution where T = 2°C and q has the consistent stationary value.

The numerical results of this method are shown on Figures 3.2-3.3 (in both cases the number of finite time elements of the collocation mesh is taken equal to N c = 100). T and T o are the values computed by the optimizer while T o,ex is the value of T o computed a posteriori by a high fidelity simulator (using a very refined upwind scheme). Figure 3.2 displays the optimal trajectories computed with w = 0.1 for two values of N (10 and 50). Good results are achieved in both cases. Refining the spatial discretization by a factor of 5 allows us to decrease the discrepancy between the value of T o predicted by the optimizer and the actual realization of the trajectory. This must however be paid by a large increase in computation time (by a factor close to 40). Furthermore, in both cases, the resolution of the problem suffers from ill-conditioning. More precisely, the Newton method that IPOPT follows requires the Jacobian of the KKT matrix to be definite positive. This is true if he constraint Jacobian is full row rank and the reduced Hessian is positive definite. If this condition is not verified, or the numerical conditioning of the matrix is too poor, IPOPT must regularize it by adding a damping term to the Hessian (this behaviour appears for N ≥ 8) 1 . The interested reader can refer to [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF] for a detailed presentation of the line-search method implemented in IPOPT. The log file of IPOPT for the resolution of the problem in the case with w = 0.1 and N = 50 can also be found in Appendix A (the definition of the various outputs of IPOPT are available in [START_REF]Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT[END_REF]).

Interestingly, in both cases, the structure of the solution is qualitatively the same and displays an apparently pseudo-periodic behaviour during the transient state. Despite its complex structure, this type of solution does not necessarily contradict the physical intuition. Indeed, low flow-rate regimes correspond to phases where the water in the tank gets heated while high flow-rate regimes lead to a quick flush of the outlet pipe (thus allowing hot water to reach faster the outlet of the pipe). This cyclic functioning could hence be described as "heat and flush".

Figure 3.3 displays the optimal trajectories computed with w = 0.001 for the same two values of N (10 and 50). Interestingly, it shows that when the regularization parameter w becomes sufficiently small, the solution of the optimization are subject to strong spurious oscillations around steady state. Figure 3.3 evidences the fact that refining the spatial discretization does not seem to solve the problem but actually worsens it.

Numerical results can still be improved by refining the time mesh of the collocation scheme while keeping the number of actual degrees of freedom on the control input constant (equality constraints are imposed on a groups of n b adjacent inputs). Figure 3.4 shows the results obtained when refining the time mesh by a factor 5. While the precision of optimizer's prediction is much improved and some of the spurious oscillations of the solution is removed, this is paid by an increase of computation time by a factor roughly equal to 5.

Comparable results are obtained using a full discretization, both w.r.t. time and space, of the PDE using the following second order accurate scheme

T j (t i+1 ) = T j (t i ) - ∆t ∆x (F j+ 1 2 (t i ) -F j-1 2 (t i ))
1 In IPOPT these properties can be assessed through the inertia of the matrix. If a regularization term is used, it is listed as lg(rg). where

F j+ 1 2 (t i ) = q(t i ) S T j-1 (t i ) + 1 2 q(t i ) S (1 - q(t i )∆t S∆x )(T j (t i ) -T j-1 (t i ))
It should be remembered that this type of finite volumes numerical schemes is only stable if the Courant-Friedrichs-Lewy [START_REF] Allaire | Numerical analysis and optimization: an introduction to mathematical modelling and numerical simulation[END_REF] condition is verified

q∆t S∆x < 1 (3.7)
As a consequence, with such an explicit scheme in time, satisfying (3.7) for a reasonably large value of N requires a very refined time discretization. In this approach, the transport PDE is directly transformed into a set of algebraic equations and we only need to apply collocations to the remaining ordinary differential equation of the tank. To allow N = 40, we discretize the PDE at 3500 time instants. To keep the comparison meaningful, equality constraints are once again imposed on adjacent inputs in order to keep the effective number of degrees of freedom of the control input equal to 100. The results are displayed on Figure 3.5. As can be seen, this approach achieves excellent accuracy but at the price of a large increase in computation time with respect to the findings of Figure 3.3.

Conclusions

Altogether, this study shows that while interesting results can be obtained using an approach based on the discretization of the transport PDE (which is an exact representation of the delay), it presents two main drawbacks:

• even for an elementary system, accurately capturing sharp transient fronts requires a fine spatial discretization which increases the size of the state and leads to large optimization problems.

Chapter 3. Direct simultaneous strategies for hydraulic delays .5: Optimal trajectories with full PDE discretization, w q = 0.001, N = 40, CPU time = 142.17 s. The "heat and flush" behaviour is also found. High accuracy of the optimizer prediction is achieved

t [s]
• refining the spatial discretization of the PDE leads to problems of numerical conditioning of the Hessian, which adversely affect the performances of the solver and lead to undesirable behaviour of the solutions (spurious oscillations).

In the remaining of the thesis, we will focus our investigations on the derivation of a methodology allowing us to exploit our knowledge of the analytical solution of the transport PDE as a delay equation in order to overcome these issues.

Chapter 4

Calculus of variations for varying delays

Calcul des variations pour des systèmes à retards variables. Dans ce chapitre, nous étudions le calcul des variations de problèmes de commande optimale de systèmes soumis à différents types de retards variables. Nous montrons en particulier que, dans le cas de retards hydrauliques apparaissant dans la commande, le problème n'est pas différentiable. Une formulation régularisée est proposée pour éviter cette difficulté et ses conditions de stationnarités sont calculées.

Introduction

Prior to developing a sort of indirect numerical method for the problem under consideration in this thesis, we carry a preliminary study on the associated calculus of variations, and focus on the derivation of stationarity conditions. The study of the calculus of variations for optimal control problems subject to time-varying delays in the dynamics is not a new subject. Following the early works of the optimal control community in the 50s and the 60s, [START_REF] Banks | Necessary conditions for control problems with variable time lags[END_REF] derived necessary stationarity conditions for a system subject to multiple fixed (with respect to the state and the control) time-varying delays in the state. [START_REF] Asher | Optimal control of systems with state-dependent time delay[END_REF] studied the case of a generic state-dependant state delay. However, since then the field has received limited attention and these works are not implemented in any available software package. Interestingly, none of the works cited above are directly treating the case of input-dependant input delays, which is our primary focus. As will be seen, this case is, in fact, more involved.

In the first part of this chapter, we will study the calculus of variations of an optimal control problem with a fixed time-varying delay1 in the state and the input. This is a fairly straightforward extension of the reordering techniques presented in [START_REF] Hughes | Variational and optimal control problems with delayed argument[END_REF] for fixed delays and will serve as an introduction to the second part of the chapter where we study the case of an input-dependant hydraulic delay in the state and the control. We first establish the conditions under which such a system is Gâteaux differentiable and outline that a problem with an explicit input-dependant input delay dependency generally does Chapter 4. Calculus of variations for varying delays not verify these conditions. We then introduce a smooth regularized version of the problem and derive its stationarity conditions. Interestingly, this latter result can be seen as a refinement of the results of the pioneer work of [START_REF] Asher | Optimal control of systems with state-dependent time delay[END_REF] to the case of an hydraulic delay, under regularization.

Notations

We denote u ∈ R p the control and x ∈ R m the state of the system.

Let L 2 (E, F ) be the set of functions of integrable square on E with values in F . With T > 0, we denote C 1 pw ([0; T ], R p ) the class of piecewise C 1 functions from [0; T ] to R p , having a finite number of jumps in their values or derivatives on their interval of definition.

Let D 1 ([0; T ], R p ) be the class of functions from [0; T ] to R p that are differentiable but whose derivative are not necessarily continuous.

We denote the function sign : R → {-1; 0; 1} that maps strictly positive (resp. negative) arguments to 1 (resp. -1) and 0 onto itself.

Let g ∈ L 2 (R, R n ) and I be an interval of R. We denote g I the restriction of g to I.

Let f be a function of a real variable, we note the one-sided limit whose s-argument defines how the t-argument is approached

lim τ → s t f (τ ) =    lim τ →t + f (τ ) if 0 ≤ s lim τ →t -f (τ ) if 0 > s

Fixed time-varying delays

Consider r a smooth, strictly increasing function such that for all t, r(t) < t. This function defines a delayed time law.

Let

L : [0; T ] × R m × R p → R and f : [0; T ] × R m × R m × R p × R p → R m be smooth functions. Take (u 0 , x 0 ) ∈ L 2 ([r(0); 0], R p ) × D 1 ([r(0); 0], R p ).
Consider the following optimal control problem having r as fixed time-varying delay

P r : min u T 0 L(t, x(t), u(t)) dt + ψ(x(T )) J r (u) s.t. ∀t ∈ [0; T ], ẋ(t) = f (t, x(t), x(r(t)), u(t), u(r(t))) x [r(0);0] = x 0 , u [r(0);0] = u 0
We seek to establish necessary stationarity conditions characterizing optimal solutions of P r . Following the classic approach of calculus of variations (see e.g. [START_REF] Gelfand | Calculus of Variations[END_REF]), these conditions are equivalent to the stationarity conditions of the augmented functional where the constraints of the dynamics have been adjoined

Jr (x, u, λ) = T 0 L(t, x(t), u(t)) -λ(t) T ẋ(t) + λ(t) T f (t, x(t),

x(r(t)), u(t), u(r(t))) dt + ψ(x(T ))

Using an integration by parts, one has

Jr (x, u, λ) = T 0 L(t, x(t), u(t)) + λ(t) T x(t) + λ(t) T f (t, x(t), x(r(t)), u(t), u(r(t))) dt -λ(T ) T x(T ) + λ(0) T x(0) + ψ(x(T )) (4.1)

Fixed time-varying delays 49

In order to compute the Gâteaux derivatives of Jr (i.e. the directional derivatives, see [START_REF] Rudin | Functional analysis[END_REF]), given any δ ∈ R, x, u, λ, let us first consider the cost variation associated with a perturbation of its first argument in a direction h. We will use the following simplified notation for greater convenience Jr (x + δh, u, λ) -Jr (x, u, λ) = Jr (x + δh) -Jr (x) Using (4.1), we have

Jr (x + δh) -Jr (x) = δ T 0 ∂L ∂x (t, x(t), u(t))h(t) + λ(t) T h(t) + λ(t) T ∂f ∂x (t, x(t), x(r(t)), u(t), u(r(t)))h(t) + λ(t) T ∂f ∂x r (t, x(t), x(r(t)), u(t), u(r(t)))h(r(t)) dt -δλ(T ) T h(T ) + δ ∂ψ ∂x (x(T ))h(T ) + o(δ)
This immediately leads to the expression of the Gâteaux derivative w.r.t. the x-argument

D h Jr (x) = T 0 ∂L ∂x (t, x(t), u(t))h(t) + λ(t) T h(t) + λ(t) T ∂f ∂x (t, x(t),

x(r(t)), u(t), u(r(t)))h(t) + λ(t) T ∂f ∂x r (t, x(t), x(r(t)), u(t), u(r(t)))h(r(t)) dt -λ(T ) T h(T ) + ∂ψ ∂x (x(T ))h(T )

This last expression is not handy for the coming derivation of stationarity conditions because the expression under the integral sign mixes the values of h at both time t and time r(t). Since h is an admissible perturbation, for all t ≤ 0, h(t) = 0. This gives a first simplification. Then, using a change of variables, one finds

T 0 λ(t) T ∂f ∂x r (t, x(t), x(r(t)), u(t), u(r(t)))h(r(t)) dt = r(T ) 0 λ(r -1 (t)) T • ∂f ∂x r (r -1 (t), x(r -1 (t)), x(t), u(r -1 (t)), u(t))(r -1 ) (t)h(t) dt
Finally, this leads to

D h Jr (x) = -λ(T ) T + ∂ψ ∂x (x(T )) h(T ) + T 0 ∂L ∂x (t, x(t), u(t)) + λ(t) T + λ(t) T ∂f ∂x (t, x(t), x(r(t)), u(t), u(r(t))) + 1 [0;r(T )] (t)(r -1 ) (t) • λ(r -1 (t)) T • ∂f ∂x r (r -1 (t), x(r -1 (t)), x(t), u(r -1 (t)), u(t)) h(t) dt (4.2)
Chapter 4. Calculus of variations for varying delays Similarly, we establish the Gâteaux derivative w.r.t. the input

D h Jr (u) = T 0 ∂L ∂u (t, x(t), u(t)) + λ(t) T ∂f ∂u (t

, x(t), x(r(t)), u(t), u(r(t))

)

+ 1 [0;r(T )] (t)(r -1 ) (t) • λ(r -1 (t)) T • ∂f ∂u r (r -1 (t), x(r -1 (t)), x(t), u(r -1 (t)), u(t)) h(t) dt (4.3)
and finally,

D h Jr (λ) = T 0 -h(t) T ẋ(t) + f (t, x(t), x(r(t)), u(t), u(r(t))) dt (4.4)
Any stationary solution (x * , u * , λ * ) of Jr is characterized by the relations

∀(h 1 , h 2 , h 3 ), D h 1 Jr (x * ) = D h 2 Jr (u * ) = D h 3 Jr (λ * ) = 0
Then, using Dubois-Reymond lemma (see [START_REF] Gelfand | Calculus of Variations[END_REF]), we classically establish the following result Theorem 4.1:

The stationarity conditions necessarily verified by locally optimal solutions of P r are given by the following two-point boundary value problem (TPBVP)

ẋ(t) = f (t, x(t), x(r(t)), u(t), u(r(t))) (4.5) x(0) = x 0 λ(t) = - ∂L ∂x (t, x(t), u(t)) T - ∂f ∂x (t

, x(t), x(r(t)), u(t), u(r(t))) T λ(t)

-

1 [0;r(T )] (t)(r -1 ) (t)• ∂f ∂x r (r -1 (t), x(r -1 (t)), x(t), u(r -1 (t)), u(t)) T λ(r -1 (t)) (4.6) λ(T ) = ∂ψ ∂x (x(T )) T 0 = ∂L ∂u (t, x(t), u(t)) T + ∂f ∂u (t, x(t), x(r(t)), u(t), u(r(t))) T λ(t) + 1 [0;r(T )] (t)(r -1 ) (t)• ∂f ∂u r (r -1 (t), x(r -1 (t)), x(t), u(r -1 (t)), u(t)) T λ(r -1 (t)) (4.7)

Input-dependant delays: a non-smooth problem

Gâteau differentiability with input-dependant delays

We now consider a more complex case when the delay depends on the input signal, according to an integral (hydraulic) law. Let φ : R p → R * + be a smooth function. Take (u 0 , x 0

) ∈ C 1 pw ([r 0 ; 0], R p )×D 1 ([r 0 ; 0], R p ), r 0 < 0 with 0 r 0 φ (u 0 (τ )) dτ = 1
4.4. Input-dependant delays: a non-smooth problem 51 Consider the optimal control problem with input-dependant delays

P 0 : min u T 0 L(t, x(t), u(t)) dt + ψ(x(T )) J 0 (u) s.t. ∀t ∈ [0; T ], ẋ(t) = f (t, x(t), x(r u (t)), u(t), u(r u (t))) x [r 0 ;0] = x 0 , u [r 0 ;0] = u 0
where r u is implicitly defined by the relation

D h r u (t) lim δ→0 r u+δh (t) -r u (t) δ = 1 φ(u(r u (t))) t ru(t) ∂φ ∂u (u(τ ))h(τ ) dτ (4.9)
Similarly, for any t ∈ [0; r u (T )[

lim δ→ s 0 r -1 u+δh (t) -r -1 u (t) δ = - 1 lim τ → s r -1 u (t) φ(u(τ )) r -1 u (t) t ∂φ ∂u (u(τ ))h(τ ) dτ (4.10)
where

s = sign -s • r -1 u (t) t ∂φ ∂u (u(τ ))h(τ ) dτ
and if u is continuous at r -1 u (t), the Gâteaux derivative is given by

D h r -1 u (t) = - 1 φ(u(r -1 u (t))) r -1 u (t) t ∂φ ∂u (u(τ ))h(τ ) dτ (4.11)
Proof. From (5.2), we have

1 = t ru(t) φ(u(τ )) dτ = t r u+δh (t) φ(u(τ ) + δh(τ )) dτ
Chapter 4. Calculus of variations for varying delays Then, from the smoothness of φ, one deduces that

t ru(t) φ(u(τ )) dτ = t r u+δh (t) φ(u(τ )) + δ ∂φ ∂u (u(τ ))h(τ ) dτ + o(δ)
It follows that

r u+δh (t) ru(t) φ(u(τ )) dτ = δ t r u+δh (t) ∂φ ∂u (u(τ ))h(τ ) dτ + o(δ) (4.12)
Since we know a priori that r 0 ≤ r u+δh (t) ≤ T , one notices that the integral in the right-hand side is expressed on a bounded domain over which its argument is bounded. We immediately deduce that

r u+δh (t) ru(t) φ(u(τ )) dτ --→ δ→0 0
Using the fact that φ is a strictly positive function by definition, we obtain the continuity of r u (t) with respect to the input

r u+δh (t) -r u (t) --→ δ→0 0
Using this result in (4.12), we have Using this, the desired results (4.8) are finally obtained by taking alternatively the limit of (4.13) when δ goes to zero from above or below. Otherwise, when (4.14) fails, we directly get 1 δ

r u+δh (t) ru(t) φ(u(τ )) dτ = o(1)
and

lim δ→0 r u+δh (t) -r u (t) δ = 0
The results regarding the variation of r -1 u (t) are established symmetrically after noticing that the definition of the delay implies that for all t ∈ [0; min(r u (t), r u+δh (t))]

1 = r -1 u (t) t φ(u(τ )) dτ = r -1 u+δh (t) t φ(u(τ ) + δh(τ )) dτ
4.4. Input-dependant delays: a non-smooth problem 53 Using this preliminary result, we can continue our analysis. For notational ease, we write

[t, x, u] v = t, x(t), x(r v (t)), u(t), u(r v (t))
Let us derive the stationary points of the augmented functional

J0 (x, u, λ) = T 0 L(x(t),u(t)) + λT (t)x(t) + λ T (t)f ([t, x, u] u ) dt + λ T (0)x(0) -λ T (T )x(T ) + ψ(x(T ))
The Gâteaux derivative of J0 w.r.t. its second argument at point (x, u, λ) in direction h is

J0 (u + δh) -J0 (u) = T 0 L(t, x(t),u(t) + δh(t)) -L(x(t), u(t)) + λ T (t) f [t, x, u + δh] u+δh ) -f [t, x, u] u
We know that there exists a finite number N of distinct time instants r 0 < t i ≤ r u (T ) i=1..N at which the control input u is not smooth. For δ small enough, u + δh has the same jumping points as u, plus those generated by δh which will all have negligible contributions. The calculus is decomposed over a mesh allowing us to cover both cases when the image of the jumps of u by the inverse perturbed delayed time law, r -1 u+δh (t), are each approached from below or above. Using this notation, we have

J0 (u + δh) -J0 (u) δ = T 0 ∂L ∂u (t, x(t), u(t))h(t) + λ(t) T ∂f ∂u [t, x, u] u • h(t) dt + min(r -1 u (t 1 ),r -1 u+δh (t 1 )) 0 ∆(t, δ) dt + N i=1 1 δ max(r -1 u (t i ),r -1 u+δh (t i )) min(r -1 u (t i ),r -1 u+δh (t i )) λ T (t)• f t, x(t), x(r u+δh (t)), u(t) + δh(t), u(r u+δh (t)) + δh(r u+δh (t)) -f t, x(t), x(r u (t)), u(t), u(r u (t)) dt + N -1 i=1 min(r -1 u (t i+1 ),r -1 u+δh (t i+1 )) max(r -1 u (t i ),r -1 u+δh (t i )) ∆(t, δ) dt + T max(r -1 u (t N ),r -1 u+δh (min(t N ,r u+δh (T )))) ∆(t, δ) dt + o(1) (4.16) with ∆(t, δ) = λ T (t) • ∂f ∂x r [t, x, u] u ẋ(r u (t)) • r u+δh (t) -r u (t) δ + ∂f ∂u r [t, x, u] u u(r u (t)) • r u+δh (t) -r u (t) δ + ∂f ∂u r [t, x, u] u • h(r u+δh (t))
where ∂f ∂xr , ∂f ∂ur designate the partial derivatives of f w.r.t. its 3 rd and 5 th arguments respectively.

Chapter 4. Calculus of variations for varying delays Using (4.10) from Proposition 4.1, we know that on a neighbourhood of δ = 0, if the upper and lower Gâteaux derivative of r -1 u (t i ) are non zero at u, we have

sign r -1 u+δh (t i ) -r -1 u (t i ) = sign -δ r -1 u (t i ) t i ∂φ ∂u (u(τ ))h(τ ) dτ
The strict monotonicity of r u and r u+δh gives that

r -1 u (t i ) ≤ t ≤ r -1 u+δh (t i ) =⇒ r u+δh (t) ≤ t i ≤ r u (t) (4.17) and r -1 u+δh (t i ) ≤ t ≤ r -1 u (t) =⇒ r u+δh (t) ≥ t i ≥ r u (t) (4.18)
Both of these inequalities (4.17) and (4.18) are instrumental for the evaluation of the integrals max(.) min(.) (.) in (4.16), by determining the arguments of f as δ goes to zero. This gives

1 δ max(r -1 u (t i ),r -1 u+δh (t i )) min(r -1 u (t i ),r -1 u+δh (t i )) λ T (t)• f t, x(t), x(r u+δh (t)), u(t) + δh(t), u(r u+δh (t)) + δh(r u+δh (t)) -f t, x(t), x(r u (t)), u(t), u(r u (t)) dt = 1 δ max(r -1 u (t i ),r -1 u+δh (t i )) min(r -1 u (t i ),r -1 u+δh (t i )) λ(t) T • f t, x(t), x(r u (t)), lim τ → r -1 u (t i ) u(τ ), lim τ → - t i u(τ ) -f t, x(t), x(r u (t)), lim τ → r -1 u (t i ) u(τ ), lim τ → t i u(τ ) dt + o(1) (4.19)
Otherwise, if the Gâteaux derivative of r -1 u (t i ) is equal to zero, we have

lim δ→0 1 δ max(r -1 u (t i ),r -1 u+δh (t i )) min(r -1 u (t i ),r -1 u+δh (t i )) λ T (t)• f t, x(t), x(r u+δh (t)), u(t) + δh(t), u(r u+δh (t)) + δh(r u+δh (t)) -f t, x(t), x(r u (t)), u(t), u(r u (t)) dt = 0
Finally, for s ∈ {-1; 1}, gathering the smooth and jump parts of the calculus and using (4.8) along with (4.10), we have lim

δ→ s 0 J0 (u + δh) -J0 (u) δ = T 0 ∂L ∂u (t, x(t), u(t))h(t) + λ T (t) ∂f ∂u [t, x, u] u h(t) + λ T (t) ∂f ∂x r [t, x, u] u • ẋ(r u (t)) φ(u(r u (t))) t ru(t) ∂φ ∂u (u(τ ))h(τ ) dτ + λ T (t) ∂f ∂u r [t, x, u] u • u(r u (t)) φ(u(r u (t))) t ru(t) ∂φ ∂u (u(τ ))h(τ ) dτ 4.4. Input-dependant delays: a non-smooth problem 55 + λ T (t) ∂f ∂u r [t, x, u] u h(r u (t)) dt + N i=1 λ(r -1 u (t i )) T • f r -1 u (t i ), x(r -1 u (t i )), x(t i ), lim τ → s (t i ) r -1 u (t i ) u(τ ), lim τ → -s (t i ) t i u(τ ) -f r -1 u (t i ), x(r -1 u (t i )), x(t i ), lim τ → s (t i ) r -1 u (t i ) u(τ ), lim τ → s (t i ) t i u(τ ) • s (t i ) lim τ → s (t i ) r -1 u (t) φ(u(τ )) r -1 u (t i ) t i ∂φ ∂u (u(τ ))h(τ ) dτ (4.20)
where

s (t) = sign -s • r -1 u (t
N i=1 λ(r -1 u (t i )) T • f r -1 u (t i ), x(r -1 u (t i )), x(t i ), lim τ → s 1 (t i ) r -1 u (t i ) u(τ ), lim τ → -s 1 (t i ) t i u(τ ) -f r -1 u (t i ), x(r -1 u (t i )), x(t i ), lim τ → s 1 (t i ) r -1 u (t i ) u(τ ), lim τ → s 1 (t i ) t i u(τ ) • s 1 (t i ) lim τ → s 1 (t i ) r -1 u (t i ) φ(u(τ )) r -1 u (t i ) t i ∂φ ∂u (u(τ ))h(τ ) dτ = N i=1 λ(r -1 u (t i )) T • f r -1 u (t i ), x(r -1 u (t i )), x(t i ), lim τ → -s 1 (t i ) r -1 u (t i ) u(τ ), lim τ → s 1 (t i ) t i u(τ ) -f r -1 u (t i ), x(r -1 u (t i )), x(t i ), lim τ → -s 1 (t i ) r -1 u (t i ) u(τ ), lim τ → -s 1 (t i ) t i u(τ ) • -s 1 (t i ) lim τ → -s 1 (t i ) r -1 u (t i ) φ(u(τ )) r -1 u (t i ) t i ∂φ ∂u (u(τ ))h(τ ) dτ (4.21)
where

s 1 (t) = sign - r -1 u (t
u(t) =      1 if t ∈ [-1; 0] 2 if t ∈ ]0; 0.5] 3 if t ∈ ]0.5; 1] and ∀t ∈ [0; 1], h(t) = 1
along with the following functional

J 0 (u) = 1 0 u(r u (t)) dt
where φ is the identity,

φ(u) = u t ru(t) u(τ ) dτ = 1
Then, we have

J 0 (u) = r -1 u (0) 0 1 dt + r -1 u (0.5) r -1 u (0) 2 dt + 1 r -1 u (0.5) 3 dt = 5 3 where r -1 u (0) = 0.5, r -1 u (0.5) = 5 6 If δ < 0, then r -1 u+δh (t) > r -1 u (t) and 1 0 u(r u+δh (t)) dt = r -1 u (0) 0 1 dt + r -1 u+δh (0) r -1 u (0) 1 dt + r -1 u (0.5) r -1 u+δh (0) 2 dt + r -1 u+δh (0.5) r -1 u (0.5) 2 dt + 1 r -1 u+δh (0.5)
3 dt and

J 0 (u + δh) -J 0 (u) = r -1 u+δh (0) r -1 u (0) (1 -2) dt + r -1 u+δh (0.5) r -1 u (0.5) (2 -3) dt + δ 1 0 h(r u (t)) dt + o(δ)
4.4. Input-dependant delays: a non-smooth problem 57

Taking s = -1 in (4.10) lim δ→ s 0 r -1 u+δh (t) -r -1 u (t) δ = - 1 lim τ → s r -1 u (t) u(τ ) r -1 u (t) t 1 dτ with s = sign -s r -1 u (t) t 1 dτ = 1 Finally lim δ→0 - J 0 (u + δh) -J 0 (u) δ = r -1 u (0) -0 2 + r -1 u (0.5) -0.5 3 + 1 2 = 1 4 + 1 9 + 1 2 = 31 36 (4.22)
Conversely, if δ > 0, then r -1 u+δh (t) < r -1 u (t) and

J 0 (u + δh) -J 0 (u) = r -1 u (0) r -1 u+δh (0) (2 -1) dt + r -1 u (0.5) r -1 u+δh (0.5) (3 -2) dt + δ 2 + o(δ)
and

lim δ→0 + J 0 (u + δh) -J 0 (u) δ = r -1 u (0) -0 3 + r -1 u (0.5) -0.5 3 + 1 2 = 1 6 + 1 9 + 1 2 = 7 9 (4.23)
There is indeed a mismatch between the left and the right limits (4.22) and (4.23).

As a consequence, P 0 is actually a non-smooth optimization problem and its optimal solutions cannot be characterized using the standard technique of imposing that all the variations of the augmented cost be equal to zero.

This result also has important practical consequences. Indeed, any standard optimization technique requiring first (or second) order regularity properties is expected to fail to properly solve P 0 .

Stationarity conditions of a regularization of the problem

To overcome the mathematical difficulty stressed by Theorem 4.2, we consider a regularized version of P 0 where the input u of the system having x as state is itself made to be the state of a pure integrator of an underlying input v.

Take (v 0 , u 0 ) ∈ L 2 ([r 0 ; 0], R p ) × D 1 ([r 0 ; 0], R p ), r 0 < 0 with 0 r 0 φ (u 0 (τ )) dτ = 1 and ∀t ∈ [r 0 ; 0], u 0 (t) = u 0 (0) + t 0 v 0 (τ ) dτ
Let P ∈ M p (R) be symmetric definite positive. The regularized optimal control problem is

P : min v T 0 L(t, x(t), u(t)) + 1 2 v(t) T P v(t) dt J(v) s.t. ∀t ∈ [0; T ], ẋ(t) = f (t, x(t), x(r u (t)), u(t), u(r u (t))) ∀t ∈ [0; T ], u(t) = v(t) x [r 0 ;0] = x 0 , u [r 0 ;0] = u 0 , v [r 0 ;0] = v 0
The stationarity conditions of P are the same as those of the free functional where the constraints have been adjoined

J(x, u, λ) = T 0 L(t, x(t), u(t)) -λ(t) T ẋ(t) -ν(t) T u(t) + ν(t) T v(t) + λ(t) T f (t, x(t), x(r u (t)), u(t), u(r u (t))) dt + ψ(x(T ))
Using two integrations by parts, one for x and one for u, we classically deduce

J(x, u, λ) = T 0 L(t, x(t), u(t)) + λ(t) T x(t) + ν(t) T u(t) + ν(t) T v(t) + λ(t) T f (t, x(t), x(r u (t)), u(t), u(r u (t))) dt -ν(T ) T u(T ) + ν(0) T u(0) + ψ(x(T ))
Consider the variation of J w.r.t. its second argument, u, which one state variable

D h J(u) = T 0 ∂L ∂u (t, x(t), u(t))h(t) + ν(t) T h(t) dt -ν(T ) T h(T ) + T 0 λ(t) T ∂f ∂x r (t, x(t), x(r u (t)), u(t), u(r u (t)))• f (t, x(t), x(r u (t)), u(t), u(r u (t))) φ(u(r u (t))) t ru(t) ∂φ ∂u (u(τ ))h(τ ) dτ dt + T 0 λ(t) T ∂f ∂u (t, x(t), x(r u (t)), u(t), u(r u (t)))h(t) dt + T 0 λ(t) T ∂f ∂u r (t, x(t), x(r u (t)), u(t), u(r u (t)))h(r u (t)) dt + T 0 λ(t) T ∂f ∂u r (t, x(t), x(r u (t)), u(t), u(r u (t)))• v(r u (t)) φ(u(r u (t))) t ru(t) ∂φ ∂u (u(τ ))h(τ ) dτ dt
Using a change of variable in the third integral and Fubini's theorem (see [START_REF] Dibenedetto | Real analysis[END_REF]) in the fourth and the fifth ones, we find that

D h J(u) = T 0 ∂L ∂u (t, x(t), u(t))h(t) + ν(t) T h(t) dt -ν(T ) T h(T ) + T 0 λ(t) T ∂f ∂u (t, x(t), x(r u (t)), u(t), u(r u (t)))h(t) dt + ru(T ) 0 λ(r -1 u (t)) T • ∂f ∂u r (r -1 u (t), x(r -1 u (t)), x(t), u(r -1 u (t)), u(t))(r -1 u ) (t)h(t) dt + T 0 r -1 u (min(τ,ru(T ))) τ λ(t) T ∂f ∂x r (t, x(t), x(r u (t)), u(t), u(r u (t)))• f (t, x(t), x(r u (t)), u(t), u(r u (t))) φ(u(r u (t))) dt ∂φ ∂u (u(τ ))h(τ ) dτ + T 0 r -1 u (min(τ,ru(T ))) τ λ(t) T ∂f ∂u r (t, x(t), x(r u (t)), u(t), u(r u (t)))• v(r u (t)) φ(u(r u (t))) dt ∂φ ∂u (u(τ ))h(τ ) dτ
Finally, Dubois-Reymond lemma gives the following necessary conditions2 , where two particular terms of interest, discussed later on, are highlighted in red

ν(t) = - ∂L ∂u (t, x(t), u(t)) T - ∂f ∂u (t, x(t), x(r u (t)), u(t), u(r u (t))) T λ(t) -1 [t 0 ;r(t 0 +T )] (t)(r -1 u ) (t)• ∂f ∂u r (r -1 u (t), x(r -1 u (t)), x(t), u(r -1 u (t)), u(t)) T • λ(r -1 (t)) - r -1 u (min(t,ru(T ))) t λ(τ ) T ∂f ∂x r (τ, x(τ ), x(r u (τ )), u(τ ), u(r u (τ )))• f (τ, x(τ ), x(r u (τ )), u(τ ), u(r u (τ ))) φ(u(r u (τ ))) dτ ∂φ ∂u (u(t)) T - r -1 u (min(t,ru(T ))) t λ(τ ) T ∂f ∂u r (τ, x(τ ), x(r u (τ )), u(τ ), u(r u (τ )))• v(r u (τ )) φ(u(r u (τ ))) dτ ∂φ ∂u (u(t)) T (4.24) ν(T ) = 0
Similarly calculating the Gâteaux derivatives with respect to x, λ and ν, we establish the following result Theorem 4.3:

The stationarity conditions of P are given by the following TPBVP ẋ(t) = f (t, x(t), x(r u (t)), u(t), u(r u (t))) (4.25)

x(0) = x 0 u(t) = v(t) (4.26) u [r 0 ;0] = u 0 λ(t) = - ∂L ∂x (t, x(t), u(t)) T - ∂f ∂x (t, x(t), x(r u (t)), u(t), u(r u (t))) T λ(t) -1 [t 0 ;ru(t 0 +T )] (t)(r -1 u ) (t)• ∂f ∂x r (r -1 u (t), x(r -1 u (t)), x(t), u(r -1 u (t)), u(t)) T • λ(r -1 (t)) (4.27) λ(T ) = ∂ψ ∂x (x(T )) T ν(t) = - ∂L ∂u (t, x(t), u(t)) T - ∂f ∂u (t, x(t), x(r u (t)), u(t), u(r u (t))) T λ(t) -1 [t 0 ;r(t 0 +T )] (t)(r -1 u ) (t)• ∂f ∂u r (r -1 u (t), x(r -1 u (t)), x(t), u(r -1 u (t)), u(t)) T • λ(r -1 (t)) - r -1 u (min(t,ru(T ))) t λ(τ ) T ∂f ∂x r (τ, x(τ ), x(r u (τ )), u(τ ), u(r u (τ )))• f (τ, x(τ ), x(r u (τ )), u(τ ), u(r u (τ ))) φ(u(r u (τ ))) dτ ∂φ ∂u (u(t)) T - r -1 u (min(t,ru(T ))) t λ(τ ) T ∂f ∂u r (τ, x(τ ), x(r u (τ )), u(τ ), u(r u (τ )))• v(r u (τ )) φ(u(r u (τ ))) dτ ∂φ ∂u (u(t)) T (4.28) ν(T ) = 0 0 = P v(t) + ν(t) (4.29)

Conclusions

In this chapter, we have derived the stationarity conditions of an optimal control problem subject to a fixed time-varying delay. The Gâteaux differentiability of the optimal control problem subject to an input-dependant hydraulic delay was investigated. We concluded that if the delayed input appears in the dynamics, the considered optimal control problem is not Gâteaux differentiable. We proposed a regularized version of this problem where the physical input is considered to be a state resulting from the pure integration of a new control input. The stationarity conditions of this problem were derived. Interestingly, these conditions only differ from the case of a fixed time delay case by the introduction of two distributed terms (emphasized in equation (4.24)) stemming from the sensitivity of the objective w.r.t. a change of delay law produced by the variation of the input.

Chapter 5

Iterative resolution algorithm for the regularization of an input-dependant hydraulic input delay

Résolution itérative du contrôle optimal régularisé de systèmes à retards hydrauliques dépendant de la commande. Dans ce chapitre, nous proposons un algorithme itératif permettant d'approcher la solution du problème de commande optimale de systèmes soumis à un retard hydraulique en résolvant une série de problèmes de commande optimal soumis avec des retards variables mais indépendants de la commande.

Une preuve de convergence de cet algorithme est établie et une application numérique est présentée.

Introduction

In this chapter we build an iterative algorithm to solve the problem of optimal control with hydraulic input-dependant input delays. Following the result of Chapter 4 showing that this problem is not Gâteaux differentiable, we consider its regularization where the input is transformed into a state resulting from the integration of a new underlying auxiliary control. In Section 5.3, we use the stationarity conditions derived in Chapter 4 to propose an iterative algorithm solving, instead of the original problem, a sequence of simpler modified problems with fixed time-varying delays. A detailed convergence analysis is carried out, showing that in a limit case, this algorithm behaves like a gradient descent. Practically, this approach is attractive as it allows to use straightforwardly the classic tools of direct simultaneous optimization to solve the intermediate problems. This is explained in Section 5.3. In Section 5.4, a numerical example is treated to illustrate the interest of this method.

Notations

We use the notations of Chapter 4. In addition, given T > 0, n ∈ N * and ∈ L 2 ([0; T ], R n ), we denote 1 and 2 the norm 1 and 2 of the function . For convenience, we recall Chapter 5. Iterative resolution algorithm that, using Cauchy-Schwarz inequality,

1 ≤ √ nT 2 Let g ∈ L 2 (R, R n
) and I be an interval of R. We denote g I the restriction of g to I.

Let us denote u ∈ R p the control and x ∈ R m the state of the system.

Optimization algorithm

Let T > 0, x 0 ∈ R m and P ∈ M p (R) be symmetric definite positive. Let φ :

R p → R * + , L : [0; T ] × R m × R p → R and f : [0; T ] × R m × R p × R p → R m be smooth functions. Take (v 0 , u 0 ) ∈ L 2 ([r 0 ; 0], R p ) × D 1 ([r 0 ; 0], R p ), r 0 < 0 with 0 r 0 φ (u 0 (τ )) dτ = 1 and ∀t ∈ [r 0 ; 0], u 0 (t) = u 0 (0) + t 0 v 0 (τ ) dτ (5.1)
Consider the following optimization problem, defined as a simplified case of the general problem covered in Theorem 4.3 (no terminal cost, no dependency on the past values of the state), without loss of generality

P : min v T 0 L(t, x(t), u(t)) + 1 2 v(t) T P v(t) dt J(v) s.t. ∀t ∈ [0; T ], ẋ(t) = f (t, x(t), u(t), u(r u (t))) ∀t ∈ [0; T ], u = v x(0) = x 0 , u [r 0 ;0] = u 0 , v [r 0 ;0] = v 0
where r u is implicitly defined by the relation

t ru(t) φ(u(τ )) dτ = 1 (5.2)
and in particular r 0 = r u (0)

Let us consider the operator P :

L 2 ([0; T ], R p ) → D 1 ([0; T ], R p ) × D 1 ([0; T ], R m ) 3 such that P(v) = (u, x, λ, ν) is defined, according to Theorem 4.3, by u(t) = v(t), u [r 0 ;0] = u 0 (5.3) ẋ(t) = f (t, x(t), u(t), u(r u (t))), x(0) = x 0 (5.4) λ(t) = - ∂L ∂x (t, x(t), u(t)) T - ∂f ∂x (t, x(t), u n (t), u(r u (t))) T λ(t) (5.5) λ(T ) = 0 ν(t) = - ∂L ∂u (t, x(t), u(t)) T - ∂f ∂u (t, x(t), u(t), u(r u (t))) T • λ(t) -1 [0;ru(T )] (t)(r -1 u ) (t)• ∂f ∂u r (r -1 u (t), x(r -1 u (t)), u(r -1 u (t)), u(t)) T • λ(r -1 u (t)) 5.3. Optimization algorithm 63 - r -1 u (min(t,ru(T ))) t λ(τ ) T • ∂f ∂u (τ, x(τ ), u(τ ), u(r u (τ ))) v(r u (τ )) φ(u(r u (τ ))) dτ ∂φ ∂u (u(t)) (5.6) ν(T ) = 0
Using these notations, the stationarity conditions of P are given by (u, x, λ, ν) = P(v)

P v + ν = 0 (5.7)
Solving P directly is difficult. Instead, we would rather solve a sequence of simpler auxiliary problems (P n ), such that for all n ≥ 1

P n+1 : min v n+1 T 0 L(t, X n+1 (t), u n+1 (t)) + 1 2 v n+1 (t) T P v n+1 (t) + S n (t)(u n+1 (t) -u n (t)) + α 2 v n+1 (t) -v n (t) 2 2 dt s.t. ∀t ∈ [0; T ], Ẋn+1 = f (t, X n+1 (t), u n+1 (t), u n+1 (r un (t))) ∀t ∈ [0; T ], un+1 = v n+1 X n+1 (0) = x 0 , u n+1 [r 0 ;0] = u 0 , v n+1 [ru(0);0] = v 0 where S n (t) = r -1 un (min(t,ru n (T ))) t λ n (τ ) T • ∂f ∂u (τ, x n (τ ), u n (τ ), u n (r un (τ ))) v n (r un (τ )) φ(u n (r un (τ ))) dτ ∂φ ∂u (u n (t))
is one of the two extra terms highlighted in Section 4.4.2 and is the sensitivity of the objective with respect to the change of the delay law caused by a change of the control input as derived from the calculus of variations. In the definition of S n and the general statement of P n , (u n , x n , λ n , ν n ) are defined as

v n → (u n , x n , λ n , ν n ) P(v n )
Throughout the rest of the discussion, the following assumptions are considered Assumption 3. L is twice continuously differentiable while f , φ are continuously differentiable, and there exists K ≥ 0 such that

∀(t, x, u) ∈ [0; T ] × R m × R p , ∇ 2 L(t, x, u) 1 ≤ K and ∀(t, x, u, u r ) ∈ [0; T ] × R m × R p × R p , ∇f (t, x, u, u r ) 1 ≤ K and ∀u ∈ R p , ∇φ(u) 1 ≤ K and, ∇ 2 L, ∇f , ∇φ are K-Lipschitz continuous.
Assumption 4. There exists

J * ∈ R such that ∀v ∈ L 2 ([0; T ]), J * ≤ J(v) λ(T ) = 0
Note the slight (but important) differences between P defined by (5.3)-(5.6) and Q. The second argument of Q is used to define the time-varying delay appearing in the righthand side of equations (5.13)-(5.14). Based on Assumption 3 and Cauchy existence and uniqueness theorem, Q is clearly defined. Given the couples of arguments (v 1 , w 1 ), (v 2 , w 2 ) and (v, w), we define

(u 1 , q 1 , x 1 , λ 1 ) = Q(v 1 , w 1 ), r 1 r q 1 (u 2 , q 2 , x 2 , λ 2 ) = Q(v 2 , w 2 ), r 2 r q 2 and (u, q, x, λ) = Q(v, w), r r q
which are used to formulate the subsequent lemma Lemma 5.1 (Lipschitz continuity of Q). The two following inequalities hold

∀t ∈ [0; T ], u 2 (t) -u 1 (t) 1 ≤ √ pt v 2 -v 1 2 (5.15) ∀t ∈ [0; T ], u(t) -u 0 1 ≤ √ pt v 2 (5.16)
There exists some positive parameters

(k 1 , k 2 , k 3 , k 4 ), (l 1 , l 2 , l 3 , l 4 ) independent of α such that ∀t ∈ [0; T ], x 2 (t) -x 1 (t) 1 ≤ k 1 v 2 -v 1 2 + k 2 (1 + v 1 2 )• (1 + w 1 2 + w 2 2 ) w 2 -w 1 2 (5.17) ∀t ∈ [0; T ], x(t) -x 0 1 ≤ k 3 + k 4 v 2 (5.18) ∀t ∈ [0; T ], λ 2 (t) -λ 1 (t) 1 ≤ l 1 (1 + v 1 2 ) v 2 -v 1 2 + l 2 • (1 + w 1 2 + w 2 2 )(1 + v 1 2 ) 2 w 2 -w 1 2 (5.19) ∀t ∈ [0; T ], λ(t) 1 ≤ l 3 + l 4 v 2 (5.20) Proof. See Appendix B.1.
The newly defined operator Q plays a key role with respect to the sequence (v n ). Indeed, the stationarity conditions of P n+1 are given by

(u n+1 , X n+1 , Λ n+1 ) = Q(v n+1 , v n ) Ṅn+1 (t) = - ∂L ∂u (t, X n+1 (t), u n+1 (t)) T - ∂f ∂u (t, X n+1 (t), u n+1 (t), u n+1 (r un (t))) T Λ n+1 (t) -1 [0;ru n (T )] (t)(r -1 un ) (t)• ∂f ∂u r (r -1 un (t), X n+1 (r -1 un (t)), u n+1 (r -1 un (t)), u n+1 (t)) T • Λ n+1 (r -1 un (t)) -S n (t) T 0 = P v n+1 + N n+1 + α(v n+1 -v n ) N n+1 (T ) = 0 (5.21) 1 α g vn , n+1 + 1 0 g vn+(v n+1 -vn)s -g vn , v n+1 -v n ds Finally J(v n+1 ) -J(v n ) ≤ - 1 α g vn 2 2 + 1 α g vn 2 n+1 2 + 1 0 g vn+(v n+1 -vn)s -g vn 2 v n+1 -v n 2 ds (5.23)
In order to go further into the convergence analysis, we need to establish a bound for n+1 2 . This is given in the following proposition

Proposition 5.1:

There exists some positive parameters (κ 1 , κ 2 , κ 3 , κ 4 ) independent of α such that such that

∀t ∈ [0; T ], N n+1 (t) -ν n (t) 1 ≤ (κ 1 + κ 2 v n 2 ) v n+1 -v n 2 (5.24) and ∀t ∈ [0; T ], ν n (t) 1 ≤ κ 3 + κ 4 v n 2 (5.25) Proof. See Appendix B.2.
Recalling (5.22), we have

v n+1 -v n 2 ≤ 1 α ( g vn 2 + P 2 v n+1 -v n 2 + N n+1 -ν n 2 )
5.3. Optimization algorithm 67 Then, using (5.24)

v n+1 -v n 2 ≤ 1 α g vn 2 + 1 α ( P 2 + κ 1 + κ 2 v n 2 ) v n+1 -v n 2
As a consequence, if

P 2 + κ 1 + κ 2 R n < α (5.26)
(which is always possible for α large enough as the left-hand side of (5.26) is independent of α). We find that

v n+1 -v n 2 ≤ 1 α -P 2 -κ 1 -κ 2 R n g vn 2 (5.27)
In particular, we deduce an a priori bound on the norm of v n+1 , using (5.25)

v n+1 2 ≤ R n + P 2 R n + κ 3 + κ 4 R n α -P 2 -κ 1 -κ 2 R n
Incidentally, this also leads to

n+1 2 ≤ P 2 + κ 1 + κ 2 R n α -P 2 -κ 1 -κ 2 R n g vn 2 (5.28)
To go further into the analysis of (5.23), we now have to prove the Lipschitz continuity of g v = P v + ν with respect to v. In order to do this, consider (v 1 , v 2 ) and the associated functions (u 1 , x 1 , λ 1 , ν 1 ) and (u 2 , x 2 , λ 2 , ν 2 ) such that (u 1 , x 1 , λ 1 , ν 1 ) P(v 1 ) and (u 2 , x 2 , λ 2 , ν 2 ) P(v 2 )

Proposition 5.2:

There exists a continuous function K : R 3 + → R + increasing with each of its arguments and independent of α such that

∀t ∈ [0; T ], ν 2 (t) -ν 1 (t) 1 ≤ K( v 0 2 , v 1 2 , v 2 2 ) v 2 -v 1 2
(5.29)

Proof. See Appendix B.3
One can then investigate further the decrease of cost formulated in (5.23). Using (5.28) and (5.29), ones gets

J(v n+1 ) -J(v n ) ≤ - 1 α (1 - P 2 + κ 1 + κ 2 R n α -P 2 -κ 1 -κ 2 R n ) g vn 2 2 + F n (α) 2 v n+1 -v n 2 2
where

F n (α) = K( v 0 2 , R n , R n + P 2 R n + κ 3 + κ 4 R n α -P 2 -κ 1 -κ 2 R n ) + P 2
Chapter 5. Iterative resolution algorithm Then, using (5.27)

J(v n+1 ) -J(v n ) ≤ - 1 α 1 - P 2 + κ 1 + κ 2 R n α -P 2 -κ 1 -κ 2 R n - αF n (α) 2(α -P 2 -κ 1 -κ 2 R n ) 2 g vn 2 2
Since F n is a decreasing function of α, there exists a value of α large enough such that 1 -

P 2 + κ 1 + κ 2 R n α -P 2 -κ 1 -κ 2 R n - αF n (α) 2(α -P 2 -κ 1 -κ 2 R n ) 2 C(α, R n ) > 0
(5.30) and J(v n+1 ) -J(v n ) < 0. In particular, this guarantees that v n+1 ∈ X n . By induction, this implies that if one picks a value α = α 1 such that (α 1 , R 1 ) satisfy (5.26) and (5.30), then for all rank n, v n ∈ X 1 and (5.26) and (5.30) hold. Then, for all rank n

∀n ∈ N * , J(v n+1 ) -J(v n ) ≤ - C(α 1 , R 1 ) α 1 g vn 2 2
This leads to

N i=0 g v i 2 2 ≤ α 1 C(α 1 , R 1 ) (J(v 0 ) -J(v n+1 ))
Finally we derive

N i=0 g v i 2 2 ≤ α 1 C(α 1 , R 1 ) (J(v 0 ) -J * )
and lim n→∞ g vn 2 = 0 which gives the conclusion.

Numerical example

We now illustrate the solution method studied in Theorem 5.1 using an example studied in [START_REF] Bresch-Pietri | Prediction-based stabilization of linear systems subject to input-dependent input delay of integral-type[END_REF]. Consider a second order unstable linear system with dynamics given by

ẍ(t) -ẋ(t) + x(t) = u(r u (t)) u(t) = v(t)
(5.31) having the following initial conditions

x(0) = 1, ẋ(0) = 0 u [r 0 ;0] = 1, v [r 0 ;0] = 0
This can equivalently be recast as

Ẋ(t) = AX(t) + Bu(r u (t)) u(t) = v(t)
5.4. Numerical example 69 where

X = x ẋ , A = 0 1 -1 1 , B = 0 1
The optimal control problem is

P : min v T 0 x(t) -x r 2 2 + w u u(t) -u r 2 2 + w v v(t) 2 2 dt s.t. Ẋ(t) = AX(t) + Bu(r u (t)) u(t) = v(t)
with T = 10, w u = 0.1, w v = 0.1 and x r = u r = 1.5. Given α = 51 , we approach iteratively a solution of P by constructing an α-admissible sequence (v n ). We pick the trivial initialization value v 1 = 0 and for all n ≥ 1 apply the following algorithm :

• given v n , compute u n and the delay law r n

• compute (x n , λ n ) and deduce S n

• solve P n+1 and obtain v n+1

The flowchart of the algorithm is presented on Figure 5.1. Practically, the resolution of P n+1 is performed using a direct collocation transcription method (see [Biegler]), with AMPL as algebraic modelling language and IPOPT 3.11.8 as NLP solver. The time horizon is divided into 100 finite elements of equal size, each of them containing 3 Radau collocation points. λ n is numerically estimated by solving Pn : min

v=vn T 0 x n (t) -x r 2 2 + w u u n (t) -u r 2 2 + w v v(t) 2 2 dt s.t. ẋn (t) = Ax n (t) + Bu n (r un (t)) un (t) = v(t)
and retrieving the adjoint of the discretized problem. The algorithm is terminated when we reach 100 iterations. More sophisticated termination criteria could be considered, for instance based on the satisfaction by the current solution of the stationarity conditions. The results are presented on Figures 5.2-5.5. Figures 5.2-5.3 report the optimal trajectory that is computed and the associated delay law. Figure 5.4 pictorially shows how this trajectory is progressively approached by the sequence. Figure 5.5 exhibits some indicators regarding the convergence properties of the algorithm : the cost J along with the relative steps size measured by log (∆v) log(

vn-v n-1 2 vn 2
) and log (∆J) log(

Jn-J n-1 2 Jn 2
) at successive iterations. As expected, the cost decreases monotonically and the linear shape of the cost decrease on the semi-log plot is evocative of a first order steepest descent-like method. The total computation time for the first 100 iterations displayed on Figure 5.5 is equal to 10.78 seconds, 8.07 seconds being actually spent in the solver.

Conclusions

In this chapter, we have proposed an iterative algorithm to solve the problem of optimal control of systems with hydraulic input-dependent input delays. A convergence proof was detailed and numerical results were given illustrating the practical interest.

A straightforward extension would be to extend the calculus of variations and the deduced iterative optimization algorithm to the case of systems with hydraulic inputdependent state delays. This case is of importance since it is instrumental in the modelling of recycling loops or cascades of reacting units.

The derivation of a method on a second order approximation of the stationarity conditions would be interesting. This would lead to an algorithm analogous to a trustregion method. where δ α > 0 is an increment parameter. It should be possible to develop more advanced selection procedures based on the theory of trust-region algorithms. However, we have not investigated these methods yet.

Further, we will consider in this chapter problems P subject to constraints. These constraints will directly be incorporated into the definition of problems P n and Pn . This is beyond the scope of Chapter 5 but, while a full theoretical analysis is challenging, a practical extension is possible, as we will show numerically.

The flowchart associated to these updates of the optimization algorithm is presented on Figure 6.1.

6.3. Relaxation to a non-regularized input for input-dependant input delays 77

Relaxation to a non-regularized input for inputdependant input delays

Approximate regularization of the control

In this section, we investigate the possibility of applying a variant of the algorithm laid out in Section 6.2 to a case where the input of the system under consideration has not been explicitly regularized into a new state. Instead, we consider an approximate regularization of the original problem still framed in terms of the u variable using a numerical approximation of its derivative as regularizing term

P : min u T 0 L(t, x(t), u(t)) + w ∆u • ψ 1 (u, t) 2 2 dt J(u) s.t. ẋ(t) = f (t, x(t), u(t), u(r u (t))) x(0) = x 0 , u [r 0 ;0] = u 0
where ψ 1 (u, .) is an approximation of the function u defined as

ψ 1 (u, t) = u(t n ) -u(t n-1 ) t n -t n-1 , t ∈ [t n-1 ; t n [ (6.2)
Hence, the sequence of problems that we solve is

P n+1 : min u n+1 T 0 L(t, X n+1 (t), u n+1 (t)) + w ∆u • ψ 1 (u n+1 , t) 2 2 + S n (t)(u n+1 (t) -u n (t)) + α 2 u n+1 (t) -u n (t) 2 2 dt s.t. Ẋn+1 (t) = f (t, X n+1 (t), u n+1 (t), u n+1 (r un (t)))
X n+1 (0) = x 0 , u n+1 [r 0 ;0] = u 0 with (according to the developments of Chapter 5)

S n (t) = r -1 un (min(t,ru n (T ))) t λ n (τ ) T • ∂f ∂u (τ, x n (τ ), u n (τ ), u n (r un (τ ))) ψ 2 (u n , r un (τ )) φ(u n (r un (τ ))) dτ ∂φ ∂u (u n (t))
where (x n , λ n ) are defined using the notations of Chapter 5 and ψ 2 (u, .) is another approximation of the function u defined as

ψ 2 (u, t) = u(t n+1 ) -u(t n-1 ) 2(t n+1 -t n-1 ) , t ∈ [t n-1 ; t n [ (6.3)
In light of our previous investigations, the problems P n are not differentiable. Practically, we will manage this issue (along with the problems related to the numerical approximation of u) by formulating the problem in such a way to enforce approximate continuity of u. Three reasons motivate this choice despite its hurdles :

• it avoids introducing new states and thus restrains the dimension of the problem

• if the problem involves state constraints, it avoids artificially increasing their index (also refereed to as relative degree, e.g. in [START_REF] Bryson | Applied Optimal Control[END_REF] or [START_REF] Isidori | Nonlinear Control Systems[END_REF]), which could lead to numerical difficulties Chapter 6. Examples Semi-log plot of cost and input step-sizes The results reported in on Figures 6.5-6.6. The cost decrease and the overall convergence of the sequence is greatly improved. Approximate convergence is obtained after 20 iterations (representing about 12 seconds of computation). An improvement by a factor of approximately 2 is hence obtained compared to the previous case. Interestingly in both cases where w ∆u = 1 or w ∆u = 10, almost half of the total CPU time used by the algorithm is actually spent in the AMPL script performing the various side tasks presented in Figure 6.1 rather than in the optimization solver itself. This means that the numerical implementation could still be greatly improved by the use of a more efficient computing environment.

Case with soft constraints

We now consider an updated scenario where the the following soft constraint is considered

y b 2 ≤ y max 2 = 0.5
Obviously, if we translate this rigorously as a constraint, the problem does not have any admissible solution since y b 2 (0) > 0.5. Instead, we are going to enforce the constraint by penalizing its infringement. Thus, using the approach presented in (6.9), we consider the Semi-log plot of cost and input step-sizes

L(t, y b (t), u(t)) = w y • y b (t) -y r 2 2 + w u • u(t) -u r 2 2 + w p • p s (y b 2 (t), y max 2 )
with w y = 2, w u = 0.15, w ∆u = 10, = 10 -7 . The results are reported in Figures 6.7-6.8. As can be seen, during the transient phase the controller imposes higher injection ratios of product 3. This allows to dilute the inescapable inlet of 2 which is already in the pipe at the initial pipe and curtails the constraint infringement. However, this also leads to increasing the overall delay and leads to a longer "flushing" phase. Approximate convergence is reached around 20 iterations (representing about 15 seconds of computation).

Case with soft and hard constraints

Additionally, a hard constraint enforcing y b 1 ≥ 0.13 is introduced. As can be seen on Figures 6.9-6.10, the hard constraints on the concentration of product 1 does not allow us to bring down the concentration of product 2 as much as in the previous case. Indeed, in the first phase of the transient (when the dead-volume has not yet been flushed out of its .9: Optimal trajectory with soft and hard constraints. The controller increases the injection ratio of 3 in the first part of the batch to lower the injection of 2 and to diminish the constraint violation. This is limited by the respect of the hard constraint on product 1 initial content), the concentration of product 1 and 2 at the outlet of the unit verifies the following algebraic relation imposed by the initial composition of the product in the pipe

∀t ∈ 0; r -1 u (0) , y 2 (t) = 4y 1 (t)
Approximate convergence is reached around 20 iterations (representing about 15 seconds of computation).

Extension to input-dependant state delays

Formally, the principle of the algorithm can directly be extended to the case of a system with state delays.

u n+1 T 0 L(t, X n+1 (t), u n+1 (t)) + S n (t)(u n+1 (t) -u n (t)) + α 2 u n+1 (t) -u n (t) 2 2 dt s.t. Ẋn+1 (t) = f (t, X n+1 (t), X n+1 (r un (t)), u n+1 (t)) X n+1 [r 0 ;0] = x 0
where

S n = r -1 un (min(t,ru n (T ))) t λ n (τ ) T ∂f ∂x r (τ, x n (τ ), x n (r un (τ )), u n (τ ))• f (τ, x n (τ ), x n (r un (τ )), u n (τ )) φ(u n (r un (τ ))) dτ ∂φ ∂u (u n (t))
and could apply straightforwardly the algorithm of Figure 6.1 to attempt to solve P.

Similarly, a slight modification of the problem under consideration in the calculus of variations accounting for a delayed term in the running cost would also allow to tackle the system studied in Chapter 3

P = min q T 0 (T (r q (t)) -T r (t)) 2 + w • (q(t) -q r (t)) 2 dt s.t. Ṫ = f (T, q) T (t ≤ 0) = T 0 q min ≤ q(t) ≤ q max
Practically, the sequence (P n+1 ) would be defined as

P n+1 : min u n+1 T 0 (T n+1 (r qn (t)) -T r (t)) 2 + w • (q n+1 (t) -q r (t)) 2 + S n (t)(q n+1 (t) -q n (t)) + α 2 (q n+1 (t) -q n (t)) 2 dt s.t. Ṫn+1 = f (T n+1 , q n+1 ) T (t ≤ 0) = T 0 q min ≤ q(t) ≤ q max where S n (t) = 2 r -1 qn (min(t,rq n (T ))) t (T (r qn (τ )) -T r (τ )) • f (T n (r qn (τ )), q n (r qn (τ )))
q n (r qn (τ )) dτ

Conclusions

In this chapter, we have presented a practical extension of the algorithm studied in Chapter 5. In particular, we showed on a numerical example that it allows to handle constrained problems that were not treated in our theoretical analysis. We studied the possibility to treat a problem subject to an hydraulic input delay when the control is not rigorously regularized. We showed that if some practical precautions are taken, the algorithm provides good numerical results. We also sketched a possible extension of the algorithm to the case of a system subject to a state delay.

Chapter 7 Conclusion

Conclusion. Dans ce chapitre, nous résumons les contributions de la thèse et proposons des axes pour poursuivre le travail de recherche.

In this thesis, we have addressed several problems related to the control and optimization of systems subject to varying delays. Early on in our work, we started by studying the influence of uncertain, varying metrology delays on the stability of a simple IMC. We showed in Chapter 2 that this can yield closed-loop instability and provided explicit guarantees on the controller gain insuring robust asymptotic convergence. This result confirms the classic fact of experience according to which the robust operation of controllers in presence of uncompensated delays requires some detuning of the gain, and thus loss of performance.

Later on, we tackled the problem of dynamic optimization of a system subject to complex but structured input-dependant delays: the so-called hydraulic delays. In Chapter 3, we showed that, while it can yield sound results on many problems, the state-of-the-art to handle this issue in direct simultaneous methods is not fully satisfactory. Indeed, it requires to discard the delay equation and to replace it by some discretization of a mathematically equivalent transport PDE. While doable, this leads to a less parsimonious representation of the problem and various numerical challenges. Based on this statement, we carried in Chapter 4 the calculus of variations of this problem to get further insight into the structure of the optimal solutions. After a lengthy derivation, this establishes the important result that the straightforward formulation one could consider for dynamic optimization of a system subject to input-dependant input delays is ill-posed in the sense that it does not yield a smooth problem. Following this, we established the stationarity conditions of a regularization of this problem. Interestingly, this work highlighted the difference between optimizing a system subject to a fixed time-varying delay or an explicitly input-dependant one. We used this to derive in Chapter 5 an iterative algorithm which only requires to solve a sequence of auxiliary problems with fixed time varying delay laws. This is of practical importance as it then allows us to use state-of-the-art simultaneous optimization methods in the resolution of each auxiliary problem. Our convergence analysis showed that, similarly to a trust region method, our algorithm becomes equivalent to a gradient descent in the limit where the allowed step-size goes to zero. Finally, in Chapter 6 we extend the scope of our algorithm to problems with input and state constraints as well as systems where the control sequence is only approximately regularized. We demonstrated where a(s) min(r 1 (s), r 2 (s)) and b(s) max(r 1 (s), r 2 (s)). Since r 1 and r 2 are strictly increasing functions, a and b also are and they are invertible. From their respective definitions, it is also clear that a(t) ≤ b(t) and a(0) = b(0) = r 0 . Then, using Fubini's theorem On the other hand, ∀i ∈ {1, 2}, using the Lipschitz continuity of φ, (5.1), (5.12) and integrating q either backward or forward, we find r i (t) = φ(q i (t)) φ(q i (r i (t))) Using Grönwall's lemma

x 2 (t) -x 1 (t) 1 ≤ φ(u 0 (0)) + K( √ pT ( w 1 2 + w 2 2 ) + √ pr 0 v 0 2 ) φ 2 min • 2K 2 T pT ( v 1 1 + v 0 1 ) w 2 -w 1 2 + 2KT pT v 2 -v 1 2 e Kt Synthetically, a conservative estimate is as follows

x 2 (t) -x 1 (t) 1 ≤ k 1 v 2 -v 1 2 + k 2 (1 + w 1 2 + w 2 2 )(1 + v 1 2 ) w 2 -w 1 2
We also have ∀t ∈ [0; T ], x(t) -x 0 1 ≤ t 0 f (τ, x(τ ), u(τ ), u(r q (τ ))) 1 dτ Using the Lipschitz continuity of f x(t) -x 0 1 ≤ t 0 f (0, x 0 , u 0 (0), u 0 (r 0 )) 1 + Kτ + K x(τ ) -x 0 1 + K u(τ ) -u 0 (0) 1 + K u(r(τ )) -u 0 (r 0 ) 1 dτ With Gröwall's lemma and (B.2), we find x(t) -x 0 1 ≤ T ( f (0, x 0 , u 0 (0), u 0 (r 0 ))

1 + KT 2 + 2K pT v 2 + K √ pr 0 v 0 2 )e KT
Finally, using the various Lipschitz continuity results established in Lemma 5.1, we find

N n+1 (t) -ν n (t) 1 ≤KT pT + k 1 + 2l 1 (1 + v n 2 ) + 2(l 3 + l 3 v n 2 )(2 pT + k 1 ) v n+1 -v n 2
This can be rewritten as

N n+1 (t) -ν n (t) 1 ≤ (κ 1 + κ 2 v n 2 ) v n+1 -v n 2
It is also straightforward to show that, for some positive constants κ 3 , κ 4 , one has -∂f ∂u (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) + 1 [0;r 2 (T )] (τ )(r -1 2 ) (τ )• ∂f ∂u r (r -1 2 (τ ), x 2 (r -1 2 (τ )), u 2 (r -1 2 (τ )), u 2 (τ )) T λ 2 (r -1 2 (τ ))

ν n (t
-1 [0;r 1 (T )] (τ )(r -1 1 ) (τ )• ∂f ∂u r (r -1 1 (τ ), x 1 (r -1 1 (τ )), u 1 (r -1 1 (τ )), u 1 (τ )) T λ 1 (r -1 1 (τ )) ∂f ∂u r (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) dτ 1

We get (a -1 (s) -b -1 (s)))

ν 2 (t) -ν 1 (t) 1 ≤ K(T -t) sup τ ∈[0;T ] ( x 2 (τ ) -x 1 (τ ) 1 + u 2 (τ ) -u 1 (τ ) 1 + 2 λ 2 (τ ) -λ 1 (τ )
• (l 3 + l 4 v 2 2 ) K 2 φ min pT (φ(u 0 (0)) + K √ pT v 2 2 ) φ min ( v 2 2 + v 0 2 )
Finally, after a few lines of calculus similar to (B.4), we get

B ≤ 2(l 3 + l 4 v 2 2 ) K 2 φ min pT (φ(u 0 (0)) + K √ pT v 2 2 ) + K √ pr 0 v 0 2 φ min ( v 2 2 + v 0 2 )• KT √ pT (φ(u 0 (0)) + K √ pT v 1 2 + K √ pr 0 v 0 2 ) φ 2 min v 2 -v 1 2
Using the same kind of computations on C, we show that ∀t ∈ [0; T ], ν 2 (t) -ν 1 (t)

1 ≤ K( v 0 2 , v 1 2 , v 2 2 ) v 2 -v 1 2
where K : R 3 + → R + is a continuous function such that for all i x i ≤ z i =⇒ K(x 1 , x 2 , x 3 ) ≤ K(z 1 , z 2 , z 3 ) which gives the conclusion Appendix C

Reporting of the numerical results

Variables Constraints CPU time In a first part, we studied the closed-loop control of a system subject to varying and uncertain metrology delays. We established new results on robust stability under explicit conditions on the controller gain. In a second part, we tackled the problem of the dynamic optimization of systems exhibiting input dependent delays due to transport phenomena in complex hydraulic architectures. We designed an iterative optimization algorithm and guaranteed its convergence through a detailed analysis.

  38). The results of these simulations are shown on Figure 2.4.

Figure 3

 3 Figure 3.2: Optimal trajectories, w = 0.1

  Figure3.5: Optimal trajectories with full PDE discretization, w q = 0.001, N = 40, CPU time = 142.17 s. The "heat and flush" behaviour is also found. High accuracy of the optimizer prediction is achieved

initialize v n , u n integrate r n , r - 1 nFigure 5 . 2 :Figure 5 Figure 5 . 4 :

 152554 Figure 5.1: Flowchart of the algorithm

Figure 6 . 3 :

 63 Figure 6.3: Optimal trajectory with minimal regularization, w ∆u = 1. The control displays a strong initial discontinuity

Figure 6 . 5 :

 65 Figure 6.5: Optimal trajectory with increased regularization of the control, w ∆u = 10

Figure 6 . 7 :

 67 Figure 6.7: Optimal trajectory with soft constraint. The controller increases the injection ratio of product 3 in the first part of the batch to lower the injection of 2 and diminish the constraint violation compared to Figure 6.5

  Figure6.9: Optimal trajectory with soft and hard constraints. The controller increases the injection ratio of 3 in the first part of the batch to lower the injection of 2 and to diminish the constraint violation. This is limited by the respect of the hard constraint on product 1

t 0 u 1 u 1

 011 (r 2 (τ )) -u 1 (r 1 (τ )) 1 dτ ≤ (r 2 (τ )) -u 1 (r 1 (τ )) 1 dτ ≤ ( sup s∈[r 0 ;a(t)]

(a - 1 u 1

 11 (s) -b -1 (s)) + sup s∈[a(t);b(t)] (t -b -1 (s))) • ( v 1 1 + v 0 1 )wherev 0 1 is used to denote v (r 2 (τ )) -u 1 (r 1 (τ )) 1 dτ ≤ sup s∈[r 0 ;a(t)]

(a - 1

 1 (s) -b -1 (s)) + a -1 (a(t)) -b -1 (a(t)) ( v 1 1 + v 0 1 ) For any s ∈ [r 0 ; a(t)] a -1 (s) -b -1 (s) = y 2 -y 1where y 1 and y 2 are uniquely defined by s = a(y 2 ) = b(y 1 )

0 u 1 1 ) 0 q 2 pT φ min w 2 -w 1 2 t 0 u 1 0 x 2

 0110220102 ≥ φ min φ(u 0 (0)) + K( √ pT w i 2 + √ pr 0 v 0 2 )a is a scalar function whose rate of change is lower bounded by the minimum of the two expressions of the previous equation. As a consequenceφ min φ(u 0 (0)) + K( √ pT ( w 1 2 + w 2 2 ) + √ pr 0 v 0 2 ) |y 2 -y 1 | ≤ |a(y 2 ) -a(y 1 )| (B.4) Then, since |a(y 2 ) -a(y 1 )| = |a(y 1 ) -b(y 1 )| = |r 2 (y 1 ) -r 1 (y 1 )|, one has t (r 2 (τ )) -u 1 (r 1 (τ )) 1 dτ ≤ 2 φ(u 0 (0)) + K( √ pT ( w 1 2 + w 2 2 ) + √ pr 0 v 0 2 ) φ min • sup s∈[0;T ] |r 2 (s) -r 1 (s)| • ( v 1 1 + v 0 2 (τ ))dτ | ≤ K t (τ ) -q 1 (τ ) 1 dτAnd then|r 2 (t) -r 1 (t)| ≤ Kt √ Then (r 2 (τ )) -u 1 (r 1 (τ )) 1 dτ ≤ 2KT pT φ(u 0 (0)) + K( √ pT ( w 1 2 + w 2 2 ) + √ pr 0 v 0 2 ) φ 2 min • ( v 1 1 + v 0 1 ) w 2 -w 1 2 (B.5)Substituting in (B.3), this leads tox 2 (t) -x 1 (t) 1 ≤ K t (τ ) -x 1 (τ ) 1 dτ + 2KT pT v 2 -v 1 2 + φ(u 0 (0)) + K( √ pT ( w 1 2 + w 2 2 ) + √ pr 0 v 0 2 ) φ 2 min • 2K 2 T pT ( v 1 1 + v 0 1 ) w 2 -w 1 2

+ r - 1 2

 1 (min(τ,r 2 (T ))) τ λ 2 (s) T • ∂f ∂u r (s, x 2 (s), u 2 (s), u 2 (r 2 (s))) v 2 (r 2 (s)) φ(u 2 (r 2 (s)x 1 (s), u 1 (s), u 1 (r 1 (s))) v 1 (r 1 (s)) φ(u 1 (r 1 (s))) ds ∂φ ∂u (u 1 (τ )) T dτ 1Hence, using a change of variable in two of the integrals above, and after a Cauchy-Schwarz inequality, one getsν 2 (t) -ν 1 (t) 1 ≤ T t ∂L ∂u (τ, x 2 (τ ), u 2 (τ )) T -∂L ∂u (τ, x 1 (τ ), u 1 (τ )) T B.3. Proof of Proposition 5.2 113 + ∂f ∂u (τ, x 2 (τ ), u 2 (τ ), u 2 (r 2 (τ ))) T λ 2 (τ ) -∂f ∂u (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) dτ 1 + T r -1 2 (min(t,r 2 (T ))) ∂f ∂u r (τ, x 2 (τ ), u 2 (τ ), u 2 (r 2 (τ ))) T λ 2 (τ ) dτ -T r -1 1 (min(t,r 1 (T ))) ∂f ∂u r (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) dτ + x 2 (s), u 2 (s), u 2 (r 2 (s))) v 2 (r 2 (s)) φ(u 2 (r 2 (s))) ds ∂φ ∂u (u 2 (τ )) T r -1 1 (min(τ,r 1 (T ))) τ λ 1 (s) T • ∂f ∂u r (s, x 1 (s), u 1 (s), u 1 (r 1 (s))) v 1 (r 1 (s)) φ(u 1 (r 1 (s))) ds ∂φ ∂u (u 1 (τ )) T dτ 1 Noting that T r -1 2 (min(t,r 2 (T ))) ∂f ∂u r (τ, x 2 (τ ), u 2 (τ ), u 2 (r 2 (τ ))) T λ 2 (τ ) dτ -T r -1 1 (min(t,r 1 (T ))) ∂f ∂u r (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) dτ 1 x 2 (τ ), u 2 (τ ), u 2 (r 2 (τ ))) T λ 2 (τ )-∂f ∂u r (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) 1 dτ + r -1 2 (min(t,r 2 (T )))r -1 1 (min(t,r 1 (T )))

1 )+ 0 u 1

 101 2K(T -t)(l 3 + l 4 v 1 2 )• sup τ ∈[0;T ] ( x 2 (τ ) -x 1 (τ ) 1 + 2 u 2 (τ ) -u 1 (τ ) 1 ) + T (r 2 (t)) -u 1 (r 1 (t)) 1 dτ + r -1 2 (min(t,r 2 (T ))) r -1 1 (min(t,r 1 (T ))) ∂f ∂u r (τ, x 1 (τ ), u 1 (τ ), u 1 (r 1 (τ ))) T λ 1 (τ ) dτ 1 A + T t r -1 2 (min(τ,r 2 (T )))r -1 1 (min(τ,r 1 (T ))) λ 2 (s)

(a - 1 0 v 2 ≤

 102 (s) -t) + sup s∈[b(t);T ] (a -1 (s) -b -1 (s)))• (l 3 + l 4 v 2 2 ) K 2 φ min T 0 v 2 (r 2 (s)) 1 dsBesides, by the Cauchy-Schwarz inequalityT (r 2 (s)) 1 ds ≤ pT pT (φ(u 0 (0)) + K √ pT v 2 2 + K √ pr 0 v 0 2 ) φ min ( v 2 2 + v 0 2 ) and since a ≤ b B ≤ (a -1 (b(t)) -b -1 (b(t)) + sup s∈[b(t);T ]

Figure

  Figure C.1: Numerical results of Chapter 3

  , "Run-to-run control with nonlinearity and delay uncertainty", in Proc. of the 11th IFAC Symposium on Dynamics and Control of Process Systems, including Biosystems, DYCOPS-CAB 2016

• C.-H. Clerget, J.-P. Grimaldi, M. Chèbre, and N. Petit, "Optimization of dynamical systems with time-varying or input-varying delays", in Proc. of the 55th IEEE Conference on Decision and Control 2016 • C.-H. Clerget, N. Petit, and L. T. Biegler "Dynamic optimization of a system with input-dependant time delays", in Proc. of the Foundations of Computer Aided Process Operations / Chemical Process Control FOCAPO/CPC 2017

  • C.-H. Clerget, J.-P. Grimaldi, M. Chèbre, and N. Petit, "Run-to-run control with nonlinearity and delay uncertainty", in Proc. of the 11th IFAC Symposium on Clerget, N. Petit, and L. T. Biegler "Dynamic optimization of a system with input-dependant time delays", in Proc. of the Foundations of Computer Aided Process Operations / Chemical Process Control FOCAPO/CPC 2017

	• M. Chèbre, N. Petit, C.-H. Clerget, and J.-P. Grimaldi "Scalable integrated solution
	for real time estimation, control and optimization of the quality of fuels manufactured
	in refineries: an industrial story", in Proc. of the IFAC 2017 World Congress
	• C.-H. Clerget, J.-P. Grimaldi, M. Chèbre, N. Petit, "An example of robust internal
	model control under variable and uncertain delay", in Journal of Process Control,
	2017 (in press)

Dynamics and Control of Process Systems, including Biosystems, DYCOPS-CAB 2016 • C.-H. Clerget, J.-P. Grimaldi, M. Chèbre, and N. Petit, "Optimization of dynamical systems with time-varying or input-varying delays", in Proc. of the 55th IEEE Conference on Decision and Control 2016 Chapter 1. Introduction • C.-H.

  )

	t	∂φ ∂u	(u(τ ))h(τ ) dτ
	Using (4.20), we can formulate the following result
	Theorem 4.2:

J0 is Gâteaux differentiable w.r.t. its second argument at point (x, u, λ) in direction h if and only if

  )

	t	∂φ ∂u	(u(τ ))h(τ ) dτ
	Remark 4.1. Note that (4.21) does not trivially hold. This fact means that when f
	is explicitly depending upon its delayed input, the augmented cost associated cannot be
	guaranteed to be differentiable with respect to the input at any given point if the function u
	is not continuous for t > 0 (interestingly, however, discontinuities in the control history
	do not raise issues) and counter examples are straightforward to build (see Remark 4.2).
	Remark 4.2. Consider for instance (u, h) such that

  ) 1 ≤ κ 3 + κ 4 v n 2

B.3 Proof of Proposition 5.2

Proof. From (5.6), one has

ν 2 (t) -ν 1 (t) 1 ≤ T t ∂L ∂u (τ, x 2 (τ ), u 2 (τ )) T -∂L ∂u (τ, x 1 (τ ), u 1 (τ )) T + ∂f ∂u (τ, x 2 (τ ), u 2 (τ ), u 2 (r 2 (τ ))) T λ 2 (τ )

  1: Numerical results of Chapter 3 Variables Constraints Total CPU time Solve CPU time Approx. conv. Table C.1: Numerical results of Chapter 5 and 6 for 100 iterationsRésuméDans cette thèse, nous avons étudié le contrôle et l'optimisation de systèmes dynamiques sujets à des retards variables. L'existence de retards, de commande ou d'état, est un problème classique en automatique, susceptible de réduire les performances du système en régime transitoire, voire de remettre en cause la stabilité de contrôleurs en boucle fermée. De tels phénomènes de retards variables jouent un rôle important dans de nombreuses applications en génie des procédés.Dans une première partie, nous avons étudié la régulation en boucle fermée d'un système soumis à des retards de métrologie variables et incertains. Nous avons établi de nouveaux résultats garantissant la stabilité robuste sous certaines conditions explicites sur le gain du contrôleur. Dans une seconde partie, nous avons abordé le problème de l'optimisation dynamique de systèmes présentant des retards variables dépendant de la commande liés à des phénomènes de transport dans des réseaux hydrauliques. Nous avons proposé un algorithme itératif d'optimisation et garanti sa convergence grâce à une analyse détaillée.

				Abstract
				This Ph.D. work studied the control and
				optimization of dynamical systems subject to
				varying time delays.
				State and control time delays are a well-
				known problem in control theory, with a
				potential to decrease performances during
				transient regimes, or even to jeopardize
				controllers closed-loop stability. Such variable
				delays play a key role in many applications in
				process industries.
	Figure 5.2	2375	2300	10.78 s	8.07 s	≈ 8 s
	Figure 6.4	7493	7600	59.21 s	31.98 s	≈ 25 s
	Figure 6.6	7493	7600	54.70 s	28.57 s	≈ 12 s
	Figure 6.8	7493	7600	72.50 s	39.21 s	≈ 15 s
	Figure 6.10	7493	8200	75.25 s	41.71 s	≈ 15 s
	Mots Clés
	Optimisation dynamique, Retards variables,
	Contrôle de procédés

In practice, it can result from the analysis of sensitivity look-up tables obtained from experiments and derivation of interpolating models.

Because r depends on t only, its definition being independent from variations of x or u, we say that the system is subject to a fixed time varying delay.

On the interval [0; r u (T )] covered by the indicator function, the function r -1 u employed in (4.24) is well defined.

2 +

This value was chosen using a trial and error approach, knowing a priori that some large enough value of α will actually provide convergence.

Remerciements

Chapter 5. Iterative resolution algorithm

Assumption 5. There exists φ min > 0 such that ∀u ∈ R, φ min ≤ φ(u) Remark 1. Assumptions 3-4 are classically considered in the optimization literature. Assumption 5 is usually considered for systems with input varying delays of hydraulic type [DBP] so that r u be bounded away from zero and the input keep on reaching the plant. Definion 5.1. Given α ≥ 0, a sequence (v n ) n∈N * is called α-admissible if for all n ≥ 2, v n is a solution (possibly local) of P n .

Let us define

the set of L 2 functions such that their J-level set is included in a ball of L 2 and note

The main result concerning the sequence (P n ) is as follows Proof. Given n ∈ N * , let us assume that v n ∈ X (which is true for n = 1 by assumption) and, by extension of (5.8), define

which is a bounded set in the sense of the L 2 norm, i.e. there exists R n > 0 such that ∀v ∈ X n , v 2 ≤ R n (5.10)

Chapter 5. Iterative resolution algorithm

From this, we directly deduce that the solutions of P n and P n+1 are related by

In turn, the cost variation between v n and v n+1 is given by

Using the adjoint state method (e.g. [START_REF] Strang | Computational Science and Engineering[END_REF]), one computes, after a few lines of calculus,

Chapter 6

Examples

Exemples. Dans ce chapitre, nous présentons les résultats de l'application d'une variante de l'algorithme présenté au Chapitre 5 sur un problème d'intérêt industriel. Ce chapitre ne propose pas d'analyses théoriques mais met en évidence la possibilité pratique d'étendre l'algorithme au-delà du cadre dans lequel nous avons établi sa convergence (en particulier en présence de contraintes).

Introduction

In this chapter, we present the numerical results of a variant of the optimization algorithm presented in Chapter 5. No further convergence analysis is performed (although this could deserve further investigations). Chiefly, this variant is developed to account for constraints.

It is applied to the dynamic optimization of a mixing system of industrial interest. A generalization to the case of systems with state delays is presented.

Desirable features for a practical implementation of the algorithm

The optimization algorithm developed in Chapter 5 is not directly suitable for practical implementation. Indeed, in this analysis, the tuning parameter α was considered constant over all steps (in particular, Theorem 5.1 established the convergence properties of the algorithm under this assumption). However, practically, keeping α constant is not a judicious choice. Usually, a desirable behaviour of the cost decrease should consist of two phases. In a first phase, the cost decreases quickly as achieving effective cost decrease at each step does not require a large value of α. At this point, α should be kept as large as possible to speed up convergence. In a second phase, the cost decrease is much slower and requires a higher α to generate a cost decrease at each step. Indeed, similarly to a trust-region method (see e.g. [START_REF] Nocedal | Numerical Optimization[END_REF]), the choice of α implicitly selects the step-size. Then, we should start with a small value of the parameter α and steadily increase it in order to insure a continued cost decrease after each iteration. As a consequence, we will allow α to vary between steps according to the following very simple transition rule

Figure 6.2: Schematic of the mixing unit

• it allows us to stay closer to the original formulation of the problem and is handy, for instance, to formulate hard hydraulic constraints on the admissible variations of u between two time steps

The approach is now illustrated on one example of practical interest.

Batch control of a mixing unit with pre-blend, the "paint" problem

In this section, we consider an industrial unit mixing three products (see [START_REF] Petit | Motion planning for two classes of nonlinear systems with delays depending on the control[END_REF] and [START_REF] Clerget | Optimization of dynamical systems with time-varying or input-varying delays[END_REF]). The products 1 and 2 are pre-blended before going through a dead-volume V pip . At the outlet, they are blended with product 3 and immediately reach a storage tank. The total flow rate, F , of the unit is fixed and we control the ratios of the different products (also referred to as recipe), u, that are injected at a given time. Our goal is to control the average composition of the outlet tank, y b . We note y the instantaneous composition of the product entering the tank and V the instantaneous volume of product in the outlet tank. A schematic view of the system can be found on Figure 6.2. Under a plug-flow assumption, the equations and constraints governing the system are given by

The optimization problems we are going to formulate encompass several elements:

• a reference recipe to follow, u r

• a target of composition of the outlet tank to track, y r

• soft inequality constraints that should be enforced but may be violated, at a cost on the solution, to avoid infeasibility. They are transcribed using a further smooth quadratic penalty term in the objective function (see [START_REF] Fiacco | Nonlinear programming: Sequential unconstrained minimization techniques[END_REF]). For example, if one wishes to enforce some inequality a ≤ b, we will add to the cost function a penalty proportional to

where is a small regularizing term

• hard inequality constraints that may not be violated by the system and that are explicitly formulated as constraints in the problem

In our simulations, we will take the following setting : F = 4 m 3 .s -1 , V 0 = 100 m 3 , V pip = 200 m 3 . Furthermore, we consider a transition from an initial state y 0 = ( 1 6 , 4 6 , 1 6 ), u(t < 0) = ( 16 , 4 6 , 1 6 ) to y r = (0.5, 0.25, 0.25), u r = (0.5, 0.25, 0.25). All ratios are subject to the following hard bound constraints 0.01 ≤ u i ≤ 1 (6.10)

Pure tracking problem

In this first scenario, we do not consider any soft or hard constraint other that the bound constraints (6.10). The optimization cost function is defined by

The resolution sequence is initialized taking α = 0.3, δ α = 0.5 and initial values y b = y 0 , u = u 0 . We first consider a case where w y = 2, w u = 0.15, w ∆u = 1. The results are presented on Figures 6.3-6.4. As is apparent from Figure 6.4, the convergence properties achieved by the algorithm are relatively poor. This is due to the strong discontinuity of the optimal input at the beginning of the time horizon that makes the problem non-differentiable. The oscillation of the solutions makes the termination of the algorithm difficult, one could however consider the sequence to have approximately converged after roughly 40 iterations (representing about 25 seconds of computation).

This problem can however be mitigated by further regularizing the control sequence. This can be alternatively achieved by increasing the regularizing cost w ∆u or by imposing bounds on the amplitude of the variations of the control between adjacent points of the collocation mesh. For this reason, we now consider the case where w ∆u = 10 and we impose ∀i ∈ 1; 3 , ∀n ≥ 0,

Chapter 7. Conclusion the good numerical performances of our formulation on a problem of industrial interest.

A full numerical benchmark between our method and the existing approaches remains, however, to be conducted. While this was not the core topic of this thesis, the case of input-dependant state delays is of great practical importance in many process applications (recycling loops and cascade of reacting units for instance) and should be treated. Interestingly, the analysis conducted in Chapter 4 shows that a problem where only delayed state values (as opposed to delayed input values) appear in the dynamics is smooth in the sense that it is Gâteaux differentiable. The results of our calculus of variations lead us to think that a straightforward extension of our algorithm is possible in this case.

An exciting problem raised by our work is the possibility of extending our approach to a second order method. This would indeed greatly improve convergence performances. It is however not clear that the problem has sufficient regularity to directly allow such an extension. Establishing or disproving this would require the computation of the problem second variation. Proof. Using Cauchy-Schwarz inequality

Appendices

and similarly

We also have

Let us define µ : t → λ(T -t). Then, integrating backwards, one gets

Using the Lipschitz continuity of ∂L ∂x , the boundedness of ∂f ∂x , (B.1) and (B.6), we find

We deduce that the norm of the adjoint state is bounded

We also have

And, reusing (B.5), we get

Again, using Grönwall's lemma, one finds

B.2 Proof of Proposition 5.1

Proof. We have

Using a change of variables in the second integral, we find

ds ∂φ ∂u (u 2 (τ ))

-λ 1 (s) T • ∂f ∂u r (s, x 1 (s), u 1 (s), u 1 (r 1 (s))) v 1 (r 1 (s)) φ(u 1 (r 1 (s))) ds ∂φ ∂u (u 1 (τ )) dτ 1

C

We have A ≤ r -1 2 (min(t,r 2 (T )))

r -1 1 (min(t,r 1 (T )))

The definition of the delay (5.2) gives us Then, using (B.8) and performing the same calculation (r -1 2 (min(t, r 2 (T ))) -r -1 1 (min(t, r 1 (T )))) ≤

Finally, we have

To treat B, we note a : [r 0 ; max(r 1 (T ), r 2 (T ))] → [0; T ] t → min(r -1 1 (min(τ, r 1 (T ))), r -1 2 (min(τ, r 2 (T )))) and b : [r 0 ; min(r 1 (T ), r 2 (T ))] → [0; T ] t → max(r -1 1 (min(τ, r 1 (T ))), r -1 2 (min(τ, r 2 (T ))))