domain: spi phys math info

text:
Sergio Pérez-Roca
email: sergio.perezroca@onera.fr

Julien Marzat
email: julien.marzat@onera.fr

Hélène Piet-Lahanier
email: helene.piet-lahanier@onera.fr

Nicolas Langlois
email: nicolas.langlois@esigelec.fr

François Farago
email: francois.farago@cnes.fr

Marco Galeotta
email: marco.galeotta@cnes.fr

Serge Le Gonidec
email: serge.le-gonidec@ariane.group

Trajectory planning and tracking via MPC for transient control of liquid-propellant rocket engines

Keywords: Liquid-propellant rocket engines, model predictive and optimisationbased control, control of constrained systems, tracking, trajectory planning, robustness

The new scenarios associated with launchers reusability force the enhancement of control algorithms for their liquid-propellant rocket engines. The transient phases of these engines are generally executed in open loop. The goal of this paper is to improve the control performance and robustness throughout the fully continuous phase of the start-up transient of a generic gas-generator cycle. The controller has to guarantee an accurate tracking in terms of combustion-chamber pressure and chambers mixture ratios, as well as to satisfy a set of hard operational constraints. The selected strategy comprises a nonlinear preprocessor and a linearised MPC (Model-Predictive Control) controller, making use of nonlinear state-space models of the engine. The former plans the reference trajectory of states and control, which is tracked by the latter. Control goals are attained with sufficient accuracy while verifying constraints within the desired throttling range. Robustness to internal parameters variations is considered in the MPC controller by means of an epigraph formulation of the minimax robust optimisation problem, where a finite set of parameter-variation scenarios is treated.

Introduction

The current context of launcher vehicles design is strongly related to the reusability feature. From the automatic control perspective, this potential need for reusable liquid-propellant rocket engines (LPRE) is translated into stricter robustness requirements, mainly forced by their multi-restart and thrust-modulation capabilities. These demanding requirements arise from the possible endogenous perturbations due to components faults or evolving parameters and from exogenous perturbations related to the more complex mission profiles envisaged by new launchers. Multivariable-control studies of main-stage LPRE have attained a short throttling envelope (70%-120%) in test benches [START_REF] Gonidec | Automatic & Control applications in the European space propulsion domain. From need expression to preparation for an uncertain future[END_REF]. In flight, the control system is generally conceived to attain the nominal operating point. One of the design capabilities of the future European PROMETHEUS engine is to throttle down to 30% of thrust [START_REF] Baiocco | Technology demonstration for reusable launchers[END_REF]. Hence, it seems necessary to enlarge the controlled operating domain. Tracking and robustness have to be maintained at those low throttle levels, where physical phenomena are more difficult to anticipate. The main control problem in these multivariable systems is the tracking of setpoints in combustion-chamber pressure and mixture ratio, whose references stem from launcher needs. Control-valves opening angles are varied so as to adjust engine's operating point while respecting a series of constraints. Most of the control approaches in the literature make use of linearised models about operating points for developing steady-state controllers, commonly based on PID techniques (such as [START_REF] Nemeth | Reusable rocket engine intelligent control system framework design, phase 2[END_REF]). In general terms, initial MIMO (Multi Input Multi Output) systems are decoupled into dominant SISO (Single Input Single Output) subsystems. Other conventional closed-loop (CL) linear strategies considered multivariable state feedback [START_REF] Zhang | State-space analysis of the dynamic characteristics of a variable thrust liquid propellant rocket engine[END_REF]. Off-line optimisation studies have also been developed [START_REF] Dai | Damage-Mitigating Control of a Reusable Rocket Engine: Part II-Formulation of an Optimal Policy[END_REF]. More advanced works identified in the literature, incorporating some nonlinear [START_REF] Lorenzo | Nonlinear control of a reusable rocket engine for life extension[END_REF], hybrid [START_REF] Musgrave | Real-time accommodation of actuator faults on a reusable rocket engine[END_REF] or robust [START_REF] Saudemont | Study of a robust control law based on H infinity for the Vulcain rocket engine -Etude d'une commande robuste a base de commande H infinity pour le moteur Vulcain[END_REF] techniques, improve certain aspects of performance and robustness. Besides, regarding the treatment of component faults, some strategies for reconfiguration control have been proposed [START_REF] Musgrave | Real-time accommodation of actuator faults on a reusable rocket engine[END_REF][START_REF] Sarotte | Actuator Fault Tolerant System For Cryogenic Combustion Bench Cooling Circuit[END_REF]. To the best of our knowledge, there are no publications which deal with not only the steady state but also the demanding transient phases at the same level of performance and robustness [START_REF] Pérez-Roca | A survey of automatic control methods for liquid-propellant rocket engines[END_REF]. Pre-defined sequences of engine operation (start-up and shutdown), are commonly performed in open loop with narrow adaptation margins. The first part of the start-up transient presents a succession of discrete events including valves openings and chambers ignitions. Once these commands have all been activated, the second part of the transient, which is completely continuous, takes place until the steady state is reached. The main reasons for carrying out open-loop (OL) control in the initial (discrete-event) phase, exposed in [START_REF] Nemeth | Reusable rocket engine intelligent control system framework design, phase 2[END_REF], are controllability and observability problems at very low mass flows. Transient control via valves becomes possible once all events have taken place. This observation has also been considered in this paper, where only the second part of the start-up transient, fully continuous, is controlled. The main goal of this work is to control the start-up transient of a pump-fed LPRE. Concretely, accurate tracking of combustion pressure (related to thrust) and mixture ratio all along the transient is targeted. At the same time, a set of hard operational constraints has to be satisfied, mainly concerning mixture ratios, turbopumps rotational speeds and valves actuators angular velocities. The type of LPRE studied is a gas-generator-cycle engine, but the methodology is intended to be generally applicable to other cycles. As in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF], the control loop presented in this paper is based on Model Predictive Control (MPC), which is accompanied by a preprocessor for reference-trajectory generation. The MPC method has been chosen as the most adequate for this type of complex systems with hard operational constraints, as introduced in the next sections. In fact, it is gaining popularity at academic and industrial levels and can take into account robustness [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] or hybrid aspects, which are relevant for future work on this topic. Regarding the generation of reference trajectories, the algorithm in this paper is also based on MPC, which is common in the literature [START_REF] Kuwata | Robust Constrained Receding Horizon Control for Trajectory Planning[END_REF][START_REF] Singh | Trajectory generation for a UAV in urban terrain, using nonlinear MPC[END_REF]. Other trajectory-planning methods have been proposed, especially in the field of robotic motion, but they are also mostly based on optimisation [START_REF] Betts | Survey of Numerical Methods for Trajectory Optimization[END_REF]. This paper is organised as follows. Section 2 serves as a recapitulation of modelling considerations. The state-space system considered in the subsequent sections is stated there. Section 3 explains the control strategy carried out, mainly based on MPC techniques. The main results are presented and analysed in Section 4. Finally, Section 5 serves as a conclusion.

Modelling

The modelling approach considered in this paper was first introduced in [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF] and further explained and updated in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF]. Basically, a set of different models is used, characterised by diverse degrees of complexity according to their purposes. In order to represent the plant to be controlled, a simulator of the engine was first developed. In this simulator, constructed component-wise, the main thermo-fluid-dynamics of LPRE elements are modelled: mass, energy and momentum conservation equations. Besides, since the application considered in this paper, representative of the Vulcain 1 engine, is based on a gas-generator (GG) cycle, turbopumps dynamics are also modelled via shaft mechanical conservation equations. In these cycles, the hot-gas flow required for spinning turbines comes from a GG, which is a small combustion chamber receiving a portion of the main propellant flow. Fig. 1 depicts the main components of the engine and lists the main acronyms used. Actuators in this system are mainly five continuously-controllable valves (VCH, VCO, VGH, VGO and VGC). There are also two discrete actuators: one binary igniter (i CC ) and one binary starter (i GG ). Nevertheless, discrete inputs are deemed active in this paper since the goal is to manage the continuouslycontrolled part ot transients (up from 1.5s after start command). Valves opening angles (α), having a nonlinear and monotone relation to sections (A), control the flows to the main combustion chamber (VCH and VCO), to the GG (VGH, VGO), and to the oxidiser turbine (VGC). Their fully open position is 90➦. The latter is the main contributor to determining mixture ratio (M R), defined as the quotient between oxidiser and fuel mass flow rates M R = ṁox / ṁfu . This ratio is a key performance indicator in LPRE and is set at three levels: at an engine's global level (M R P I ), accounting for pumped propellants; in the combustion chamber (M R CC ) and in the GG (M R GG ). The simulator was rewritten as a nonlinear state-space (NLSS) model by joining components equations symbolically. As explained in [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF][START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF], throughout that translation several simplifications, comprising physical assumptions and mathematical rewriting, were carried out until obtaining the here-called simplified NLSS such that ẋ = f s (x, u). The number of states is n = 12 and m = 5 is the number of control inputs. The state vector x contains the two turbopumps speeds ω H and ω O , the four pressures in the system (p CC of combustion chamber, p GG of the GG, p LT H for hydrogen-turbine inlet cavity and p V GC for oxygen-turbine inlet cavity) and six mass flows, including the ones streaming through control valves ( ṁV CH , ṁV CO , ṁV GH , ṁV GO and ṁV GC ) and the one streaming through the hydrogen-turbine inlet pipe ṁLT H .

x = [ω H ω O p CC p GG p LT H p V GC ṁLT H ṁV CH ṁV CO ṁV GH ṁV GO ṁV GC ] T . (1)
The states with greater tracking importance are incorporated into a reduced state vector x z :

x z = [p CC ṁV CH ṁV CO ṁV GH ṁV GO ] T . (2)
The control input u comprises the sections of the five control valves:

u = [A V CH A V CO A V GH A V GO A V GC ] T . (3)
All equations, states and control have been rendered non-dimensional with respect to the nominal equilibrium values. The form of the dynamic system ẋ = f s (x, u) is the following, in which a i , b i , ..., h i ∈ R + are internal parameters and w t is an exogenous input corresponding to starter mass flow:

