# Model predictive control for helium compression station: simulation and experimental results

X.H. Pham, F. Bonne, Mazen Alamir, P. Bonnay, A. Attard

► To cite this version:

X.H. Pham, F. Bonne, Mazen Alamir, P. Bonnay, A. Attard. Model predictive control for helium compression station: simulation and experimental results. Cryogenic Engineering Conference (CEC) 2019, Jul 2019, Connecticut, United States. IOP Conference Series: Materials Science and Engineering, 755 (1), pp.012075, 2020, 10.1088/1757-899X/755/1/012075. cea-04811311

HAL Id: cea-04811311

https://cea.hal.science/cea-04811311v1

Submitted on 2 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Model Predictive Control for Helium Compression Station: Simulation and experimental Results

X.H.Pham, F.Bonne, M.Alamir, P.Bonnay and A.Attard

UMR-E 9004 CEA/UJF-Grenoble 1

INAC, SBT, 17 rue des Martyrs, 38054 Grenoble, France

E-mail: francois.bonne@cea

**Abstract.** This paper deals with multivariable on-line model predictive control (MPC) for helium Warm Compression Stations (WCS). During WCS operation, control algorithms must ensure that the operational constraints are respected. These constraints can be imposed by the system itself (valves open from 0 to 100%, compressor maximum electrical current and pressures, ...), or imposed by the users (valves which must remain closed or open to a minimum other than 0, pressures which should not be too low or too high). The MPC controller takes into account the constraints and set points into one optimization problem, which makes it the ideal candidate to control the WCS. The papers presents experimental results obtained on the SBT WCS, showing that the WCS is running safely while taking into account the constraints. The experimental tests shows that using MPC leads to high stability and fast disturbance rejection such as those induced by a turbine or a compressor stop, which is a key-aspect in the case of large-scale cryogenic refrigeration. The proposed control scheme can be used to achieve precise control of pressures in normal operation or to avoid reaching a stopping criteria (such as excessive pressures) under high disturbances (such as the pulsed heat load expected to take place in future fusion reactors).

## 1. Introduction

Large superconducting tokamak devices produce significant pulsed heat loads on magnets, due to huge eddy currents encountered in the magnetic system, to AC losses and to neutron flux radiations coming from the plasma. Such high pulsed loads disturb the cryogenic plant that is cooling the magnets, and make it necessary to use appropriate control strategies. The aim is to maintain the stability of the overall process subject to the variable thermal load and to satisfy operational and safety constraints (turbine operational temperature range, maximum capacity of the helium tank, compressor suction and discharge pressures, etc.).

Currently, technological solutions (such as thermal buffers, as described in [12], by-pass valves, etc.) are studied to smooth the effect of the thermal disturbance on the cryoplant and to avoid the over-dimensioning of the process. These solutions have to be combined with specific control algorithms, resulting in optimally designed closed-loop systems that can operate near their maximum capacity without the need for too conservative security margins.

With the appearance of several international projects (such as Large Hadron Collider), modeling and control of cryogenic plant have become one of the most important emerging interest. In particular, several dynamic simulators have been proposed by [6, 13, 5, 9, 15] for operator training, dimensioning and/or control design. Thanks to a better dynamic modelling of the

<figure><img src="image_2.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
underlying process, advanced control schemes have been proposed which were often dedicated to a particular key objective. For instance, scalar model predictive control (MPC) of the helium bath temperature at 1.8 K using a Joule-Thomson expansion valve has been investigated in [16]. In [14], the bath pressure control is studied, while in [7, 10], the high pressure level is monitored in order to control the bath level. In [8], a study of an optimal multivariable control of a refrigerator is conducted, considering pulsed heat loads.

In this paper, we propose a schema to control of warm compression stations (WCS), using model predictive control (MPC), a kind of model-based synthesized controller that manages both constraints and optimality by formulation [11]. This paper is organized as follows: section 2 introduces the WCS along with its modeling, while section 3 presents the proposed control schema and its elements. Section 4 is focused on synthesis of an observer for the controller, while section 5 presents the formulation of the MPC as well as experimental result. Section 6 concludes the paper and gives ideas for future work.

