Integration of Prognostics and Control of an Oil/CO2 Subsea Separation System

Lucas Ferreira Bernardino 1 , André Felipe Ferreira de Souza 1, Argimiro Resende Secchi 1,* , Maurício Bezerra de Souza Jr. 1,2 and Anne Barros 3 1 Programa de Engenharia Química, UFRJ, Rio de Janeiro 21941-450, Brazil; lbernardino@peq.coppe.ufrj.br (L.F.B.); andref@peq.coppe.ufrj.br (A.F.F.d.S.); arge@peq.coppe.ufrj.br (A.R.S.) 2 Escola de Química, UFRJ, Rio de Janeiro 21941-909, Brazil; mbsj@eq.ufrj.br 3 Department of Mechanical and Industrial Engineering, NTNU, 7491 Trondheim, Norway; anne.barros@ntnu.no * Correspondence: arge@peq.coppe.ufrj.br; Tel.: +55-21-25628307


Introduction
In the context of process systems engineering, plant automation extends control possibilities and allows for operations in places virtually inaccessible to humans. The subsea environment is a classic example of an inaccessible place in which process operation is performed, mostly related to oil and gas exploitation in deep and ultradeep waters. Another issue that arises with low accessibility is maintenance planning, which needs to be done in order to prevent process interventions and production downtime, but also needs to take into account the difficulty of performing such maintenance. For this, reliability analysis aims to quantify rate or probability of casualties, in order to aid in the decision-making process [1].
In the early industrial periods, maintenance was done only at process breakdown ("run-to-failure maintenance"). This behavior often leads to major financial losses, so this kind of maintenance has mostly given way to preventive maintenance, in which parts are repaired or replaced based on historical data. As competition in the industries became tighter, the time between repairs became a critical variable. If this time interval is too high, the chances of a process breakdown rise. If it is too low, the introduced prognostic is based on a new, original stochastic dynamic degradation model, and a particle filter, which performs online state estimation based on the model. Finally, the application focuses on a new HISEP TM process of high CO 2 percentage separation. There are no published works in the open literature of HAC of this process, as it is a technology that is being developed right now.

Methodology
In this section, the studied process and the methods related to its control and prognostics are detailed.

Process Description
The process flow diagram (PFD) for the process considered in this work is displayed in Figure 1. This PFD is based on a process conception developed by Passarelli [10] and further analyzed by Souza [12] and Souza et al. [13].  [20].
In this process, the crude oil in reservoir conditions is submitted to heat exchange with a hot stream from the process, which provides part of the heat necessary for separation. This stream is then sent to a separator drum, in which the remaining heat necessary for separation is provided. The CO 2 rich stream, withdrawn at the top outlet, is cooled, elevating its density so that it becomes adequate for pump operation. The pump elevates the fluid pressure, allowing its reinjection, but it also significantly raises the fluid temperature. This temperature elevation enables heat integration between this stream and the crude oil extracted from the reservoir. After heat integration, the CO 2 rich stream is sent to the reinjection section. The separation drum bottom stream is then sent to topside processing. In addition, to stabilize the process, a proportional integral (PI) controller is implemented to control the flash drum liquid content by manipulating the bottom valve opening.
The dynamic model for the system studied in this work and the process variables are detailed in Appendix A. Simulations were carried out in Python 3.7 [22]. The CO 2 Therm package developed by Souza [12] was interfaced to Python using the Boost:Python library [21]. Differential-algebraic equation system integration was performed using the Implicit Differential-Algebraic (IDA) algorithm from the Assimulo package [23].

Setpoint Tracking Control System
A discrete linear state-space dynamic model, with a sampling time of 10 s, was identified. The interested reader may find all the details of the linear identification in [20]. The identified model was used as the internal model for the model predictive control (MPC) framework, in which the controlled variables were selected as y = [F in , P f lash , T f lash , T 2 ], and the manipulated variables were selected as u = [x v,t , x v , Q f lash , Q cool , W]. This choice of controlled and manipulated variables was made according to a previous work [13], in which unconstrained MPC was compared with feedback control.
The implemented MPC framework comprises a discrete Kalman filter (KF), which was used to perform state estimation and model correction, and an optimizer used to calculate the optimal control actions sequence.
For the KF implementation, the temporal and measurement update equations were implemented as in Equations (1) and (2), respectively [24].
The dynamic optimization problem for the control actions calculations was formulated as in Equation (3). This optimization problem was solved using the scipy.optimize.minimize routine, using the sequential least squares programming (SLSQP) method (tolerance f tol = 10 −6 ) [25]. Tuning parameters for the control system were selected based on the literature for the MPC [26] and simulation tests. The used parameters are presented in Table 1.
subject to: u min k+j ≤ u k+j ≤ u max k+j , j = 0, · · · , N c − 1 −∆u max k+j ≤ ∆u k+j ≤ ∆u max k+j , j = 0, · · · , N c − 1 x k+j = Fx k+j−1 + Gu k+j−1 + w k , j = 1, · · · , N p y k+j = Cx k+j + Du k+j , j = 1, · · · , N p (3) With this setup, closed-loop simulations were performed to evaluate the performance of the controller with the proposed tuning. All dynamic simulations had the nominal steady state, described in Appendix A, as the initial condition, and at t = 0, a change in the flash drum pressure setpoint of +0.5 MPa was implemented.