ẋ1 = -a 1 x 5 x 1 + b 1 x 5 + (-c 1 x 8 -d 1 x 10 ) x 1 + e 1 x 8 2 + f 1 x 8 x 10 , (4) ẋ2
= -a 2 x 2 2 -b 2 x 6 x 2 + c 2 x 6 -(d 2 x 9 + e 2 x 11 ) x 2 + f 2 x 9 2 + g 2 x 9 x 11 , (5) ẋ3
= a 3 x 8 + b 3 x 9 -c 3 x 3 , (6)
ẋ4 = a 4 x 10 + b 4 x 11 -c 4 (x 7 + x 12 ) + d 4 w t , (7) ẋ5
= a 5 x 7 -b 5 x 5 , (8) ẋ6
= a 6 x 12 -b 6 x 6 , (9) ẋ7
= a 7 (x 4 -x 5 ) - b 7 x 7 2 x 4 , (10) ẋ8
= a 8 x 1 2 -b 8 x 8 2 -c 8 x 10 x 8 -d 8 x 3 + e 8 u 1 2 -f 8 x 8 2 (g 8 u 1 + h 8 ) u 1 , (11) ẋ9
= a 9 x 2 2 -b 9 x 9 x 2 -c 9 x 9 2 -d 9 x 3 + e 9 u 2 - f 9 x 9 2 u 2 , (12) ẋ10
= x 4 -x 6 - a 12 x 12 2 x 4 u 5 2 u 5 . (14) ẋ12
It is noticeable that the system presents numerous nonlinearities and is non-affine with respect to control. Besides, full-state estimation is considered as perfect in this paper. This assumption is realistic for ω and p. Nevertheless, measuring some mass flows is generally inconvenient in practice. Hence, discussion on estimation will be relevant in future work.

Controller design

The synthesised controller is composed of two main blocks. First, a preprocessing block serves to generate state and control trajectories for arriving to a constant steady-state reference, x r , up from some given initial conditions. The main component is the MPC block, which receives those full-state and control reference trajectories and drives the system to them while meeting constraints. At the end of the start-up transient, the state should track x r , with a particular interest in having a small tracking error in x z . Simultaneously, a set of hard constraints on x and u has to be satisfied throughout the transient. This second goal is major during these operation phases in order to avoid excessive temperatures, p or ω along the operation. The entire control loop is depicted in Fig. 2.

Fig. 2. Control-loop diagram

The remaining elements in that loop are the following. On the right-hand side the rocket engine is simulated (complex simulator), at an integration time step of 10 -5 s in order to capture fast dynamics and to be robust to numerical stiffness. The inputs of that simulator and of the state-space model considered for control are valve sections u. However, the actuators model (internal valve actuators) is expressed in terms of α. Consequently, a conversion block defining static and monotone nonlinear relations is required. The MPC controller provides valve sections that are then converted into angles. Valve actuators are considered as a separate element since they consist in an internal servo-loop, in which the angular position of the valve is regulated via a hydraulic or electrical actuator, modelled as a second-order system. This is a simplified modelling assumption since internal phenomena such as hysteresis and solid friction are not taken into account.

Preprocessor block: final-reference and trajectory generation

The preprocessing block serves as an off-line reference generator for the MPC controller. First of all, the end reference has to be constructed because the set of reference commands derived from launcher needs (p CC,r , M R P I,r , M R CC,r and M R GG,r ) is not sufficient to provide a full-state target equilibrium point x r to the engine controller. This reference generation is presented in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF], in which differential equations ( 4) to [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF] are equated to zero. A steady-state reference in terms of control inputs u r is also obtained in this calculation, since it is required for defining trajectories and for linearising the system. Once these end targets (x r ,u r ) have been computed, a reference start-up trajectory (X t ,U t ) from some given x 0 and u 0 can be built. X and U are defined in general as the series of x and u at each time step k along a horizon N (valid throughout the whole paper):

X = [x 1 , ..., x k , ..., x Np ] T , U = [u 1 , ..., u k , ..., u Nu ] T . ( 16
)
To do so, an optimisation-and model-predictive-based scheme is used. It can be regarded as an OL finite-horizon MPC scheme in which the prediction horizon is set to cover the duration of the start-up buildup transient, considered between 1.5s, the end of the discrete sequential phase, and 2.5s, the time when it is aimed at attaining the reference in the studied engine. This algorithm is based on the minimisation of a classical quadratic cost function J OL , defined as:

J OL (X t , U t ) =   N p,OL k=1 ∆x T k Q OL ∆x k + N u,OL k=1 ∆u T k R OL ∆u k )   ∆t, (17)
where ∆x k = x kx r and ∆u k = u ku r are the variables to minimise, that is to say, the distances with respect to the reference equilibrium point. Q OL = I n and R OL are the weight matrices associated to states and control respectively. Diagonal terms in R OL are set to 10 10 so as to minimise control action. N p,OL and N u,OL are states and control prediction horizons, which in this case are taken equal to the horizon (1s) over the discretisation time ∆t = 10ms, the latter selected to comply with controller's computational restrictions. Concerning the dynamics considered to predict the behaviour at each time step k, as announced in Section 2, the system is highly nonlinear. Hence, neglecting nonlinear dynamics can lead to relevant prediction errors, especially at points far from the equilibrium one. However, the main repercussion of imposing nonlinear constraints in optimisation problems is generally the loss of convexity of the optimised function and hence the increase in resolution complexity. The compromise chosen in this paper, related to the specific behaviour of the system, is the inclusion of nonlinear dynamic constraints until the system approaches its reference values to within 90%. This coincides approximately with the first half of the transient, where modelling errors of linearisation would be relevant if linear dynamics were used. Concretely, the aforementioned f s is discretised via an Euler implicit scheme:

x k+1 = x k + f s (x k+1 , u k+1 )∆t, k ∈ [1, N p,OL 90% ]. ( 18
)
This scheme has been selected since it is the most stable among the first-order integration methods, required for lowering the complexity of the optimisation by reducing the interdependencies between decision variables. Once the pressure p CC (x 3 or x z,1 ) attains its reference value, x z states are forced to be equal to the end-reference values and only linear dynamics is imposed, in order to simplify the optimisation. A few time steps before the start of those end-state constraints, a smooth transition between nonlinear and linear dynamics is set via maximum slope constraints. Linear dynamics (A c , B c ) stems from the linearisation of f s about (x r , u r ), which are then discretised via zero-order hold at ∆t (A d , B d ):

∆x k+1 = A d (x r , u r )∆x k + B d (x r , u r )∆u k , k ∈ [N p,OL 90% + 1, N p,OL + 1]. ( 19
)
Having defined the different dynamics, the optimisation algorithm which is executed once for the whole horizon (OL trajectory planning) under constraints is the following:

min Xt,Ut J OL (X t , U t ) (20) s.t. X t ∈ X, U t ∈ U A ineq [X t U t ] T ≤ b ineq x k+1 ≤ x k + f s (x k+1 , u k+1 )∆t + ε, ∀k ∈ [1, N p,OL 90% ] x k+1 ≥ x k + f s (x k+1 , u k+1 )∆t -ε, ∀k ∈ [1, N p,OL 90% ] A eq [∆X t ∆U t ] T = b eq (including x Np+1 = x r ).
X and U are the allowable sets for states and control (compact subsets of R n(N p,OL +1) and R mN u,OL respectively). The set U for the first control u 1 is specially constrained to comply with actuators capacity as in [START_REF] Luo | Model predictive dynamic control allocation with actuator dynamics[END_REF]:

u 1 ∈ [max(U , u 0 -umax ∆t), min(U , u 0 + umax ∆t)], ( 21
)
where umax is the maximum sectional velocity of valves.

Regarding the rest of constraints, (20) contains linear inequality constraints (defined by A ineq and b ineq ), for satisfying M R and actuators sectional-velocity bounds, as well as for defining a monotonically increasing pressure buildup. Nonlinear dynamic constraints are not defined as strict equality constraints, but are treated as inequalities with a small margin ε = 10 -2 × 1 n×1 (non-dimensional) so as to simplify the computation of a feasible optimal solution. Linear dynamics [START_REF] Sarotte | Actuator Fault Tolerant System For Cryogenic Combustion Bench Cooling Circuit[END_REF], initial conditions and end-state reaching are considered in the equality constraints (defined by A eq and b eq ). The end-state hard constraint forces the trajectory to precisely finish at the desired point. Fig. 3 depicts the set of generated trajectories for the different operating points (70%, 100%, 120% of thrust) in terms of combustion-chamber pressures. The interior-point optimisation software IPOPT [START_REF] Waechter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] is used to solve [START_REF] Saudemont | Study of a robust control law based on H infinity for the Vulcain rocket engine -Etude d'une commande robuste a base de commande H infinity pour le moteur Vulcain[END_REF]. Due to the inclusion of nonlinear constraints in this quadratic-cost optimisation problem, the solution might not be global. Initial conditions are fixed according to the end of the discrete part of the transient (until t = 1.5s). The first control inputs (VCO and VCH), related to the main chamber, are pretty constrained in terms of operation, thereby imposing straight trajectories.

MPC algorithm for trajectory tracking

MPC predicts the future behaviour of the modelled system along a horizon, and optimises control inputs according to a cost function related to a reference trajectory or to an end state. In this paper, the strategy is to track the predefined trajectory (X t , U t ). A similar algorithm, directly tracking an end state, was presented in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF]. However, tracking a trajectory enables a more precise control of transients, which is advantageous in this application. Dynamics is predicted in this controller in a linearised way, so as to enable a real-time control computation. The discrete-time matrices A d and B d are evaluated about each step in the trajectory, improving the prediction representativeness with respect to the use of single end-state matrices.

∆x k+1 = A d (x t,k , u t,k )∆x k + B d (x t,k , u t,k )∆u k . (22)
Thus, the goal of the controller is to drive X and U to X t and U t . The matrix A d is always stable along the trajectory, which is a particularity of these GGcycle LPRE models. The approach selected here incorporates the quasi-infinite horizon (QIH) notions from [START_REF] Chen | A Quasi-Infinite Horizon Nonlinear Model Predictive Control Scheme with Guaranteed Stability[END_REF]. That reference presents proofs for guaranteed stability and end-state reachability of MPC via the consideration of a terminal region and a fictitious local controller K acting at the end of the state prediction horizon, N p + 1. The main role of that controller appears in the computation of the P matrix of a Lyapunov function V (x) = ∆x T P ∆x, by solving the following Lyapunov equation:

(A K (x t,Np+1 , u t,Nu ) + κI) T P + P (A K (x t,Np+1 , u t,Nu ) + κI) = -Q K -K T R K K. (23)
In [START_REF] Zhang | State-space analysis of the dynamic characteristics of a variable thrust liquid propellant rocket engine[END_REF], the system considered is the composition of the end linear system with a simple LQR feedback controller,