**2. The SBT Warm Compression Station**

The DSBT's warm compression station will be introduced and modelled. The first subsection presents the objectives and key figures of the SBT's warm compression station while section two presents how to set-up the simplified model used to generate the MPC controller.

**2.1. DSBT's warm compression station**

Figure 1 gives an illustration of the unit. The low and high pressures, respectively, $L$ and $H$ are the variables to be controlled. These pressures are disturbed by the flowrates, namely $m^L$ and $m^H$ coming from or to the cold box (not presented). The station has five actuators to achieve the control objective: two compressors ($C1$ is driven by variable frequency, $C2$ is constantly driven by 50 Hz) and three control valves. $V1$ is the bypass valve, used in the case where the cold box is not connected (to bypass the flow generated by the compressors) or if the flowrates $m^H$ and $m^L$ are lower than the minimum fluid flow handleable by the compressors. The valves $V2$ and $V3$ are respectively the charge and discharge valve, used to add and withdraw gas from the process. The SBT compression station is capable of handling a 72 g/s flow rate. It is the most energy consuming part of the refrigerator with the electrical power requirement of 330 kW.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Synoptic view of the SBT's compression station. Arrows represent the flow direction. The block between the $V2$ and $V3$ represents the capacity buffer of $10 m^3$</figcaption></figure>
## 2.2. Modelling to synthesize and simulate the MPC controller

To find the mathematical model the compression station in order to validate the controller designed in the following sections, the Simcryogenics library for Matlab/Simulink/Simscape [4] developed by DSBT of CEA is used. The dynamical and algebraic equations that are implemented in the library (especially to model the WCS) can be found in [3]. The validation model for the controller made with Simcryogenics can be found in [2].

According to [3], the evolution of pressures is expressed by equation (1), in which $m^L$ and $m^H$ respectively denote the flowrates going to and from the cold box. These two flows are considered to be the source of process disturbances. $m_{C1}$ and $m_{C2}$ respectively are the flowrates handled by the first and second compressor, while $m_{V1}$ represents the flow rate passing through the bypass valve V1. $m_{V2}$ and $m_{V3}$ depict the flows generated by the charge valve (V2) and discharge valve (V3). These five variables are the control (manipulated) variables of the process. $K^H$ and $K^L$ are inertia terms linked to the pipes volumes.

$$ \begin{bmatrix} K^H & 0 \\ 0 & K^L \end{bmatrix} \underbrace{\begin{pmatrix} \dot{H} \\ \dot{L} \end{pmatrix}}_{\dot{x}(t)} = \begin{bmatrix} 1 & -1 \\ 1 & -1 \\ -1 & 1 \\ 0 & 1 \\ -1 & 0 \end{bmatrix}^T \underbrace{\begin{pmatrix} m_{C1} \\ m_{C2} \\ m_{V1} \\ m_{V2} \\ m_{V3} \end{pmatrix}}_{u(t-\tau)} + \underbrace{\begin{bmatrix} -1 & 0 \\ 0 & +1 \end{bmatrix}}_{w(t)} \begin{pmatrix} m^H \\ m^L \end{pmatrix} \quad (1) $$

The state space model described by equation (1) is non-linear since the flows $m_{V1}$, $m_{V2}$, $m_{V3}$, $m_{C1}$ and $m_{C2}$ are non-linear in pressure $H$ and $L$. The control effort $u$ is also delayed by a time constant $\tau$ due to numerous control systems. To synthesize a linear MPC controller, a discrete-time linear time invariant model of the following form is to be used:

$$ z(k+1) = Gz(k) + Hw(k) \quad (2) $$

To get such model, equation (1) will be linearized and time discretized around an operating point of interest, namely: $u = u_0$, $w = w_0$, $x = x_0$ where $u_0$, $w_0$ and $x_0$ are chosen for $x$ to be in steady state condition, i.e. $||\dot{x}|| = 0$. It leads to the following expression:

