X H Pham

M Alamir
email: mazen.alamir@grenoble-inp.fr

F Bonne

P Bonnay

Revisiting a fixed-point hierarchical control design for cryogenic refrigerators under constraints, nonlinearities and real-time considerations

Keywords: Hierarchical MPC, Cryogenic refrigerator, Fixed-point iteration, Real-time, Modular Control Design, Constrained Control

This paper presents an extension of a recently proposed hierarchical control framework applied to linearization-based unconstrained control of a cryogenic refrigerator. The extension concerns the validation of the framework in the presence of both nonlinear models and constraints. It is also shown that real-time handling of these features requires a specific complexity reduction technique. This technique performs the distribution of the optimization process over cyclically changed decision variables aiming at limiting the number of iterations per updating period. Numerical simulations are proposed in order to show the impact of the parameters choices and to assess the real-time implementability of the proposed framework.

Introduction

Cryogenic refrigerators are necessary to cool down the superconducting devices in many applications including nuclear fusion reactors as well as particle accelerators [START_REF] Henry | Analysis of the ITER cryoplant operational modes[END_REF][START_REF] Claudet | Economics of large helium cryogenic systems: experience from recent projects at CERN[END_REF]. In fact, theses refrigerators are composed of many subsystems that span over large areas. For many reasons, it is common to make changes in these subsystems (PID re-tuning, valve changing) without any consistent assessment of how these changes might affect the overall system behavior. Several centralized architectures have been proposed and experimented on a experimental station at CEA 1 /IRIG 2 /DSBT 3 . However, these control methods appear to be non scalable and even inappropriate from an operational point of view. This is in particular true when some actuators need to be changed or some sub-systems need to be turned off/on. Hence, the need for scalable and modular design that smoothly accommodates for changes in local sub-systems level triggered interest in distributed/hierarchical design investigation leading to the framework proposed in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] which the present work extends.

In the last decades, non-centralized control architectures [START_REF] Scattolini | Architectures for distributed and hierarchical model predictive control -a review[END_REF][START_REF] Negenborn | On 35 approaches for distributed MPC made easy[END_REF][START_REF] Barcelli | Hierarchical multirate control design for constrained linear systems[END_REF][START_REF] Doan | A Hierarchical MPC Approach with Guaranteed Feasibility for Dynamically Coupled Linear Systems[END_REF][START_REF] Doan | A distributed optimizationbased approach for hierarchical MPC of large-scale systems with coupled dynamics and constraints[END_REF][START_REF] Stewart | Hierarchical cooperative distributed model predictive control[END_REF] attracted attention within the control community due to their ability to break the problem into smaller tractable ones inducing modularity, security, and easiness of implementation. These approaches can be split into two categories: Distributed and Hierarchical.

Distributed architectures involve subsystems exchanging information with adjacent neighbors. In contrast, in hierarchical frameworks, the subsystems communicate only with a coordinator that uses the information received from the subsystems to compute new evaluation of the coupling signals that are then sent back to the subsystems in order to enhance the overall quality of the system behavior. The survey book [START_REF] Negenborn | On 35 approaches for distributed MPC made easy[END_REF] listed more than 35 different approaches to such frameworks.

This paper is an extension of the framework proposed in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] where the reader can find the precise positioning concerning the landscape of available approaches as well as the first application to a cryogenic refrigerator. Nevertheless, the following items recall the main points that enable to position the framework studied in this paper and [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] in the landscape of related contributions:

Many studies have focused on the optimality and/or the fulfillment of constraints by assuming the closed-loop behavior to be stable [START_REF] Barcelli | Hierarchical multirate control design for constrained linear systems[END_REF][START_REF] Doan | A Hierarchical MPC Approach with Guaranteed Feasibility for Dynamically Coupled Linear Systems[END_REF][START_REF] Doan | A distributed optimizationbased approach for hierarchical MPC of large-scale systems with coupled dynamics and constraints[END_REF]. This strong assumption is not required for the fixed-point iterations proposed in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] and which is revisited in the present paper. Although this method is not proven to be stable when constraints and nonlinearities are present, it enforces at least local closed-loop stability. The examples provided in this paper partially assess its effectiveness even when nonlinearities and constraints are present.

In several studies, the same assumption is re-introduced by stipulating a weak coupling between subsystems. More precisely, it is assumed that the coupling signals do not lead to the destabilizing effect.

In many works, the structure of the coupling between the subsystems is constrained. For instance, in [START_REF] Stewart | Hierarchical cooperative distributed model predictive control[END_REF], it is assumed that there is only a coupling through control input actions. In [START_REF] Picasso | An MPC approach to the design of two-layer hierarchical control systems[END_REF], a hierarchical structure is proposed where the higher level is coupled with the control input provided by the lower level. Such structural assumptions on the coupling are not used here since they might highly compromise the scalability and the scope of application. Instead, in order to overcome the challenge of strong coupling in the scheme proposed in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF], the coordinator takes appropriate explicit actions in order to stabilize the iterations with local subsystems.

The present work extends the results of [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] in three directions:

This paper proposes an extended set of numerical experiments in which the feedback used in the subsystems level is a constrained and/or Nonlinear Model Predictive Control NMPC. More precisely, constrained NMPC is used for the Brayton cycle while constrained linear MPC is used for the Joule-Thomson cycle (see Fig. 1(a)).

Due to the non negligible computation time associated with constrained nonlinear framework which has to be performed several times in a single sampling period (due to the fixed-point iterations), it is shown that standard use of the state-of-the-art numerical software [START_REF] Andersson | Casadi: a software framework for nonlinear optimization and optimal control[END_REF] leads to a nonimplementable scheme. The paper proposes a complexity reduction technique based on the ideas proposed in [START_REF] Alamir | A framework for real-time implementation of low-dimensional parameterized NMPC[END_REF] in order to derive a real-time implementable version of the framework.

As far as the comparison with a decentralized framework is concerned, while in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] the local controllers (at the subsystems level) assume nominal values of the coupling signals, in this paper, local observers are implemented which highly improve the performance of the decentralized setting rendering the comparison fairer and more challenging for the hierarchical settings.

The paper is organized as follows: Section 2 describes the system and states the control problem. Section 3 recalls the hierarchical framework and its formulation as stated in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] in a more general nonlinear and constrained setting. Section 4 proposes a method to distribute the optimization over time towards achieving real-time implementability. The simulation results and analysis are given in section 5 while section 6 concludes the paper and gives hints for future investigations.

Problem description

Brief Description of the Cryogenic Refrigerator

This section briefly describes the cryogenic refrigerator that is studied in this paper 4 . This plant is used for the purpose of conducting experiments (testing cryogenic components, study of super-fluid helium, etc) aiming at improving the design and operation of cryogenic refrigerators. Fig. 1(a) shows the decomposition of the overall process into four interconnected subsystems. More precisely, the system can be viewed as the interconnection of four subsystems: Joule-Thomson (J-T) cycle (S 1 ), Brayton cycle (S 2 ), LN2 precooling (S 3 ), and warm compression station (WCS) (S 4 ).

