Fixed-point Based Hierarchical MPC Control Design For a Cryogenic Refrigerator

TL;DR

This paper introduces a fixed-point based hierarchical model predictive control framework for cryogenic refrigerators, enabling flexible, stable, and smooth control of interconnected thermodynamic cycles without redesigning local controllers.

Contribution

It presents a novel hierarchical MPC scheme that handles destabilizing interconnections and allows easy behavior modifications by adjusting coordinator parameters.

Findings

Successfully stabilizes cryogenic refrigerator cycles

Allows behavior changes without local controller modifications

Enables smooth operator handover

Abstract

In this paper, a simple and general hierarchical control framework is proposed and validated through the interconnection of the Joule-Thompson and the Brayton cycle stages of a cryogenic refrigerator. The proposed framework enables to handle the case of destabilizing interconnections through state and/or control signals (which is the case of the cryogenic refrigerator example). Moreover, it offers the possibility to simply change the behavior of the overall system (depending on the context) by only changing the coordinator problem's parameters without changing the set of local controllers used by subsystems which is a common industrial requirement regarding industrial control architectures. Finally, the proposed scheme enables a smooth operator handover on a specific subsystem and/or actuator.

Equations

$${}_{155}^0$$

$$u_1:=\left[\begin{array}{c}C{V}_{155}-C{V}_{155}^0\\ NC{R}_{22}^{(a)}-NC{R}_{22}^{(a,0)}\end{array}\right];u_2:=C{V}_{156}-C{V}_{156}^0$$

$$y_1:=\left[\begin{array}{c}Lt{b}_{131}-Lt{b}_{131}^0\\ Tt{b}_{108}-Tt{b}_{108}^0\end{array}\right]x;y_2:=Tt{b}_{130}-Tt{b}_{130}^0$$

$$w=NC{R}_{22}^{(w)}-NC{R}_{22}^{(w,0)}$$

$$U_1^0:=\left[\begin{array}{c}C{V}_{155}^0\\ NC{R}_{22}^{(a,0)}\end{array}\right],U_2^0:=C{V}_{156}^0,Y_1^0:=\left[\begin{array}{c}Lt{b}_{131}^0\\ Tt{b}_{108}^0\end{array}\right],Y_2^0:=Tt{b}_{130}^0$$

$$x_1^{+}$$

$$x_2^{+}$$

$$v_1$$

$$v_2$$

$$y_1$$

$$y_2$$

$$v_1:=(P_h,T_h,P_c{)}^T\in {\mathbb{R}}^3$$

$$v_2:=(M_h,T_c,M_c{)}^T\in {\mathbb{R}}^3$$

$$x_s^{+}$$

$$v_1$$

$$v_2$$

$$y_s$$

$${\mathbf{v}}_1(k+i)$$

$${\mathbf{v}}_2(k+i)$$

$${\mathbf{v}}_1={\mathbf{g}}_1({\mathbf{u}}_2,x_2(k),{\mathbf{v}}_2)$$

$${\mathbf{v}}_2={\mathbf{g}}_2({\mathbf{u}}_1,x_1(k),{\mathbf{v}}_1)$$

$$J_c(\mathbf{u},x(k)):=\sum _{s=1}^2\left[\sum _{i=1}^N{\left(\frac{i}{N}\right)}^q\left[{\Vert y_s(k+i)-r_s^d\Vert }_{Q_c^{(s)}}^2+{\Vert U_s^0+u_s(k+i)\Vert }_{R_c^{(s)}}^2\right]\right]$$

$$J_c^{ext}(\mathbf{u},\mathbf{v},r^d,x(k))=\sum _{s=1}^2{J}_1({\mathbf{u}}_s,r_s^d,x_s(k)|{\mathbf{v}}_s)$$

$$\displaystyle:=\sum_{s=1}^{2}\Bigl{[}\sum_{i=1}^{N}\left(\dfrac{i}{N}\right)^{q}\left[\left\|h_{s}(k+i)-r_{s}^{d}\right\|_{Q_{c}^{(s)}}^{2}+\left\|U_{s}^{0}+u_{s}(k+i)\right\|_{R_{c}^{(s)}}^{2}\right]\Bigr{]}$$

$$h_{s}(\sigma):=h_{s}\Bigl{(}\bm{X}_{s}(\sigma,\bm{u}_{s},x_{s}(k)|\bm{v}_{s}),\bm{u}_{s}(\sigma),\bm{v}_{s}(\sigma),0\Bigr{)}$$

$${\mathbf{v}}_1={\mathbf{g}}_1({\mathbf{u}}_2,x_2(k),{\mathbf{v}}_2)$$

$${\mathbf{v}}_2={\mathbf{g}}_2({\mathbf{u}}_1,x_1(k),{\mathbf{v}}_1)$$