$$ \tilde{x}(k+1) = A\tilde{x}(k) + B\tilde{u}(k-1) + F\tilde{w}(k) \quad (3) $$

where $\tilde{x}$, $\tilde{u}$ and $\tilde{w}$ depict the deviations from the linearization point of the variables $x$, $u$ and $w$. The sampling period has been chosen to be exactly the observed delay on actuators $\tau$. Hence, the following state space representation is constructed in order to be in the equation (2) form to be used for the stage after:

$$ \underbrace{\begin{pmatrix} \tilde{x}(k+1) \\ \tilde{u}^d(k+1) \end{pmatrix}}_{\tilde{x}_o(k+1)} = \underbrace{\begin{bmatrix} A & I \\ 0 & 0 \end{bmatrix}}_{A_o} \underbrace{\begin{pmatrix} \tilde{x}(k) \\ \tilde{u}^d(k) \end{pmatrix}}_{\tilde{x}_o(k)} + \underbrace{\begin{bmatrix} 0 \\ B \end{bmatrix}}_{B_o} \tilde{u}(k) + \underbrace{\begin{bmatrix} F \\ 0 \end{bmatrix}}_{F_o} \tilde{w}(k) \quad (4) $$

where $\tilde{u}^d$ depicts the control effort to be applied at the next sampling period. Since $\tilde{u}^d(k+1) = B\tilde{u}(k)$, we have $\tilde{u}^d(k) = B\tilde{u}(k-1)$ which corresponds to the delayed control effort. The model is now of the standard form proposed by equation (2). This state representation will be used throughout this paper.
### 3. Proposed control scheme

Essentially, the control process is illustrated schematically in figure 2. The pressures are the controlled outputs of the system. Since the disturbances w are non-measured variables, an observer is essential to estimate its values which are useful for prediction stage of MPC and $u^{ref}$ calculation (see section 4). In order for the control effort $u_{flow}$ given by MPC to be applied to the system, it is converted into appropriate units for the actuators (Hz for compressors and % for valves). The MPC is capable of taking into account real time constraints by default or given by an operator. These constraints are defined in physical unit and then are converted into to flow rate unit $g/s$ for the MPC. To summarize, the control architecture contains :

(i) a MPC controller

(ii) an observer (state, disturbance + actuator error)

(iii) a converter from flow rates to actuators values

(iv) a converter of actuators values to flow rates

(v) a reference calculation for the control input

<figure><img src="image_4.png" /><figcaption>Figure 2. Complete architecture of MPC control of the WCS pressures.</figcaption></figure>

### 4. Disturbance Observer

One may notice that the steady state of the system can be obtained by $u = B_o^{-1}F_o w$ deduced from equation (4). It makes it necessary to estimate the disturbances values to find a reference $u^{ref}$ for the control action $u$ that stabilizes the pressures. Since it has five variables but two equations, the minimal norm is used to solve for $u^{ref}$ under positivity constraints. The values of the disturbances are also used for the prediction stage of the MPC. These two reasons reinforce the need for the disturbance observer which is presented in this section

#### 4.1. Synthesis

To synthesize the observer gain, the model is extended by including the unmeasured disturbances $m^L$ and $m^H$ as components of the extended state vector. Hence, the new system is defined as;

$$x_e(k+1) = \underbrace{\begin{bmatrix} A_o & F \\ 0 & I \end{bmatrix}}_{A_e} \underbrace{\begin{pmatrix} \tilde{x}_o(k) \\ \tilde{w}(k) \end{pmatrix}}_{x_e(k)} + \underbrace{\begin{bmatrix} B_o \\ 0 \end{bmatrix}}_{B_e} \tilde{u}(k), \quad y_e(k) = \underbrace{\begin{bmatrix} I & 0 \end{bmatrix}}_{C_e} x_e \quad (5)$$