A cryogenic refrigerator provides cooling power by using a cryogenic turbine (the one following the valve CV 156 in S 2 of Fig. 1(a)) to extract thermal energy from the fluid and by exchanging heat power through a series of heat exchangers (represented by NEF x in Fig. 1(a)). The main objective is to reject the disturbance represented by the heat power induced by the operation of the experimental facility. In our experimental facility, this disturbance is represented by the heat source denoted by NCR 22 (see S 1 on Fig. 1(a)). A compressor is used to close the thermodynamic cycle in the so-called warm compression station (S 4 ). The role of this subsystem is to maintain the pressures (denoted by H p (high pressure) and B p (low pressure) in Fig. 1(a)) within some tightly imposed bounds which correspond to an optimal operation regime of the overall system.

Concerning the dynamic model of the cryogenic refrigerator, a library has been developed in [START_REF] Bonne | Simcryogenics: a library to simulate and optimize cryoplant and cryodistribution dynamics[END_REF] in which nonlinear models of all the cryogenic components necessary for system modeling (valves, heat exchangers, Helium bath, turbines, compressor) are implemented and connected before a linearization is performed around the operational point of the refrigerator. All the details concerning this library can be found in [START_REF] Bonne | Simcryogenics: a library to simulate and optimize cryoplant and cryodistribution dynamics[END_REF] and the references therein. This paper extends the work of [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] (by addressing the constrained nonlinear setting as well as the induced real-time implementation issues) while keeping the framework in which the warm area (subsystem S 3 and S 4 ) is assumed to be properly controlled. Therefore, this paper focuses on the control of the so-called cold zone, including subsystems S 1 and S 2 . The corresponding control problem is described in the next subsection. In order to clearly state the control problem, the inputs, outputs and constraints are first defined in the next section.

The Manipulated Inputs

There are three control inputs which are CV 155 , NCR 22 belonging to Joule-Thomson cycle (S 1 ) and ∆P 156 which is a part of the Brayton cycle (S 2 ). These actuators are defined below:

1. CV 155 ∈ [20%, 90%]: This valve is situated at the inlet of the helium bath. 2. NCR 22 : This actuator is located inside the helium bath (S 1 ). Note that this heating resistance is also used to simulate the heat pulses coming from the operation of physical experiments served by refrigerators. More precisely, the total power of the resistance is split into two components NCR 22 = NCR It should be noted that the valve CV 156 is used to control the pressure drop ∆P 156 between its inlet and outlet pressure. To do this, the local NMPC of S 2 computes and sends an appropriate value of the pressure drop ∆P 156 to the PID controller, which acts on the opening position of the valve CV 156 (Fig. 1(a)).

The Regulated Outputs

There are three regulated outputs (see Figure 1(a). for the notation) :

1. Ltb 131 : The helium liquid level (%) that must be controlled to ensure that some thermal charges deposited inside the phase separator (e.g. used to cool super-critical helium at liquid helium temperature to be ready for the final customer) are immersed with liquid helium. The set-point is chosen by the operator. In the usual operation, it is set at Ltb sp 131 = 60.5%. 2. Ttb 108 : The temperature at the inlet of the J-T valve must be tightly controlled in order to ensure the efficiency of the liquefaction of the helium.

3. Ttb 130 : Since the cryogenic turbine is a sensitive component, the temperature at its outlet must be regulated to avoid the risk of solid droplet forming at the outlet, potentially destructive for the turbine.

The Measured Outputs

The design of the local controllers needs the states of the subsystems to be reconstructed. This is done through dynamic observers. In order to synthesize the observers for these two subsystems, the measured outputs need to be specified. The observer of J-T cycle uses the measurement of Ltb 131 and Ttb 108 while the observer of the Brayton cycle uses the measurements of Ttb 130 .

The nominal operation point

The nominal operation point is the usual operating point of the refrigerator. This operation point is determined by the steady values of the inputs, outputs and heating loads described above, namely:

CV 0 155 , NCR (a,0) 22 , ∆P 0 156 , Ltb 0 131 , Ttb 0 130 , NCR (w,0) 22 (1)
The nominal values are shown in Table 1. These values are used to define the following deviation variables:

u 1 := CV 155 -CV 0 155 ; NCR (a) 22 -NCR (a,0) 22 , (2)
u 2 := ∆P 156 -∆P 0 156 , (3)
y 1 := Ltb 131 -Ltb 0 131 ; Ttb 108 -Ttb 0 108 , (4)
y 2 := Ttb 130 -Ttb 0 130 , (5)
w := NCR

In the sequel, the following condensed notation is used to refer to the nominal values of the controls and outputs to be regulated or constrained for the two subsystems involved in the study:

U 0 1 := CV 0 155 NCR (a,0) 22 , U 0 2 := ∆P 0 156 , (7)
Y 0 1 := Ltb 0 131 Ttb 0 108 , Y 0 2 = Ttb 0 130 (8)
Note that the operation point of the heating power is the sum of NCR