$$\displaystyle\min_{\bm{u}_{1},\bm{u}_{2},\bm{v}_{1},\bm{v}_{2}}\sum_{s=1}^{2}J_{s}(\bm{u}_{s},r_{s}^{d},x_{s}(k)|\bm{v}_{s})\quad\mbox{\rm under}\quad\left\{\begin{array}[]{l}\bm{v}_{1}=\bm{g}_{1}(\bm{u}_{2},x_{2}(k),\bm{v}_{2})\\ \bm{v}_{2}=\bm{g}_{2}(\bm{u}_{1},x_{1}(k),\bm{v}_{1})\end{array}\right.$$

$${\mathbf{v}}_1^{(0)},{\mathbf{v}}_2^{(0)}$$

$$J_s({\mathbf{u}}_s,r_s^d,x_s(k)|{\mathbf{v}}_s^{(0)})$$

$${\mathbf{u}}_s^{\ast }(r_s^d,x_s(k)|{\mathbf{v}}_s^{(0)}),J_s^{\ast }(r_s^d,x_s(k)|{\mathbf{v}}_s^{(0)})$$

$$\displaystyle\mbox{\rm subsystem $S_{1}$ send to the coordinator}\quad\rightarrow\quad\hat{\bm{v}}_{2}^{(1)}:=\bm{g}_{2}(\bm{u}_{1}^{*},x_{1}(k),\bm{v}_{1}^{(0)})$$

$$\displaystyle\mbox{\rm subsystem $S_{2}$ send to the coordinator}\quad\rightarrow\quad\hat{\bm{v}}_{1}^{(1)}:=\bm{g}_{1}(\bm{u}_{2}^{*},x_{2}(k),\bm{v}_{2}^{(0)})$$

$${\mathbf{v}}_s^{(1)}:=(1-\beta ){\mathbf{v}}_s^{(0)}+\beta {\hat{\mathbf{v}}}_s^{(1)}$$

$$J_s^{MPC}({\mathbf{u}}_s,r_s,x_s(k)|{\mathbf{v}}_s^{(\sigma )}):=\sum _{i=1}^N\Vert x_s(k+i)-x_s^d(r_s)){\Vert }_{Q_s}^2+\Vert {\mathbf{u}}_s(k+i)-u_s^d(r_s)){\Vert }_{R_s}^2$$

$${\mathbf{u}}_s^{opt}(r_s,x_s(k)|{\mathbf{v}}_s^{(\sigma )})$$

$$\displaystyle\mbox{\rm$S_{1}$ send to the coordinator}\quad\rightarrow\quad\hat{\bm{v}}_{2}^{(\sigma+1)}:=\bm{g}_{2}(\bm{u}_{1}^{opt}(r_{1}),x_{1}(k),\bm{v}_{1}^{(\sigma)})$$

Full text

††journal: Journal of Process Control

1 Introduction

Nowadays it is widely admitted that non centralized control architectures are to be preferred, when possible, to fully centralized ones [15, 9]. Many reasons for this can be invoked including breaking large problems into small tractable ones, robustness, modularity, privacy preservation and easiness of maintenance.

Non centralized architectures can be mainly divided into distributed and hierarchical. In the former, no coordinator is used and communication takes place between adjacent neighbors while in the latter, subsystems communicate exclusively with a coordinator which attempts to coordinate the subsystems behavior in order to achieve a global task or to optimize some globally defined performance index.

Excellent reviews of the state of the art regarding these architectures are given in [15, 9] as far as Model Predictive Control (MPC) design is concerned. These reviews show a highly dynamic research areas. In particular, the survey book [9] enumerates no less than 35 distinguishable approaches. This makes hard for any new proposed solution to claim full novelty. Instead, let us summarize the items that might help positioning the hierarchical solution proposed in this paper w.r.t existing works on hierarchical control architectures111As far as process control is concerned, hierarchical architectures are preferred as there is always some global goal to achieve such as global efficiency and/or product quality [12, 16]. Distributed architectures are more adapted to manage the behavior of a set of similar agents (robots, vehicles, drones, etc.) although some distributed schemes can also be used in process control.:

• $\checkmark$

Many works focus on enforcing the optimality and/or constraints fulfillment by assuming the existence of a fully decentralized stabilizing feedback [1, 7, 6]. This strong assumption is not required here nor is it valid for the cryogenic process under interest. On the other hand, our stability result does not theoretically hold in the presence of constraints although it can be technically applied with potential success. This obviously makes the two families of solutions non comparable as on one hand, many processes does not fit the fully centralized stabilizability assumption while on the other hand, many processes need constraints to be explicitly handled. Now it is a fact that the situation in the process industry leaves dozens if not hundreds of non coordinated PID-like controllers that do not address the constraints anyway nor they guarantee any kind of stability even in the presence of relevant scenarios.

• $\checkmark$

In some works, a dedicated assumptions regarding the nature of the coupling signals are introduced in order to avoid strong destabilizing coupling effects. For instance, the framework of [17, 10] assumes only coupling through control input actions while in [11] a two layer hierarchical structure is studied where the the higher layer system is slow with control input given by the lower level. Such assumptions circumvent the major difficulties in strongly coupled systems where a typical fixed-point iterations lies beneath with a potential risk of non convergence during the iterations between the coordinator and the subsystems. This mechanism is clearly highlighted in the paper and its convergence is analyzed.

• $\checkmark$

A large part of the literature borrows the viewpoint according to which non centralized frameworks are mainly useful to break down the complexity/dimensionality of the centralized problem. In this perspective, the full knowledge of the system is still assumed in the coordinator level but the computation is broken into several parallel parts that are achieved by the subsystems. This approach does not address an important industrial concern according to which, it should be possible to change one subsystem locally (switch to a new generation of valves or change the logic of the local controller from PID to MPC for instance) without being obliged to propagate the consequences of this change to every line of algorithm at the remaining local controllers and/or the coordinator’s level.

• $\checkmark$