A K = A c (x t,Np+1 , u t,Nu )+B c (x t,Np+1 , u t,Nu )K, κ ∈ R + (satisfying κ < -λ max (A K )) and Q K and R K are positive definite sym- metric matrices Q K ∈ R n×n , R K ∈ R m×m .
As explained in the aforementioned references, the first use of the computed P ∈ R n×n is to add a terminal-region penalty term in the MPC cost (24). It is used in constraints too, which are explained later. Besides, an integral action is added for increasing tracking precision on x z (in a similar manner to [START_REF] Santos | On-line implementation of nonlinear MPC: an experimental case study[END_REF]). Integrator decision variables are denoted by z (contained in Z) and are also penalised in the cost through a weight matrix S ∈ R nz×nz , whose diagonal is [1, 0.1, 0.1, 0.1, 0.1]. Thus, the MPC cost J to minimise consists in:

J(X, U, Z) =   Np k=1 (x k -x t,k ) T Q(x k -x t,k ) + Nu k=1 (u k -u t,k ) T R(u k -u t,k ) + Np k=1 z T k Sz k   ∆t + (x Np+1 -x t,Np+1 ) T P (x Np+1 -x t,Np+1 ). ( 24
)
Basically, the integral and terminal costs are added to the traditional quadratic cost on states and controls, with a prediction horizon N p = 10 steps (0.1s) and a control horizon N u = 5. Implicitly, the last control u Nu is used for k ≥ N u . Since these horizons are shorter than those used in trajectory generation (Section 3.1), the last step does not necessarily correspond to the end reference point. Longer horizons (more costly) did not enhance tracking or constraints satisfaction. Q and R are positive-definite symmetric weighting matrices Q ∈ R n×n , R ∈ R m×m , whose diagonals have been tuned off-line via Kriging-based black-box optimisation as in [START_REF] Marzat | Automatic tuning via Kriging-based optimization of methods for fault detection and isolation[END_REF]. The criterion for that weight definition involves the minimisation of static error and overshoot in simulations. Moreover, some robust considerations have been implemented. The minimisation of the previous J under constraints is not intrinsically robust. In fact, robustness to parameters and initial conditions variations and to modelling error is of interest in this application. In robust MPC approaches, the minimax optimisation, which minimises the worst-case scenario, is usually employed. A generic expression of this problem with uncertain dynamics is the following [START_REF] Loefberg | Minimax Approaches to Robust Model Predictive Control[END_REF], in which ∆ represents uncertainty in dynamic matrices

A d (•, •, ∆), B d (•, •, ∆): min U max ∆ J(X, U) s.t. X ∈ X ∀∆ ∈ ∆ c U ∈ U ∀∆ ∈ ∆ c . (25)
Nevertheless, solving (25) for large compact uncertainty sets ∆ c is too computationally costly for this application. Thus, only a finite set of uncertain scenarios are considered (inspired from [START_REF] Calafiore | Robust Model Predictive Control via Scenario Optimization[END_REF]); and an equivalent formulation to the minimax in [START_REF] Loefberg | Minimax Approaches to Robust Model Predictive Control[END_REF] is solved. Precisely, a scalar γ ∈ R + is minimised via an epigraph formulation. That γ constrains the previous cost J defined in (24) evaluated at some perturbed states propagations X i :

X i = [x i,1 , ..., x i,k , ..., x i,Np+1 ] T , i ∈ I ∆x i,k+1 = A d (x t,k , u t,k , ∆ i,k )∆x i,k + B d (x t,k , u t,k , ∆ i,k )∆u k , k ∈ [0, N p + 1], (26)
where ∆ i,k are certain selected internal parameter variations belonging to ∆ = {∆ i,k , i ∈ I, k ∈ [0, N p + 1]}. I is a finite set which indexes the perturbed cases. The epigraph formulation shifts robustness considerations into the constraints list, thereby not requiring maximisation. Only a smooth convex nonlinear programme (NLP) is solved. That minimisation problem, in which decision variables are extended to account for all X i , is stated below:

min Xi,U,Zi,γ γ (27)
s.t. J(X i , U, Z i ) ≤ γ ∀i ∈ I X i ∈ X, U ∈ U ∀i ∈ I A ineq [X i U] T ≤ b ineq ∀i ∈ I A i,eq [X i U] T = b i,eq ∀i ∈ I (x i,Np+1 -x t,Np+1 ) T P i (x i,Np+1 -x t,Np+1 ) ≤ α P ∀i ∈ I z i,k+1 = z i,k + ∆tK I (x z,i,k -x z,t,k ) ∀i ∈ I, k ∈ [0, N p ].
As in Section 3.1, the set U for the first control u M P C ≡ u 1 is specially constrained as in [START_REF] Singh | Trajectory generation for a UAV in urban terrain, using nonlinear MPC[END_REF], but in this closed-loop case u 0 is the previous-step control (warm start is performed). In (27) there are also equality constraints (defined by A i,eq and b i,eq ) for the different linear dynamics along the whole horizon (26) and the same linear inequality constraints as in [START_REF] Saudemont | Study of a robust control law based on H infinity for the Vulcain rocket engine -Etude d'une commande robuste a base de commande H infinity pour le moteur Vulcain[END_REF] applied to all X i . Concerning the terminal region, a constant α P refers to the neighbourhood in which the Lyapunov term of J is constrained in a nonlinear way (further details on the QIH method in [START_REF] Chen | A Quasi-Infinite Horizon Nonlinear Model Predictive Control Scheme with Guaranteed Stability[END_REF]). A specific P i is used in each perturbed case since matrices are different. The series of ∆ i , where I = {1, 2, 3}, are taken constant for all k and represent the following parameter-variation cases:

-i = 1: nominal parameters.

-i = 2: the main varying parameters are increased by 10%: interdependent tanks pressures and temperatures. -i = 3: the main varying parameters are decreased by 10%: interdependent tanks pressures and temperatures.

These cases consist in the worst ones which the controller must face in terms of tanks parameters variations, corresponding to the bounds of ∆. It can be checked in simulations that these would be the extreme dynamics-perturbing cases related to the compact uncertainty set ∆ c , for which a costly minimax would be necessary. Since the resulting U computed in ( 27) is forced to comply with all these uncertainty scenarios and all propagated perturbed states must verify all constraints, the robustness of the controller is increased. The last set of constraints in (27) consists in the integrator dynamics [START_REF] Santos | On-line implementation of nonlinear MPC: an experimental case study[END_REF], where K I is a gain matrix computed off-line in the same manner as Q and R.

Analysis of results

The interior-point optimisation software IPOPT [START_REF] Waechter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] has been employed due to the smoothness and convexity of the NLP (27). Simulations of the aforementioned control loop within the MATLAB environment are executed from 1.5s until 3s after the engine-start command, which is the time window in which continuous control is applied in the start-up transient. Mixture-ratios initial values commence very far from the allowable area, owing to the low initial mass flows that alter the definition of quotients. In fact, chambers are still not physically ignited during the first instants (even if igniters have been activated). Hence, M R are not representative quantities at the very beginning of the start-up. Fig. 5 comprises the results of p CC tracking at three non-dimensional pressure levels: p CC,r = 1 (nominal), p CC,r = 0.7 (minimum for this engine) and p CC,r = 1.2 (maximum). At all three points, the final reference mixture ratios remain equal M R CC,r = 6, M R GG,r = 1 and M R P I,r = 5.25, as usual during start transients. M R tracking in the nominal case is shown in Fig. 6 (with end-reference curves for enhancing visualisation). Tracking is accomplished with acceptable accuracy in p CC for all cases (under 0.12% in nominal, under 0.9% in off-nominal) and in M R (under 0.2% in nominal, under 3% in off-nominal). Simultaneously, constraints are respected up from the time when it is considered feasible and acceptable to respect them in practice (1.9s). The presence of a slight overshoot (2.6% in nominal) before the final convergence is generated by the nonlinear influence of the GG-starter input mass flow, which is not accounted for in the linearised model. It is more elevated in the minimum case since the relative impact of the starter is greater. Planned trajectories from Fig. 3 do not present overshoot in order to counteract this aspect of real engine behaviour. Fig. 7 shows the computed control as well as the planned trajectories. Some differences are present, especially during the first 0.5s in α V GC , due to the more realistic initial mixture ratios in the simulator, which lie outside the bounds considered in trajectory generation. Even though a nonlinear model is used in the planning of those first instants, modelling errors are still relevant due to simplifications [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF][START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF].

The controller achieves that tracking performance even if random initial conditions (after the sequential transient phase) are set, while constraints-satisfaction times vary by some hundredths of seconds. Computational times in MATLAB are generally ten times longer than real time, not ruling out a future real-machine implementation.

Comparison with open loop and alternative MPC controller

Table 1 contains the comparison in terms of some performance indicators between this trajectory-tracking MPC (T.MPC) proposal, the end-state tracking (E.MPC, with different robustness considerations) presented in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF], and OL simulations. The nominal OL corresponds to engine's original input command, which is precisely tuned for the standard case, as traditionally done in flight-ready engines. The minimum and maximum OL commands have been calculated via the preprocessor explained in Section 3.1, maintaining the same control as in the reference u r along the transient. It seems clear that the strategy in this paper results in a higher level of performance than OL and that it reduces overshoot with respect to the alternative approach in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF]. Indeed, the real benefit of this CL MPC control is more noticeable at off-nominal points, where it is more problematic to accomplish multivariable tracking at high performance while meeting constraints during the transient. Furthermore, even if the nominal pressure and M R are targeted, the precomputed u r would not always drive the system in OL to the same state x r . As explained in the next subsection, robustness is also an advantage.

Robustness analysis

In Section 3.2, the different perturbed cases have been listed. The system has to be robust to those scenarios in which some internal parameters can vary in real operation. These scenarios have been simulated by intentionally altering those parameters in the simulator. The performance indicators obtained with the different approaches for the nominal thrust are shown in Table 2. Results point to the moderately greater robustness to parameters variations of the approach presented in this paper in comparison to [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF]; and considerably greater with respect to OL. Comparing both CL approaches, the enhancements of T.MPC in terms of static errors and settling time are not elevated. The noticeable reduction of overshoot is the major advantage.

Conclusion