It is important to notice in order to keep the same operation point NCR which leads to the reduction of the amount of real heat disturbance power that can be rejected. This is the reason why NCR (a,0) 22 is chosen at the minimal value, which it can adopt in real-life configurations with the constraint (NCR In what follows, the bold-faced notation p denotes the profile of a vector variable p over a prediction horizon of length N . More precisely:

p = [p T (k), . . . , p T (k + N -1)] T ∈ R N •np
where n p is the dimension of p.

The dynamic model

Recall that all the details regarding the physical models are given in [START_REF] Bonne | Simcryogenics: a library to simulate and optimize cryoplant and cryodistribution dynamics[END_REF] and references therein. Only the formal and compact representation is given here in order to explain the hierarchical control framework. Using the sub-index s ∈ {1, 2} to refer to the subsystem S s , the nonlinear model governing the two subsystems defined, in terms of the deviation variable mentioned above, takes the following form:

x + s = f s (x s , u s , w s , v s ) y s = h s (x s , u s , v s ) (10)
v s ′ = g s ′ (x s , u s , v s )
with s ′ = s, s, s ′ ∈ {1, 2}. Note that the model of subsystem S 1 is fully linearized, namely f 1 (•), g 1 (•) and h 1 (•) are linear functions which is compatible with the linear nature of underlying phenomenon. On the contrary, a nonlinear model is necessary to describe the dynamics of S 2 . The notations x s , y s , v s and w s are defined as follows:

• x 1 ∈ R 20 and x 2 ∈ R 20 express the deviations from the steady state values of subsystem S 1 and S 2 respectively.

• y 1 ∈ R 2 and y 2 ∈ R express the deviations of the regulated outputs of each subsystem.

• v 1 ∈ R 3 and v 2 ∈ R 3 are coupling signals by which the two systems interact dynamically. Indeed as indicated in the last equation of ( 10), v 1 depends on (x 2 , u 2 , v 2 ) while v 2 depends on (x 1 , u 1 , v 1 ). More precisely, as illustrated in Fig. 1(a), the coupling variable v 1 is given by

v 1 := (P h , T h , P c ) T ∈ R 3 (11)
where P h and T h are the pressure and temperature (deviations) at the downstream inlet of the heat exchanger NEF 1 , whereas P c stands for the pressure at the upstream outlet of the same exchanger. Similarly, the coupling signal v 2 is given by:

v 2 := (M h , T c , M c ) T ∈ R 3 (12)
in which M h stands for the mass flow rate at the downstream inlet of the heat exchanger NEF 1 while T c and M c are the temperature and the mass flow rate at the upstream outlet of the same heat exchanger.

For any initial state x(k) and any control profile u defined over some prediction horizon of length N , the corresponding nominal [disturbance-free] trajectories, denoted by x(•, u, x(k)|v), which are obtained by applying the recursive formula [START_REF] Picasso | An MPC approach to the design of two-layer hierarchical control systems[END_REF], lead to the following so-called coherence constraints defined for all i ∈ {0, . . . , N -1}:

v 1 (k + i) = g 1 (x 2 (k + i, u 2 , x 2 (k)|v 2 ), u 2 (k + i), v 2 (k + i)), v 2 (k + i) = g 2 (x 1 (k + i, u 1 , x 1 (k)|v 1 ), u 1 (k + i), v 1 (k + i)),
These constraints simply express the fact that the signal profiles v := (v 1 , v 2 ) are compatible with the system's nominal coupled equations. This can be shortly written by introducing the following straightforward notation:

v 1 = g 1 (u 2 , x 2 (k), v 2 ) ( 13
)
v 2 = g 2 (u 1 , x 1 (k), v 1 ) (14)
Indeed, the r.h.s of ( 13) and ( 14) can be viewed as functions of x s (k), u s and v s (for s ∈ 1, 2), since the trajectories x s depend on x s (k), u s and v s (for s ∈ 1, 2).

Overview of the the hierarchical control architecture

In this section, the concept of coordinator that is invoked in the hierarchical design is introduced. The role of the coordinator is to control the system through two operation modes which are described hereafter:

1. In the first mode, the objective is to regulate the system around the nominal point x = 0 in spite of the unmeasured disturbances (disturbancerejection mode). This is the main objective of the cryogenic refrigerator.

2. In the second mode, the coordinator can temporarily drive the system to a different steady-state corresponding to a new set-point y = 0. For instance, the operator might decide to change the level Ltb 131 of liquid helium in the bath or the temperature Ttb 108 . This corresponds to a change in the corresponding set-points.

These two modes can be taken into account by using different set-points and different weighting matrices in the centralized cost function, namely:

J c (u, x(k), r d |v) = 2 s=1 J s (u s , x s (k), r d s , |v s ) (15)
where

J s (u s , x s (k), r d s , |v s ) := N i=1 r d s -h s (k + i) 2 Q (s) c + u s (k + i) + U 0 s 2 R (s) c
where h s (k + •) represents the output profile given by:

h s (σ) = h s (x s (σ, u s , x s (k)|v s ), u s (σ), v s (σ)) (16) Note that J s (u s , x s (k), r d s |v s ), s ∈ {1, 2} stand for the local costs that are com- puted for given v s . Q (1) c ∈ R 2×2 , Q (2) c ∈ R, R (1) c ∈ R (2×2) and R (2) c
∈ R + are positive semidefinite weighting matrices to penalize the distance of outputs to set-points, and to penalize the usage of the control action. Note that the output set-points r d 1 ∈ R 2 and r d 2 ∈ R represent the desired values of the deviation from the real set-points r d s + Y 0 s .

Remark 2.1. It is important to underline the difference between the central cost used by the coordinator and the stability oriented costs used by the local MPC controllers. Indeed, stability-oriented MPC formulations (local MPC) need to incorporate stabilizing terms that are not directly linked to the economic cost (central) and thus optimize a cost that is not exactly equal to the economic cost.

The coordinator role focuses on the economic cost using the stable results coming from the local subsystems.

At this stage, we have all we need to state the hierarchical control requirements:

Hierarchical control requirements.

A key feature in the proposed hierarchical framework is the concept of auxiliary set-points. Indeed, while r d stands for the desired set-points at the central level, there is a need to define the so-called auxiliary set-points. These are the set-points that are to be used in the local level so that the central cost function (15) is minimized. These auxiliary set-points have no reasons to be equal to the central (original) set-points as the expressions and the weighting used in the local and the central levels are different, as mentioned above.

A hierarchical control scheme should be defined where both local MPC and NMPC (designed for S 1 and S 2 respectively), receive the appropriate optimal auxiliary set-points r opt 1 and r opt 2 that minimize the central cost (15). Moreover, the following operational conditions have to be satisfied:

1. S 1 and S 2 communicate only with the coordinator.

2. The coordinator has no knowledge regarding the mathematical models of S 1 and S 2 nor the details of their local controllers' settings.

3. The modes mentioned before are handled by changing only the desired set-points r d and the weighting matrices Q c and R c . This should leave the local controller formulations and design unchanged.

Therefore, the coordinator computes the optimal auxiliary set-points r opt in order to minimize the central cost J c corresponding to the desired set-point r d . In order to do this, the hierarchical control algorithm illustrated in Fig. 2 is proposed. More precisely, the coordinator constructs a grid of auxiliary set-points G. Then, the fixed-point iteration described in section 3 is executed in order to evaluate the central cost J c (r) at all candidate auxiliary set-points r ∈ G.

Using these evaluations, a quadratic approximation Ĵc of the central cost J c is obtained.

Having this map at its disposal, the coordinator computes the candidate auxiliary set-points r opt c by solving a quadratic optimization problem inside a continuously updated trust-region R. A new round of set-point iterations is used to check whether this presumably optimal value does effectively improve the central cost. Depending on the answer, the trust-region size used to define the grid of auxiliary set-points in the next sampling instant is updated (size is increased in case of success and decreased otherwise). Depending on the case, the optimal value r opt (k) to be used is either the one just found, namely r opt c or the previous value r opt (k -1) adopted in the previous sampling period.

All these steps are successively explained in the remainder of the paper before a set of numerical investigations is proposed to illustrate the efficiency of the proposed framework.

Recalls on the fixed-point iteration framework

Let us first recall the fixed-point iteration based algorithm proposed in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] which is used to evaluate the central cost for any candidate auxiliary set-point r. Then the way the loop is closed by finding the best auxiliary set-point is explained.

Computing the central cost for a given auxiliary set-point: The Fixed-Point Iteration

In this section, a global set-point

r d = (r d 1 , r d 2 )
is given together with a value r of the auxiliary set-point at some sampling instant k. The task of the coordinator is to evaluate the value of the central cost J c (15) if the auxiliary set-point r = (r 1 , r 2 ) is adopted by the local controllers. This has to be done in spite of the absence of knowledge regarding the mathematical model, the current states or control design used at the local subsystems level.