Finally, putting aside the always questionable novelty assertion with regards to the huge literature on the topics, it is a fact that too many proposed schemes are validated through toy examples such as steered tank reactor or coupled water tanks, etc222A nice counter example, among a few others, where a relevant and challenging example is handled can be found in [10] where hierarchical framework is applied to the Barcelona’s drinking water network.. The example addressed in this contribution is quite significant: two coupled systems with $10$ and $14$ states respectively, each one acts on the other through $3$ coupling signals including destabilizing coupling effect. All these characteristics make the example a nice benchmark for existing and future hierarchical control architectures and algorithms. To this respect, the paper can be viewed as an introduction to such a benchmark together with a first successful hierarchical control architecture that might be compared to future candidate hierarchical design.

This paper is organized as follows: First of all, the system under study is presented and the control objective is explained in Section 2. This presentation enables the general setting and the proposed solution sketched in Section 3 to be better understood. Section 4 gives the results obtained by means of the proposed hierarchical framework on the cryogenic refrigerator. Finally, Section 5 gives further discussion regarding the proposed solution and how it can be used to smoothly handle the case of operator handover at several levels. Finally, the paper ends by Section 6 which summarizes the paper and gives some hints for further investigations.

2 Problem description

2.1 Presentation of the cryogenic refrigerator

Cryogenic refrigerators are used in the experimental facilities containing supra-conducting circuits in order to provide cooling power to cool them down. Figure 1 shows a sketch of the 400W@1.8 K experimental refrigerator (in the 400W@4.4 K configuration) at CEA333Commissariat à l’Energie Atomique /INAC444Institut Nanoscience et Cryogénie/SBT555Service des Basses Températures [14]. This figure shows the system to be controlled which is composed of subsystems $S_{1}$ (The Joule-Thomson cycle) and $S_{2}$ (The Brayton cycle).

Generally speaking, a cryogenic refrigerator implements a thermodynamic cycle which produces cooling power by extracting work from a fluid by means of cryogenic turbine (the one following the valve CV156 in $S_{2}$ of Figure 1) and by exchanging heat power through a series of heat exchangers (denoted by NEFx in Figure 1). The main objective is to absorb the disturbing heat power generated by the operation of the experimental facility which is represented in the 400W experimental facility by the heat source denoted by NCR22 in system $S_{1}$ of Figure 1.

Note that in order to close the thermodynamic cycle, a compressor is used in the so-called warm zone (composed of subsystems $S_{3}$ and $S_{4}$ on Figure 1) with the objective to maintain pressures (denoted by $H_{p}$, $B_{p}$ in Figure 1) at some prescribed level which is presumed to adapted to a well-conditioning of the overall system. For more details regarding the principle and the modeling of cryogenic plants, interested reader can refer to [5, 4].

In this contribution, it is assumed that the warm zone (systems $S_{3}$ and $S_{4}$) is appropriately controlled (see for instance [3]) and we focus 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 section.

2.2 The control problem of the cold zone $S_{1}$-$S_{2}$

Using the notation of Figure 1, the control problem of the cold zone including $S_{1}$ and $S_{2}$ is described hereafter by successively describing the control inputs, the regulated outputs, the dynamical model and the control objective:

2.2.1 The control inputs

There are three control inputs, two of them belong to the Joule-Thomson cycle ($S_{1}$) and one belongs to the Brayton cycle ($S_{2}$), more precisely (see Figure 1):

CV155: the position ($\in[0,100]$) of the valve situated at the input of the Helium bath (Figure 1). This valve is referred to as the JT-valve. 2. 2.

NCR${}_{22}^{(1)}$ : the heating power (W) of the resistance inside the Helium bath. Recall that this same actuator is used to simulate heat pulses coming from the operation of the physical experimental facility served by the cryogenic refrigerator. This explains the upper index used in NCR${}^{(1)}_{22}$ as the whole power is split into two components NCR22=NCR${}^{(a)}_{22}$+NCR${}^{(w)}_{22}$, the first of which is used as a control input inside $S_{1}$ while the second is used as a disturbance signal. Note that although an observer can be built to estimate the disturbance NCR${}^{(w)}_{22}$ (see [2] for instance), this disturbance is considered here as an unmeasured disturbance because of the difficulties to predict its future behavior in an MPC design666Indeed, the temporal characteristics of the heat pulse depend on the experiment being conducted on the experimental facility. This is required to be unknown by the cryogenic refrigerator.. 3. 3.

CV156: the position ($\in[0,100]$) of the valve situated at the inlet of the cryogenic turbine of $S_{2}$.

2.2.2 The regulated outputs

Three output variables need to be regulated, as many as control inputs in each subsystem. More precisely:

Ltb131: the level (%) in the liquid helium bath. This variable represents the thermal storage that is immediately available for cooling purposes. 2. 2.

Ttb108: the temperature ($K$) at the inlet of the JT-valve in subsystem $S_{1}$ (see Figure 1). The tight regulation of this variable is crucial in order to maintain the efficiency of the process. 3. 3.

Ttb130: the temperature ($K$) at the outlet of the cryogenic turbine. Maintaining this temperature above a lower limit is mandatory in order to avoid the formation of liquid droplets that can break the turbine.

2.3 The nominal operation point

A nominal operation point is defined which is considered to be optimal given the sizing of the refrigerator. This operational point is defined by the stationary values of the inputs, outputs and disturbance described above, namely:

[TABLE]

These values enable the deviations in the inputs, outputs and disturbances to be defined according to:

[TABLE]

[TABLE]

and

[TABLE]

In what follows, the following notation is used to denote the values at the operation point:

[TABLE]

2.3.1 The dynamic model

Using the index $s\in\{1,2\}$ to refer to the subsystem $S_{s}$, the linearized model of the two subsystems, defined around the operation point described above takes the form:

[TABLE]

in which

• $\checkmark$

$x_{1}\in\mathbb{R}^{10}$ and $x_{2}\in\mathbb{R}^{14}$ are the states of subsystems $S_{1}$ and $S_{2}$ respectively.

• $\checkmark$

$y_{1}\in\mathbb{R}^{2}$ and $y_{2}\in\mathbb{R}$ are the deviations of regulated outputs described above.

• $\checkmark$

$v_{1}\in\mathbb{R}^{3}$ and $v_{2}\in\mathbb{R}^{3}$ are coupling signals that make the two systems interact dynamically as $v_{1}$ depends on $(x_{2},u_{2})$ while $v_{2}$ depends on $(x_{1},u_{1})$. As a matter of fact $v_{1}$ and $v_{2}$ both depend on the whole state $x$ and the control $u$ as $v_{1}$ depends also on $v_{2}$ and vice-versa [see (8) and (9)].

More precisely, using the notation of Figure 1, the coupling variable $v_{1}$ is given by

[TABLE]

where $P_{h}$ and $T_{h}$ are the pressure and temperature (deviations) at the downstream inlet of the heat exchanger NEF1 while $P_{c}$ stands for the pressure (deviation) at the upstream outlet of the same heat exchanger.

Similarly, the coupling variable $v_{2}$ is given by:

[TABLE]

where $M_{h}$ is the mass flow rate (deviation) at the downstream inlet of the heat exchanger NEF1 while $T_{c}$ and $M_{c}$ stand respectively for the temperature and the mass flow rate (deviations) at the upstream outlet of the same heat exchanger.

From these equations, it comes clearly that the two subsystems $S_{1}$ and $S_{2}$ are strongly coupled since the dynamics of $S_{1}$ defined by (6) depends on the variable $v_{1}$ which depends through (8) on $x_{2}$, $u_{2}$. The same can be said about the dynamics of $S_{2}$ which depends on $S_{1}$ through the coupling variable $v_{2}$ which depends through (9) on $x_{1}$ and $u_{1}$.

In the sequel, the equations (6)-(11) are sometimes shortly rewritten in the following compact form when linearity is not explicitly involved:

[TABLE]

Note that for any initial state $x(k)$ and any control profile $\bm{u}$ defined over some prediction horizon of length $N$, the corresponding nominal [disturbance-free] trajectories $\bm{X}(\cdot,\bm{u},x(k)|\bm{v})$ that are obtained by integrating the dynamics (14)-(16) enables the following so-called coherence constraints to be defined:

[TABLE]

which are obviously two conditions on the coupling variable profiles $\bm{v}:=(\bm{v}_{1},\bm{v}_{2})$ to be compatible with the system’s nominal (disturbance-free) coupled equations. This condition can be shortly written as follows:

[TABLE]

with obvious appropriate notation. These last conditions are referred to as the coherence constraints. IN what follows, the obvious notation $x=(x_{1}^{T},x_{2}^{T})^{T}$, $y=(y_{1}^{T},y_{2}^{T})^{T}$, $v=(v_{1}^{T},v_{2}^{T})^{T}$ are used.

2.3.2 The control objective

There are two modes that define the control objective:

In the first mode which is the by-default mode, the control has to regulate the system $S_{1}\cup S_{2}$ around $x=0$ (the nominal steady state operation point invoked above) despite of the presence of non measured disturbance. This is a disturbance rejection mode representing the raison d’être of the cryogenic refrigerator. 2. 2.

In the second mode, the operator should be able to temporary steer the system to a different steady state that corresponds to some new set-point $y\neq 0$. For instance, the operator would like to change the level of the liquid helium in the bath (new value of the set-point on Ltb131) or on the critical temperature Ttb108.

These two modes can be accounted for by using different set-points and different weighting matrices in the following centralized cost function:

[TABLE]

where $\bm{u}$ is a candidate sequence of control inputs to be applied starting from the initial state $x(k)$ over the prediction horizon $[k,k+N]$. Note that the output set-points $r_{s}^{d}$, for $s=1,2$ represent the desired values of the deviation on the output corresponding to the set-point $r_{s}^{d}+Y_{s}^{0}$. Finally, $Q_{c}^{(1)}\in\mathbb{R}^{2\times 2}$, $Q_{c}^{(2)}\in\mathbb{R}_{+}$, $R_{c}^{(1)}\in\mathbb{R}^{2\times 2}$ and $R_{c}^{(2)}\in\mathbb{R}^{+}$ are non negative weighting matrices.

Remark 1

Note that the time-dependent weighting term $(i/N)^{q}$ for some $q\in\mathbb{N}$ enables to put higher weight on the tail of the prediction horizon by taking high value of $q$. Note that taking $q=0$ recover the standard time-invariant stage cost.

Remark 2

It is worth underlining that in the cost function (20), the total value of the control inputs are penalized and not the deviation between the input and the steady input that is typically used in a standard Model Predictive Control formulation. This is because in common MPC formulations where the control is penalized through a term of the form $\|u_{s}(k+i)-u_{s}^{d}\|$, the stage cost at the desired position is always [math] which obviously does not represent the difference in the real cost that might correspond to each different set-point. Another difference lies in the fact that in order to guarantee the stability of MPC scheme, the whole state should be penalized using a term of the form $\|x_{s}(k+i)-x_{s}^{d}(r_{s}^{d}))\|$ (see [8] for more details). This again leads to a cost function that does not rigorously represent the economic/performance criterion one would intuitively like to minimize. We will see that such formulations which are necessary for the stability of the resulting MPC will be used in the formulation of the local MPC controllers while the realistic and relevant cost function (20) is used to formulate the centralized problem’s cost. Such considerations are intimately linked to the concept of economic Model Predictive Control [13] although the latter is generally studied in a non hierarchical framework.

