Plantwide control of an oil production network

In this paper, we consider Real-Time Optimization (RTO) and control of an oil production system. We follow a systematic plantwide control procedure. The process consists of two gas-lift oil wells connected to a pipeline-riser system, and a separator at the topside platform. When the gas injection rates are low, the desired steady flow regime may become unstable and change to slug flow due to the casing-heading phenomenon. Therefore, a regulatory control layer is required to stabilize the desired two-phase flow regime. To this end, we propose a new control structure using two pressure measurements, one at the well-head and one at the annulus. For the optimization layer, we compare the performance of nonlinear Economic Model Predictive Control (EMPC), dynamic Feedback-RTO (FRTO) and Self-Optimizing Control (SOC). Based on dynamic simulations using the realistic OLGA simulator, we find that SOC is the most practical approach. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Modeling, estimation, control, and optimization methodologies are becoming exceedingly important in the upstream petroleum industries. Concepts from Advanced Process Control (APC) are being developed and deployed in offshore oil and gas production ( Campos et al., 2015 ). There are however numerous remaining challenges ( Foss, 2012 ).
One standard approach for control and optimization of multiinput multi-output processes is centralized model-based control ( e.g. , Nonlinear Model Predictive Control, NMPC) which simultaneously uses all the inputs and outputs of the system ( Engell, 2007 ). In theory, such a control scheme can optimally handle the dynamic interactions between different input/output pairings, and provide inputs for optimal operation of the system. However, the success of this solution depends on obtaining a good process model and the ability to update it. In addition, optimization using detailed dynamic models (with hundreds of state variables) is computationally demanding and often not suitable for real-time applications. Campos et al. (2015) says that many numerical issues need to be addressed before dynamic optimizers can be widely used in the offshore oil and gas production. Instead, fast local controllers can on steady-state models. These models are usually detailed physical models, but they may also contain empirical nonlinear equations. Such models can be described, for example, by piecewisedefined numeric functions in the optimization problem formulation ( Gunnerud and Foss, 2010 ). The fastest regulatory layer typically controls levels, flow rates, and pressures. There is also a growing interest in introducing a slower secondary advanced (supervisory) control layer (APC), for example using model predictive control (MPC), for coordination purposes and for taking into account the constraints and interactions ( Foss, 2012 ).
The plantwide control structure design method presented by Skogestad (2004) is divided into a Top-down analysis and Bottomup design ( Skogestad, 2004 , see Table 1). The Top-down analysis starts with the definition of operational objectives, then the identification of the manipulated variables and degrees of freedom (DOF), optimization and finally selecting the primary variables for control (CV 1 in Fig. 1 ). In contrast, the Bottom-up design starts with the stabilizing control layer, then supervisory controller, and ends up with providing the integration with the optimization layer. The idea is that the control layer should contribute to the optimization.
For the top-down step a key decision is the selection of economic controlled variables (CV 1 ). For the bottom-up step a key desicion is the selection of stabilizing controlled variables (CV 2 ). Often the CV 1 variables are active constraints and they are sometimes moved into the fast regulatory layer (as part of CV 2 ).
Many papers on oil and gas production optimization focus only on the optimization layer, assuming a perfect regulatory control  ( Krishnamoorthy et al., 2016 ). However, gas-lifted wells and multiphase risers may become dynamically unstable, in particular at the optimum economic point ( Di Meglio et al., 2012 ). The usual remedy to these instabilities is to use more lift-gas, which may be costly or not available. Moreover, injecting more gas increases the friction pressure loss which reduces the production rate. Another passive solution is to choke a downstream valve, which may lead to sub-optimal operation due to an increased back-pressure. Stabilizing the desired non-slug flow regime by feedback control is shown as an optimal solution to gas-lifted well instabilities ( Dalsmo et al., 2002 ). Both robustness and tracking performance are necessary for the stabilizing regulatory control layer. The regulatory control layer must remain stable despite the variations in the process gain due to change of the operating point.
The primary objective of this study is to apply the plantwide control design procedure on the gas-lift production system. We compare three online optimizing control strategies, namely, 1) Economic NMPC, 2) Self-Optimizing Control, and 3) Direct input adaptation using a new model-based feedback-RTO ( Krishnamoorthy et al., 2019 ). Moreover, we study the effect of unmeasured disturbances in reservoir pressure and gas-oil ratio. We use the Olga simulator as the "real" plant and use a simplified dynamic ODE model for state estimation and optimization. This approach allows us to study the effect of modeling errors.
The article is organized as follows. In Section 2 , we introduce the oil production network. The top-down analysis is explained in Section 3 . The bottom-up design procedure is presented in Section 4 where the design and the robustness properties of the regulatory control layer are presented. In Section 5 , we consider the optimization layer including the Economic NMPC, selfoptimizing control and the dynamic Feedback-RTO. Dynamic simulations are presented in Section 6 . After discussing the findings and challenges of this work in Section 7 , the concluding remarks are given in Section 8 .