Initialization

Note that the central optimization problem can be viewed as a problem in the extended decision variable (u, v) provided that one adds the coherence constraints ( 13)-( 14) on v, namely:

min u1,u2,v1,v2 2 s=1 J s (u s , x s (k), r d s |v s ) (17) subject to v 1 = g 1 (u 2 , x 2 (k), v 2 ) v 2 = g 2 (u 1 , x 1 (k), v 1 ) (18)
Note that the two individual costs are now conceptually decoupled for any given choice of coupling signal profiles v 1 and v 2 . The coupling now appears in the equality constraints (18).

In [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF], a method has been proposed to estimate the central cost defined by (17) using a fixed-point iteration in which, given a pair of set-points r = (r 1 , r 2 ), the coordinator starts with some initial guesses regarding the coupling signals:

v (σ) 1 , v (σ) 2
; σ = 0 (19)

These current guesses are sent to the subsystems S s , s ∈ {1, 2} so that each subsystem can compute the corresponding optimal control profile (should the coupling signal profile be correct), namely:

u opt s (r s , x s (k)|v (σ) s ) (20)
In order to compute u opt s , given the candidate auxiliary set-points r s and coupling profile v (σ) s , the local controllers solve in parallel and independently their local stability-oriented optimization problems, namely:

P s : min us N i=1 i N q x s (k + i) -x sp s (r s ) 2 Qs (21) + u s (k + i) -u sp s (r s ) 2

Rs

subject to :

x s (k + i) = x s (k + i, u s , x s (k)|v (σ) s ) u s ≤ u s ≤ u s
where (x sp s (r s ), u sp s (r s )) is the steady pair that is consistent with the candidate set-point r s and the given exogenous signal v (σ) .

x s (k + i) = x s (k + i, u s , x s (k)|v (σ)
s ) is the predicted state of the subsystem S s at instant k + i given the control profile u s , the initial state x s (k) and the presumed coupling profile v

(σ) s . Q 1 ∈ R 20×20 , Q 2 ∈ R 20×20 , R 1 ∈ R 2×2 and R 2 ∈ R
are the weighting matrices on the states and the control inputs of subsystems S 1 and S 2 , respectively.

Remark 3.1. As mentionned in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF], the time-dependent weighting term (i/N ) q for some q ∈ N allows to put high weight on the tail of the prediction horizon by choosing non-zero values q.

Note that the steady pair (x sp s (r s ), u sp s (r s )) corresponding to the candidate set-point r s , s ∈ {1, 2}, is computed by solving the optimization problem below:

P sp s : min x sp s ,u sp s y sp s -r s 2 Q sp s + u sp s R sp s ( 22
)
subject to :

y sp s = h s (x sp s , u sp s , v end s ) (23)
x sp s = f s (x sp s , u sp s , w s , v end s ) (24)
where

v end s := v (σ) s (k + N -1) is the last element of the exogenous profile v (σ) s , whereas Q sp 1 ∈ R 2×2 , Q sp 2 ∈ R, R sp 1 ∈ R 2×2 and R sp 2 ∈
R are respectively the weighting matrices on output and input. Note that (24) is the stationary condition associated to the set-point r s .

Once the control input profiles u opt s , s ∈ {1, 2}, are obtained, each subsystem computes the corresponding coupling profile and sends it to the coordinator :

S 1 sends to coordinator → v(σ+1) 2 := g 2 (u opt 1 , x 1 (k), v (σ)
1 ) S 2 sends to coordinator → v(σ+1)

1 := g 1 (u opt 2 , x 2 (k), v (σ) 2 )
In order to force the convergence of the iteration, a filtering dynamics has been proposed in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF] which is defined by:

v (σ+1) = F (v (σ) , v(σ+1) ) = (I -Π) • v (σ) + Π • v(σ+1) with v (σ) = v (σ) 1 v (σ) 2 (25)
The computation of the Π matrix and the resulting proof of convergence (in the linear unconstrained case) are described in [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF].

The filtered estimates v (σ+1) s

, are sent by the coordinator to the subsystems. Based on these values, the subsystems are able to compute the new optimal control profiles u opt s (r s , x s (k)|v (σ+1) s ) and corresponding profile of v(σ+2) , which is then sent back to the coordinator for a new iteration. It is important to note that it is only when the iteration of the fixed point converges to some v (∞) (r) that the coherence constraints ( 13)-( 14) become satisfied and the quadruplet:

u opt 1 (r 1 , x 1 (k)|v (∞) 1 (r)), u opt 2 (r 2 , x 2 (k)|v (∞) 2 (r)), v (∞) 1 (r), v (∞) 2 (r) (26)
can be viewed an admissible sub-optimal solution to the constrained central problem (17)-(18). Practically, the iterations are stopped if one of two termination criteria is reached:

σ ≥ σ max or ǫ := max(|v (σ+1) -v (σ) |) ≤ ǫ max (27)
Then, each subsystem computes its contribution to the central cost (using the central r d ), namely:

J s (u opt s (r s ), x s (k), r d s |v (∞) s (r))
and send it to the coordinator. After receiving these evaluations, the coordinator computes the central cost (15), which is now considered as a function of the auxiliary set-point r := (r 1 , r 2 ), r d and v (∞) (r):

J c (r|r d , v (∞) (r)) := 2 s=1 J s u opt s (r s ), x s (k), r d s |v (∞) s (r) (28)
Note that the estimation process of the central cost (left-hand side of (28)) by the coordinator, as described above, does not involve the knowledge of the state x s (k) nor the knowledge of the optimal control u opt s (r s ) but the estimation of the local costs and the coupling signals sent by the subsystems. Indeed, the estimation of local costs depends on the coupling signals profiles v (∞) (r) which encompasses all the information inaccessible to the coordinator and which are transmitted to the latter upon the convergence of the fixed point iterations. The previous discussion is summarized in Fig. 3 and Algorithm 1.

Fixed-point:

v (σ+1) = F (v (σ) , v(σ+1) ) v(σ+1) 2 , J 1 (u opt 1 (r 1 ), x 1 (k), r d 1 |v (σ) 1 ) v(σ+1) 1 , J 2 (u opt 2 (r 2 ), x 2 (k), r d 2 |v (σ) 2 ) r 1 , r d 1 r 2 , r d 2 v (σ) 1 v (σ) 2 J 1 (•, x 1 (k), r d 1 |v (∞) 1 ) J 2 (•, x 2 (k), r d 2 |v (∞) 2
)

J MP C 1 (•, x 1 (k), r 1 |v (σ) 1 ) J N MP C 2 (•, x 2 (k), r 2 |v (σ) 2 ) Compute u opt 1 (r 1 , x 1 (k)|v (σ) 1 ) Compute u opt 2 (r 2 , x 2 (k)|v (σ) 2 ) Compute J 1 (u opt 1 (r 1 ), x 1 (k), r d 1 |v (σ) 1 ) Compute J 2 (u opt 2 (r 2 ), x 2 (k), r d 2 |v (σ) 2 ) v(σ+1) 2 ← g 2 (u opt 1 (r 1 ), x 1 (k), v (σ) 1 ) v(σ+1) 1 ← g 1 (u opt 2 (r 2 ), x 2 (k), v (σ)
2 )