Prognostics Module
The implemented prognostics module comprises a stochastic dynamic model, which represents the pump degradation process as a function of its operating point, and a sequential importance resampling (SIR) particle filter, which performs online state estimation based on the degradation model [27].
The degradation state evolution (λ) was modeled inspired by Paris' law for crack propagation [28] and is given by Equation (4), in which the p.d.f. corresponds to a Gamma(k, θ) distribution as defined by Equation (5), and W k represents the pump operating power at time k.
The considered observation for the degradation system is the degradation state, corrupted by white noise, according to Equation (6). Table 2 presents the parameters for the degradation model. For this set of parameters, the operational degradation limit was fixed as λ lim = 1.0.  The particle filter was implemented with the same model and parameter set as the considered degradation process, except from the standard deviation, which was considered to be σ = 0.001. The algorithm was based on Speekenbrink [27], and is further described in Appendix B. The number of particles was N part = 100, the resampling rate was c = 0.5, and particles at time k = 0 were sampled from N (0.001, 10 −10 ). All statistical operations were performed using the scipy.stats algorithms.
Simulations of the degradation process were performed considering a planned horizon for the manipulation of the pump power W. Using the filter results, a prediction of RUL was made at each measurement. To evaluate prediction performance, this was compared to the actual end-of-life time in the simulation.

Health-Aware Control System
By combining the setpoint tracking control system and the prognostics module, an HAC structure was implemented. It consists of an MPC, with the objective function written as in Equation (7), with w H AC being the weighting factor between RUL extension and control objectives, and RUL being the average of RUL between the propagated particles. Tuning parameters were kept as in Table 1.
The model proposed in Equation (4), however, cannot be directly used in a deterministic optimization problem, due to its randomness. Furthermore, as the degradation model is represented in discrete time, RUL has no direct correlation in continuous space. To address these issues, the strategy used to calculate RUL during the optimization procedure is as given in Algorithm 1.

Algorithm 1: RUL calculation
k represents the quantile associated with the k-th evolution of the i-th particle, W end represents pump power at the end of the control horizon, and Q(·) represents the inverse cumulative distribution function associated with the distribution in Equation (4), depending on the pump power which varies along the optimization. The set of µ (i) k is sampled from a standard uniform distribution, and kept constant in the optimization course. As a result of the difference between timescales, W end was considered as the dominant effect in the degradation pattern, and the variation during the control horizon can be neglected.
For simplicity, the MPC internal model was considered to be representative of the process, being used as the process model. The used MPC tuning was the same as described in Table 1, except for the control horizon, which was N c = 2. For higher control horizons, the optimization problem is ill-posed.
Closed-loop simulations were performed to evaluate the influence of the number of particles (N part ) and the value of w H AC , in Equation (7), over the controller performance.

Setpoint Tracking Control System
In order to attain the main objective, which is to develop an HAC tool, first it is necessary to evaluate the system control separately. Guaranteeing that the controller is well tuned is fundamental before adding complexity to the control framework. Furthermore, these results serve as a benchmark for the proposed HAC framework. Figure 2 presents the closed-loop simulation results for a setpoint change of the flash drum pressure. System settling occurred at around 600 s, which justifies the chosen prediction horizon of 750 s. The implemented controller produced no offset, even though there is a mismatch between the controller internal model and the plant. This result is due to the KF-based model update translated in the variable w. For comparison, the same simulation without the model update, i.e., w k = 0 ∀k, is also displayed, in which an offset can be observed. This highlights the importance of feedback, not only in state estimate, but in the process internal model in order to attain a good performance in the control strategies. The adopted model update, although simple, was enough to represent the process nonlinearities that come with the steady state change.