Oil production network
The oil production network is simulated using the Olga simulator and is represented in Fig. 2 . Olga is the industry-standard tool for simulating dynamic multiphase flow ( Schlumberger, 2018 ). In this work, the Olga model consists of two wells operated by gaslift. The oil wells feed a common pipeline-riser connected to a top- We consider fourteen measurements: • Pressure measurements (7) at the riser top, at the inlet of the pipeline, at the wellheads, at the top of each annulus and in the separator • Mass flow rate measurements (6) gas and liquid rates at the two wellheads and at the riser top • Level measurement (1) in the separator The two vertical wells are assumed to be geometrically identical with a tubing and annulus length of 2048 m. The inner diameter of the tubing is 0.124 m, and the annulus is modeled by a cylindrical (not annular) pipe with 0.2 m diameter. The roughness of the pipes is set to 4.5E-5 m. The reservoir temperature is assumed to be 108 • C. The well inflow relation is assumed to be linear ( i.e. , w res = PI (P res − P bh ) ) with a Productivity Index (PI) of 0.247 kg/s/bar. Nominally, the feed from the reservoir is oil for the present case study, that is, the produced gas-oil ratio (GOR) and water-cut are assumed to be negligible at nominal conditions. The two reservoir pressures are considered to be different; the nominal values are 160 bar for well A and 170 bar for well B.
The pipeline length is 4300 m, where the last 2300 m has a negative inclination of 1 • . The pipeline goes to a riser with height 300 m. The pipeline and riser have a diameter of 0.2 m and roughness of 2.8E-5 m. The separator operates at a constant pressure of 5 bar.
In Olga, the fluid properties can be specified by a black-oil model or as PVT Tables. We use PVT tables generated by PVTSim©. For a given feed composition, these tables contain all fluid properties such as viscosity, gas density, oil density, and gas mass fraction as functions of pressure and temperature. The viscosity of the oil considered in this paper ranges from 0.2 to 1 cP, which is not sufficient to classify it as a heavy oil. In Fig. 3 , we show the oil density and gas mass fraction as a function of pressure at 88 • C. The produced fluid is at 'undersaturated condition' and does not have free gas in the range of bottom-hole pressures considered. With these fluid conditions and the low reservoir pressure, the wells consid-  Aamo et al., 2005 ). ered in this work are not naturally flowing . Therefore, gas-lift is required to assist the production.
The two gas-lifted oil wells operate in the casing-heading instability region when the injected gas rates are low. As shown in Fig. 4 , the theoretically optimal steady-state operating point is in this region ( Aamo et al., 2005 ). The casing-heading instability causes slugging flow regimes which are characterized by large flow and pressure oscillations in the form of a limit cycle behavior. The slugging causes safety problems and inadequate separation of oil and gas. Moreover, as seen from Fig. 4 , the average production under slugging is less than the production with non-slug flow. To avoid slugging, stabilizing control is required. The stabilization is performed by low-level controllers in the regulatory layer.

Models
Three different models of the oil production network are used in this work as described next.

Reference model in the Olga simulator
The Olga model of the oil production network is used as the "real" process in the subsequent dynamic simulations. When the development of Olga was initiated in the1980's at NTNU and Sintef, one of the motivations was to study the dynamics of slow flow regimes for offshore oilfields ( Bendiksen et al., 1991 ). The simulator includes empirical correlations which have been fine-tuned by data from the oil-fields in the North Sea.
However, for the user, the simulator is a black box and the user cannot see the model equations. This means it is difficult to use the commercially available Olga simulator for production optimization.

Spline surrogate model
With the Olga simulator, it is possible to find the global optimal operating point and optimal cost. To this end, we generated data points from simulations using Olga and constructed a surrogate spline model Grimstad et al., 2016 ). The surrogate spline model consists of piecewise third-order or fourth-order polynomials. These models have the necessary accuracy and smoothness (two times differentiable) for global optimization solvers, and are then used to find the optimal operating points for different disturbance scenarios. The optimal cost values are used as a benchmark to compare the different optimizing control methodologies.