Figure 3: Schematic view of the fixed-point iteration at instant k [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF]. The states x 1 (k), x 2 (k) and the set-points rs and r d s are considered frozen. Note that subsystem S 2 uses NMPC as a local controller.

Optimizing the auxiliary set-points

In the previous section, it has been shown how the fixed-point iterations help the coordinator computing the central cost for a given auxiliary set-point r. Recall, however, that the role of the coordinator is to optimize the choice of the auxiliary set-points so that the central cost can be minimized. This section explains how the coordinator can use successive evaluations of the central cost Algorithm 1: Fixed-point-iteration-based-algorithm for a given setpoint vector r

Initialization: σ ← 0, v (0) s ← 0, ǫ ← ∞; while (σ ≤ σ max )
(σ+1) s -v (σ) s |) s=1..2 ; σ ← σ + 1 end for s ← 1, ..., 2 do
Subsystem s computes h s then J s , and sends J s to coordinator; end Coordinator computes J c (r) ← 2 s=1 J s (r s ) for different candidate auxiliary set-points to build a quadratic approximation of the central cost (as a function of the auxiliary set-points at the current sampling time k) in order to derive a candidate optimal auxiliary set-point r opt .

Approximating the central cost

Using the fixed-point iteration, the coordinator can compute for each auxiliary set-point r = (r T 1 , r T 2 ) T the corresponding value of the central cost:

J c (r|r d , v ∞ (r)) (29)
after convergence of the fixed-point iteration.

The central problem in the coordination layer can now be defined as follows:

r opt = argmin r J c (r |r d , v ∞ (r)) (30)
In order to solve (30), the central cost (15) will be approximated by a quadratic function, namely:

Ĵc = 1 2 r T Qr + f T r + c (31)
where Q ∈ R nr ×nr , f ∈ R nr and c ∈ R, with n r being the dimension of vector r.

These unknown parameters can be identified if the coordinator disposes of the values of the central cost at, at least, (n r + 1)(n r + 2)/2 different auxiliary set-points. The remainder of this section is devoted to explaining the way this is done by the coordinator. Note that this is a single possibility among many other possibilities of optimizing a black-box given function through different evaluations of its values at a set of possible points within its domain of definition. This is linked to the general task of derivative-free optimization.

In [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF], the grid of auxiliary set-points G(k) is constructed using a fixed grid around the original set-point r d . In the present version, a moving grid around the suboptimal solution found at the last instant. Moreover, the size ρ of the trust-region is modified, which will be described in section 3.2.2, depending on the relevance of quadratic approximation. At each sampling instant k, the grid G(k) of auxiliary set-points for the evaluation (at the sampling instant k) of the central cost is defined around the previous optimal value r opt (k -1) as follows:

G(k) := Pr r opt (k -1) + ∆(ρ(k -1)), R (32)
where

• R is an admissible set of r, namely:

R = {r | r min ≤ r ≤ r max } (33)
where r min , r max ∈ R nr are a priori defined bounds on possible values of the set-points.

• for a discrete subset D ⊂ R nr , the notation Pr (D, R) denotes the discrete set obtained by projecting all the elements of D on the hypercube R.

• ρ ∈ R + is positive real (size of the trust region where the quadratic approximation is presumably relevant).

• ∆(ρ) is a discrete set of displacements in R nr defined around 0 with distances that are proportional to ρ so that r opt (k -1) + ∆(ρ) represent the set of different auxiliary set-points around the previous optimal value to be visited and where the cost is to be evaluated). More precisely, the subset ∆(ρ) ⊂ R nr is defined by (m is supposed to be odd):

∆(ρ) := - (m -1) 2 ρ, . . . , -ρ, 0, ρ, . . . , (m -1) 2 ρ nr (34)
The trust region size ρ is updated at each instant k, which will be described later. Recall that the identifiability of the quadratic form coefficients is possible provided that m nr ≥ (n r + 1)(n r + 2)/2.

Based on the above definitions, the grid G(k) is constructed by using (32), the evaluation of the central cost J c (•) at every set-point r ∈ G(k) is performed by using the fixed-point methodology introduced in section 3. The values J c (r (j) ), j = 1, . . . , n ev ≤ m nr enable to compute the parameters of the quadratic form:

(Coordinator) min

Q,f,c nev j=1 J c (r (j) ) - 1 2 r (j) 2 Q + f T r (j) + c (35)
Once Q, f and c are available, a candidate optimal set-point r opt c (that minimizes the quadratic approximation) can be computed. Note, however, that since the central cost is not necessarily quadratic, this candidate optimal cost does not necessarily induce a decrease in the central cost. This can happen when the trust-region parameter ρ is too large for the quadratic approximation to be relevant. In such case, the size ρ should be reduced. This mechanism is discussed in the next section.

Trust region updating law of ρ

As mentioned previously, the parameter ρ defines the size of the neighborhood of the current desired set-point r d over which the better value is computed based on the current quadratic approximation of the cost function. On one hand, ρ must be sufficiently high to ensure a rapid decrease of the cost value. On the other hand, small values of ρ might be required in order for the quadratic approximation to be relevant. Hence, ρ should be updated accordingly: ρ is increased if the quadratic approximation induces a decrease of the cost function while ρ is decreased otherwise.

Concretely, the following quadratic problem is first solved to obtain the candidate value r c (k)

r opt c (k) = argmin r∈P(k) Ĵc (r) (36)
where P(k) is given by:

P(k) := Conv{Pr (r d (k) + ∆(ρ(k -1)), R)} (37)
Once the candidate r opt c (k) is obtained, the corresponding cost is computed by launching the algorithm 1 to obtain J c (r opt c ). The quadratic approximation is said relevant if it meets the condition below:

J c (r opt c ) < min{J c (r (j) ) | r (j) ∈ G(k)} (38)
Therefore, the trust-region size ρ is updated according to:

ρ(k) := β + • ρ(k -1) if (38) is satisfied β -• ρ(k -1) otherwise ( 39
)
where β + ≥ 1 and β -∈ (0, 1) denote respectively the expansion and the contraction factors. Finally, the updating law for r opt is given by:

r opt (k) := r opt c (k) if (38) is satisfied r opt (k -1) otherwise ( 40
)
where r opt (k -1) is the solution found at the previous instant k -1. The so adopted set-point r opt (k) is then sent to the subsystems with an end-of-iterations flag, which allows the subsystems to compute their corresponding control profiles. Finally, according to the MPC definition, the first action in each profile, namely:

u s (k) := [I n u s , O n u s , . . . , O n u s ]u opt s (r opt s (k), x s (k)|v (σmax) s ) ( 41
)
is applied by subsystem S s during the sampling period [k, k + 1].