The control of the transient phases of liquid-propellant rocket engines is normally carried out in open loop due to its highly nonlinear and hybrid behaviour. In this paper, the control of the fully continuous part of the start-up of a gasgenerator-cycle LPRE has been treated. Continuous valves are the considered actuators for regulating pressure in the main chamber and mass-flow mixture ratios. The tracking of these quantities, states in the model used, is accomplished via an MPC controller, which at the same time verifies a set of hard operational constraints. Concretely, what the controller tracks are pre-generated trajectories of states and control. These are computed off-line in a preprocessor which takes into account a full-state end reference constructed according to launcher commands, via a nonlinear state-space model of the engine. The linearised MPC controller with integral action achieves trajectory tracking with sufficient accuracy and constraints are satisfied during the transient. Robustness considerations are also included in the controller so as to face a set of internal-parameter variation scenarios, plausible in real engine operation. The computationally costly nested minimax optimisation, commonly used in robust MPC approaches, has been rewritten as the minimisation of a scalar cost. The controller demonstrates robustness against the predefined worst-case scenarios, since performance is not degraded. In future work, other ways of posing this robustness consideration globally, for instance via the estimation of perturbations, will be investigated. The use of discrete actuators in control may also be a relevant factor to consider. A more extensive validation study with respect to parametric variations will be carried out.

Fig. 1 .

 1 Fig. 1. Vulcain 1 simplified flow plan with engine acronyms

= a 10 x 1 2 -b 10 x 8 2 -c 10 x 10 x 8 -d 10 x 10 2 -e 10 x 4 + f 10 u 3 - g 10 x 10 2 u 3 ,a 11 x 2 2 -b 11 x 9 x 2 -c 11 x 9 2 -d 11 x 11 2 -e 11 x 4 + f 11 u 4 - g 11 x 11 2 u 4 ,

 2282432322224424

Fig. 3 .

 3 Fig. 3. Trajectory planning results in pCC for pCC,r = 1 (nominal), pCC,r = 0.7 (minimum) and pCC,r = 1.2 (maximum)

Fig. 4 .

 4 Fig. 4. Trajectory planning results in u for pCC,r = 1 (nominal). Maximum opening angles are 90➦.

Fig. 5 .

 5 Fig. 5. Tracking results in pCC for pCC,r = 1 (nominal), pCC,r = 0.7 (minimum) and pCC,r = 1.2 (maximum)

Fig. 6 .

 6 Fig. 6. Tracking results in M R for pCC,r = 1 (nominal)

Fig. 7 .

 7 Fig. 7. Control results for pCC,r = 1 (nominal)

Table 1 .

 1 Performance-indicators comparison between this MPC proposal, end-state MPC[START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF] and OL at the three selected operating points

	Operating		Nominal		Minimum	Maximum
	point
	Indicator	OL E.MPC T.MPC OL E.MPC T.MPC OL E.MPC T.MPC
	Settling time	2.8 2.51	2.54 2.67 2.55	2.57	2.69 2.53	2.51
	(pCC,r ± 1%) [s]
	Overshoot (%	6.31 5.04	2.6	15.1 11.46	11.29 3.34 4.04	2.57
	in pCC )
	Constraints	1.81	1.8	1.82 1.83 1.76	1.78 1.77 1.81	1.84
	verification [s]

Table 2 .

 2 Performance-indicators comparison between this MPC proposal, end-state MPC[START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF] and OL at the perturbed scenarios

	Perturbed		Nominal	Increased by 10% Decreased by 10%
	cases
	Indicator	OL E.MPC T.MPC OL E.MPC T.MPC OL E.MPC T.MPC
	Settling time	2.8 2.51	2.54	2.7	2.6	2.55 2.66 2.52	2.53
	(pCC,r ± 1%) [s]
	Overshoot (%	6.31 5.04	2.6	5.72 4.79	2.67 5.3	5.45	2.6
	in pCC )
	Constraints	1.81	1.8	1.82 1.79 1.79	1.82	1.8	1.8	1.82
	verification [s]
	pCC static er-	0.25 0.26	0.12 0.35 0.36	0.22 0.64 0.17	0.04
	ror (%)
	M RCC static	0.17 0.01	0.13	3.07 0.19	0.3	3.48 0.17	0.05
	error (%)
	M RGG static	1.39 0.05	0.12	0.6 0.0008	0.09 0.83	0.1	0.16
	error (%)
	M RP I static	1.43	0.3	0.19 3.24 0.12	0.01 3.65 0.48	0.37
	error (%)

ACD2019, 028, v5 (final): 'Trajectory planning and tracking via MPC for transient control . . .

ACD2019, 028, v5 (final): 'Trajectory planning and tracking via MPC for transient control . . .

Transient control of liquid-propellant rocket engines

Sergio Pérez-Roca et al.

ACD2019, 028, v5 (final): 'Trajectory planning and tracking via for transient control . . .

ACD2019, 028, v5 'Trajectory planning and tracking via MPC for transient control . . .


text:
Sergio Pérez-Roca
email: sergio.perez_roca@onera.fr

Julien Marzat
email: julien.marzat@onera.fr

Hélène Piet-Lahanier
email: helene.piet-lahanier@onera.fr

Nicolas Langlois
email: nicolas.langlois@esigelec.fr

Marco Galeotta
email: marco.galeotta@cnes.fr

François Farago
email: francois.farago@cnes.fr

Serge Le Gonidec
email: serge.le-gonidec@ariane.group

Robust Transient Control of Reusable Liquid-Propellant Rocket Engines

Keywords: Liquid-propellant rocket engines, Model Predictive Control, transient phases, robustness, reusability

The current trend towards a more affordable access to space is generally materialising in reusable launchers and engines. From the control perspective, these reusable liquid-propellant rocket engines (LPRE) imply more demanding robustness requirements than expendable ones, mainly because of their multi-restart and thrustmodulation capabilities. Classically, the control system handles LPRE operation at a finite set of predefined points. That approach reduces their throttability domain to a restricted interval in which they are designed to be safe in nominal conditions. Moreover, the operation of their transient phases, which have a great impact on the duration of engine life, is not robust to the possible engine evolution. Hence, the goal of this work is to develop a control loop which is adapted to the whole set of operating phases, transient and steady-state, and which is robust to internal parametric variations. Several blocks have been assembled to constitute the control loop: engine simulation, reference generation and several controllers. First, simulators representative of the gas-generator-cycle (GG) Vulcain 1 and PROMETHEUS engines were built. The purely thermodynamic modelling of the cycle was subsequently adapted to the control framework, obtaining a nonlinear state-space model. The available actuators are continuously controllable valves, binary igniters and binary starters. These actuators are related to discrete events in transient phases. Regarding the start-up operation, the igniter, starter and valves are activated during the first seconds. Up from the end of those activations, the whole system behaves in a fully continuous way. Hence, a different control strategy is proposed for each sub-phase. For the first and discrete sub-phase, a discrete optimisation of events timing is proposed, in which the time differences between events are adapted according to operation criteria and constraints. This trajectory planning, still under implementation, is to be performed off-line. The subsequent continuous sub-phase is feedback controlled to track pre-computed reference trajectories. Apart from the start-up, throttling scenarios also present a dedicated end-state-tracking algorithm. A model-based control method, Model Predictive Control, has been applied in a linearised manner with robustness guarantees to all these scenarios, in which a set of hard state and control constraints must be respected. Tracking of pressure (thrust) and mixture-ratio operating points within the design envelope is achieved in simulation along the continuous sub-phase while respecting constraints. Robustness to variations of the parameters, which are checked to be predominant according to analyses, is also demonstrated.

Introduction

The current context of launcher vehicles design is strongly linked to their reusability. From the automatic control perspective, this potential need for reusable liquid-propellant rocket engines (LPRE) implies stricter robustness requirements, related to their multi-restart and thrust-modulation capabilities. These demanding requirements stem from the possible internal perturbations caused by components faults or evolving parameters and from exogenous perturbations linked to the more complex mission profiles conceived for new launchers.

Multivariable-control developments of main-stage LPRE have attained a short throttling envelope (70%-120%) in test benches [START_REF] Gonidec | Automatic & Control applications in the European space propulsion domain. From need expression to preparation for an uncertain future[END_REF]. In real engines, the control system is generally designed to achieve the nominal operating point. One of the new features of the future European PROMETHEUS engine is the ability to throttle down to 30% of thrust [START_REF] Baiocco | Technology demonstration for reusable launchers[END_REF]. Consequently, an enlarged controlled operating domain has to be feasible. Tracking and robustness must be kept at those low throttle levels, where physical phenomena are harder to anticipate. The main control goal in these multivariable systems consists in tracking reference points in combustion-chamber pressure and mixture ratio, which come from a high-level order according to launcher needs. Controlvalves opening angles are changed so as to adjust engine's operating point while verifying a set of constraints. The majority of control studies in the literature employs linearised models about operating points for synthesising steady-state controllers, mostly based on PID (proportional, integral, derivative) techniques (such as [START_REF] Nemeth | Reusable rocket engine intelligent control system framework design, phase 2[END_REF]). In most papers, initial MIMO (Multi Input Multi Output) systems are decoupled into dominant SISO (Single Input Single Output) subsystems. Off-line optimisation studies have also been carried out [START_REF] Dai | Damage-Mitigating Control of a Reusable Rocket Engine: Part II-Formulation of an Optimal Policy[END_REF]. Other works in the literature, incorporating more advanced techniques such as nonlinear [START_REF] Lorenzo | Nonlinear control of a reusable rocket engine for life extension[END_REF], hybrid [START_REF] Musgrave | Real-time accommodation of actuator faults on a reusable rocket engine[END_REF] or robust [START_REF] Saudemont | Study of a robust control law based on H infinity for the Vulcain rocket engine -Etude d'une commande robuste a base de commande H infinity pour le moteur Vulcain[END_REF] ones, reinforce certain indicators of performance and robustness. However, there is an absence of publications which concern not only the steady state but also the demanding transient phases at the same level of performance and robustness, as reviewed in [START_REF] Pérez-Roca | A survey of automatic control methods for liquid-propellant rocket engines[END_REF]. There is also a lack of method comparisons on a common benchmark, even simulated. Besides, only narrow throttling domains are feasible. Sequential transient phases of engine operation (start-up and shutdown), are generally executed in open loop with narrow correction margins. The first sub-phase of the start-up transient is determined by a series of discrete activations. After all these commands, the second sub-phase of the transient, which is fully continuous, takes place until the steady state is reached. Open-loop (OL) control is normally applied due to the controllability and observability issues at very low mass flows [START_REF] Nemeth | Reusable rocket engine intelligent control system framework design, phase 2[END_REF]. In this paper, different control approaches are proposed for each sub-phase. Throttling transient scenarios are also treated differently. The main objective of this work is to control the transient phases of pump-fed LPRE. Concretely, tracking of combustion pressure (linked to thrust) and mixture ratio all along the transients is targeted. Simultaneously, a set of hard operational constraints has to be respected, related to mixture ratios, turbopumps rotational speeds and valves actuators angular velocities. The LPRE cycle studied is the gas-generatorcycle engine, but the method here is conceived to be applicable to other cycles. Control approaches in this paper are based on Model Predictive Control (MPC) techniques. MPC has been selected as the most adequate for this type of complex systems with hard constraints. Indeed, it is gaining popularity at academic and industrial levels and can incorporate robustness [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF] or hybrid aspects, which are important for future work on this topic. This paper is organised as follows. In Section 2, the modelling approach for simulation and control is described. In Section 3, the derived models are analysed, especially from the point of view of sensitivity to parameters. Section 4 describes the different control strategies developed, mainly based on MPC techniques. The main results, concerning the selected scenarios, are presented and discussed in Section 5. Finally, Section 6 serves as a conclusion.