Simplified dynamic model
A simplified nonlinear dynamic model of the process is used for the state and parameter estimation and for solving the dynamic optimization problem (NMPC). The state equations of the model are the mass balances in different parts of the network. The gas-lift well model is similar to the one used by Jahanshahi et al. (2012) with some modifications, and the pipelineriser model is as presented by Jahanshahi and Skogestad (2014) .
The dynamic model can be represented as a set of ordinary differential equations: where x represents the dynamic states, i.e. , the mass of liquid and gas in the pipeline segments, u represents the manipulated variables (MV), i.e. , the valve openings and the mass injection of liftgas, and d represent the disturbances, e.g. accounting for variations in the reservoir pressure and GOR. The outputs y includes variables of particular interests, such as pressure measurements and constrained outputs.
We fine-tuned the simplified dynamic models to match the Olga simulator. The valve coefficients were used as the tuning parameters for this purpose. The model tuning purpose is to match both the dynamic and steady-state behavior. The dynamic behavior consists of the critical valve opening at which the slugging starts as we open the valve gradually, and the frequency of the slugging just after the critical valve opening. The steady-state behavior concerns the pressure and flow rate values of the stable (stabilized) process. The model tuning criteria and procedure are described by Jahanshahi and Skogestad (2014) :

Definition of operational objectives
The primary goal is to optimize an economic objective which in most cases corresponds to recover as much oil as possible. The objective can be seen from a long-term perspective, where the reservoir dynamics play an essential role ( Jansen et al., 2009 ), or from a short-term perspective (commonly known as daily production optimization), where the reservoir is considered to be static ( Kosmidis et al., 2005;Gunnerud and Foss, 2010;Codas et al., 2012;Grimstad et al., 2016 ). This work considers the short-term perspective and therefore disregards the reservoir dynamics. The short-term optimization of gas-lift wells is often constrained by the maximum amount of available lift-gas ( Camponogara et al., 2009 ).
In this work, we minimize the operation cost, J = injected gas cost -produced oil value. This is the negative of the economic profit. The constraints include the maximum injection rates and the minimum valve pressure drops. The latter are for controllability purposes. The constraints are controlled in the regulatory layer as explained in the bottom-up design, Section 4.2 . Ideally, the constraints should be controlled at their corresponding physical bounds, but a back-off is introduced to ensure feasibility and stable operation.
More precisely, the steady-state optimization problem is: Cost as a function of gas injection rates of two wells, generated by surrogate models.
are the two gas injection rates and three valve openings. α o , α g are oil and gas prices ($/kg). Eq. (2d) is the constraint on the total available gaslift, (2e) are the constraints on selected outputs ( e.g. valve pressure drops), and (2f) are the constraints on the five decision variables.

Manipulated variables and degrees of freedom
As mentioned, the degrees of freedom for the steady-state optimization are the three valve openings (v w A , v w B , v p ) and the two gas injections rates ( w inj, A , w inj, B ). There are actually five valves, but it is assumed that two downstream valves are used to keep constant separator level and pressure. The degrees of freedom resulting from the three valves are replaced by the valve pressure drop setpoints (CV 1s ). This transformation does not have any impact on the control performance, but it is essential for optimization purposes due to the uncertainty of the valve models.

Optimal operation
The steady-state real-time optimizer uses the surrogate spline model to represent the gas-lift system. The two injection rates are the unconstrained steady-state degrees of freedom for the optimization.
To obtain the surrogate model, we used the built-in parametric study tool in Olga, and we applied 0.01 kg/s steps in the two gas injection rates which gave 420 points for each disturbance scenario. Fig. 5 shows a plot of the cost as a function of the two gas injection rates with other parameters at their nominal values.
In addition to the nominal operating point, we also optimize the process for disturbances in the two reservoir pressures (D1, D2) and their gas-oil (mass) ratios (D3, D4). The nominal reservoir pressures are 160 and 170 bars, and the nominal gas-oil mass ratios are 0. Table 1 summarizes the optimal inputs and corresponding cost for various disturbances.

Active constraints
The primary economic controlled variables (CV 1 ) consist of the active constraints plus self-optimizing controlled variables. In general, it is always economically optimal to open the valves as much as possible. A maximum valve opening is equivalent to a minimum pressure drop of the valve. That is, the three valve pressure drop constraints on the three valves are always active and these are thus selected as primary controlled variables (CV 1 ). As a result, only the two gas injection rates are unconstrained degrees of freedom for the steady-state optimization.