The system is considered observable since the observability matrix of the pair $(A_e, C_e)$ is full rank. Hence, the gain of a standard optimal Kalman state observer can be calculated by solving
the standard Ricatti equation in infinite time. The linear correction gain $L_e$ is used on the extended linear model in equation (5) to finally construct the observer expressed as follows :

$$ \hat{x}_e(k + 1) = (A_e - L_e C_e)\hat{x}_e(k) + B_e \tilde{u}(k) + L_e y_e(k) \quad (6) $$

Hence, the estimated disturbance is given by:

$$ \hat{w}(k) = [0 \ I] \hat{x}_e(k) \quad (7) $$

## 4.2. Simulation results

Several scenarios have been conducted to validate the observer. The weighting matrices Q and R of the Matlab command has been chosen to be diagonal matrices. These weights are given in table 1.

<table><caption>Table 1: Assigned weights to the measured variables $y$ and to the extended state $x_e$ used by the observer respectively corresponding to the diagonal of the R and Q matrices used in the Matlab command `dlqr(·)`. It can be noticed that $m^H$ and $m^L$ have the highest weights since it is the variables to be estimated.</caption><thead><tr><th>Var.</th><th>Physical meaning</th><th>Weight</th></tr></thead><tbody><tr><td>y<sup>1</sup></td><td>High pressure <i>H</i></td><td>250</td></tr><tr><td>y<sup>2</sup></td><td>Low pressure <i>L</i></td><td>200</td></tr></tbody></table><table><thead><tr><th>Var.</th><th>Physical meaning</th><th>Weight</th></tr></thead><tbody><tr><td>ˆx<sup>1</sup><sub>e</sub></td><td>High pressure <i>H</i></td><td>0.01</td></tr><tr><td>ˆx<sup>2</sup><sub>e</sub></td><td>Low pressure <i>L</i></td><td>1</td></tr><tr><td>ˆx<sup>3</sup><sub>e</sub></td><td><i>u</i><sup>d1</sup></td><td>1</td></tr><tr><td>ˆx<sup>4</sup><sub>e</sub></td><td><i>u</i><sup>d2</sup></td><td>1</td></tr><tr><td>ˆx<sup>5</sup><sub>e</sub></td><td>Disturbance <i>m</i><sup>H</sup></td><td>20</td></tr><tr><td>ˆx<sup>6</sup><sub>e</sub></td><td>Disturbance <i>m</i><sup>L</sup></td><td>100</td></tr></tbody></table>

Since the system is not stable, the observer should be tested in closed-loop as depicted in figure. 3. Figure 4 and figure 5 show the result of 2 scenarios corresponding to the test with the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Synoptic view of the observer connection with the system to control and the reference calculator. The observer output is the estimation of extended state.</figcaption></figure>

linear and the nonlinear model, respectively. From the results it is clear that the performance of the observer is validated. In the scenario two, the estimation quality falls short of that achieved in the scenario one with the linear model. However, the pressures are stabilized by the $u^{ref}$ calculated from $w$ estimated. It is shown that this estimation error is due to non-linearity and uncertainties but the overall stability is obtained. It only remains static error that needs to be compensated by the MPC presented in the next section.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Observation results with a linear model</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Observation results with the nonlinear model</figcaption></figure>

## 5. MPC design

In the previous section, the observer is synthesized to estimate the disturbance used to calculate $u^{ref}$. This control action shows the capability to stabilize the pressures but not to pilot it to the reference. This section presents a multi-variable constrained controller which allows control of the pressures by taking into account pre-defined constraints.

### 5.1. Model predictive control

According to MPC methodology, an optimal sequence of control moves is computed in order to minimize the cost function, leading to the following optimization problem:

$$ \min \quad J = \sum_{i=k}^{k+Np-1} |x(i+1) - x^{ref}(i+1)|_Q^2 + |u(i) - u^{ref}(i)|_R^2 \quad (8) $$

s.t