Unfortunately, because of the presence of nonlinearity and constraints, the computation of the n ev necessary evaluations for central cost approximation might require a computation time that goes beyond the available sampling time T s .

The following section proposes a method to reduce the computation time with a rather little impact on the quality of the resulting closed-loop performance.

Distributing the optimization over time

Since using constrained nonlinear MPC induces a significant increase in the computation time, it might be impossible to compute a solution r opt (following the scheme of the previous section) for the next sampling period in the presence of limited computational resources. In order to overcome this potential issue, this section proposes a technique inspired by [START_REF] Alamir | A framework for real-time implementation of low-dimensional parameterized NMPC[END_REF] which is based on the idea of distributing the optimization over time. To facilitate the following explanations, the notation k and k + 1 are used to refer to instants kτ u and (k + 1)τ u with τ u being the control updating period, namely, the time during which the computation of a new optimal open-loop sequence is recomputed to implement the MPC feedback. Note that τ u is not necessarily equal to the sampling time T s . The process described in this section will be executed during the updating period [k, k + 1] as long as the computation time does not exceed τ u .

Recall that the approximation of the cost function J c (r) needs the evaluation of J c at n ev ≥ (n r + 1)(n r + 2)/2 values of the auxiliary set-points. By reducing the number of degrees of freedom (DOF) of vector r to be improved from n r to n z < n r , only (n z + 1)(n z + 2)/2 realizations would be needed, which accordingly leads to a decrease of the computation burden per updating period.

More precisely, a change in the decision variable is cyclically operated by defin-ing a reduced dimensional parameterization of r of the form:

r = M r + Dz (42)
where r ∈ R nr , M ∈ R nr ×nr and D ∈ R nr ×nz . Moreover, the transformation matrices M and D are changed in a cyclic way in order to explore all the degrees of freedom of r after a finite number of successive iterations. This is explained in a more detailed way in the remainder of this section.

At the beginning of each updating period k, the optimization problem to be solved is given by:

z ⋆ (k) = argmin z Ĵc (M (j k ) r ⋆ (k -1) + D (j k ) z) (43)
where the transformation defined by the matrices M (j k ) and D (j k ) is defined in order to assign some components of the vector r to be equal to the corresponding components of the previous solution r ⋆ (k -1) while leaving as degrees of freedom the n z remaining components that define the reduced dimensional decision variable z. Note that the definition of the transformation matrices depends on the updating instant k through the upper index j k which is a cyclic variable defined by:

j k = (j k-1 + 1) mod n r (44)
In the numerical investigation, the following two configurations are tested in order to illustrate the proposed methodology:

4.1. Configuration 1: n z = 1, n r = 3 M (0) =   0 0 0 0 1 0 0 0 1   D (0) =   1 0 0   (45)
M (1) =   1 0 0 0 0 0 0 0 1   D (1) =   0 1 0   (46)
M (2) =   1 0 0 0 1 0 0 0 0   D (2) =   0 0 1   (47) 4.2. Configuration 2: n z = 2, n r = 3 M (0) =   0 0 0 0 0 0 0 0 1   D (0) =   1 0 0 1 0 0   (48)
M (1) =   0 0 0 0 1 0 0 0 0   D (1) =   1 0 0 0 0 1   (49)
M (2) =   1 0 0 0 0 0 0 0 0   D (2) =   0 0 1 0 0 1   (50)
Note that the same methodology explained before regarding the definition of the grid of points is adopted with r and n r respectively replaced by z and n z . The only difference is that the number of degrees of freedom to be considered at the beginning of each updating period is reduced, and the significance of the degrees of freedom in terms of the components of r changes at each updating period.

When a sub-optimal solution z ⋆ (k) to ( 43) is obtained (after the allowed number of iterations), the corresponding candidate sub-optimal auxiliary set-point

r ⋆ c (k) is given by r ⋆ c (k) = M (j k ) r ⋆ (k -1) + D (j k ) z ⋆ (k) (51)
This candidate value is then used to update the size of the trust region in a similar way as explained above. The method can be simply sketched by Algorithm 2 for a given updating cycle involving n d iterations. More precisely, for-loop in Algorithm 2 allows to perform n d iterations within the updating period. Indeed, if the computation time does not exceed the updating period [k, k + 1], the whole process mentioned in this section can be repeated in order to improve the sub-optimal candidate auxiliary set-point r ⋆ c (k).

Simulation results

In this section, some numerical simulations are proposed in order to illustrate the different concepts and solutions introduced in the paper.

The simulation parameters

First of all, a cost index is necessary in order to evaluate and compare the performances in terms of closed-loop costs associated to the different framework settings. The commonly used closed-loop central cost will be adopted, namely:

J sim = 1 N sim 2 s=1 Nsim i=1 y(i) -y d Qc,s + u(i) Rc,s (52)
Algorithm 2: Pseudo code for the distributed-in-time optimization for l ← 1, ..., n d do Coordinator defines a grid of auxiliary set-points G(M (j k ) r ⋆ (k -1) + D (j k ) ∆(ρ(k -1)), R); Coordinator evaluates the cost function for each element r in the grid G(M (j k ) r ⋆ (k -1) + D (j k ) ∆(ρ(k -1)), R); Coordinator computes the quadratic approximation Ĵc (z) of J(z); Coordinator finds z ⋆ (k) by solving (43) Coordinator computes the candidate auxiliary set-point r ⋆ c (k) according to (51); Coordinator updates ρ and r ⋆ (k) using ( 39) and (40). end Coordinator sends r ⋆ (k) to the subsystems.

where N sim is the length of the simulation (in terms of sampling periods T s = 5s), Q c,s and R c,s are the penalty matrices used in the central cost (15). There are two sets of weighting matrices for the central cost (15), which will be used in the two distinct modes describes earlier (see the beginning of section 2.2), namely:

Mode 1: For disturbance rejecting scenario:

Q c,1 = 10 4 0 0 10 4 , R c,1 = 0 0 0 1 (53) Q c,2 = 10 6 , R c,2 = 0.1 (54)
Mode 2: For set-point tracking scenario:

Q c,1 = 10 6 0 0 0.1 , R c,1 = 0 0 0 0 (55) Q c,2 = 10 4 , R c,2 = 0 (56)
Where the second mode is dedicated to the tight regulation of the first output of the first subsystem, namely the liquid helium level in the bath (Ltb 131 ).

The penalty matrices for the local MPCs are fixed regardless of the scenario to the following values:

Q 1 = C T 1 • 10 0 0 1 • C 1 + 10 • I 20×20 R 1 = 1 0 0 1 Q 2 = C T 2 • 10 3 • C 2 + 10 • I 20×20 R 2 = 1
where C 1 is the matrix involved in the expression of the regulated output of the subsystem S 1 :