Based on the above definitions and notation, the control problem can be stated as follows:

Problem statement

Define a hierarchical control scheme in which two local MPC controllers, defined for $S_{1}$ and $S_{2}$ respectively, receive appropriate set-points $r_{1}$ and $r_{2}$ (send by a coordinator) that minimize the centralized cost function (20) defined for some desired set-point $r^{d}$. Moreover, the following conditions should be satisfied:

$S_{1}$ and $S_{2}$ exclusively communicate with the coordinator 2. 2.

The coordinator ignores the details of the mathematical models inside $S_{1}$ and $S_{2}$ and the details of their controllers 3. 3.

The different modes described above are handled by simply changing the desired set-point $r^{d}$ and/or the weighting matrices $Q_{c}^{(s)}$ and $R_{c}^{(s)}$ of the centralized problem without changing the tuning of the local controllers.

As it is explained later on, there is generally no reason to have $r=r^{d}$ as the local costs and the central sub-costs are quite different. The optimal local set-points $r_{s}$, $s=1,2$ are to be computed by the coordinator in order to minimize the central cost corresponding to the original set-point $r^{d}$.

In the next section, a solution is proposed for the above hierarchical control problem in the general case before it is validated on the cryogenic refrigerator in Section 4. It is worth underlying that while the presentation of the next section considers only two subsystems for the sake of clarity and in order to simplify the notation, the framework is obviously applicable almost without change to the case of a networked of several coupled subsystems. Only the notation would be heavier which we prefer to avoid here.

3 General setting and proposed solution

In this section, it is first shown (section 3.1) that a fixed point-like iteration can be defined so that if convergence arises, the coordinator can get a hierarchical estimation of the central cost for a candidate value $r=(r_{1},r_{2})$ of the set-point to be sent to the subsystems. For the sake of clarity of exposition, it is first assumed that this convergence occurs systematically and the way the loop is closed to get a hierarchical feedback control addressing the hierarchical control problem stated in the previous section is shown (Section 3.2). The section is then ended by the analysis of the convergence of the fixed point iteration (Section 3.3).

3.1 Hierarchical estimation of the central cost by fixed point negotiation

Let us consider two subsystems $S_{1}$ and $S_{2}$ which are described by the coupled system of equations given by (14)-(16). Consider also the hierarchical control problem stated at the end of the previous section with the central cost (20). During this subsection, the states $x_{s}(k)$, the central set-point $r_{s}^{d}$ and the auxiliary individual set-points $r_{s}$ are supposed to be given and frozen.

Note first of all that the central optimization problem (20) can be redefined by considering the extended vector of degrees of freedom $(\bm{u},\bm{v})$ as follows:

[TABLE]

where $h_{s}(k+\cdot)$ represents the nominal output profile defined by (17):

[TABLE]

provided that the coupling signal profiles $\bm{v}=(\bm{v}_{1},\bm{v}_{2})$ satisfy the coherence constraints (18)-(19) stated in the previous section:

[TABLE]

Note that in the above cost (21), the two individual costs are now conceptually decoupled for any given choice of the coupling signal profiles $\bm{v}_{1}$ and $\bm{v}_{2}$. The coupling appears now in the equality constraints (24)-(25).

Consequently, the central optimization problem can be rewritten in the following form:

[TABLE]

The last form (28) of the central optimization problem is already amenable to hierarchical distribution although it is not the one that will be finally used. Indeed, the coordinator can begin with some initial guesses:

[TABLE]

of the coupling variables profiles. It sends $\bm{v}_{s}^{(0)}$ to subsystem $S_{s}$. Each subsystem is now able to compute its own contribution to the central cost, namely:

[TABLE]

for any candidate control profile $\bm{u}_{s}$. Therefore, each subsystem can solve its own optimization problem and get its own optimal profile and cost:

[TABLE]

Moreover, each subsystem sends what would be the corresponding coupling profile it would apply in this case and send it to the coordinator:

[TABLE]

where $\bm{u}_{s}^{*}:=\bm{u}_{s}^{*}(r_{s}^{d},x_{s}(k)|\bm{v}_{s}^{(0)})$.

Having the new estimations $\hat{\bm{v}}_{s}^{(1)}$, the coordinator elaborates a filtered version of the coupling profiles:

[TABLE]

and therefore, a fixed-point round of iterations can take place.

The reason why the fixed-point iteration will not be defined as explained above lies in the fact that the centralized cost formulation, when used to built up an MPC feedback design, does not necessarily lead to a stable behavior of the system. Indeed, as it is noticed earlier, the stability of a finite horizon MPC schemes requires some conditions on the formulation of the cost function (penalizing the whole state, zeros cost at the desired state, positive definite penalty on the control action and so on). Now modifying the central cost in order to meet these requirements may lead to an irrelevant cost from the economic/performance point of view.

Instead, stabilizing formulations are used in the local level in order to enforce stability while only the local set-points $r_{s}$ sent by the coordinator to the subsystems are kept as decision variables for the coordinator. This is explained in the remainder of the current section.