$$ \begin{array}{ll} u_{min} \le u(k + i) \le u_{max}, & i = 0, 1, \dots, Np - 1 \\ \Delta u_{min} \le \Delta u(k + i) \le \Delta u_{max} & i = 0, 1, \dots, Np - 1 \\ y_{min} \le y(k + i) \le y(max) & i = 0, 1, \dots, Np \end{array} $$

where $u(k + i)$ and $x(k + i)$ refers respectively to the i-step ahead control action and predicted state. $u_{min}$, $u_{max}$, $\Delta u_{min}$ and $\Delta u_{max}$ are respectively the lower and upper constrain of control action and of its increment. $Np$ is the length of the prediction horizon. From equation (4), the i-step ahead prediction of the state vector, can be easily obtained, as described below:

$$ x_{k+i} = A^i x_k + \sum_{m=0}^{i-1} A^m Bu(k-m+i-1) + A^m Fw(k-m+i-1) \quad (9) $$

Similarly, the prediction of output vector can be expressed in terms of control sequence $u_{k+i}$ and initial state $x_k$.

$$ y(k + i) = CA^i x(k) + \sum_{m=1}^{i} (CA^{m-1}Bu(k - m + i) + CA^{m-1}Fw(k - m + i)) \quad (10) $$
with (i = 0... Np). Hence, the optimization problem can be expressed in terms of initial state x(k), current disturbance w(k) and the control sequence p = [u(k), u(k+1), ..., u(k+Np-1)]T into the standard quadratic form that a solver can deal with:

$$ J(x_k, w_k, p) = p^T H p + 2p^T f(x_k, w_k) \\ \text{s.t.} \quad \underline{p} \le p \le \bar{p} \\ \Gamma p - \gamma(x_k, w_k) \le 0 \quad (11) $$

This type of problem can be solved by a quadratic problem solver to find the optimal control sequence p* that minimize J. The problem has been solved with the algorithm inspired from [1], which will be detailed in another paper. After solving the optimization problem, only the first element of the control sequence is actually applied to the system, and the new optimal control sequence is calculated at the next sampling period.

## 5.2. Experimental result

The control scheme is tested on the Warm Compression Station of the DSBT in CEA. Figure 6 shows the scenario conducted. The scenario is described as below:

* $t_1$: the set point of H is changed 16 bar from to 14 bar.

* $t_2$: offset of V1 is changed from 0% to 20% to simulate disturbances.

* $t_3$: offset of V1 is reset to 0%.

* $t_4$: valve V3 is blocked.

* $t_5$: offset of valve V2 is set to 20% to simulate disturbances.

* $t_6$: offset of V2 is reset to 0%.

* $t_7$: the V3 is released.

<figure><img src="image_8.png" /><figcaption>Figure 6. Pressures, actuators and estimated disturbances : experimental results</figcaption></figure>

It is possible to see on figure 6 that the system is well controlled by the MPC controller associated to a disturbance observer. The disturbances (here simulated by valves offset) is rejected rapidly. The set point can be changed and it is followed by the process. The actuators can be taken into manual mode (or have minumum and maximum applied to them). On this results, when the valve V3 is blocked by the user at 0%, it can no longer hold the pressure to their set points if a disturbance is coming from the cold box low pressure flow: this is why the high pressure increases while the V2 valve disturbs the plant. Once the V3 valve is released, the plant can be controlled efficiently again.
**6. Conclusion**

This paper has proposed a new approach to control warm compression station by using an
multi-variable constrained controller. As it was expected from the simulation, the performance
is encouraging. The disturbances are rejected and the pressures are piloted to the set points.
The controller is capable of taking in to account the real time constraints (such as actuator
saturation, pressure limits,...). During normal operation, control objective is achieved by driving
the compressor frequency instead of using the bypass valve. It can lead to significant energy
saving, depending on the flow margin taken on bypass valves. Furthermore, the controller is
calculated with PLC, making the proposed control scheme available for the industry. Our next
focus will be on setting constraints on the compressor current in order to make the high pressure
going down automatically in the case of compressor current overdrive.

**7. References**