Modelling approach

The modelling approach considered in these studies was first described in [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF] and revisited in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF][START_REF] Pérez-Roca | Trajectory planning and tracking via MPC for transient control of liquid-propellant rocket engines[END_REF]. Here it is again revisited and updated for considering different engine case-studies. Firstly it is relevant to clarify that several model structures are employed in the different blocks of the control loop in this paper. Concerning the plant to which the control is applied, a simulator of the real plant was constructed in the first place, instead of testing on a real engine or on a hardware-in-the-loop test bench. Indeed, a simple, dynamic and efficient way of modelling generic LPRE was sought, instead of using more accurate programmes or computations. An easy integration into Simulink was also preferred so as to easily test different control methods. Along these lines, a new Simulink library of rocket-engine components has been developed to build a simulator of the gasgenerator-cycle engine. It has been named T-RETM, Toolbox for Rocket-Engine Transient Modelling. This simulator, whose structure is built componentwise, contains the basic thermo-fluid-dynamics and mechanics of LPRE elements: mass, energy and momentum conservation equations. Simplified 0-D models of resistive components (valves, pipes, turbopumps) and capacitive ones (cavities, combustion chambers) have been developed. Subsequently, all components are joined in Simulink according to the selected engine. The engines considered in this paper, representative of Vulcain 1 and PROMETHEUS, present a GG (gas generator) cycle. The former consists in a LOX/LH 2 (liquid oxygen as oxidiser, liquid hydrogen as fuel) engine, while the latter presents the LOX/LCH 4 mixture (liquid oxygen, liquid methane). In Fig. 1 the Vulcain 1 cycle and main components are depicted and the main acronyms are summarised, while Fig. 2 is devoted to PROMETHEUS. The hot-gas flow necessary to drive turbines comes from a GG, a small combustion chamber that receives a small portion of the main propellant flow. The actuators considered in this paper are four continuously controllable valves (VCF, VCO, VGF, VGO) and a fifth one in the case of Vulcain 1 (VGC). Apart from those, there are three discrete actuators: two binary chamber igniters (𝑖 𝐶𝐶 , 𝑖 𝐺𝐺 ) and one binary starter ( 𝑖 𝑠𝑡𝑎 ). In Vulcain 1 𝑖 𝐺𝐺 ≡ 𝑖 𝑠𝑡𝑎 can be assumed according to the sequence. That GG starter injects hot gas into the cavity during less than 1.5s so as to start driving turbines. Vulcain 1 requires two turbopump shafts while PROMETHEUS only needs one, accounting for the density differences between propellants. This consists in the main contribution to the start-up transient phases of LPRE, because once turbines start rotating, pumps can provide more flow to chambers, which increases combustion pressures and temperatures. These increases also lead to greater shaft speeds until a steady-state is achieved, at around three to four seconds after start. Indeed, the activations of the aforementioned actuators, both the initial opening of valves and the ignitions, constitute the so-called start-up sequence. That consists in the first and discrete sub-phase of the start-up transient. Once all actuators are active, a continuous sub-phase takes place, where valves are continuously adjusted so as to reach a desired steady state, which is still far when the discrete sub-phase ends. Valves angles (α), which present a nonlinear but direct relation to sections (A), control the flows to the main combustion chamber (VCF and VCO), to the GG (VGF, VGO), and to the oxidiser turbine (VGC, in Vulcain 1). The latter is the main contributor to adjusting mixture ratio (MR), defined as the quotient between oxidiser (O) and fuel (F) mass flow rates 𝑀𝑅 = 𝑚̇𝑂 𝑚̇𝐹 . This ratio, a major behaviour indicator in LPRE, is established at three levels: at an engine's global level (𝑀𝑅 𝑃𝐼 ), taking pumped propellants into account; in the combustion chamber (𝑀𝑅 𝐶𝐶 ) and in the GG (𝑀𝑅 𝐺𝐺 ). Since the beginning it was aimed at capturing in models the transient behaviour of the engine with the available engine parameters. No identification techniques from real data or precise multi-physical simulation platforms were available to the project. Thus, the developed simulator was then translated into a nonlinear statespace model by joining components equations symbolically. That is to say, a global expression of system's dynamic behaviour as a function of its state variables, internal parameters and control inputs was targeted.

This process was performed via the symbolic mathematical environment Maple, taking into account the internal definition of each LPRE component and of the global system (interconnection of components, as shown in Fig. 1 and2). At this stage, having already performed some simplifications with respect to the initial simulator, the model is referred to as complex NLSS (nonlinear statespace) or 𝑓 𝑐 (𝒙, 𝒖). However, this model presented an excessive complexity for control design, since many ODEs (ordinary differential equation) extended over several pages. Thus, it was further reduced until attaining the here-called simplified NLSS such that 𝒙̇= 𝑓 𝑠 (𝒙, 𝒖), more tractable for the derivation of control laws (more details in [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF]).

Fig. 2. PROMETHEUS flow plan considered in model

The number of states is n (12 in Vulcain 1 (V.) and 9 in PROMETHEUS (P.)) and m is the number of control inputs (5 in V. and 4 in P.). Here, the state vector 𝒙, of both NLSS, comprises turbopumps speeds, the several pressures in the system (combustion chamber, GG, and pre-turbine cavities), and mass flows streaming through valves and pre-turbine pipes. Thus, in V., the state vector is defined as: 𝒖 𝒅 = [𝑖 𝐶𝐶 𝑖 𝐺𝐺 𝑖 𝑠𝑡𝑎 ] 𝑇 . Besides, all equations, states and control are rendered non-dimensional with respect to the nominal equilibrium values. An overview of the Vulcain 1 statespace system can be found in [START_REF] Pérez-Roca | Trajectory planning and tracking via MPC for transient control of liquid-propellant rocket engines[END_REF]. Here, the form of the PROMETHEUS dynamic system 𝒙̇= 𝑓 𝑠 (𝒙, 𝒖) is presented, where 𝑎 𝑖 , 𝑏 𝑖 , . . . , 𝑘 𝑖 ∈ ℝ are internalparameters conglomerates and 𝑤 𝑡 is an exogenous input corresponding to starter mass flow: .