More precisely, for each subsystem $s$, given a set-point $r_{s}$ (not necessarily equal to the central set-point $r_{s}^{d}$) together with a current estimation of the coupling profile777$\sigma=0$ was used in the previous discussion $\bm{v}^{(\sigma)}$, a local, standard well defined and stability inducing cost function is defined:

[TABLE]

where

• $\checkmark$

$(x_{s}^{d}(r_{s}),u_{s}^{d}(r_{s}))$ is the steady pair that is compatible with the output set-point $r_{s}$.

• $\checkmark$

$x_{s}(k+i):=\bm{X}_{s}(k+i,\bm{u}_{s},x_{s}(k)|\bm{v}_{s}^{(\sigma)})$ is the nominal predicted state at instant $k+i$ given the initial state $x_{s}(k)$, the control profile $\bm{u}_{s}$ and the presumed coupling profile $\bm{v}_{s}^{(\sigma)}$.

The optimal solution of this well-posed optimization problem in the decision variables $\bm{u}_{s}$ is denoted by

[TABLE]

and it is now this sequence that is used to construct the fixed point iteration instead of $\bm{u}^{*}(r_{s}^{d},x_{s}(k)|\bm{v}_{s}^{(\sigma)})$ as it is suggested earlier [see (29) and the development that followed].

Therefore, rephrasing (34)-(35) with $\bm{u}_{s}^{opt}$ instead of $\bm{u}_{s}^{*}$ and $\bm{v}_{s}^{(\sigma)}$ instead of $\bm{v}_{s}^{(0)}$ leads to:

[TABLE]

where $\bm{u}_{s}^{opt}(r_{s}):=\bm{u}_{s}^{opt}(r_{s},x_{s}(k)|\bm{v}_{s}^{(\sigma)})$.

Having the new estimations $\hat{\bm{v}}_{s}^{(\sigma+1)}$, the coordinator elaborates a filtered version of the coupling profiles:

[TABLE]

and therefore, a fixed point-like round of iterations can take place. The whole scheme is sketched in Figure 3.

Note that in order for the coordinator to compute the value of the central cost function for a given pair $(r_{1},r_{2})$ of local set-points, each subsystem evaluates its contribution to the central cost (depending on the central set-points $r_{s}^{d}$), at the optimal profile $\bm{u}_{s}^{opt}(r_{s})$, namely:

[TABLE]

and sends it to the coordinator. Upon receiving these evaluations, the coordinator can compute the estimate the value of the central cost for the set-point $r:=(r_{1},r_{2})$:

[TABLE]

Remark 3

It is important to underline here that the estimation so obtained by the coordinator does not involve the knowledge of the states $x_{s}(k)$ nor that of the optimal sequence of control $\bm{u}^{(opt})$, these quantity are only known locally in the subsystem which sends the resulting corresponding value of $J_{s}$.

It is important to underline that the above estimation is irrelevant unless the fixed point iteration converges toward some fixed point $\bm{v}^{(\infty)}$ since the coherence constraint is not fulfilled otherwise. Said equivalently, if one has:

[TABLE]

then the quadruplet:

[TABLE]

is an admissible sub-optimal solution to the constrained central optimization problem (28). Therefore, only in this case, one can consider that $J(r|r^{d},\bm{v}^{(\infty)})$ is the true value of the central cost when the control profiles $\bm{u}_{s}^{opt}$ are applied by the subsystems.

In order to smoothly introduce the remaining part of the control strategy, the analysis of the conditions under which the above defined fixed-point iteration converges is delayed to Section 3.3. Meanwhile, the next section explains how to built the overall hierarchical feedback in case this convergence unconditionally holds.

3.2 Closing the loop

For each value of the auxiliary set-point $r=(r_{1},r_{2})$ the coordinator settles to define the fixed-point round of negotiation with the subsystems, the coordinator obtains after convergence the value of the central cost for this set-point vector $r\in\mathbb{R}^{n_{r}}$, namely

[TABLE]

which is defined by (38). Now since the above function is quadratic in $r$, the coordinator needs to evaluate its value over some grid including more than or equal to $(n_{r}+1)(n_{r}+2)/2$ different nodes in order to reconstruct the central cost as a quadratic function of $r$, namely:

[TABLE]

where $Q\in\mathbb{R}^{n_{r}\times n_{r}}$, $f\in\mathbb{R}^{n_{r}}$ and $c\in\mathbb{R}$ depends obviously on $r^{d}$ and $x(k)$ in a way that is analytically unknown to the coordinator.

This can obviously be done by taking a regular grid around $r^{d}$ such as the one defined by:

[TABLE]

for some $\Delta>0$ and some integer $m$ such that $m^{n_{r}}\geq(n_{r}+1)(n_{r}+2)/2$. The coefficients $Q$, $f$ and $c$ of the quadratic form can therefore be computed by solving the least squares problem:

[TABLE]

Once $Q$, $f$ are available, the optimal auxiliary set-point can be computed according to:

[TABLE]

This optimal set-point is therefore sent to the subsystems with an end-of-iterations flag which is interpreted by the subsystem as the definitive optimal set-point at instant $k$. Based on this information and in accordance with the receding-horizon principle, the first action in the corresponding optimal sequences, namely

[TABLE]

is applied by subsystem $S_{s}$ during the sampling period $[k,k+1]$.

Note that this implicitly assumes that the computation time for the $m^{n_{r}}$ iterations is negligible when compared to the sampling period of the control loop. Note also that the scheme should include an anticipative action in the sense that the computation of $u_{s}(k)$ should be done during the previous sampling period $[k-1,k]$ based on each subsystem estimation of its future state at instant $k$. These are standard implementation tricks that have been skipped here in order to focus on the main message.