Prognostics Module
Using a representative model of the degradation process, online state estimation was performed using a particle filter. Despite the fact that a simple tool such as the Kalman filter was used for the process control problem, the prognostics problem requires a more sophisticated tool. This is due to the absence of a representative linear model of the degradation process, and the importance of keeping the statistical information of the model.
In the first set of simulations, a constant pump power corresponding to its nominal value (W(t) = 4.246 MW) was considered. Figures 3 and 4 present the results of the SIR filter implementation for the a priori and a posteriori distributions, respectively. It can be seen that the distributions remain well conditioned throughout the simulation due to the resampling step. Furthermore, as expected, the measurement update effectively narrows the sample distribution a posteriori, as the a priori information is combined with the measurement likelihood function, leading to smaller confidence regions for the degradation state. True state SIR estimate Figure 3. Degradation state online estimation-a priori distribution (dot size represents particle weight).
All parameters from the stochastic model were implemented in the particle filter as their actual values, meaning that no influence of modeling error was considered in this implementation. Only the likelihood function standard deviation was changed to σ = 0.001, solely due to numerical issues related to the very narrow actual likelihood function. This led to a more permissive filter, which means that a posteriori distributions are wider than necessary to describe the state.
Using the particle filter framework, RUL predictions can be performed. However, as RUL prediction generally involves a high number of state transitions, leading to high-dimensional statistical problems, it is known that a distribution reconstruction using importance sampling is nearly infeasible. The chosen strategy, then, was to retrieve information about the distributions by standard Monte Carlo sampling. As a discrete-time state transition model was used, predicted RUL belongs to a discrete set. Thus, the very number of occurrences of each realization was considered as representative of the probability. Results of RUL prediction for the SIR filter are presented in Figure 5, in which predicted RUL is displayed in terms of the number of discrete model evolution steps. It can be seen that SIR filter RUL distribution narrows around the true RUL as time passes. This is due to the resampling step, which redistributes particles in the most likely state values.  A simulation with time-varying pump power was also performed (W(t) = 10.615 × (5 − t)/5 MW, t in years). Figures 6 and 7 show results for the sequential importance resampling (SIR) filter in the cases of a priori and a posteriori estimates, respectively. It can be seen that, for the a priori estimates, the distribution starts to degenerate when lower pump power is applied. This is due to the lack of degradation in this region. Nevertheless, a posteriori weights are normalized and this effect is counteracted.   Regarding the RUL prediction for this case, results for the SIR filter are given in Figure 8. Even though RUL prediction in early states is biased, as new measurements are incorporated, the predictions progressively become more accurate. This result encourages the use of particle filters as auxiliary tools in RUL prediction, even in the presence of modeling errors.

Health-Aware Control System
The HAC objective function described in Equation (7) aims to combine the health extension and control objectives in the same optimization problem. As the control problem has four independent controlled variables and five independent manipulated variables, there are enough degrees of freedom to find an operating point which satisfies the control objectives with the maximum equipment RUL. It is evident from the models that one should aim for the maximum top valve opening (x v,t ) to minimize pump effort and thus extend its lifetime.
The influence of w H AC over the HAC strategy is explored in the simulations presented in Figures 9 and 10, with lower and higher values of w H AC , respectively. As a result of the randomness of sampling µ (i) k at each control step, the RUL estimate (RUL) fluctuates at each control action calculation, which changes the optimization problem at each time step and thus makes the attainment of a steady state impossible. As expected, the variable x v,t is brought to its allowed maximum, and W is minimized to some extent in all cases. In the case with lower w H AC , even though the control objectives are not ignored, an offset is produced. As the operation reaches a point where health and control objectives compete against each other, the obtained solution is a compromise between these objectives. In the case with higher w H AC , the optimization problem prioritized the RUL extension objective over the control objective to the point where the system diverges from the reference trajectory.
When stochastic problems are solved by sampling strategies such as particle filters, the number of samples is a key tuning parameter to attain an acceptable performance. As such, strategies to reduce the noise associated with sampling were investigated, and the results are summarized in Figure 11. By increasing the number of particles from 3 to 20, some reduction in the oscillation amplitude of RUL is observed. This means that the number of particles needs to be raised some orders of magnitude to attain a representative mean value that stays constant with the resetting of quantile samples.
As the elevation of the number of particles proved to be not very effective in stopping oscillatory behavior, a different strategy was adopted. Instead of resetting the values of µ (i) k at each control step, these values were kept constant throughout the closed-loop simulation. For this approach, the attainment of a noise-free steady-state value for the manipulated variables is evidenced, at the expense of a biased RUL.     For comparison purposes, the results presented in Figure 2, which correspond to an HAC controller with w H AC = 0, are complemented with the corresponding RUL presented in Figure 11. The difference between the steady RUL of HAC and MPC highlights the RUL extension capability of the HAC strategy. In all of the presented HAC simulations, it was not possible to eliminate the offset when the health extension objective conflicts with the reference tracking objective. This suggests that a different approach is necessary to attain both objectives.