Primary controlled variables CV 1 for self-Optimizing control
In this section, we apply the self optimizing control ( Skogestad, 2004 ) to select the two associated controlled variables (CV 1 ) to be kept at constant setpoint. These controlled variables are regulated by the two unconstrained DOFs (gas injection rates). The objective is to achieve an acceptable loss with constant setpoints when disturbances occur without re-optimizing the process for the disturbances. The setpoint values are chosen as the nominal optimal values. The loss relative to the ideal (reoptimized) costs for each disturbance are shown for seven candidate controlled variables (CV 1 ) in Table 2 . Here, the loss is the difference between the cost at the steady-state and the optimal cost ( L = J ss − J opt ). The loss has the same unit as the cost ($/s). The optimal cost values are given in Table 1 . The losses in Table 2 are obtained from simulations in the Olga simulator without any simplification. We also considered the open-loop case (Alternative 0) where the two gas injection rates are kept constant.
We conclude that the well-head gas flow rates (Alternative 7) are the best controlled variables (CV 1 ) for self-optimizing control. Note that the well-head gas flow rate is the sum of the injection rate, the produced gas from the reservoir, and the gas flashed from the petroleum at the lower pressure of the well-head compared to the reservoir pressure. Controlling the well-head liquid flow rates (Alternative 6) gives the worst result, even worse than the openloop case (Alternative 0). Constant GOR at the well-head also gives small losses (Alternative 1).

Production rate
The location of the throughput manipulator (TPM) is not considered in this work, because it is optimal to minimize the production cost, and the production rate is set indirectly by the optimization. The optimal production rate is affected by the prices of injected gas and produced oil, as well as the different disturbances. However, in general, the throughput manipulator should be located close to the bottleneck. For example, it should be located close to the receiving facilities if the maximum capacity of processing units (e.g., separation) is the bottleneck.

Selection for regulatory controlled variables CV 2
From a controllability point of view, the bottom-hole pressure is the preferred Controlled Variable (CV 2 ) for stabilizing control of gas-lift oil wells ( Jahanshahi et al., 2012 ). However, the bottomhole pressure is often not available, and if it is available, it has long sampling time intervals only suitable for monitoring purposes. Bottom-hole sensors can also easily fail and are generally not replaced in case of failure. Aamo et al. (2005) proposed to apply a nonlinear observer to estimate the bottom-hole pressure from a combination of top-side measurements. These topside measurements are the tubing head pressure, casing head pressure, and the multi-phase fluid density. Jahanshahi et al. (2012) have carried out a controllability analysis on unstable gas lift oil wells, and have concluded that a combination of topside measurements results in a robust stabilizing control system. Codas et al. (2016) considered combining the tubing head pressure (CV 1 ) and the casing head pressure (CV 2 ) by applying cascade control. Both CVs are located at the seabed and are relatively robust and easy to measure. Fig. 6 shows such a control structure (CS1) where the annulus pressure (CV 2 ) is controlled in the inner loop (slave) and the tubing pressure (CV 1 ) is controlled by the outer loop (master). Here, the production choke valve opening is used as the Manipulated Variable (MV) for the regulatory control, and the tubing pressure setpoints (CV 1s ) are available as a DOFs for the optimization layers.
Alternatively, the pressure drops over the two wellhead valves (CV 1 ) can be controlled by the master controller. This is reasonable because, as noted, the minimum pressure drop is an active constraint and should be selected as CV 1s . Such a control structure is shown as Control Structure 2 (CS2) in Fig. 6 . A similar control structure is applied to control the pipeline-rise subsystem (see Fig. 7 ).
The master loops of the cascade controllers are here classified as part of the supervisory control. The master loop helps for the stabilization by avoiding drift from the controller design point ( i.e. keeping the system in the linear region), and at the same time track the optimal setpoints given by the optimization layer. The tracking performance of the master control loops is important for optimal operation.

Setpoints for active constraints
In this work, we have used 2 bar pressure as the minimum pressure drop P for the wellhead valves, and 0.7 bar for the riser choke.
For optimal operation, it is desired to minimize the pressure drop in the network. Therefore, these P constraints will always be active, and should be selected as primary controlled variables (CV 1 ). The proposed control structure (CS2 in Fig. 6 ) allows for controlling the pressure drops at their constraints.