𝒙 𝑽 = [𝜔
𝑥̇1 = (𝑎 1 𝑥 1 2 + 𝑏 1 𝑥 1 𝑥 4 + 𝑐 1 𝑥 4 + 𝑥 1 (𝑑 1 𝑥 6 + 𝑒 1 𝑥 7 + 𝑓 1 𝑥 8 + 𝑔 1 𝑥 9 ) + ℎ 1 𝑥 6 2 +
It is clear that this system representation, even if it has been simplified, presents numerous nonlinearities and is non-affine with respect to control.

Models analysis

The analysis of the previous models is necessary in order to select and build the most appropriate control strategy. Several aspects can be analysed, ranging from the effect of simplifications, stability, measurability and observability, controllability, sensitivity to parameter variations, etc. The different simplifications carried out obviously increase the modelling error slightly. It is specially present in mass flows, which can present a mismatch of 10 to 25% at each step of simplification (simulator, 𝑓 𝑐 , 𝑓 𝑠 , and linearised models). Errors in the rest of states remain below 10% at each step.

Stability is a general property of these GG LPRE models, which describe a system that compensates itself when a perturbation appears in the cycle due to its high coupling [START_REF] Pérez-Roca | A survey of automatic control methods for liquid-propellant rocket engines[END_REF]. The possible destabilising phenomena in LPRE, such as combustion instabilities or frequency couplings, are not modelled in these state-space models for the sake of simplicity. Hence, when analysing linearised models about transient trajectories, all eigenvalues present negative real parts. The state is assumed to be completely measurable in the real engine. This is a realistic assumption for 𝜔 and 𝑝. However, measuring some mass flows would be problematic in terms of engine design. Mass flows are normally not measured in LPRE, but estimated through pressure, temperature and volumetric flow measurements. This estimation process is assumed perfect in this paper, but consists in a future block to add to the loop. In fact, during the discrete part of the transient, some observability issues can appear due to the low volumetric flows characteristic of that phase [START_REF] Nemeth | Reusable rocket engine intelligent control system framework design, phase 2[END_REF]. This, apart from other issues highlighted in Section 4, hinders the realisation of closed-loop (CL) control during the discrete sub-phase of transients. Controllability can be verified after linearising 𝑓 𝑠 about the whole continuous start-up trajectory. In other words, all states can be controlled via valves during the continuous sub-phase. However, as commented for observability, controllability is not ensured during the discrete sub-phase, where there are valves that have to remain closed in order to execute a safe sequence. As a consequence, not all states can be controlled. Indeed, the first valves to open are CC-injection ones, which are considerably less influential in the control of the system than GG ones, which open at the end of the sequence, after ignitions and starter activation. The previous models establish a complex but deterministic relation between control inputs and system states. No stochastic considerations are explicitly included, which could arise from the fact that some parameters or dynamics are not perfectly known a priori. This is an issue when modelling this kind of complex thermodynamic systems, since some of the parameters used in the differential equations are estimated from test data or tuned in simulations. In order to analyse the effect of parameters on the behaviour of the engine, a series of sensitivity analyses have been performed. In the following, only Vulcain 1 will be used as case study due to the availability of a more consolidated set of internal parameters in contrast to PROMETHEUS, still under development. The list of considered engine parameters that might vary during operation or that are estimated for modelling (and hence not 100% certain) is summarised in Table 1. In order to study the influence of each individual parameter on the transient behaviour of the engine, the Sobol sensitivity analysis [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] has been applied. This method is appropriate for large, highly-coupled nonlinear systems such as this one. In addition, it is a global sensitivity-analysis method, which allows to identify the relative contributions of each parameter independently as well as the combined variations to the variance in the output. Here, it is of interest to study the variations in terms of a global performance criterion J Sobol , defined as the weighted addition of several simulated indicators: The effects of variations in the rest of parameters can be neglected according to the sensitivity analysis.

J Sobol =

Control approaches

The control goals mainly concern reaching a desired end-state or following a predefined trajectory while complying with a set of hard constraints on 𝒙 and 𝒖. This second goal is very important during transient phases, in order to avoid excessive mixture ratios (and hence temperatures), pressures or rotational speeds. Reusability requirements also point to this goal. The proposed structure of the whole control loop is depicted in Fig. 3.

Fig. 3. Complete control loop, with areas of future work

To the right there is the plant block, where the LPRE simulator is placed. It simulates the behaviour of the engine via integration of the system of ODE as a function of its state, parameters and control inputs. The provided control inputs come from the controller, which computes the orders sent to valves (actuators block). These actuators are modelled as a separate entity via simplified second-order models, requiring an angular input instead of sections, which are considered in the controller. In fact, there is a nonlinear monotone relation between angles and sections, and hence it is more convenient to avoid that additional equation in the statespace models used in the controller. This controller receives the full state measurement from the simulator, since estimation is still considered perfect (future work).

It also receives a reference to follow, which consists in pre-generated trajectories. Indeed, that reference is generated in a preprocessor connected to the controller.

Preprocessor

The role of the preprocessor is to generate in the first place a target reference (𝒙 𝒓, 𝒖 𝒓 ) in terms of the full state and control vectors. In fact, the orders provided by the launcher to the engine control system only concern the final desired 𝑝 𝐶𝐶,𝑟 , 𝑀𝑅 𝐶𝐶,𝑟 , 𝑀𝑅 𝐺𝐺,𝑟 and 𝑀𝑅 𝑃𝐼,𝑟 . In order to obtain the whole 𝒙 used in the state-space models in this work, a calculation has to be performed. This is explained in [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF]. Basically, system's ODE are equated to zero in order to obtain an equilibrium point that coincides with the provided input elements. It is paramount to obtain this full-state and control equilibrium point because the control approach in this study is fully model-based, as explained in the next subsections.

Once these (𝒙 𝒓, 𝒖 𝒓 ) are computed, if the controller has to deal with the continuous start-up transient, trajectories (𝑿 𝒕, 𝑼 𝒕 ) for states and control are also pregenerated off-line according to the initial conditions and to the desired 𝒙 𝒓 . 𝑿 and 𝑼 are defined in general as the series of 𝒙 and 𝒖 at each time step k along a horizon N (valid throughout the whole paper): 𝑿 = [𝒙 1 , … , 𝒙 𝑘 , … , 𝒙 𝑁 ] 𝑇 𝑼 = [𝒖 1 , … , 𝒖 𝑘 , … , 𝒖 𝑁 ] 𝑇 This trajectory planning is explained in [START_REF] Pérez-Roca | Trajectory planning and tracking via MPC for transient control of liquid-propellant rocket engines[END_REF], where an optimisation-and model-predictive-based method is proposed. This computation method can be regarded as an OL finite-horizon MPC scheme, typically used in trajectory planning [START_REF] Betts | Survey of Numerical Methods for Trajectory Optimization[END_REF]. A simple cost based on the current quadratic error with respect to (𝒙 𝒓, 𝒖 𝒓 ) is considered. The NLSS 𝑓 𝑠 is used for defining the nonlinear-dynamics constraints, more relevant during the initial build-up phases. Close to the steady state, linearised models are used, since modelling error is reduced. In addition, all important operating constraints can be imposed in the planning, which in the end provides a reference control trajectory that is safer than the constant OL openings. Thus, the preprocessor serves to define the references to be tracked by the controller, either an end state or a trajectory.

MPC Controllers

Controllers in this study are based on MPC. This family of techniques is an optimisation-based control approach, where the control law is recomputed at each time step 𝑘 with the goal of minimising a cost function. Basically, control inputs are computed according to the minimisation of the cost while respecting all dynamic and operating constraints along a discretised time horizon, in which the behaviour of the system is predicted with models. Only the first computed control of the horizon, corresponding to the first time step, is factually applied to the plant, since at the next step the whole process is repeated. In this paper, different controllers are proposed according to the scenario that the engine is facing. As explained in previous sections, one of the goals of this work is to control start-up transients in CL. These transients contain a discrete and a continuous phase, which are treated differently here. The latter was considered in a first place so as to simplify the discrete elements of models, which are already active in that phase. Apart from the start-up, operating-point changes or throttling scenarios (fully continuous) are also targeted.

Continuous MPC controllers

Regarding the controllers for fully continuous scenarios, the tracking goal is different whether throttling or startup is to be controlled. Throttling: an operating-point change is performed here by tracking a new desired equilibrium (𝒙 𝒓, 𝒖 𝒓 ), without specifying a reference trajectory. The algorithm for this kind of scenarios can be taken from the previous paper [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF], where it was originally presented as an also valid option for controlling the continuous start-up. The linear MPC algorithm considered there presented linear dynamic constraints in the form: Δ𝒙 𝑘+1 = 𝐴 𝑑 (𝒙 𝑟 , 𝒖 𝑟 )Δ𝒙 𝑘 + 𝐵 𝑑 (𝒙 𝑟 , 𝒖 𝑟 )Δ𝒖 𝑘 ,

(2) where 𝐴 𝑑 and 𝐵 𝑑 are the zero-hold-discretised and linearised state-space matrices. In (2) they are always evaluated at the end point to track. The goal is to find the control Δ𝒖 = 𝒖 -𝒖 𝑟 = 0 that drives the system to Δ𝒙 = 𝒙 -𝒙 𝑟 = 0. Hence, the cost to be minimised is then:

𝐽 𝑒𝑛𝑑 (𝑿, 𝑼, 𝒁) = (∑ Δ𝒙 𝑘 𝑇 𝑄Δ𝒙 𝑘 𝑁 𝑝 𝑘=1 + ∑ Δ𝒖 𝑘 𝑇 𝑅Δ𝒖 𝑘 𝑁 𝑢 𝑘=1 + ∑ 𝒛 𝑘 𝑇 𝑆𝒛 𝑘 𝑁 𝑝 𝑘=1 )) Δ𝑡 + Δ𝒙 𝑁 𝑝 +1 𝑇 𝑃Δ𝒙 𝑁 𝑝 +1 , (3)
where Q and R are weight matrices for state and control respectively. 𝑁 𝑝 and 𝑁 𝑢 are the states and control prediction horizons. S is related to the addition of integral error states, represented by 𝒛, and the last term is related to an end-state penalty term, as in [START_REF] Chen | A Quasi-Infinite Horizon Nonlinear Model Predictive Control Scheme with Guaranteed Stability[END_REF]. The MPC optimisation programme to solve at each time step is given by: min 𝑿 𝑖 ,𝑼,𝒁 𝑖 ,𝛾 𝛾 (10) The problem solved here is a minimisation of a scalar 𝛾 ∈ ℝ + which constrains the cost evaluated at a set of perturbed states 𝑿 𝑖 . This consists in robust approach to MPC, based on [START_REF] Calafiore | Robust Model Predictive Control via Scenario Optimization[END_REF][START_REF] Loefberg | Minimax Approaches to Robust Model Predictive Control[END_REF], where a control 𝒖 is computed so as to face different perturbed dynamic propagations at the same time. Perturbations consist in internal parameter variations 𝛥 𝑖,𝑘 inducing different 𝐴 𝑑 and 𝐵 𝑑 . 𝐼 is the set of perturbed scenarios cases, in this case an amount of three plus the nominal scenario. In this sense, (2) is reformulated as: Δ𝒙 𝑖,𝑘+1 = 𝐴 𝑑 (𝒙 𝑟 , 𝒖 𝑟 , 𝛥 𝑖,𝑘 )Δ𝒙 𝑖,𝑘 + 𝐵 𝑑 (𝒙 𝑟 , 𝒖 𝑟 , 𝛥 𝑖,𝑘 )Δ𝒖 𝑘 . In [START_REF] Musgrave | Real-time accommodation of actuator faults on a reusable rocket engine[END_REF] the bounds of states 𝑋 and control 𝑈 are imposed, concerning the maximum thermo-fluid-dynamic values and valve opening limits. Inequality [START_REF] Saudemont | Study of a robust control law based on H infinity for the Vulcain rocket engine -Etude d'une commande robuste a base de commande H infinity pour le moteur Vulcain[END_REF] serves to establish the minimum and maximum mixture-ratio constraints, vital to avoid high temperatures or extinction. In addition, actuators speed constraints are also included there. In [START_REF] Pérez-Roca | A survey of automatic control methods for liquid-propellant rocket engines[END_REF], linear dynamic constraints between states and control are required. Inequality (9) ensures end-state reachability [START_REF] Chen | A Quasi-Infinite Horizon Nonlinear Model Predictive Control Scheme with Guaranteed Stability[END_REF] and (10) describes integrator's dynamics [START_REF] Santos | On-line implementation of nonlinear MPC: an experimental case study[END_REF], where a gain 𝐾 𝐼 is defined so as to reduce error in 𝒙 𝑧 . And all the constraints (5)- [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF] have to be verified for all 𝑖. Start-up control: if the start-up transient is to be executed, another approach is proposed, in which a set of planned trajectories (𝑿 𝒕, 𝑼 𝒕 ) are tracked, and not only a final point. Trajectory tracking is more convenient in this case where the system evolves in a highly nonlinear way and where trajectories can be pre-computed without tight computational-time limits. In a quick throttling request, there might not be enough time to generate adequate trajectories (around one minute in MATLAB), but only an end reference (about 3s). Since this is also a continuous MPC controller, only the continuous phase of the start-up, which starts after the opening of the GG-injection valves, is controlled. In this case, the NLSS 𝑓 𝑠 model is linearised about trajectories at each 𝑘, which reduces the prediction error: Δ𝒙 𝑖,𝑘+1 = 𝐴 𝑑 (𝒙 𝑡,𝑘 , 𝒖 𝑡,𝑘 , 𝛥 𝑖,𝑘 )Δ𝒙 𝑖,𝑘 + 𝐵 𝑑 (𝒙 𝑡,𝑘 , 𝒖 𝑡,𝑘 , 𝛥 𝑖,𝑘 )Δ𝒖 𝑘 . The cost function ( 3) is modified so as to penalise deviations from trajectories: )) Δ𝑡 + (𝒙 𝑁 𝑝 +1 -𝒙 𝑡,𝑁 𝑝 +1 ) 𝑇 𝑃(𝒙 𝑁 𝑝 +1 -𝒙 𝑡,𝑁 𝑝 +1 ). The minimisation problem in this case presents the same robust structure as in ( 4)- [START_REF] Pérez-Roca | Derivation and Analysis of a State-Space Model for Transient Control of Liquid-Propellant Rocket Engines[END_REF], but now it is related to trajectory tracking. This implies the cost function is substituted by 𝐽 𝑡𝑟𝑎𝑗 and that deviations from trajectories are considered instead of final references. This method is presented in detail in [START_REF] Pérez-Roca | Trajectory planning and tracking via MPC for transient control of liquid-propellant rocket engines[END_REF].