This completely define the hierarchical control algorithm. In the next section, the analysis of the convergence of the fixed point iteration is done.

3.3 Convergence of the fixed-point iteration

In order to analyze the convergence of the fixed point iteration, we begin by establishing the dynamics that governs the successive iterates $\bm{v}^{(\sigma)}$. To do this, note that at each iteration, subsystem $S_{s}$ solves a finite horizon unconstrained quadratic optimization problem in which the only exogenous information is represented by the current state $x_{s}(k)$, the auxiliary set-point $r_{s}$ and the exogenous coupling signal $\bm{v}_{s}^{(\sigma)}$, therefore, the optimal control profile $\bm{u}^{opt}_{s}$ computed by subsystem $S_{s}$ takes the following form:

[TABLE]

with straightforward definition of the gain matrices.

On the other hand, the updated value $\hat{\bm{v}}_{1}^{(\sigma+1)}$ [resp. $\hat{\bm{v}}_{1}^{(\sigma+1)}$] is nothing but output of a linear system with initial state $x_{1}(k)$ [resp. $x_{2}(k)$], input profiles given by (47) [resp. (48)] and affected by coupling signal $\bm{v}_{1}^{(\sigma)}$ [resp. $\bm{v}_{2}^{(\sigma)}$]. This can be written as follows:

[TABLE]

withstraightforward definition of the gain matrices. Combining (47)-(50) obviously leads to:

[TABLE]

where

[TABLE]

and finally, using (51) in the updating rule (36) one obtains the dynamics that governs the fixed point iteration:

[TABLE]

This clearly indicates that the convergence of the fixed-point iteration is conditioned by the spectrum radius of the matrix $Z(\beta)$ defined above. More precisely:

[TABLE]

where $\rho(Z)$ is the spectrum radius of the matrix:

[TABLE]

The above arguments enables to state the following result:

**Fixed-point Convergence certification

**

The matrices $\{\bar{M}_{s}^{(v)}\}_{s\in\{1,2\}}$ summarize the static information the coordinator needs to have in order to be able to check and to certify the convergence (if any) of the fixed-point iteration that is in the heart of the proposed hierarchical framework. Having these matrices, the coordinator might appropriately choose the filtering parameter in order to enhance the convergence.

Remark 4

Note that the matrices $\bar{M}_{s}^{(v)}$ represent a high level condensed information that does not depend on the state of subsystem $S_{s}$, neither it depends on the number of actuators or the precise control law settings that are used in the local level.

Remark 5

The filtering law (36) is obviously over simplified. More elaborated filters can be used in case the simple rule does not allow convergence for any possible value of $\beta\in(0,1)$. This is not elaborated here for the sake of simplicity and since the simple rule is sufficient for the process under study.

4 Simulation results

In this section, the hierarchical MPC framework proposed in the previous section is applied to address the hierarchical control problem of the cryogenic refrigerator. More precisely, the following validation topics are successively addressed.

First of all, the convergence of the fixed-point iteration is analyzed for different values of the filtering parameter $\beta$. 2. 2.

Then the by-default regulation scenario is simulated and the performance of the closed-loop system is shown. 3. 3.

Two typical scenarios of set-point changes on the outputs Ltb131 and Ttb108 are successively simulated under the proposed control framework. 4. 4.

The behavior of the fully decentralized control (assuming zero coupling signals by each local controller) is shown in order to underline the benefit from the use of the hierarchical structure.

Before addressing the issues enumerated above, some parameters that are used throughout the simulations are first described in the following section.

4.1 The simulation parameters

The sampling time for the simulation is taken equal to $\tau=5$ seconds. This is the currently used sampling time at CEA/INAC/SBT. Moreover, it is the one that has been used in all previous studies [3, 2, 4]. The prediction horizon is taken equal to $N_{p}\tau$ where $N_{p}=100$. This corresponds to roughly 8 minutes.

The local MPC controllers are based on the local cost function (33) in which the following trivial weighting matrices are used:

[TABLE]

Regarding the weighting matrices used in the central cost, the exponent $q=20$ is used in the time-varying weighting settings; the control weighting matrices were taking systematically equal to $R_{c}^{(s)}=0$, $s=1,2$ in order to focus on the regulation performance. Regarding the output regulation weighting matrices $Q_{c}^{(s)}$, three different settings are used depending on the control mode, namely:

In the by default disturbance rejection mode, the following settings is used:

[TABLE] 2. 2.

When the second mode is activated in order to steer the level of the helium bath to some desired value, the following setting is used:

[TABLE] 3. 3.

Finally when the second mode is activated in order to steer the temperature Ttb108 to some desired value, the following settings is used:

[TABLE]

Regarding the quantities involved in the identification of the quadratic cost by the coordinator (see Section 3.2), we have $n_{r}=3$ and $m=3$ and $\Delta=1$ are used.

4.2 Fixed-point convergence analysis

Recall that the convergence condition (56) for the fixed point iteration does not depends on the settings of the central cost. It rather depends on the local settings and the filtering parameter $\beta$. Figure 4 shows the evolution of the spectrum radius $\rho(Z(\beta))$ involved in the stability condition (56) as a function of the filtering coefficient $\beta$. This figure clearly shows that the choice $\beta=0.5$ is appropriate in terms of the contraction ratio.