[1] F. Bonne, M. Alamir, and P. Bonnay. Experimental investigation of control updating period monitoring in industrial plc-based fast MPC: application to the constrained control of a cryogenic refrigerator. *CoRR*, abs/1406.6281, 2014.

[2] F. Bonne, M. Alamir, and P. Bonnay. Control of warm compression stations using model predictive control: Simulation and experimental results. *IOP Conference Series: Materials Science and Engineering*, 171:012135, feb 2017.

[3] F. Bonne, M. Alamir, P. Bonnay, and B. Bradu. Model based multivariable controller for large scale compression stations. design and experimental validation on the lhc 18kw cryorefrigerator. *AIP Conference Proceedings*, 1573(1):1610-1617, 2014.

[4] F. Bonne, M. Alamir, C. Hoa, P. Bonnay, M. Bon-Mardion, and L. Monteiro. A simulink library of cryogenic components to automatically generate control schemes for large cryorefrigerators. *IOP Conference Series: Materials Science and Engineering*, 101(1):012171, 2015.

[5] B. Bradu, P. Gayet, and S.-I. Niculescu. A process and control simulator for large scale cryogenic plants. *Control Engineering Practice*, 17(12):1388 – 1397, 2009. Special Section: The 2007 {IFAC} Symposium on Advances in Automotive Control.

[6] B. Bradu, Ph. Gayet, and S.I. Niculescu. Modeling, simulation and control of large scale cryogenic systems. In *Proceedings of the IFAC world congress, Séoul*, pages 13265–13270, 2008.

[7] B. Bradu, Ph. Gayet, and S.I. Niculescu. Control optimization of a lhc 18 kw cryoplant warm compression station using dynamic simulations. *AIP Conference Proceedings*, 1218(1):1619-1626, 2010.

[8] F. Clavel, M. Alamir, P. Bonnay, A. Barraud, G. Bornard, and C. Deschildre. Multivariable control architecture for a cryogenic test facility under high pulsed loads: Model derivation, control design and experimental validation. *Journal of Process Control*, 21(7):1030–1039, 2011.

[9] C. Deschildre, A. Barraud, P. Bonnay, P. Briend, A. Girard, J.M. Poncet, P. Roussel, and S.E. Sequeira. Dynamic simulation of an helium refrigerator. In *Proceedings of the Cryogenic Engineering Conference*, 2008.

[10] V. Ganni and P. Knudsen. Optimal design and operation of the helium refrigerator system using the ganni cycle. In *Proceedings of the Cryogenic Engineering Conference*, 2009.

[11] C. E. Garcia, D. M. Prett, and M. Morari. Model predictive control: Theory and practice - a survey. *Automatica*, 25(3):335–348, 1989.

[12] C. Hoa, M. Bon-Mardion, P. Bonnay, P. Charvin, J.-N. Cheynel, B. Lagier, F. Michel, L. Monteiro, J.-M. Poncet, P. Roussel, B. Rousset, and R. Vallcorba-Carbonell. Investigations of pulsed heat loads on a forced flow supercritical helium loop - part a: Experimental set up. *Cryogenics*, 52(7-9):340–348, 2012.

[13] R. Maekawa, K. Ooba, M. Nobutoki, and T. Mito. Dynamic simulation of the helium refrigerator/liquefier for lhd. *Cryogenics*, 45(3):199 – 211, 2005.

[14] R. Maekawa, S. Takami, and M. Nobutoki. Adaptation of advanced control to the helium liquefier with c-prest. In *Proceedings of the International Cryogenic Engineering Conference*, 2008.

[15] C. Regier, J. Pieper, and E. Matias. Dynamic modeling of a liquid helium cryostat at the canadian light source. *Cryogenics*, 50(2):118 – 125, 2010.

[16] E. Blanco Vinuela, J. Casas Cubillos, and C. de Prada Moraga. Linear model based predictive control of the LHC 1.8 K cryogenic loop. In *Proceedings of the 1999 Cryogenic Engineering Conference*, 1999.