Hybrid MPC controllers

Discrete sub-phase of start-up: concerning the discrete sub-phase of the start-up transient, a computation logic is proposed here. The problem to solve in this hybrid scenario is more complex than in the purely continuous ones. Indeed, the discrete events, which consist in valve openings, ignitions and starter activation, alter the dynamic behaviour of the engine. During the sequence, each valve (CC and GG injection ones) is forced to remain closed until a certain instant, when it can start influencing the system via 𝒖. Ignition and starter activations are modelled here as discrete inputs which activate a set of terms in differential equations. All these events are executed at some precise instants which are precisely tuned for nominal conditions. In this paper a strategy for optimising the time differences between events is proposed. Furthermore, the sections of the valves which are already open are also computed. The goal is to robustify the sequence a priori and to gain controllability of valves, which otherwise are simply opened to pre-defined degrees. The nominal order of events is maintained, since it consists in a safe succession of actuations according to fuel and oxidiser properties. An optimisation problem based on MPC principles is again envisaged. An OL finite-horizon scheme, precomputed off-line as for trajectory generation, seems more appropriate than an on-line CL one. In on-line conditions, little observability and controllability are present during this phase, as explained in Section 3. In addition, computational times tend to be long due to the necessary inclusion of nonlinear dynamic constraints. Indeed, in this frame, the effects of events can be expressed via constraints. But these constraints need to be nonlinear because they must include a dependency on the additional decision variables 𝛕, which incorporate the optimal time differences between events. In the case of Vulcain 1, it is: 𝝉 = [𝜏 𝑉𝐶𝐹 𝜏 𝑉𝐶𝑂 𝜏 𝑖 𝐶𝐶 𝜏 𝑖 𝐺𝐺 𝜏 𝑉𝐺𝐹 ]. The definition of these intervals is graphically shown in Fig. 4. The function discrete_tree provides the value of discrete inputs 𝑖 according to the current time intervals and instant. If the current instant 𝑡 𝑘 surpasses the sum of the respective intervals, the corresponding event is activated. The vector 𝑼 𝑜𝑓𝑓 reflects the additional constraints which have to be verified on control if valves are forced to be closed according to the timing. The system of equations 𝑓 𝑠 adapts according to the discrete inputs, and an implicit backward Euler scheme is proposed. In this manner, the changing dynamics are expressed via these implicit nonlinear constraints. With these constraints, the following minimisation problem is to be solved, where increments are with respect to (𝒙 𝑟,𝑑 , 𝒖

Bounds are still considered in [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]. In [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF], in contrast to the previous [START_REF] Saudemont | Study of a robust control law based on H infinity for the Vulcain rocket engine -Etude d'une commande robuste a base de commande H infinity pour le moteur Vulcain[END_REF], only actuators limits are constrained, since mixture ratios are not meaningful if some valves are closed. Equation ( 15) fixes the initial conditions of the problem, but not dynamics, which is constrained in [START_REF] Calafiore | Robust Model Predictive Control via Scenario Optimization[END_REF] by running [START_REF] Pérez-Roca | An MPC Approach to Transient Control of Liquid-Propellant Rocket Engines[END_REF]. The sum of all time differences must not exceed the beginning of the continuous phase 𝑡 𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑜𝑢𝑠 [START_REF] Chen | A Quasi-Infinite Horizon Nonlinear Model Predictive Control Scheme with Guaranteed Stability[END_REF]. In [START_REF] Loefberg | Minimax Approaches to Robust Model Predictive Control[END_REF] the corresponding valves are forced to be closed, depending on 𝝉. The obtained 𝑿 𝒅 , 𝑼 𝒅 and 𝝉 would serve as reference trajectories for a robustified on-line tracking of the discrete sub-phase.

The controller performing that tracking would have to take into account the same model-structure switches determined by the planned timing 𝝉. That implies a varying number of available control inputs. Seeing the controllability and observability scarcities during this phase (mentioned in Section 3), only the tracking of 𝑼 𝒅 trajectories would be conceivable. That is to say, only valve sections would be susceptible of being modified according to that reference, since not all the states in the system would be controllable nor observable. This approach remains to be numerically evaluated with extensive simulations.

Results and discussion

In order to solve the defined optimisation problems, the interior-point optimisation software IPOPT [START_REF] Waechter | On the implementation of an interior-point filter linesearch algorithm for large-scale nonlinear programming[END_REF] has been employed within MATLAB. A general time step of Δt = 10ms has been used in all cases, due to engine computer constraints. Computational times in MATLAB are about ten times longer than real time, not ruling out a future real-machine implementation.

Continuous control results

Regarding the results of continuous-control scenarios, different cases have been successfully simulated on the Vulcain 1 engine. Firstly, continuous start-up control (trajectory tracking) results for different levels of desired end CC pressure (nominal 100%, 70% and 120%) are depicted in Fig. 5. The adjustable level of pressure or thrust at the end of the start-up is a required capability of new engines. Mixture ratios are kept constant to MR CC,r = 6, MR GG,r = 1, MR PI,r = 5.25.

Tracking is achieved with acceptable accuracy in p CC for all cases (under 0.07% in nominal, under 0.95% in off-nominal) and in MR (under 0.32% in nominal, under 3.2% in off-nominal). When testing a throttling scenario, the resolution of ( 4)-( 10) yields the results shown in Fig. 6. A throttle-down operation from 100% to 70% of thrust has been selected as a representative case. It is noticeable that the operation shift is safely performed (no constraints violated) within 0.34s. As explained in Section 4, robustness considerations are included in the algorithms. Since some internal parameters can vary during operation, and others present uncertainty, the variations of the predominant parameters in Section 3 are also taken into account. Some robustness tests are presented in [START_REF] Pérez-Roca | Trajectory planning and tracking via MPC for transient control of liquid-propellant rocket engines[END_REF], in which variations of tanks pressure are well mitigated. The considered perturbed cases in that paper, are also considered here (+10% and -10% in p tanks ), since they are the most probable. The inclusion of the selected list of additional influential parameters in a third perturbed scenario also mitigates their variations in the simulator, as shown in the overall performance indicators summary in Table 2. That third scenario is referred to as the worst case, since it consists in the worst combination of those alterations. Results point to a considerable performance improvement in all indicators with respect to the original OL, especially in perturbed scenarios. Robustness to parametric variations is noticeable in those cases.

Conclusions

The evolving design requirements of reusable launchers and their associated liquid-propellant rocket engines impose more demanding robustness requirements on control systems. Transient phases have traditionally been carried out in OL, which does not ensure robustness to possible parameter variations in the engine. These variations are more probable in reusable engines owing to their longer lives and executions. In this work, approaches to control the transient phases of a GG-cycle LPRE have been proposed. The main goal is to track combustion-chamber pressure and mixture ratios while respecting engine constraints. By making use of state-space models of these engines, strategies based on Model Predictive Control (MPC) have been developed for the different transient scenarios that the engine faces. An end-state-tracking algorithm for throttling operations, as well as a trajectory-tracking one for the start-up transient, have been synthesised with successful results. Robustness to the set of more influential parameters in the model is considered in the controller and demonstrated.

Perspectives

Several improvement areas can be mentioned. Full-state measurements are assumed perfect, which in reality has to be solved via mass-flow estimators. Discrete inputs to the system, involved in the start-up sequence, are also to be controlled in order to robustify that sequence. The aforementioned algorithm (11)-( 18) is a proposal which mimics in a relatively precise way the physics of the complex problem which is controlling the discrete subphase of the LPRE start-up. Due to the numerous nonlinearities, coupling and implicitness of constraints, the optimisation problem becomes highly non-convex. The software IPOPT is not well suited for these computations. Thus, this approach remains to be numerically evaluated with extensive simulations. General non-convex optimisation solvers are under study, as well as other ways of expressing the problem. Furthermore, more realistic simulations, such as hardware-in-the-loop ones, are also to be carried out.

Fig. 1 .

 1 Fig. 1. Vulcain 1 flow plan considered in model

  𝐽 𝑡𝑟𝑎𝑗 (𝑿, 𝑼, 𝒁) = (∑ (𝒙 𝑘 -𝒙 𝑡,𝑘 ) 𝑇 𝑄(𝒙 𝑘 -𝒙 𝑡,𝑘 ) 𝑁 𝑝 𝑘=1 + ∑ (𝒖 𝑘 -𝒖 𝑡,𝑘 ) 𝑇 𝑅(𝒖 𝑘 -𝒖 𝑡,𝑘 )

Fig. 4 .

 4 Fig. 4. Representation of considered time intervals during discrete sequenceThe opening of the last valve (VGO) is considered at a fixed time, when the continuous phase starts. The goal is to attain a reference state 𝒙 𝑟,𝑑 at that instant. Hence, a simple cost is used, only penalising the difference between the final step and that reference. The implicit dependencies on 𝝉, expressed as nonlinear constraints, are built in the following way: Algorithm 1. Definition of hybrid nonlinear dynamic constraints 𝐼𝑛𝑖𝑡𝑖𝑎𝑙𝑖𝑠𝑒 𝑿 𝑁𝐿 𝑎𝑛𝑑 𝑼 𝑜𝑓𝑓 𝑎𝑠 𝑣𝑜𝑖𝑑 𝑚𝑎𝑡𝑟𝑖𝑐𝑒𝑠 𝒇𝒐𝒓 𝑘 = 1, 2, … , 𝑁 𝑝 (11) 𝑡 𝑘 = 𝑡 0 + 𝑘Δ𝑡; [𝒖 𝑘,𝑜𝑓𝑓 , 𝑖 𝐶𝐶 , 𝑖 𝐺𝐺 , 𝑖 𝑠𝑡𝑎 ] = 𝑑𝑖𝑠𝑐𝑟𝑒𝑡𝑒_𝑡𝑟𝑒𝑒(𝝉, 𝑡 𝑘 ); 𝒙 𝑘+1 = 𝑓 𝑠 (𝒙 𝑘+1 , 𝒖 𝑘+1 , 𝑖 𝐶𝐶 , 𝑖 𝐺𝐺 , 𝑖 𝑠𝑡𝑎 )Δ𝑡 + 𝒙 𝑘 ; 𝐴𝑝𝑝𝑒𝑛𝑑 𝒙 𝑘+1 𝑡𝑜 𝑿 𝑁𝐿 𝐴𝑝𝑝𝑒𝑛𝑑 𝒖 𝑘,𝑜𝑓𝑓 𝑡𝑜 𝑼 𝑜𝑓𝑓 𝒆𝒏𝒅

Fig. 5 .Fig. 6 .

 56 Fig. 5. Start-up control: p CC tracking At the same time, constraints are respected once mixture ratios become meaningful. Further performance indicators can be consulted in Table2.

  𝐻 𝜔 𝑂 𝑝 𝐶𝐶 𝑝 𝐺𝐺 𝑝 𝐿𝑇𝐻 𝑝 𝑉𝐺𝐶 𝑚̇𝐿 𝑇𝐻 𝑚̇𝑉 𝐶𝐹 𝑚̇𝑉 𝐶𝑂 𝑚̇𝑉 𝐺𝐹 𝑚̇𝑉 𝐺𝑂 𝑚 𝑉𝐺𝐶 ] 𝑇 , and for P. it is: 𝒙 𝑷 = [𝜔 𝑝 𝐶𝐶 𝑝 𝐺𝐺 𝑝 𝐿𝑇 𝑚̇𝐿 𝑇 𝑚̇𝑉 𝐶𝐹 𝑚̇𝑉 𝐶𝑂 𝑚̇𝑉 𝐺𝐹 𝑚̇𝑉 𝐺𝑂 ] 𝑇 , The states with higher tracking importance are incorporated into a reduced state vector 𝒙 𝒛 : 𝒙 𝒛 = [ 𝑝 𝐶𝐶 𝑚̇𝑉 𝐶𝐹 𝑚̇𝑉 𝐶𝑂 𝑚̇𝑉 𝐺𝐹 𝑚̇𝑉 𝐺𝑂 ] 𝑇 . The control input 𝒖 = [𝒖 𝒄 𝒖 𝒅 ] 𝑻 contains control inputs of continuous and discrete nature, which renders the model hybrid from the control perspective. The sections of control valves are the continuous ones: 𝒖 𝒄 = [𝐴 𝑉𝐶𝐹 𝐴 𝑉𝐶𝑂 𝐴 𝑉𝐺𝐹 𝐴 𝑉𝐺𝑂 (𝐴 𝑉𝐺𝐶 𝑖𝑛 𝑉. )] 𝑇 , while igniters and starter activations consist in the discrete ones:

  𝑖 1 𝑥 7 2 + 𝑗 1 𝑥 6 𝑥 8 + 𝑘 1 𝑥 7 𝑥 9 )𝑢 7 , 𝑥̇2 = (𝑎 2 𝑥 6 + 𝑏 2 𝑥 7 )𝑢 5 + 𝑐 2 𝑥 6 + 𝑑 2 𝑥 7 + 𝑒 2 𝑥 2 , 𝑥̇3 = (𝑎 3 𝑥 8 + 𝑏 3 𝑥 9 )𝑢 6 + 𝑐 3 𝑥 8 + 𝑑 3 𝑥 9 + 𝑒 3 𝑥 5 + 𝑓 3 𝑤 𝑡 , 𝑥̇4 = 𝑎 4 𝑥 4 + 𝑏 4 𝑥 5 , 𝑥̇5 = 𝑎 5 (𝑥 3 -𝑥 4 ) + 𝑏 5 𝑥 5 2 ,

	𝑥̇6 = 𝑥̇7 = (𝑎 7 𝑥 1 (𝑎 6 𝑥 1 2 +𝑏 6 𝑥 6 2 +𝑐 6 𝑥 6 𝑥 8 +𝑑 6 𝑥 2 +𝑒 6 )𝑢 1 2 +𝑓 6 𝑥 6 2 (𝑔 6 𝑢 1 +ℎ 6 )𝑢 1 2 + 𝑏 7 𝑥 7 2 + 𝑐 7 𝑥 1 𝑥 7 + 𝑑 7 𝑥 2 + 𝑒 7 )𝑢 2 + , 𝑥̇8 = (𝑎 8 𝑥 1 2 + 𝑏 8 𝑥 6 2 + 𝑐 8 𝑥 8 2 + 𝑑 8 𝑥 6 𝑥 8 + 𝑒 8 𝑥 3 + 𝑓 8 )𝑢 3 + 2 𝑔 8 𝑥 8 , 𝑢 3 𝑥̇9 = (𝑎 9 𝑥 1 2 + 𝑏 9 𝑥 7 2 + 𝑐 9 𝑥 9 2 + 𝑑 9 𝑥 1 𝑥 7 + 𝑒 9 𝑥 3 + 𝑓 9 )𝑢 4 + (1) 2 𝑓 7 𝑥 7 , 𝑢 2 𝑔 9 𝑥 9 2
	𝑢 4