Figure 5 shows the evolution of the maximum error between two successive iterates during the fixed-point iterations. $50$ random initial values $\bm{v}^{(0)}$ are simulated. From this figure it comes out that $400$ iterations generally lead to an error less than $10^{-5}$ which is assumed to be sufficiently small to stop the iterations. Therefore, this number of iterations is systematically used in the sequel although a context-dependent varying number can be adopted by monitoring the error between two successive iterates. The fact is that a sampling period of $5$ seconds is already largely higher than the time needed for $400$ iterations ($\approx 0.31$ sec) and hence, adopting a stopping criteria does not change the result for our application.

4.3 Closed-loop simulations

Figure 6 shows the behavior of the closed-loop process in response to a periodic heat pulse disturbance. This disturbance emulates the heat disturbance that comes from the experimental facilities using the cooling power to cool down the supra-conducting circuits. The level of heat pulses used in this scenario is considered to be the sizing level for the experimental cryogenic refrigerator. Note that the computation time needed for the coordinator to deliver the optimal auxiliary set-points $r^{opt}$ never exceed 310 milliseconds. Moreover the convergence of the fixed-point iteration is assessed through the last plots where it can be observed that it is roughly around $10^{-5}$.

Figure 7 shows the behavior of the closed-loop process when mode 2 is activated in order to steer the level of liquid in the helium bath to some desired values. Note how due to the specific settings of the weighting matrices of the central cost, the regulation of the temperature Ttb108 is softened in favor of better regulation of the level. This can be shown by examining the difference between the central desired values $r^{d}$ (in axis-like red lines on the plots) and the auxiliary set-points $r^{opt}$ (in dashed blue lines on the plots).

Figure 8 shows the behavior of the closed-loop process when mode 2 is activated in order to steer the the temperature Ttb108 of the JT valve to some desired values. Note how due to the specific settings of the weighting matrices of the central cost in this case, the regulation of the level Ltb131 in the helium bath of subsystem $S_{1}$ is softened in favor of better regulation of the JT temperature. This can be shown by examining the difference between the central desired values $r^{d}$ (in axis-like red lines on the plots) and the auxiliary set-points $r^{opt}$ (in dashed blue lines on the plots).

It is worth underlying that the number of iteration ($=400$) is probably over pessimistic and a precision of $10^{-5}$ on the fixed-point iterations error is not necessary. Simulations proving this fact are not produced here as the plots are roughly the same as the ones presented. This means that the computation time can be divided by $2$ without noticeable performance losses.

Finally, figure 9 shows the behavior of the fully decentralized controller, namely, when each local controller considers that the coupling signal are vanishing. This scenario underlines the destabilizing character of the coupling signals involved in the cryogenic refrigerator. It is worth emphasizing that in this context, the last two plots (cpu time and convergence error) are not relevant.

5 Further discussion: handling operator handover

In industrial context, it is mandatory that the control architecture allows operator handover on the controlled process. These are predefined modes where the operator informs the control system that he (she) wishes to partially short-circuit the control logic in order to directly deliver his (her) own decision.

Depending on the level of the decision the operator wishes to deliver, two handover levels can be defined:

In the first, the operator decides to deliver one (or more) auxiliary set-point, say $r_{j}$ instead of this set-point being computed by the coordinator of the hierarchical architecture. Handling this case is extremely easy as all the coordinator needs to do is to replace the unconstrained expression (45) of the optimal $r^{opt}$ by the constrained one given by:

[TABLE]

Nothing else has to be changed in the control architecture as explained in the previous sections. 2. 2.

In the second level, the operator decides to take the direct control of some actuator value, say the first control input $u_{11}$ of subsystem $S_{1}$ (see Figure 10). Assume that before this decision the coupling signals between subsystems $S_{1}$ and $S_{2}$ where given by (see Figure 10):

[TABLE]

Now when the operator decides to be directly involved, it comes obviously that the operator should be viewed as an additional new subsystem, say $S_{3}$ which acts on subsystem $S_{1}$ though the signal $u_{11}$888Another way to view the situation is to consider that the operator is embedded inside subsystem $S_{2}$ which now acts on subsystem $S_{1}$ from which the control input has been removed.. Therefore, the operator handover at this level can be handled by two changes (see Figure 10.(b)) :

1. (a)

First the subsystem $S_{1}$ looses the control input $u_{11}$, 2. (b)

$u_{11}$ is transformed into a new component of the coupling signal $v_{1}$, that is:

[TABLE]

Now this change sets obviously a new hierarchical framework that needs to be assessed and tuned using the same tools and concepts described above. In particular, it can be revealed that for some handover modes, the stability of the system might be threatened and/or the tuning of the filtering and the local weights need to be different. The contribution of this paper gives the theoretical background to rationally address this crucial industrial question.

6 Conclusion and further investigations

In this paper, a complete framework is proposed for the design of a hierarchical control structure that addresses the issues of modularity, privacy and distributed computation that arise in the industrial process control context. The framework is presented in the case where MPC control design is used in the local level. However it goes without saying that the same framework can be used in the more common case where local PID controllers are used at the local level. Indeed, when PID controllers (or any other linear controllers) are used in the local level, the subsystems can still send the predictions of their coupling signals if the exogenous signals acting on them are given by the coordinator. This enables the convergence of the fixed-point iteration to be studied and the appropriate tuning changes to be undertaken if necessary. This enables the use of the proposed scheme in coordinating the overall consequences of the network of local controllers.

Natural extension of the present work concerns constraints handling as well as experimental validation on the real process. Using the available nonlinear knowledge-based model is also under investigation.