Conclusions
In this work, motivated by the recently developed process of subsea CO 2 separation, process modeling, control strategies, and equipment prognostics were assessed.
The control strategies employed in this work were the PI controller, whose main role was to stabilize plant dynamic behavior, in parallel with an MPC controller, with a linear internal model identified from the PI-stabilized plant, and a KF, used to correct state estimates and model bias. This strategy was successful, highlighting the importance of filter-based model correction in cases of model-plant mismatch. A pump wear model was proposed, which reunited most of the desired characteristics in an equipment wear model (monotonicity, time independence, dependence with important operational variables). Using the stochastic processes mathematical formulation, particle filters were successfully used to estimate states and predict the remaining useful lifetime. The development of this kind of model is challenging mainly because of the scarcity of quality data regarding all of the relevant variables. This development was however very encouraging because of the positive impact it can have on process reliability and optimization.
The main achievements of this work were then consolidated in the deployment of an HAC tool to solve the proposed case study. For the case with w H AC = 10 (Figure 12), the tool was able to extend the lifetime of the process by roughly 4 months, compared to the same operation with a conventional MPC (Figure 11) when the process is still under healthy conditions. Even though consistent results were obtained, some limitations of the proposed method were noted, most of them concerning the contraposition of health and control objectives.
The unstable behavior seen with the high RUL extension priority is a known issue which has been explored in the development of economic MPCs. For this class of control problem, stability can be guaranteed under some assumptions and modifications of the optimization problem [29].
As any optimization problem with multiple conflicting objectives, the most accurate way to express it mathematically is as a multiobjective optimization problem. The objective function formulated in Eq. (7) can be seen as simply the multiobjective problem solved using the weighted sum approach. Therefore, the treatment of this generic problem using other tools, such as the goal attainment method [30], can result in more reasonable control policies. Although this problem was solved in an one-layer control framework, the complexity of the methods necessary to describe the degradation phenomenon suggests that this issue could in future investigative works be addressed in a higher control layer and a systematic comparison could be performed. With this, a higher number of particles could be used, resulting in not only a more reliable RUL estimation, but in an estimate of the RUL p.d.f. itself. This is expected to enable the analyst to work with the constraints related to the confidence level of the degradation process, resulting in a more robust decision-making.

Funding:
This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES)-Finance Code 001 and by the Norwegian Research Council through the INTPART programme.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Process Model
In this section, the model of the subsea CO 2 separation system considered in this work is described. It is considered that the reservoir contains CO 2 , water, methane, and heavy hydrocarbons modeled as pseudocomponents F 1 , F 2 , F 3 , and F 4 . The interested reader is encouraged to consult Bernardino [20] for a more thorough process description, and Souza [12] for the description of the thermodynamic package CO 2 Therm.

Appendix A.1. Valves
A scheme which emphasizes variables of the valve model is presented on Figure A1. For this model, it was assumed that the flow through a valve is an isenthalpic process, and that the molar flow has a linear correlation with the valve opening, as in in which A represents the valve cross-section area, x v represents the valve opening, and K v represents the valve constant.

Appendix A.2. Heat Exchanger
The heat exchanger used for process integration was discretized into theoretical stages, as depicted in Figure A2.   Figure A2. Heat exchanger discretization scheme, emphasizing model variables.
By performing energy balances at each theoretical stage, the following equations apply: with n st being the number of theoretical stages.

Appendix A.3. Flash Drum
The flash drum modeled in this work is schematized in Figure A3. The equation system that models this equipment is written as The model for this equipment comprises the component material balances and the energy balance. The thermodynamic package CO 2 Therm [12] supplies the algorithm Flash3P, that performs (P, T) specified flash calculations. However, as the drum internal energy and volume are specified, algebraic constraints are necessary to ensure consistency.
The drum level is obtained from the liquid volumetric fraction at drum conditions, given by The scheme for the submersible pump employed in this virtual plant is shown in Figure A4. The model for this equipment is composed of an energy balance (Bernoulli equation) and the assumption of isentropic flow, given respectively by

. PI Controller
The PI controller was implemented using the following formulation, in which the integral term is calculated through an auxiliary differential equation in time: The controller was tuned using rules proposed by Skogestad [31] based on the identified transfer function of the pair x v,b x h f lash , with τ c = 10 s. Appendix A. 6

. Model Parameters
The parameters and variables presented in Tables A1 and A2 were used throughout all the simulations. Using these, the obtained nominal steady state is described by the variables in Table A3.

Appendix B. Particle Filter Algorithm
The particle filter algorithm implemented in this work is given by Algorithm A1. In this algorithm, X i k and w i k represent the i-th particle and its correspondent weight at time k, y k represents system measurement at time k, ψ j represents probability distribution of j-th parameter conditioned to the state and previous parameters, f represents the state evolution law, g represents the likelihood function, N part is the number of particles, and c is the filter resampling rate. The resampling algorithm used was the systematic resampling, described in Speekenbrink [27].