Regulatory control layer
A block diagram of the regulatory controllers and the disturbances are shown in Fig. 8 . The regulatory layer must keep the process stable during the transition along the optimal path provided by the supervisory controllers, and reject disturbances on a fast time scale. The reservoir pressure and gas-oil-ratio (GOR) are uncertain parameters which are considered as unknown disturbances in this work. Besides, the MVs used by the layers above are disturbances to the regulatory layer. For example, the gas injection rates used for as the DOFs for optimization are disturbances to the regulatory controllers.

Robustness
The main consideration for tuning the regulatory controllers is robustness. The robustness is evaluated at two operating points, at the initial conditions and at the steady-state optimal point. The Table 2 Loss compared to optimal (J opt ) for four disturbance scenarios. The well-head gas mass flow rates (Alternative 7) are selected for self-optimizing control.  gas injection rates to the wells are initially 1 kg/s, with the three valves 50% open. PI tunings, gain margin and phase margin for six pressure controllers are shown for the two operating points in Table 3 . The gain margins are larger than 2 for all cases. That is, if the process gain increases by a factor 2 or decreases to half due to nonlinearity, the system remains stable. and T = I − S as functions of frequency at the optimal operating point. Here, S is the transfer function from an output disturbance to CV 1 , and T is the transfer function from CV 1s to CV 1 . These are for the master controllers. The peaks of sensitivity transfer func- tions (M s -value) are not larger than 1.4 which indicates good robustness ( Skogestad and Grimholt, 2012 ).

Economic nonlinear model predictive control
In this section, we apply nonlinear economic model predictive control (EMPC) as defined by the following optimization problem.
subject to where Z + denotes the set of non-negative integers. The first term in (3a) sums the economic cost J (2a) over the prediction horizon { t k , t k + n p } where g (x k , u k ) = J , and the second summation term penalizes large control movements to regularize the problem and make the control action smooth. Note that this second term does not have any steady-state effect, because at steady-state the input change u k is zero. The weight matrix R is a tuning parameter which is chosen as a diagonal matrix. The continuous-time model (1a) is discretized using the direct single shooting method ( c.f. , Biegler, 2010 ). Thus, the equality constraints (3b) modeling the dynamic process are directly substituted in the objective function (3a) and inequality constraints (3c) by the NLP solver. Input constraints and rate of input constraints are imposed in (3d) and (3e) , respectively. An extended Kalman filter is also implemented to give full state-feedback for the EMPC implementation.

Steady-state gradient control (feedback RTO)
Recently, there have been several developments on using transient measurements along with steady-state RTO, such as works by Gao and Engell (2016) , Krishnamoorthy et al. (2018) , François and Bonvin (2013) and Krishnamoorthy et al. (2019) to name a few. However, the methods proposed by Gao and Engell (2016) ; Krishnamoorthy et al. (2018) ; François and Bonvin (2013) require explicitly solving a steady-state numerical optimization problem at each time step, unlike the method in Krishnamoorthy et al. (2019) , where the optimization is achieved via feedback control. Therefore in this paper, we consider the recently proposed dynamic feedback RTO approach ( Krishnamoorthy et al., 2019 ), which is based on estimating the steady-state gradient using process measurements y meas and a nonlinear model. A state estimation scheme is used to estimate the states ˆ x and the unmeasured disturbances ˆ d . In this paper, for the sake of demonstration, we use an augmented extended Kalman filter (EKF) for combined state and parameter estimation, see Simon (2006) for a detailed description.
Once the states and unmeasured disturbances are estimated to get an updated nonlinear model (1a) , the model is linearized to obtain a local linear dynamic model from the inputs u to the objective function J . This linear model is given in state-space form by the matrices A, B, C and D .
The steady-state gradients can then be estimated as follows ( Garcia and Morari, 1981;Krishnamoorthy et al., 2019 ). The process can be driven to its optimum by controlling the estimated steady-state gradient to a constant setpoint of zero using any feedback controller, for example a PI controller. The idea is illustrated in Fig. 10 . Note that the steady-state gradient is obtained from a dynamic model and not from the steady-state model as would be the conventional RTO approach ( François et al., 2012 ). The use of a dynamic model means that we can use the transient measurements, and thus avoid the wait time for the process to reach steady-state, as required for conventional RTO.

Self-optimizing control (SOC)
We have in Section 3.5 found that the well head gas rate is a very promising self-optimizing variable (see Table 2 ). Therefore, in addition to economic MPC and the feedback RTO approach, we also study the performance of self-optimizing control, by keeping the wellhead gas rate constant at its nominal optimal value.

Simulation results
In this section, we compare the three alternative optimization approaches by dynamic simulations.