y 1 (k) = C 1 • x 1 (k) + D 1 • u 1 (k) + D v1 • v 1 (k) (57)
Whereas, the matrix C 2 is the partial derivative of h 2 (x 2 , u 2 , z 2 ) in [START_REF] Picasso | An MPC approach to the design of two-layer hierarchical control systems[END_REF] with respect to x 2 at an operating point (x op 2 , u op 2 , v op 2 ). Note that in each of the local MPC settings, a penalty on the whole state is used in order to enforce the stability of the local closed-loop.

The prediction horizon N is chosen to have a length of N T s (where N = 100 and T s = 5 seconds). This corresponds roughly to 8 minutes. This setting is currently used at CEA/IRIG/DSBT and is also the one that has been used in many previous studies involving MPC control design. As for the definition of the stage costs, the exponents q 1 = 1 and q 2 = 1 are used for the cost functions used in the local MPC design [see ( 21)].

In order to estimate the states to be used in the local MPC implementation, each subsystem uses an appropriate observer whose synthesis is beyond the scope of this paper. The observer not only estimates the states x s but also estimates the exogenous input v s (extended observer). For the decentralized settings (without coordinator), the exogenous signal v s (k) is supposed to remain constant over the prediction horizon. For the updating rules of the trust region size, the parameters β -, β + are set to 0.7 and 1.25, respectively.

In the next sections, the following aspects will be illustrated through the numerical simulations:

The benefit from using nonlinear models To start, Fig. 4(a) shows the relation between the stationary values of Ttb 130 and ∆P 156 . From this figure, it appears clearly that the visited interval of values of ∆P 156 involves rather high nonlinearities which explains the benefit from using Nonlinear MPC. Indeed, a comparison between the coordination with Linear MPC and Nonlinear MPC used in the Brayton cycle is presented (Fig. 5). It is shown that the temperature Ttb 130 is better controlled when using the nonlinear models which results in an obvious improvement of the cost. This validates the first purpose of this contribution (need for explicit handling of the nonlinearities). Therefore, only the Nonlinear MPC design in the Brayton cycle is used in the forthcoming investigations. However, the related computation time of S 2 exceeds the time limit (subplot [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Doan | A Hierarchical MPC Approach with Guaranteed Feasibility for Dynamically Coupled Linear Systems[END_REF] of Fig. 5), which means that the simulated feedback actions cannot be really implemented. In what follows, the simulation results aim at showing the benefit of using the hierarchical control framework (compared to observer-based decentralized framework) and the effectiveness of the distributed-in-time optimization heuristic, proposed in the previous section, in addressing the computation time and real-time implementability issue. The first row presents the outputs, and the second one presents the inputs of the system. The subplot [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Doan | A distributed optimizationbased approach for hierarchical MPC of large-scale systems with coupled dynamics and constraints[END_REF] illustrates the periodic disturbance applied to the helium bath. The subplot [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Doan | A Hierarchical MPC Approach with Guaranteed Feasibility for Dynamically Coupled Linear Systems[END_REF] shows the computation time when using the Nonlinear MPC. Note also that the computation time limit is not considered for this simulation (the solver takes the time needed to perform the assigned computations). The subplot [START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Stewart | Hierarchical cooperative distributed model predictive control[END_REF] depicts the evolution of the cost index J sim in time.

The benefit from hierarchical design

In this section, the disturbance rejection mode is simulated. The closed-loop cost of the proposed hierarchical framework is compared to the one obtained under the extended observer-based decentralized approach. More precisely, two scenarios are simulated in which the limitation on the available computation time is respectively enforced or not:

Without constraints on the computation time Fig. 6 shows simulations without limitation in the computation time (the decided number of evaluations is assumed to be possible to achieve within the sampling period T s ). This means that the corresponding u opt s (k) is implemented even its computation time exceeds the available time within the sampling period. With the use of the proposed framework, the regulation of the temperature Ttb 130 is highly enforced compared to its behavior in the absence of coordination. One can also note that the level Ltb 131 seems to loosely track the desired set-point which is expected since in this disturbance rejection mode 1, the focus is on the temperature Ttb 130 as suggested by the choice of the penalty matrices [see ( 53) and ( 54)]. This can be simply checked by examining the last plot that compares the closed-loop index evolutions for the two examined settings. The comparison with the observer-based decentralized approach shows the advantage of using the hierarchical scheme should the computation be possible to perform in the allowed sampling period T s .

Figure 6: Closed-loop behavior with coordination (solid blue line) and without coordination (dash-dot green line). The first row presents the outputs, and the second one presents the inputs of the system. The hierarchical control method gives a better cost J sim than the decentralized method (closed-loop cost decreased by 56%). Note that the cost under linear MPC (Fig. 5) corresponds to a reduction of only 12% compared to the decentralized control without coordination). Note also that the computation time limit is not considered in this simulation.

Enforcing the cpu-constraint (τ u ≤ T s ) on the available computation time

In this section, the constraint induced by the limited available computation time is taken into account in order to underline the benefit from using the proposed distributed-in-time optimization framework. In order to take into account the cost increase induced by the impossibility to compute the solution r opt within the allocated time, Fig. 7 shows the closed-loop behavior when the subsystems apply the previous control u opt s (k -1) each time the computation time exceeds the available computation time T s . Indeed, in this case, the master cannot dispose of the needed information in order to update the approximation of the cost function, which is needed to update the value of the auxiliary set-point and the associated coupling signals that are needed to compute the updated control to be applied. for all the settings mentioned previously. This table clearly shows the cost increase that is induced by the limitation on the available computation time. Indeed, the cost increases from 44% to 70.64% of the baseline represented by the observer-based decentralized scheme.

In the next section, the possibility to partially recover the optimal performance through the distributed-in-time optimization scheme proposed in section 4 is investigated.

Impact of the distributed-in-time setting's parameters

Different configurations of the distributed-in-time optimization parameters (n z , n d and τ u ) are simulated, and the corresponding closed-loop costs J norm sim (normalized cost index) are reported in order to give a flavour of the impact of each choice on the results. More precisely, the testing scenario of a periodic heating disturbance is simulated again with the distributed-in-time optimization framework being implemented in the hierarchical framework. Fig. 8 shows the comparison of the behavior of the process when using the hierarchical control combined with optimization distribution and the one without distribution under time limit constraint. Fig. 9 shows the computation time for different configuration of n z and n d . Finally, table 3 shows the associated closed-loop cost for the different configurations of the distributed-in-time optimization scheme. It comes out that the setting corresponding to n z = 2, n d = 1 enables to get closer to the ideal cost index drop (44%) corresponding to the non constrained computation time simulation while being fully real-time compatible. Checking modularity: controlling the system by only tuning the central cost's definition

One of the claims of this contribution concerns the possibility to keep the local controllers unchanged (in terms of penalty) while changing the penalties of the economic cost (by the coordinator) in order to achieve different behaviors of the closed-loop system. Here, it is assumed that the operator need to change the set-point of the helium liquid level Ltb 131 , for example, in order to embed a test facility below the liquid level. Thus, the reference tracking scenario is simulated. More precisely, the closed-loop behavior is compared under the two different modes defined in section 5.1. This is done in order to illustrate the fact that the results can be affected in the desired direction by only modifying the central cost's definition while keeping the local controllers unchanged.

For the tracking set-point mode, Fig. 10 shows the comparison of the behavior of the process between using hierarchical control combined with optimization distribution and using decentralized control.

Conclusion and future works

This paper extends a recently proposed hierarchical control design framework with application to a the control of cryogenic refrigerator. It is particularly shown that the proposed framework enables the incorporation of constraints and nonlinearity in the involved models of subsystems. Moreover, a dedicated distributed-in-time optimization scheme is proposed and validated. This mechanism enables to recover a great part of the loss of optimality that might be induced by the impossibility to achieve the associated computation within the available time.

The undergoing investigations aim at pushing forward the splitting of the system into subsystems. In particular, the heat exchangers typically introduce many states with mainly linear relationship except at their boundary conditions. It should therefore be possible to split the subsystems they belong to into two subsystems of two categories: linear high dimensional ones and nonlinear small dimensional ones. However, this implies that the framework should be adapted to the case where some of the subsystems do not have controlled inputs and/or regulated outputs. Some preliminary results are quite encouraging and might be proposed in future communication.

Fig. 1 (Figure 1 :

 11 Figure 1: Block diagram of the 400@1.8K experimental refrigerator (in the 400W@4.4K configuration) at CEA/IRIG/DSBT [3]. Fig.1(a) shows the system can be decomposed into four parts: the Warm compression station (WCS), Precooling stage, Brayton cycle and Joule-Thomson cycle. Fig. 1(b) illustrates the ideal entropy -temperature diagram, which is implemented by the 400@1.8K refrigerator.

  first component is the control input inside S 1 and the second one is used as a disturbance signal. The value of NCR (a) 22 is in the range of [0, 100] W. 3. ∆P 156 ∈ [0, 10] bar: The pressure drop between the inlet pressure and outlet pressure of the valve CV 156 .

  Coordinator constructs a grid G of auxiliary set-points Fixed-point Iteration Coordinator evaluates J c for each set-point r ∈ G Optimization 1. Approximate J c by Ĵc 2. Compute r opt = min r∈R Ĵc Finalization Coordinator sends r opt to the subsystems for them to compute the control to be applied on the plant. Communication between the coordinator and the subsystems.

Figure 2 :

 2 Figure 2: block diagram of the hierarchical control algorithm.

Figure 4 :

 4 Figure 4: Fig. 4(a) shows the evolution of the stationary value of Ttb 130 as a function of ∆P 155 . Fig. 4(b) shows the histograms of the computation times of S 1 and S 2 during the simulated scenario.

Fig. 4 (

 4 Fig. 4(b) compares the computation time of S 1 and S 2 showing that since the subsystem S 1 only involves linear MPC problem, the related computation time is very small compared to the time needed to solve general non quadratic optimization problem that is associated to subsystem S 2 . The same can be said regarding the computations done by the coordinator, in view of the fact that it only performs basic operations.

Figure 5 :

 5 Figure 5: Comparison of closed-loop behavior under the proposed hierarchical framework with Linear MPC (dash-dot green line) and Nonlinear MPC (solid blue line) used in modeling the Brayton cycle.The first row presents the outputs, and the second one presents the inputs of the system. The subplot[START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Doan | A distributed optimizationbased approach for hierarchical MPC of large-scale systems with coupled dynamics and constraints[END_REF] illustrates the periodic disturbance applied to the helium bath. The subplot[START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Doan | A Hierarchical MPC Approach with Guaranteed Feasibility for Dynamically Coupled Linear Systems[END_REF] shows the computation time when using the Nonlinear MPC. Note also that the computation time limit is not considered for this simulation (the solver takes the time needed to perform the assigned computations). The subplot[START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Stewart | Hierarchical cooperative distributed model predictive control[END_REF] depicts the evolution of the cost index J sim in time.

Figure 7 :

 7 Figure7: Closed-loop cost comparison between the case where the limitation of the available computation time Ts is taken into account in the hierarchical framework's implementation (dash-dot green line) and the ideal case where this constraint is not considered (solid blue line). The first row presents the outputs, and the second one presents the inputs of the system. Note in particular how the control of the liquid helium level Ltb 131 is visibly deteriorated.

Figure 8 :

 8 Figure 8: Comparison of the closed-loop behavior under the proposed hierarchical framework with two different settings: 1) with distributed-in-time optimization being implemented and 2) without distributed-in-time implementation under available time limitation constraint. The first row presents the outputs and the second one presents the inputs of the system. The parameters nz = 2, n d = 1 and τu = 5s are used in the distributed-in-time optimization framework. subplot[START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Alamir | Fixed-point based hierarchical MPC control design for a cryogenic refrigerator[END_REF][START_REF] Doan | A Hierarchical MPC Approach with Guaranteed Feasibility for Dynamically Coupled Linear Systems[END_REF] shows the computation time associated to this configuration which now meet the constraint on the available computation time.

Figure 9 :

 9 Figure 9: Histogram of computation time of different choices of nz and n d . Note that the computation time of each configuration is always lower than the updating time τu.

Figure 10 :

 10 Figure 10: Set-point tracking scenario: closed-loop responses under coordination, using distributed-in-time optimization in two different mode for the centralized cost on one hand and without coordination on the other hand. The first row presents the outputs, and the second one presents the inputs of the system. The set-point on Ltb 131 is increased. Two configurations of Qc and Rc of Mode 1 and Mode 2 are tested. Mode 2 (corresponding to higher penalty on Ltb 131 deviations) allows better reference tracking while mode 1 which is dedicated to disturbance rejection and not especially to track set-point on the level. With the set of parameters nz = 2, n d = 1 and τu = 5s. Note that both hierarchical design with distributed optimization are real-time compatible.

Table 1

 1 shows the nominal values of the inputs and outputs described above.

	Variable	Value	Unit	Variable	Value	Unit
	CV 0 155 ∆P 0 156 NCR 0 22	55 1.5 310	% bar W	Ttb 0 130 Ttb 0 108 Ltb 0 131	12.31 5.37 60.5	K K %
	P Hw	16	bar	P Cw	1.05	bar

Table 1 :

 1 Steady state values of the inputs and outputs. P Hw and P Cw represent respectively the high and the low pressures controlled by the warm compression station.

Table 2 :

 2 The normalized cost index of for the different settings mentioned above.

		J norm sim
	Without Coordination	100 %
	Coordination with Linear MPC	87.57 %
	Coordination with Nonlinear MPC (Enforcing computation time limit)	70.64 %
	Coordination with Nonlinear MPC (Without computation time limit)	44.21%
	J norm sim

Table 3 :

 3 The normalized cost J norm sim for different configurations nz, n d and τu of the distributed-in-time optimization. n z n d τ u t max J norm

	sim

The system is located at CEA/IRIG/DSBT, Grenoble