Table 1 .

 1 |err pCC | + |err MRCC | + |err MRGG | + |err MRPI | + 0.001 • |A pCC | + 0.01 • |os pCC | ,where err are static steady-state errors in the main operating quantities (contained in 𝒙 𝒛 ), os is overshoot and A pCC is the surface between the perturbed p CC startup transient curve and the nominal one in open-loop (OL). The sampling of parameters variations according to Table1has been carried out via the Latin Hypercube sampling, recommended for calculating the Sobol global indices. In this fashion, a set of 1000 parametervariation combinations is provided to a Kriging-based Sobol algorithm based on[START_REF] Iooss | A review on global sensitivity analysis methods[END_REF]. This algorithm creates a Kriging model based on the provided input-output data, which avoids the execution of costly simulations. After running 1 million Monte Carlo evaluations of that model, the most influential parameters, related to the higher Sobol global indices, can be identified. These are AthCC, AthT and p tanks , whose indices present an order of magnitude of 0.1 (normalised with respect to one), while ResI, ResCR and η GG indices are about 0.01. AthCC and AthT are considered due to possible interengine discrepancy or varying effective hydraulic sections during operation. Oscillations in p tanks are also possible during operation, while ResI, ResCR and η GG are not certainly known, they are estimated for modelling. This shortlist of parameters, to which the model is more sensitive, is taken into account in the control algorithm (Section 4). List of Vulcain 1 engine parameters considered

	in sensitivity analysis
	Parameter Definition	Considered
			variation
			range [%]
	p tanks	Tanks pressures,	±10
		coupled with inlet
		densities and
		temperatures
	ResLC	LC lines fluidic	±10
		resistances
	ResI	Injectors fluidic	±10
		resistances
	ResCR	Cooling-circuit fluidic	±10
		resistance
	IneCR	Cooling-circuit fluidic	±10
		inertia
	η CC	CC combustion	±5
		efficiency
	η GG	GG combustion	±5
		efficiency
	ResLTH	LTH line fluidic	±10
		resistance
	IneLTH	LTH line fluidic inertia ±10
	Vcav	Pre-turbine cavities	±10
		volume
	AthCC	CC throat section	±5
	AthT	Turbine inlet sections	±5
	AthLE	Turbine outlet sections ±5

  𝒛 𝑖,𝑘+1 = 𝒛 𝑖,𝑘 + Δ𝑡 𝐾 𝐼 Δ𝒙 𝑧,𝑖,𝑘 ∀ 𝑖 ∈ 𝐼, 𝑘 ∈ [0, 𝑁 𝑝 ].

				4)
	𝑠. 𝑡. 𝐽 𝑒𝑛𝑑 (𝑿 𝑖 , 𝑼, 𝒁 𝑖 ) ≤ 𝛾	∀ 𝑖 ∈ 𝐼	(5)
	∀ 𝑖 ∈ 𝐼 𝐴 𝑖𝑛𝑒𝑞 [𝑿 𝑖 𝑼] 𝑇 ≤ 𝒃 𝑖𝑛𝑒𝑞 ∀ 𝑖 ∈ 𝐼 𝑿 𝑖 ∈ 𝑋, 𝑼 ∈ 𝑈	(6) (7)
	𝐴 𝑖,𝑒𝑞 [𝑿 𝑖 𝑼] 𝑇 = 𝒃 𝑖,𝑒𝑞	∀ 𝑖 ∈ 𝐼	(8)
	Δ𝒙 𝑖,𝑁 𝑝 +1 𝑇	𝑃 𝑖 Δ𝒙 𝑖,𝑁 𝑝 +1 ≤ 𝛼 𝑃 ∀ 𝑖 ∈ 𝐼	(9)

  𝑟,𝑑 ): min 𝑿 𝒅 ,𝑼 𝒅 ,𝝉 Δ𝒙 𝑁𝑝+1 𝑠. 𝑡. 𝑿 𝒅 ∈ 𝑋, 𝑼 𝒅 ∈ 𝑈 (13) 𝐴 𝑖𝑛𝑒𝑞,𝑑 [𝑼 𝒅 ] ≤ 𝒃 𝑖𝑛𝑒𝑞,𝑑 (14) 𝐴 𝑒𝑞,𝑑 [𝑿 𝒅 𝑼 𝒅 ] 𝑇 = 𝒃 𝑒𝑞,𝑑

	𝑇	𝑄Δ𝒙 𝑁 𝑝 +1	(12)
			(15)
	∑ 𝝉 ≤ 𝑡 𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑜𝑢𝑠	(16)
	|𝑿 𝒅 -𝑿 𝑁𝐿 | = 𝟎	(17)
	𝑼 𝑜𝑓𝑓 = 𝟎

Table 2 .

 2 Table 2 collects results for start-up and throttling control in OL, CL and in perturbed scenarios. Control performance indicators results

		Nominal	Worst case	Nominal	Worst case
	Scenario	100% Start-	100%	Throttling	Throttling
		up		Start-up	10070%	10070%
	Indicator OL	CL	OL	CL	OL	CL	OL	CL
	Settling
	time (𝑝 𝐶𝐶,𝑟 ±	2.76 2.54	-	2.51	-	0.34	-	-
	1%) [𝑠]
	Over-
	shoot (%	6.29 2
	in 𝑝 𝐶𝐶 )

.84 4.28 2.67 3.19 2.29 5.24 2.57

	𝑝 𝐶𝐶 static error (%)	0.21 0

.064 1.75 0.48 3.19 0.498 5.24 1.12

	𝑀𝑅 𝐶𝐶
	static	0.18 0.32 1.25 0.44 2.72 0.18 4.61 0.82
	error (%)
	𝑀𝑅 𝐺𝐺
	static	1.41 0.069 1.62 1.21 1.37 0.05 2.68 0.58
	error (%)
	𝑀𝑅 𝑃𝐼
	static	1.41 0.022 2.97 0.056 2.98 0.29	4.9 1.38
	error (%)