Implementation and computation time
The Olga simulator is treated as the "real" process. The controllers (including the optimization) are implemented in Python using the simplified dynamic model of the process. The simplified dynamic model of the process is implemented in Modelica. The Modelica compiler generates a functional mock-up unit (FMU), which is a standard model component that can be shared with other applications. The resulting model was imported to CasADi ( Andersson et al., 2019 ) which includes efficient automatic differentiation techniques. NMPC and EKF were implemented using the CasADi verion 2.0.0. The communication between the Olga Simulator and the controllers was done by OPC Data Access where the Olga OPC Server is a built-in module of the simulator and the OPC client is coded in Python. We have published all of the models (Olga, Modelica, Spline models) and python codes for optimization and control as a Github repository ( Jahanshahi and Codas, 2020 ).
We have used the Single Shooting formulation in this work, and Ipopt ( Wächter and Biegler, 2006 ) with the linear solver 'Mumps' was used to solve the optimization problem. The embedded DAE solver IDAS from the SUNDIALS library ( Hindmarsh et al., 2005 ) was used to integrate the dynamic system over the prediction horizon and to compute exact derivatives. Using IDAS, we were able to choose the max step time of the integration to meet integration error tolerances. Because of stiffness, high nonlinearity, and inclusion of regulatory controls in the dynamic system, we used 2 s as the max step time of the IDAS solver to promote accuracy and prevent numerical errors, which lead to failure of the optimization. Around 90% of the CPU time was used for NLP function evaluations.
The full discretization (the Collocation Method, Biegler, 2010 ) is usually faster than embedding a DAE solver with sensitivities such as IDAS. However, our dynamical system is stiff, and for accuracy, it requires very refined timesteps on some intervals. While IDAS checks the accuracy and adjusts the timesteps online to meet the tolerances, the full discretization alternative requires fixing the timestep offline, i.e., before the problem is solved. Thus, we preferred to use IDAS at the expense of a higher computational burden because it is more robust as it deals with stiffness online.
The sampling interval of the measurement updates for the state estimation using EKF is 10 s, but the control interval of the EMPC is set to 20 min. Both the prediction and control horizons of EMPC are set to 12 h. The dynamic model used for the EKF and the optimization includes 21 state variables and five inputs. Therefore, for a single shooting formulation of the EMPC, there are 5 × 12 × 3 = 180 optimization variables. We tried shorter prediction/control horizons ( e.g. , 8 h), but the process did not settle on the steady-state optimal for short prediction horizons. 12h was the shortest prediction horizon with successful results.
Choosing shorter control intervals is beneficial, and it leads to a smoother response of the control system. However, the computation time of the EMPC optimization, which is dependant on the number of the optimization variables, must be less than the chosen control interval. We found that 20 min control interval with 12 h prediction/control horizon (180 optimization variables) is a good trade-off between the smooth response and the computation time. With the settings we chose, each EMPC optimization takes about 18-20 min of computation time to solve. We used an HP ZBook workstation (with 16 GB RAM and Intel Core i78850H CPU, 6 cores/12 threads running at 4.1 GHz) to run the EMPC simulations.
The computation time of the Feedback-RTO and self-optimizing control are less than the EKF sampling interval. As a result, the feedback-RTO and self-optimizing control can update the optimal control settings (injection rates and pressure setpoints) at each sampling interval (10 s). Time scale separation is necessary for the regulatory controllers to settle to the new optimal setpoints before we perform a new re-optimization ( Foss, 2012;Saputelli et al., 2006 ). One should notice that the time-scale separation concerns the closed-loop time constant of the upper control layer, not the sampling time interval of the data acquisition and the control signal. In general, a more frequent MV signal update leads to a smoother operation when the closed-loop time constant (adjusted by the controller gain) is kept constant.

Optimization at nominal conditions
Figs. 11 and 12 compare the performance of EMPC, selfoptimizing control and dynamic Feedback-RTO when taking the system from an initial operating point to the optimal steady-state with no disturbance.
We fine-tuned the simplified dynamic model so that its optimal operating point is very close to that of the Olga model. Therefore, the dynamic optimization converges to the same values as the steady-state RTO based on the surrogate model. Fig. 11 shows the two gas injection rates and Fig. 12 shows the cost and two gradients. The gradients for EMPC and self-optimizing control are shown for comparison purposes, as the gradients are used for control only by Feedback-RTO. For self-optimizing control, we show in Fig. 13 the well-head gas flow rates which are the self-optimizing CV 1 s. The optimal setpoints for the self-optimizing control are provided by the steady-state optimizer explained in Section 3.3 . As expected, the estimated gradients for the self-  optimizing control approach to zero as the gas flow rates settle to their optimal setpoints.
Figs. 14-16 show the optimal setpoints and the performance of the low-level pressure controllers for EMPC, self-optimizing control and Feedback-RTO, respectively. As discussed earlier, the Feedback-RTO updates the optimal settings more frequently, and the process response is smoother than with EMPC. The regulatory controllers for the self-optimizing control and Feedback-RTO show only one overshoot when the optimization is turned on (at t = 1 h), whereas there is one overshoot for every 20 min for EMPC ( Fig. 14 ). The overshoots are because of the inverse response of the process to the step changes of the pressure setpoints and the gas injection rates with 20-minute intervals. The process with the P outputs is non-minimum phase .
For the self-optimizing control, we have a decentralized control layer with five CVs and five MVs. The two unconstrained DOFs (gas injection rates) are used to control the self-optimizing CV 1 s (gas flow rates at the well-heads shown in Fig. 13 ). The remaining DOFs are the three CV 1s setpoints used to control the three P (CV 1 ) on their constraints ( Fig. 15 ).
The three optimizing control methods are next compared for two disturbances scenarios as follows.

Optimization in presence of disturbance in reservoir pressure
The first disturbance scenario is a 5 bar decrease in the reservoir pressure of both wells (D1, D2). The disturbances are ramped down within 10 h. We do not use step changes because they are not typical for a real process, and also because a large step change crashes the numerical simulation. Fig. 17 show the cost and loss from the ideal for the disturbance in the reservoir pressures. As expected, reduced reservoir pressures have a negative effect on the production and increases the cost. Fig. 18 shows the gas injections for the reservoir pressure disturbances. Injecting gas to less productive wells does not pay off in terms of economy. Therefore, the decrease in the reservoir pressure causes the optimizer to decrease the gas injection rate slightly. However, it depends on the specified gas price in the cost function.

Optimization in presence of disturbance in mass gas-oil ratio
The second disturbance scenario (D3, D4) is a 3% increase in the mass gas-oil ratio of the fluid coming from the reservoir. The disturbances are ramped up within 10 h. Fig. 19 show the cost and loss from the ideal for the GOR disturbances. As expected, increased gas-oil ratio decreases (improves) the cost.  Fig. 20 show the gas injections for the gas-oil ratio disturbances. We observe that the increase in the gas-oil ratio has a more significant effect on the manipulated variables than the reservoir pressure disturbance when we compare Figs. 18 to 20 . The optimizer decreases the gas injection significantly to compensate for the extra gas coming from the reservoir due to the gas-oil ration increase. From a fluid mechanics point of view, the excess gas causes an increase in the fluid velocity and the friction in the well, which leads to more pressure drop, and consequently, the larger pressure drop would decrease the production rate.

Performance of EKF for parameter estimation
Both the EMPC and the Feedback-RTO rely on the joint state and parameter estimation by EKF. The disturbances in the reservoir pressure (D1, D2) and the gas oil ratio (D3, D4) are estimated as parameters. Figs. 21 and 22 show the performances of EKF for estimating the disturbances when the well-head pressures and flow rates are used as the available measurements. The accuracy of this estimation depends on the model used by the EKF. As shown in the Figs. 21 and 22 , there are almost 10% estimation errors compared to the actual values of the Olga simulator.

Comparison of three optimizing control methods
The losses are compared in Table 4 . For self-optimizing control, we use the variable that is found in Table 2 , namely, the gas mass flow rates at wellheads. The loss values in Table 4 show that for the disturbances in the reservoir pressures (D1, D2), none of the controllers have any significant advantage over the open-loop (gas injections kept at their nominal optimal values). However, this is dependent on the regulatory layer control design. The losses with the self-optimizing control are the lowest overall, and the loss values of the Feedback-RTO are lower than those for the Economic NMPC for the disturbances in the gas-oil ratios (D3, D4). Since the EMPC and Feedback-RTO are based on the same dynamic model and the same EKF, this is a surprising result that is discussed in the next section. Regardless of the dynamics, we expected they settle to the same steady-state.

Effect of regulatory controllers on optimization
For the disturbance in the reservoir pressure (D1, D2), we do not obtain any significant benefit from the optimizing control over the open-loop situation ( i.e. , the gas injection rates are kept at their nominal optimal values). However, it depends on the control structure of the regulatory layer. In this work, we chose to control the pressure drop of the wellhead valves because their constraints are active, and in practice, they will be controlled on a constant setpoint. Hence, it becomes easier to apply the self-optimizing control. If we control the wellhead pressures (CS1) instead of the pressure drops (CS2), we will get poor results for the disturbances in the reservoir pressures. Obviously, it is not optimal to keep the wellhead pressure constant when the reservoir pressure decreases.
The Feedback-RTO gives lower losses than the EMPC ( Table 4 ) for gas-oil ratio disturbances (D3, D4), although they are designed based on the same dynamic model. There are two reasons for better performance of the Feedback-RTO, which reinforce each other. First, the Feedback-RTO is faster to compensate for the disturbances by adjusting gas injections ( Fig. 20 ), because it acts on the system every T c = 10 s compared to T c = 1200 s of the EMPC. Second, the dynamic response of the regulatory controllers for the Feedback-RTO is smoother than EMPC (see Fig. 23 ). For the EMPC, the responses of the process to the step changes in gas injection rates cause spikes in the control variables. The spikes decrease the control performance and hinder the controllers from tracking the optimal setpoints. Therefore, the limitation of the regulatory control performance affects the optimization too. To verify this argument, we increased the control interval of the EMPC from 20 min to 1 h so that the regulatory controllers have larger spikes and hence a worse performance (see Fig. 23 ). As a result, the loss values of EMPC increased by about 38%, as seen in Fig. 24 .

Importance of robust regulatory control
Dealing with an unstable plant is a challenge for long-term dynamic optimization. Explicitly considering the regulatory controllers in the model for dynamic optimization seems necessary, for accuracy and completeness. However, this approach requires proper tuning and possibly gain-scheduling for robustness. If the  system becomes unstable, the optimization problem cannot be solved because of invalid Jacobians. This is related to the fact that when the system becomes unstable, the inputs saturate that is u = 0 , and note that elements of the Jacobians matrix are approximately ( y i )/( u j ).  Excluding the regulatory dynamics from the optimization is possible only when the two control layers operate at completely different time scales, such that the regulatory dynamics do not affect the optimization. However, this is not the case for our system.

Consideration of active constraints in MPC
As explained above, reducing the control interval of EMPC increases the number of the optimization variables. We know that three of the input constraints (related to pressure drop setpoints) are active at the optimal point. We used this information for the disturbance rejection simulations and removed these three variables from the optimization problem by replacing the inequality constraints by equality constraints. This reduces the number of optimization variables from 180 to 72, saving a significant amount of computation time.

Computation time for steady-state models
To obtain the data for the surrogate spline models, we run each simulation for 20 h to reach the steady state. It takes about 15 min to finish each simulation on a laptop with a quad-core CPU running at 3.7 GHz. It took about five days to run 420 simulations and generate the surface shown in Fig. 5 .
By using surrogate models, it is possible to incorporate the commercially available simulation models ( e.g. , Olga, LedaFlow, K-Spice) into the optimization problem formulation. However, constructing the surrogate models based on simulation data involves extensive offline computations to obtain the steady-state data. Nevertheless, this approach is becoming a viable solution with the availability of faster computers and using cloud computing with high (parallel) processing power.

Conclusion
To our knowledge, this is the first publication on the complete control structure design (including both the regulatory and the optimization layer) applied to an oil production network and tested on the Olga simulator.
We compared three approaches for optimization layer. The dynamic Feedback-RTO and self-optimizing control are able to steer the operation to the optimal point more smoothly compared to economic MPC (EMPC). In case of unknown disturbances, selfoptimizing control results in the lowest loss compared to EMPC and Feedback-RTO.
We conclude that the self-optimizing control is the most practical method for our case study. However, this depends on the right choice of the controlled variables and the control structures throughout the control hierarchy. In this way, we found the gas mass flow rate at the well-head as the best controlled variable for self-optimizing control. Also, controlling the pressure drop over the valves ( i.e. supervisory layer) is necessary to control the process on the active constraints.

Future works
As shown by the simulations, the results of the EMPC significantly improves if optimization acts on the system with shorter time intervals. However, shorter control intervals of the EMPC gives a larger number of optimization variables when the control horizon is kept constant and requires more computation to solve the problem. On the contrary, the computation time must be shorter than the control intervals. We suggest implementing methods to speed up solving the optimization problem, such as implementing a warm-start strategy and using a rolling horizon. Also, parallel processing using multi-threaded programming techniques can reduce computation time significantly. For example, the EMPC and the EKF should run on two separate threads.

Declaration of Competing Interest
None.