Control pairings of a deoiling membrane crossflow filtration process based on theoretical and experimental results

Journal of Process Control 81 (2019) 98–111

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

Control pairings of a deoiling membrane crossflow filtration process based on theoretical and experimental results Kasper L. Jepsen, Simon Pedersen, Zhenyu Yang ∗ Department of Energy Technology, Aalborg University, Esbjerg, Denmark

a r t i c l e

i n f o

Article history: Received 19 November 2018 Received in revised form 16 May 2019 Accepted 17 May 2019 Keywords: Relative gain array Membrane Control parings Separation Multi-phase Oil and gas

a b s t r a c t Despite membrane filtration being extensively used in many industries, it is not widely used for offshore produced water treatment because of the large installation footprint required. Commonly, PID-controllers are deployed to maintain the operating conditions that minimize fouling, but the design of the controllers is rarely highlighted. The interaction between the deployed controllers can cause undesired performance and even instability. In this study, the relative gain array method was used to find the set of control pairings that minimized the interaction between the control loops. Experimental results showed that by accounting for cross couplings in the design phase of the control system, the controllers performance in terms of reference tracking can be improved significantly. By improving the reference tracking and disturbance rejection capabilities of the filtration system, the selected operating point which minimizes fouling growth, can be maintained consistently subject to disturbances in terms of varying feed flow rate. Better tracking of the operating point can reduce membrane resistance, and consequently reducing the operating cost. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction In the Danish offshore oil and gas sector large quantities of oil are extracted, but with a 75% average water cut, the enormous amounts of produced water (PW) must be treated before discharge or reinjection [1,2]. To fulfill the Danish governmental regulations for 2017/2018, the Oil-in-Water (OiW) concentration in the discharge must be below 30 mg/l and the annual total oil discharge below 222 tonnes [3,4]. To ensure compliance with potentially stricter regulations in the future, new technologies must be considered. Membrane filtration is a technology that has proved to reduce OiW concentration far below the capabilities of the typically applied deoiling hydrocyclones [5–8]. However, fouling of the membranes significantly reduces the flow capacity [9,10]. Membrane filtration is rarely integrated into the conventional offshore treatment train because of the relatively large CAPEX and OPEX, which is caused by the required installation footprint [11]. Increasing the attractiveness of membrane filtration demands a reduction in CAPEX and OPEX. This can be achieved by reducing or removing fouling by optimizing cleaning methods, membrane materials, chemical additives, and process control [12].

∗ Corresponding author. E-mail addresses: [email protected] (K.L. Jepsen), [email protected] (S. Pedersen), [email protected] (Z. Yang). https://doi.org/10.1016/j.jprocont.2019.05.009 0959-1524/© 2019 Elsevier Ltd. All rights reserved.

The operating point, in terms of temperature, permeate flow rate, and crossflow (CF) significantly affect the fouling growth [13–15]. Several studies have investigated the relationship between operating points and fouling growth, to find the optimal operating point [8,16,17]. However, little attention has been given to the dynamic controllers deployed to maintain those operating conditions. Because of the relationship between operating conditions and fouling, the design of the dynamic controllers is essential. For membrane processes, the PID control technique is the most frequently deployed control scheme, but it is often deployed as a necessity for sufficient reference tracking rather than a topic of optimization [18–23]. In general, the focus within the field of membrane filtration is often dealing with added chemicals [24–26], materials [27,28], or set-point optimization [19,29–31]. Thus, tuning of the dynamic controllers and control pairings are rarely discussed. More advanced control strategies, such as Model Predictive Control (MPC), are only deployed for membranes in a few studies: (i) MPC is used to improve the system’s ability to maintain the selected operating conditions, but the study is purely based on simulations and lacks experimental results [13]; (ii) MPC is used to schedule switching to avoid water hammering [32]. Theoretically, MPC can provide superior reference tracking, compared to PID, but its reliance on an accurate model causes increased sensitivity to model deviations, therefore model updates are required if system changes are made. Furthermore, as MIMO control strategies

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

are rarely deployed within the field, this study aims to emphasize the potential improvements by accounting for interactions without changing the control strategy. The popularity of the PID control structure, compared to more advanced methods, is likely a result of simplicity, transparency, and availability of simple tuning rules [13]. Most PID design methods do not account for cross couplings and interactions between control loops. When interactions and cross couplings are ignored in the control design phase, the controllers may individually perform well. However, the interactions between control loops can cause performance degradation and ultimately instability, consequently increasing the fouling growth. As controller interaction can degrade the reference tracking and disturbance rejection performances of the controllers, this study highlights the benefit of accounting for interactions in the design phase, despite the controllers being individually tuned to perform identically. Furthermore, the work in [4] is extended by deploying frequency dependent relative gain array (RGA) to analyze and select the control pairings that minimizes control loop interaction [33]. In Section 2, the pilot plant used for the experiments and validation is described and the model structure from [4] is modified and extended. The coefficients of the described model and dynamic features are identified and analyzed in Section 3. Section 4 compares the fouling behavior subject to oscillating and constant operating conditions, and concludes that oscillating operating conditions can potentially increase fouling. To improve control and thereby ensureing constant operating conditions, different control pairings are investigated in Section 5. Here the RGA method is used to identify the control pairings that minimizes control loop interaction, thus improving reference tracking and disturbance rejection’s performances. Based on the RGA analysis, a potential best and a worse control paring candidate are compared in Section 6. The comparison highlighted the large differences in reference tracking and disturbance rejection performances caused by different control pairings. Finally, Section 7 concludes that by accounting for the cross coupling in the design of the controllers, fouling can potentially be reduced.

99

Fig. 1. Considered membrane filtration configuration, where WP1 is a centrifugal pump, V are control valves, MFM are magnetic flow meters, and PT are pressure transmitters with the exception of PT8 which is a virtual measurement.

2. Physical and mathematical description The pilot plant and its operation are described, and a model of the process system is formulated based on the pilot plant.

Fig. 2. The model structure, where the flow rates Q and the pressure difference PWP1 are defined according to Fig. 1. U are the control reference for the valves and pump, KVf is the conductance of the fouling layer, Pˆ Vx is the position of the valve, and H is the transfer functions for the different model elements. The subscript x denotes valve number.

2.1. Pilot plant The laboratory scaled pilot plant is described in [4] and sketched in Fig. 1. During filtration operations, a part of the flow permeates the membrane as purified water and the remaining flow is a concentrated form of the feed, containing a higher OiW concentration. The concentrated flow is either recirculated as CF or rejected from the CF-loop to keep the CF OiW concentration low. The reject can either be passed to another filtration unit for further processing or returned to an upstream separation stage. For the experiments in this study, batch mixing is used, and after separation both oil and water are returned to the mixing tank to be reused and remixed, which enables continual experiments. To ensure that the OiW mixture is kept well mixed, two large mechanical stirrers (Milton Roy mixing) are deployed. The membrane filtration system consists of isolation and control valves, centrifugal pumps, pressure transmitters, temperature transmitters, flow meters, and membrane filters. The membrane filters are ceramic-based and manufactured by LiqTech, with a nominal pore size of 0.04 ␮m. For experimental purposes, the pilot plant is over-actuated, i.e. the system has more actuators and degrees of freedom than controlled outputs.

2.2. Process model To carry out the interaction analysis, a model of the pilot plant must be developed. The considered system has nonlinearities, especially the pump and valves are known to behave in a nonlinear manner. The valves are modeled as a Hammerstein Wiener model, while the pump is modeled as a linear parameter-varying (LPV) transfer function (TF), where the coefficients are described in terms of scheduling parameters, a structure proposed in [34]. Based on previous work, V2 can be excluded from the RGA analysis as the orifice of the valve is too large to dynamically control the small reject flow [4]. However, as V2 is required to create sufficient reject flow resistance for proper operation, V2 is kept in the static analysis to facilitate estimation of the steady state operating point. The complete model structure is illustrated in Fig. 2, where it is emphasized that the relationship between V2 input and position is static.

100

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

2.2.1. Statics The static part of the model is identical to the model described in [4], but will be summarized here. The flow balances at the intersections are defined as:

where PWP1 (QWP1 , UWP1 ) will be described and identified Section 3.2. Combining Eqs. (3) and (6) produced the following relationship between PWP1 and P: P =

PT 1 + PWP1 + PT 6 − PT 4 2

QWP1 = Qf + Qcon ,

(1a)

Qcon = QCF − Qrej ,

(1b)

QCF = QWP1 − Qpm ,

(1c)

yc = f (z, d),

where QWP1 is the flow rate through WP1 , Qf is the feed flow rate, Qcon is the concentrated flow rate, QCF is the volumetric flow rate of the CF, Qrej is the reject flow rate, and Qpm is the permeate flow rate, all of which are in L/s and defined as in Fig. 1. The flow rate through the different valves are described as being linearly proportional to the square root of the pressure drop over the valve:

QWP1 = fs (z, d),

QCF = Qrej =



PT 6 − PT 1 · KV V 1 (Pˆ V 1 ),



Qpm =

PT 1 − PT 2 · KV V 2 (Pˆ V 2 ),



PT 4 − PT 3 · KV V 3 (Pˆ V 3 ),

PT 5 + PT 6 , 2

(2c)

(3b)

(4b)

where KVCF is the flow conductance for the CF channel and KVm is the flow conductance through the membrane wall, which is described as two restrictions in series, caused by fouling and membrane porosity: 1 KV f

,

(5)

where KVm|c is the flow conductance for a clean membrane, and KVf is the flow conductance for the fouling layer. WP1 is modelled as a pressure source depending on pump speed and the flow rate through the pump: PWP1 (QWP1 , UWP1 ) + PT 1 = PT 5 ,

d = [KV f , Qf ]T ,

(9b) T

(9d)

Qpm = KV m · P.

+

(9a)

u = [UV 1 , UV 2 , UV 3 , UWP1 ]T .

(4a)

1

T

z = [Pˆ V 1 , Pˆ V 2 , Pˆ V 3 , PWP1 ] ,

(2b)

QCF = KV CF · (PT 5 − PT 6 ),

1 KV m|c

where z, d, yc , and u are defined as:

(9c)

The cross and permeate flow channel are less restricted and prone to flow turbulence. Thus, the flow rate is assumed to be linearly proportional to the pressure drop:

KV m =

(8)

yc = [Qpm , QCF ] ,

(3a)

P = PT 8 − PT 4 .

Eqs. (1)–(6) can be solved to provide the functions f and fs , defined as:

(2a)

where PTx is the pressure at the location indicated in Fig 1, KVVx is the valve conductance as a function of normalized valve position, Pˆ Vx , where the subscript x denotes the valve number. The functions are identified in Section 3.3. The pressure measurement PT8 does not exist on the pilot plant, as it is impractical to place a pressure transmitter inside the CF channels. Because the flow in the CF channel is gradually reduced as the permeate permeates the membrane, the reduction in P across the length of the membrane is nonlinear, on the contrary to a linear pressure gradient across the CF channel if no permeated is produced. The nonlinearity caused by the permeate flow is relatively small, especially as QCF is up to 46 times larger than the Qpm [35]. As such, Qpm is found insignificant compared to QCF when calculating the pressure drop over the CF channel thus the pressure at PT8 can be approximated mean of PT5 and PT6 , and the transmembrane pressure, P, is the pressure difference between PT8 and PT4 : PT 8 =

(7)

(6)

where z is the input to the static model, u is the controllable inputs, d is the disturbance input, yc is the output from the static function. The model functions are divided into two functions, f for the considered outputs and fs for the scheduling variable used by the LPV TF HWP1 . Commonly the membranes can be controlled by either maintaining a constant P or a constant Qpm [36]. For this study, the constant Qpm control mode, where both QCF and Qpm are controlled, is selected. 2.2.2. Dynamics Different from the previous work in [4], the model is extended to include the dynamic features of the plant. Based on a preliminary study, it was observed that the Qpm measurement followed the P measurement is only delayed by the internal flow meter filter. Thus indicating that the valve is significantly slower than the hydrodynamics and the relationship between P and Qpm can be assumed instantaneous. Furthermore, since identical valves and flow meters are used throughout the pilot plant, the hydro and plant dynamics are assumed to be instantaneous for the purposes of control. By assuming instantaneous hydro and plant dynamics, the dynamics which are neglected can potentially be captured by the actuator models as the identification process is based on pressure and flow measurements. The preferred control solution must ensure disturbances in Qf and KVf do not affect the considered outputs. Thus, the slow dynamics for fouling growth (settling time is several hours/days) is not included in the model, but rather assumed to be a disturbance [16,37,38]. 3. Model identification To finalize the model, several constants and functions must be identified. The process is divided into two parts: (i) identification of model constants (KVCF and KVm|c ), and (ii) identification of the dynamics. 3.1. Model coefficients KVm|c is calculated based on manufacturer provided specifications and KVCF is estimated based on a series of experiments at different flow rates (0.06 L/s–1.89 L/s) and pressures (0.1 bar–0.68 bar). The identified coefficients and variance are shown in Table 1. The variance for KVCF over the considered estimation range is small. Thus, the assumed linear relationship between QCF and pressure over the CF channel is concluded to be acceptable.

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

101

Table 1 Model coefficients. Variable

Value

Variance

KVm KVCF

0.833 L · s−1 · bar−1 3.437 L · s−1 · bar−1

N/Aa 0.078 b

a b

Variance not provided by the manufacturer. Variance between the estimated KVCF at different pressures and flow rates.

Fig. 4. Estimated TF coefficients for WP1 , where red is a negative step response and blue is a positive step response, with a step magnitude of 0.2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 3. Examples of dynamic pump responses at different flow conductances, with a step input in UWP1 at 2s and with a magnitude of 20% of the maximum pump speed. WP1 is tested in the CF loop with V2 and V3 closed and V1 at different positions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3.2. Circulation pump The static model formulated in Section 2.2 requires a model describing the dynamic relationship between the normalized pump speed reference and pressure over the pump, as such the aim is to identify the TF relationship: HWP1 (s) =

PWP1 (s) . UWP1 (s)

(10)

However, because of nonlinearities, the TF coefficients are dependent on operating conditions [39]. Fig. 3 show examples of the dynamic behavior of WP1 subject to different pump speeds and flow conductances. The small overshoot observed under some operating conditions combined with the relatively fast and steep dynamic behavior indicate that the dynamics of WP1 is similar to a second order system and dependent on operating conditions. Furthermore, comparing the red and green responses also indicates the dynamic behavior to be directional dependent. Based on these observations, the dynamics of WP1 should be captured by: HWP1 (QWP1 , UWP1 , s) =

DC WP1 ( · )ωn2 ( · ) s2

+ 2( · )ωn ( · )s + ωn2 ( · )

,

(11)

where DCWP1 (·), (·), and ωn (·) are the DC gain, natural frequency, and damping ratio, all of which are functions of QWP1 and UWP1 . Because of the suspected asymmetric behavior (directional dependency) the TF coefficients must be individually identified for both positive and negative step directions. Furthermore, as the dynamics are dependent on operating conditions, the TF coefficients must be identified for different operating conditions. Commonly step inputs are considered poor excitation signals, but to facilitate system identification for different operating conditions and step directions, the step response approach is used. The identified coefficients as a function of operating conditions are shown in Fig. 4, where red and blue are negative and positive step responses. Based on the estimated coefficients, the dynamic behavior of the pump is indeed asymmetric and depended on operating conditions. Preferably, the coefficients should be independent of step direction, while this is the case for the DC gain, it is not the case for ωn and . The directional difference is caused by the design

Fig. 5. Estimated polynomial functions describing the TF coefficients for WP1 as a function of operating conditions. The data points used are the directional mean of the data presented in Fig. 4.

of the pump impeller which is designed for acceleration and not deacceleration of the water. ωn has the highest directional difference and differs at most by a factor of 2.8 and a mean of 2 across the positive and negative direction. It is possible to calculate the RGA for each direction. However, the directional dependent RGA can produce conflicting preferred control pairings depending on the considered direction, in which case a compromise between the directional results must be made. For simplification, it was chosen to make a compromise between the positive and negative direction in an early stage by taking the average over the direction for both ωn and . However, as the system gain is only affected by ωn and  at frequencies around (depending on ) and above ωn , the directional simplification causes uncertainties in the RGA values around and above ωn . The uncertainties introduced by the directional simplification should be considered in the control pairing phase. To describe the relationship between the operating conditions and TF coefficients, three second-degree polynomial functions where chosen, as lower or higher polynomial degree resulted in poor fitness or overfitting. The estimated polynomial functions and the adjusted coefficient of determination (Rˆ 2 ) values are shown in Fig. 5, where especially DCWP1 provides a very good estimation of the DC gain. The fitness of the remaining functions, ωn and ,

102

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

Fig. 6. Two dynamic responses of WP1 . The input is a step at 2s with a magnitude of 20% of the maximum pump speed. The operating point for the transfer function HWP1 is QWP1 = 0.9L/s and UWP1 = 70% pump speed.

Fig. 8. Estimated valve flow conductance as a function of valve position. Table 2 Identified valve coefficients and the range of flow rates through and pressures over the valve during the identification process.

Fig. 7. Estimated dynamic valve behavior.

are less precise, which is likely linked to; (1) insufficient data, (2) inadequate excitation of step inputs, (3) the polynomial model’s inability to describe the relationship, (4) the signal to noise ratio being reduced at the operating limit. Furthermore, the internal controller of the pump may cause abrupt behavioral changes in some operating range. The polynomial functions are generated based on data from the entire operating range of the pump, minimizing the risk that the polynomial functions are used for extrapolation. As the polynomial functions are only used offline, the function values are manually checked to ensure that they are within the expected range before use. Fig. 6 shows step responses for WP1 based on experiments and simulations using the estimated polynomial functions. For each simulation, a single TF is chosen based on the average flow rate and pump speed for that case, to void the complications related to simulation of an LPV TF. While the directional difference can be observed, the simulated response, using the directional mean, is an acceptable compromise between model accuracy and simplicity.

Valve

aVx

bVx

cVx

Flow range (L/s)

Pressure range (bar)

V1 V2 V3

4.0 4.2 7.3

1.2 1.7 0.17

31.5 40.0 0.2

0.4–1.3 0.06–0.77 0.04–0.08

0.15–1.8 1.79–1.81 0.05–1.5

The steady state relationships between PVx and KVVx are different for the each valve, as illustrated in Fig. 8. The estimated KVV3 is actually a combination of the conductance for both V3 and MFM1 connected in series. As the total conductance of two series connected conductances are calculated as shown in Eq. (5), the total conductance cannot exceed the smallest of the flow conductances, which explains why opening V3 beyond 50% has no effect on the total conductance. Consequently, the size of V3 can be reduced without any reduction in the combined flow conductance, while improving the control resolution in the valid operating range. 3.3.2. Dynamic relationship The dynamic relationship for V1 and V3 will be identified in this section, as V2 is excluded as a control valve for the reason mentioned in Section 2.2. The dynamic relationship between the input and the actual valve position is difficult to measure, as the valves have no position feedback. To estimate the dynamic relationship, the static functions in Eq. (12) are inverted;

3.3. Control valves

Pˆ Vx (KV Vx ) = KV −1 Vx (PVx ).

In this section, the model previously used in [4] is extended with dynamic features to facilitate frequency dependent RGA analysis. The general transient performance of the valves is shown in Fig. 7, where a small delay and a settling time of around 7 s can be observed. The valve model is divided into two parts: Firstly, the static relationship between normalized valve position, PVx , and valve conductance, KVVx . Secondly, the dynamic relationship between UVx , and PVx .

The inverted function, can be used to estimate valve position based on the observed KVVx , whereas KVVx can be estimated using:

3.3.1. Static relationship KVVx is estimated based on pressure and flow measurements according to Eq. (2), and as steady state is achieved before estimating KVVx , the set-point in valve position is assumed to be the actual position. However, the hysteresis (dead band) of the valve introduces an error of ±1% of the full operational range. To describe the static steady state relationship the following relationship is used: KV Vx (PVx (t)) =

1 1 + e−aVx · (PVx (t)−bVx ) · cVx

,

(12)

where aVx , bVx , and cVx are coefficients. The estimated functions are shown in Fig. 8, whereas the coefficients of the functions and the range of flow rates through and pressures over the valve used for the identification process, are presented in Table 2.

KV Vx =

QVx



,

(13)

(14)

PVx

where QVx and PVx are the flow rate through and pressure over Vx , respectively. By using Eq. (14), which is valid in steady state, the hydrodynamics propagate through the static functions to be captured by the dynamic relationship. Based on the observed behavior in Fig. 7, the dynamic relationship between valve input and estimated valve position, Pˆ Vx , was chosen to be modelled as a first order system with time delay, as described in: Pˆ Vx (s) DC Vx −s · t d. = e Vx s + 1 UVx (s)

(15)

The estimated  Vx and DCVx for different valve openings are shown in Fig. 9, while the time delay td for each valve was observed to be approximately 0.25 s across different operating conditions. Ideally the estimated DC gains for the valves should be close to 1 at all openings, assuming that the estimated static relationships are accurate. However, this is not the case for valve openings below 25%. The inconsistent DC gains at low openings are likely caused by the small steps in UVx combined with the valve hysteresis. For the estimated

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

103

Table 3 The maximum and minimum open loop BWs and DC gains for the system linearized over the considered operating range. System

BW (Hz)

DC gain (dB)

System order

WP1 to QCF V1 to QCF V3 to QCF WP1 to Qpm V1 to Qpm V3 to Qpm

[0.48, 1.91] 0.11 N/A [0.48, 1.91] 0.11 0.11

[−4.27, 23.74] [−3.47, 13.98] −∞ [−24.27, −9.90] [−46.15, −28.15] [−55.69, −32.20]

2 1 0a 2 1 1

a

The input has no effect on the output, hence a gain of 0 or −∞dB.

Fig. 9. Estimated TF coefficients for different valves and opening degrees.

time constant,  Vx , the value is nearly constant across the different valves and opening degrees. Despite the inconsistent DC gains at low openings, the values for both DCVx and  Vx is assumed constant across the different pressures and flow rates covered in this study (see Table 2). The identified static models are included into the static function f, while the dynamic parts are added in accordance with Fig. 2.

3.4. Combined model structure Fig. 11. Constant and oscillating operating conditions.

The static and dynamic models can be combined as illustrated in Fig. 2. Most of the model is described in terms of linear dynamic and static nonlinearities, in accordance with the Wiener model structure [40]. However, the dynamics of WP1 are described with an LPV TF. To describe the properties of the combined nonlinear relationships between the inputs and outputs, the system is linearized at different operating conditions within the considered range and the resulting envelope (upper and lower bounds) of the Bode plot is shown in Fig. 10. As the valves behave similarly with respect to the frequency, V3 is excluded to increase readability. The highest and lowest bandwidths (BW) and DC gains for the open loop linearized systems, as well as the system orders are summarized in Table 3. According to the model, V3 has no effect on QCF , whereas, in reality, the valve affects QCF to a small degree. This is caused by a combination of the model structure and assumptions; (1) assuming constant Qf , (2) concentrate flow conductance is neglected, and (3) the Qpm is considered to permeate the membrane before entering the CF channel, whereas in reality the permeation occurs along the length of the CF channel.

4. Fouling analysis To study how fouling is affected by oscillating operating conditions, an experiment is conducted on the pilot plant which compares fouling growth subject to constant and oscillating operating conditions. For the experiment two cases are considered: Constant: The dynamic controllers are used to maintain a steady CF velocity (CFV) and P. Oscillating: The reference to the dynamic controller maintaining the P is altered between 0.4 bar and 0.8 bar every 10 s. J (permeate flow per square meter membrane surface), CFV, and P from the two experiments are shown in Fig. 11, where a zoomed view shows the oscillating P and J introduced by manipulating the reference for the P controller. To ensure that the constant and oscillating cases are comparable, the OiW concentration should be identical for each experiment.

Fig. 10. Envelope Bode plot of the open loop systems, where the red and blue dashed lines show the BW for the considered operating range. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

104

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

Table 4 Summary of the oscillating and constant case, where OiW concentration is measured before each experiment.

OiW1b OiW2b Final Jc CFV P Jc a b c

Constant

Oscillating

Differencea

3213 PPM 2911 PPM 118.7 L h−1 m−2 1.9941 m/s 0.5972 bar 1658 L/m2

2960 PPM 2854 PPM 98.5 L h−1 m−2 1.9943 m/s 0.5977 bar 1528 L/m2

8.6% 2% 20.5% −0.01% −0.08% 8.49%

∂yi  ∂uj u

CLG =

∂yi  ∂uj y

k =0,k

= / j

k =0,k

= / i

(17)

For the considered system, the RGA matrix is formulated as:

te

J(t)dt,

OLG =



However, as it is difficult to recreate the exact same OiW concentration, two OiW monitors (Turner 4100XDC) are used to measure the concentration levels. Furthermore, the OiW concentration for the constant case was intentionally higher to avoid favoring the constant case in terms of OiW concentration. Despite the fact that the oscillating case initially had a higher Qpm , which could be linked to the difference in OiW concentration, the constant case clearly accumulated less fouling during the 10 h of filtration. To evaluate membrane performance, the total produced permeate per square meter membrane surface, Jc , is used and defined as:



The RGA values are a measure of interaction, defined as the open loop gain (OLG) divided by the close loop gain (CLG):



Calculated by (constant/oscillating · 100) − 100. The feed concentration is measured using two different Turner 4100XDC. Final J is the average over the last 6 min of the experiment.

Jc =

5.1. Theory

(16)

0

where te is the duration of the experiments. The OiW concentration of the feed, CFV, P, J, and Jc from the two experiments are summarized in Table 4. The table shows that for the constant case an overall 8.49% more permeate was produced and after 10 h of filtration J was 20.5% higher compared to the oscillating case. The difference in fouling between the oscillating and constant case is likely linked to the critical flux hypothesis theorizing that below a critical J, little to no fouling accumulates and above fouling growth rapidly increases [41,42]. The nonlinearity is observed experimentally in [43,44] and is the likely cause for the observed difference. In addition, oscillating CF was observed to decrease fouling growth in [45], in contrast to oscillating TMP/flux increasing fouling growth. While the results concern a different variable, it confirms that fouling can indeed be affected by oscillating operating conditions. Overall oscillating P and Qpm significantly affect the fouling growth, which emphasizes the potential loss in performance caused by poor control.

where 1,1 is the interaction between V1 and Qpm , and 2,2 is the interaction between V3 and Qcf , etc. If the OLG and CLG are nearly identical the RGA value is close to 1 and the output is mainly affected by the considered input, which are preferred. To summarize, the rearranged (preferred pairings on the diagonal of the RGA matrix) system should be the closest match to the identity matrix, to ensure minimum interaction across all the control loops. The RGA values can be very large and negative. RGA numbers above 5 should be avoided as they indicate controllability problems and strong interactions. Negative RGA values are caused by a change in sign when applying feedback control, indicating the existences of RHP-zero in; the considered loop, the combined system, or in the open loop combined system. The presence of a RHP-zero can cause issues as the interaction between control loops or a disabled controller can cause system instability. In short, negative values must be avoided if Decentralized Integral Controllability (DIC) is desired [47,48]. If a system has DIC and is combined with integral control action, the system remains stable even if controllers are disabled, the input saturates, or the controllers are detuned (the feedback gain k is replaced with k, where { ∈ R|0 ≤  ≤ 1}) [47,49–51]. Thus, if possible, pairing on negative RGA elements should be avoided. For systems considering more inputs than outputs, input scaling should be carried out before the analysis [48]. For this study, all inputs are normalized such that the actuator ranges are ∈ [0, 1]. In general, the RGA matrix is evaluated at steady state (s = 0) and all other frequencies are ignored. Furthermore, the RGA method is often applied to linear systems hiding potential problems at different operating conditions [52]. For this study, the RGA analysis is extended to cover both the frequency domain and different operating points. To achieve this, the system gain matrix can be calculated as a function of both frequency and operating point:



G(z0 , d0 , ω) =

∂f (z, d0 )  · diag(H0 (jω)), ∂z z=z

(18)

0

5. RGA analysis

where T

Based on the observed fouling behavior, oscillating operating conditions can increase fouling growth and operating cost. For welldefined process systems where the entire process is controlled and few disturbances occur, the performance may not improve noticeably by control pairings selection. However, for offshore systems where slugging is a common problem [12], reference tracking and disturbance rejection can be essential for minimizing fouling. For decentralized control design, the RGA method can be used for pairing of inputs and outputs, such that interactions between control loops are minimized and reference tracking and disturbance rejection are improved [46].

H0 (jω) = [HV 1 (jω), HV 3 (jω), HWP1 (QWP1,0 , u0 , jω)] ,

(19)

ω is the considered frequency, z0 , d0 , u0 , and QWP1,0 are the steady state operating conditions chosen in z, d, u, and QWP1 , respectively. The RGA matrix can for square systems be calculated according to: (G(z0 , d0 , ω)) = G(z0 , d0 , ω) × (G(z0 , d0 , ω)

−1 T

) ,

(20)

where × denotes the Schur product (element by element multiplication). For the considered membrane system two outputs and three inputs are investigated, hence the system is non-square, as

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111 Table 5 Designed operating point. Description Feed flow rate Permeate flow rate CF rate Conductance of the fouling layer Reject flow rate

5.3. Analytic results Parameter Qf Qpm QCF KVf Qrej

Value −1

0.11 L · s 0.1 L · s−1 1.317 L · s−1 or 2 m · s−1 0.113 L · bar−1 · s−1 0.01 L · s−1

such pseudo-inverse can be used [48]. The RGA values, , can be separated into vectors: Qpm = [1; 1, . . ., 3] QCF = [2; 1, . . ., 3]

105

The RGA results are presented in Figs. 12–17, where Figs. 12 and 13 show the steady state RGA values at different operation conditions for the two considered combinations of actuators, V1 and WP1 or V3 and WP1 , respectively. Figs. 14 and 15 extends the RGA analysis into the frequency domain to highlight potential pairings issues at higher frequencies. Because the frequency dependent RGA plot can be difficult to visually interpret, the data are dimensionally reduced according to: Qpm = max(SQpm ) − max(SQpm {max(SQpm )}), QCF = max(SQCF ) − max(SQCF {max(SQCF )}),

(21)

where Qpm and QCF are the RGA values for the two outputs, and the notation [2 ; 1, . . ., 3] is the 1st to 3th element in the second row. 5.2. Scenario design Ideally the plant should be operated at the designed operating point in terms of Qf , Qpm , QCF , Qrej , and KVf . Based on recommendation from the membrane manufacturer the designed operating point used for the RGA analysis is estimated for steady state and listed in Table 5. However, when applying membrane filtration for offshore OiW separation there are a few parameters which are directly uncontrollable; KVf caused by fouling, and the Qf which is determined by the upstream processes. As these parameters are uncontrollable, the RGA values are analyzed over a reasonable range in these parameters to ensure the chosen control pairings are suitable for the entire range. On the pilot plant, KVf was observed to be between 0.028 and 0.4 L s−1 bar−1 , depending on the current fouling state. For the considered operating mode QCF and Qpm are controlled, hence only two actuators are required to control the system. Multiple actuators could be coordinated to control the system, but as the focus remains to find the control parings for a decentralized control scheme, each considered valve is combined with the pump; V1 , WP1 and V3 , WP1 . For each of these actuator combinations, the RGA values are calculated to determine which control pairings are best suited to control the system. The valves not present in the considered combination (e.g. V2 and V3 for the V1 , WP1 combination) are manually adjusted by trial and error to identify a fixed opening degree that provides a wide control range. The input to the considered combination can then be calculated subject to the model constraints and the operating conditions. Varying KVf : All but KVf are fixed according to Table 5, where the valve and pump are used to maintain the equilibrium in the range; KVf ∈ [0.01, 0.5]. Corresponding to a Qpm between 0.01 and 0.5 L s−1 bar−1 approximately, where the highest observed clean membrane Qpm was 0.4 L s−1 bar−1 . Varying Qf : As the Qf varies the Qpm and Qrej are adjusted according to: Qf ∈ [0, 0.35]

(22a)

Qpm = min{Qf · 0.7, 0.1}

(22b)

Qrej = Qf − Qpm

(22c)

As a result, the permeate is 70% of Qf , but maximum 0.1 L/s, whereas the remaining flow is directed through the reject outlet. QCF and KVf remain fixed according to Table 5. The valve and pump inputs are calculated to maintain the equilibrium point.

(23)

where SQpm = {Qpm },

(24a)

SQCF = {QCF }.

(24b)

Qpm and QCF are an expression for the distance between the two largest RGA elements for the specific output. As all the RGA values are real valued numbers and between 0 and 1, larger distance corresponds to less control loop interaction. The resulting surface plots in Figs. 16 and 17 have a single surface, where the magnitude is the distance between the two largest RGA values for that operating point, and where the colors indicate which actuator has the RGA value closest to one (the preferred actuator). V1 scenario (Figs. 12, 14, and 16): Crossflow; V1 and WP1 are the only actuators which affect QCF . At the designed operating point the preferred actuator is V1 and it is only at higher Qf that the preferred actuator becomes to WP1 . At higher frequencies the RGA values remain the same across different Qf (Figs. 14c and 16c) and across different KVf the distance between the two largest RGA elements is increased, further decoupling V1 and WP1 (Figs. 14a and 16a). Permeate; WP1 is the best candidate for controlling Qpm . In general, WP1 is preferred over the entire considered range in both KVf and Qf . It is only at high Qf where the preferred pairings change. With varying KVf , V1 and WP1 are further decoupled as the frequency is increased (Figs. 14b and 16b), and with varying Qf the decoupling remains unaffected as the frequency increases (Figs. 14d and 16d). V3 scenario (Figs. Figure 13, and Figure 17): Crossflow; V1 and WP1 are both candidates for controlling the CFV. At steady state, V1 is preferred with both varying Qf and KVf . However, for varying KVf and increasing frequency, the decoupling performance is significantly reduced, especially with a KVf above 0.25 L · s−1 · bar−1 (Figure 15 and Figure 17Fig. 15a and 17a). Furthermore, for higher Qf and frequencies the pairing that minimizes interaction changes from V1 to WP1 , which would have been unnoticed if the frequency domain was not included in the analysis (Figs. 15c and 17c). Permeate; V3 and WP1 are the preferred candidates for controlling Qpm . This paring is problematic as WP1 is preferred at the defined operating point, but at higher KVf and both lower and higher Qf , V3 is preferred, indicating that V3 is better suited for controller Qpm in most of the operating range. Across the frequency domain and varying KVf , WP1 becomes slightly better decoupled at the operating point but remains a poor choice for higher KVf (Figs. 15b and 17b). For the varying Qf case, higher frequencies improves the range where WP1 is suitable for controlling Qpm , but the range remains relatively small (Figs. 15d and 17d). Based on the results, the preferred control pairings according to the RGA analysis are; V1 to control QCF and WP1 to control Qpm . However, the results also show that at high Qf the preferred parings changes. This effect is a result of the feed pressure increasing and

106

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

Fig. 12. RGA values when using V1 and WP1 to control the process given the different Qf and KVf , where the remaining valves are fixed at UV2 = 30% and UV3 = 50%.

Fig. 13. RGA values when using V3 and WP1 to control the process given the different Qf and KVf , where the remaining valves are fixed at UV1 = 40% and UV2 = 30%.

Fig. 14. RGA values across different frequencies, KVf , and Qf , where the manipulated variables are UV1 and UWP1 and where the colors represent: V1 is blue, V3 is yellow, WP1 is purple, and the gray mesh grid is the preferred operating points. The remaining valves are fixed at UV2 = 30% and UV3 = 50%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 15. RGA values across different frequencies, KVf , and Qf , where the manipulated variables are UV3 and UWP1 and where the colors represent: V1 is blue, V3 is yellow, WP1 is purple, and the gray mesh grid is the preferred operating points. The remaining valves are fixed at UV1 = 40% and UV2 = 30%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

107

Fig. 16. Qpm and QCF values across different frequencies, KVf , and Qf , where the manipulated variables are UV1 and UWP1 and where the colors represent: V1 is blue, V3 is yellow, WP1 is purple, and the gray mesh grid is the preferred operating points. The remaining valves are fixed at UV2 = 30% and UV3 = 50%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 17. Qpm and QCF values across different frequencies, KVf , and Qf , where the manipulated variables are UV3 and UWP1 and where the colors represent: V1 is blue, V3 is yellow, WP1 is purple, and the gray mesh grid is the preferred operating points. The remaining valves are fixed at UV1 = 40% and UV2 = 30%. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the pumps inability to reduce the pressure. Furthermore, the directional simplification has little effect on the chosen control pairings as the RGA values are well separated, see Fig. 14 .

Table 6 Upper and lower bounds for the closed loop system BW (Hz), and the controllers designed. Bandwidth (Hz)

6. Validation and results To validate the results from the RGA analysis, the preferred control pairing (V1 and QCF , WP1 and Qpm ) are tested and compared to the opposite pairing (unpreferred), while the remaining actuators, V2 and V3 are fixed. To test both the preferred and unpreferred parings, a total of four controllers were developed:

Controllers

QCF

Qpm

QCF

WP1

[0.034, 0.1]

[0.008, 0.044]

0.086 s

V1

[0.006, 0.044]

[0.006, 0.056]

0.018 +

Qpm



1.1 s 0.049 s

− 2.4 +

6.48 s



• Preferred: – CV1 ,QCF – CWP 1 ,Qpm • Unpreferred: – CV1 ,Qpm – CWP 1 ,QCF Where the notation CWP 1 ,QCF is the controller using WP1 to control QCF . The closed-loop (CL) TF can for each controller be defined as: CL(s) =

C(s)H(s) 1 + C(s)H(s)

(25)

where H and C are the transfer functions for the plant and the controller, respectively. The combined CL systems will be denoted as CLV1 ,QCF , CLWP 1 ,QCF , CLV1 ,Qpm , and CLWP 1 ,Qpm . To make a fair comparison, the CL system controlling the same output, e.g. CLV1 ,QCF and CLWP 1 ,QCF should ideally have identical system responses as the objective is to highlight that controllers, tuned to perform individually identical, can be greatly affected by poorly selected control pairings. To achieve this, a set PI-controllers is tuned using the robust response time method from Matlab with tuning goals of a CL bandwidth (BW) of 0.03 Hz and phase margin of 60deg. The method deploys the plant model to tune each controller to get nearly identical CL characteristics ignoring the cross coupling between the controllers. Because the cross coupling is ignored in the tuning phase, the observed difference in controller performance is mainly caused by interactions between the controllers, hence the performance difference is caused by control pairing selection. Table 6

Fig. 18. Simulated CL unit step response for the four controllers, CLV1 ,QCF , CLWP 1 ,QCF , CLV1 ,Qpm , and CLWP 1 ,Qpm . The controllers in each graph, respectively, are designed to give similar output responses.

shows the four designed controller and the respective range of the CL BW, where the CL BW range is defined as the minimum and maximum CL BW for the system linearized across the defined operating range. For the controllers using V1 , a PI controller could meet the tuning goals previously defined, whereas, for the controllers using WP1 , the P-term was selected as zero as only the I-term was necessary to meet the design specification. The simulated response of each controller is shown in Fig. 18, and the envelopes of the CL systems’ bode plots are shown in Fig. 19. The behavior of the CL systems are almost identical, with a small difference in the BW and frequency response within the BW. Thus,

108

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

Fig. 19. Evenlope of the closed loop system’s bode plot. The upper and lower bounds for the BW are marked with dashed lines.

Table 7 Operating conditions for the experimental validation. Description

Parameter

Value

Feed flow rate Permeate flow rate CF rate Conductance of the fouling layer Reject flow rate

Qf Qpm QCF KVf Qrej

0.04 L · s−1 0.036 L · s−1 1.317 L · s−1 or 2 m · s−1 Uncontrollable 0.004 L · s−1

a comparison between the preferred and unpreferred control pairings are assumed to be reasonable. The two set of controllers will be compared based on simulations and experiments. 6.1. Simulations The process model is implemented in Simulink and simulated with both the preferred and unpreferred controllers. The simulation results in Fig. 20 shows that the preferred control pairings are stable and provide good reference tracking, subject to disturbances. However, for the unpreferred control pairing, the system becomes unstable, despite the individual CL being stable. Based on simulation it was observed that the controllers individually performed according to design specification, but once both controllers were in operation the system became unstable, indicating that the interaction between the controllers causes instability, confirming results from the analysis. The controllers are tuned to achieve identical rise time, which amplifies the interaction between the control loops, consequently causing the observed instability. Stability could likely be achieved by time-scale separation of the controllers (selecting different CL BW). However, as the objective is to highlight the possible performance degradation caused by poor control pairings, time-scale separation would conceal the potential magnitude of the problem. 6.2. Experiments For the experiments, the membranes were significant fouled and unable to operate under the previous defined operating point, consequently, the operating conditions are changed to the values in Table 7. According to the RGA analysis, the updated operating point should be of no consequences to the preferred pairings. To replicate the constant Qf assumed in the analysis, an additional controller is deployed to control Qf using an upstream supply pump. The controller maintained Qf even when the considered controllers compensated for disturbances. Furthermore, the setpoint for Qf is periodically manipulated to observe the transient performances of the preferred and unpreferred controllers, subject to process distur-

bance. It would have been preferred to push the system in terms of varying KVf , thus possible with backwashing, it is nearly impossible to replicate the same fouling conditions for the two sets of controllers making the results in incomparable. Furthermore, fouling naturally accumulates very slowly and does not excite the controllers to a degree where the dynamic responses of the controllers can be observed. Contrary to the simulation, the experimental results in Fig. 21 show that the unpreferred pairings do not cause instability. It is conceivable that the unmodelled interaction between the plant and the Qf controller can cause a stabilizing effect. In addition, model assumptions and inaccuracies can also contribute to the deviation between simulation and experiments. For the preferred pairings, time is required for the system to settle after initialization. Once the system has settled, disturbance in Qf does not affect Qpm or QCF significantly. For the unpreferred control pairings, the Qf controller is performing differently compared to the preferred case. The difference is caused by the interaction between the unpreferred controllers and the controller maintaining Qf , as disabling the unpreferred controllers reduced the difference to acceptable levels. The strong interaction between Qf and unpreferred controllers indicates that not only is the preferred pairings causing less interaction between the control loops, it also affects the upstream processes less. For the unpreferred control pairings, compared to the preferred pairings, the setpoints are poorly maintained. To confirm that the observed issues with the unpreferred pairing are caused by control loop interaction, each controller is individually disabled, and the results are shown in Fig. 22. In each case, the performance, in terms of reference tacking, of the active controllers is improved significantly, which confirms that the interaction between the unpreferred controllers is the source of the poor reference tracking in Fig. 21b. Based on the experimental results, the RGA analysis successfully predicted which input/output pairings that resulted in less interaction between the active controllers. The analysis highlights how essential it is to consider inputoutput pairings as the interaction between the controllers, that individually perform acceptably, can cause significant degradation in controller tracking and disturbance rejection. Alternatively, MIMO control schemes and different controller tuning techniques could provide superior reference tracking compared to the methods and techniques used throughout this work, but the aim is to highlight the potential performance different caused mainly by control pairing selection and thereby highlighting the need to consider alternative pairings when SISO control schemes are chosen.

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

109

Fig. 20. Comparison of the preferred and unpreferred control pairings at identical conditions. The dashed green lines are the setpoint to the controllers. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 21. Experimental results for both the preferred and unpreferred control pairing, where the dashed green lines are the references. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 22. Experimental results for the unpreferred control pairing with a single controller disabled, where the dashed green lines are the references. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

7. Conclusion This study highlighted the need to consider control pairings. The different control pairings were analyzed with the RGA method to identify the control pairings that resulted in the least amount of interaction. The interaction analysis was carried out over a range of different operating conditions and frequencies, that the system will likely experience during operation in an offshore environment. Fur-

thermore, the best control pairings according to the RGA analysis, was experimentally tested and compared to a possible worst paring, where it was clear that both the simulations and experiments supported the RGA results. Experimental results showed that constant conditions resulted in 20.5% higher flux after 10 hours of filtration compared to oscillating conditions. Thus, accounting for the interactions to improve the tracking and disturbance rejec-

110

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

tion’s performances of the controllers is essential to minimize fouling. If the system must be controlled using a decentralized control strategy and have the least amount of interaction, the best paring scenario is to pair V1 with QCF and WP1 with Qpm . Logically, it could be argued that the pump, configured as a circulation pump, creates flow and should be used to control the CF. However, results showed that the interaction between the control loops caused poor reference tracking performance. As such, this study highlighted the need to analysis control pairings, especially if decentralized control is to be deployed. It was preferred to include V2 as a possible control valve. However, the orifice of the valve was chosen too large for the reject flow rate in the system, which required the valve to operate in the per mille range causing numerical model issues. In reality, the valve was unable to operate at these openings due to the hysteresis. As process systems are rarely completely decoupled, MIMO control strategies could utilize the couplings instead of avoiding them, to further improve the control systems disturbances reject capabilities. However, for the commonly deployed PID control schemes the control pairing proved to be crucial. Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

References [1] Danish Energy Agency, Oil and Gas Production in Denmark, 2014, Available from: https://ens.dk/sites/ens.dk/files/OlieGas/oil and gas in denmark 2014 pdf (accessed 26.06.18) [Online]. [2] OSPAR-Commission, North Sea Manual on Maritime Oil Pollution Offences, 2012, Available from: https://www.bonnagreement.org/site/assets/files/ 3948/north-sea-manual-on-maritime-oil-pollution-offences.pdf (accessed 26.06.18) [Online]. [3] Danish Environmental Protection Agency, General Authorization for Maersk Oil and Gas, 2016, Available from: http://mst.dk/media/92144/20161221ann-generel-udledningstilladelse-for-maersk-olie-og-gas-2017-18.pdf [Online]. [4] K.L. Jepsen, L. Hansen, Z. Yang, Control parings of a de-oiling membrane process, IFAC-PapersOnLine 51 (8) (2018) 126–131, http://dx.doi.org/10. 1016/j.ifacol.2018.06.366. [5] S.M. Seyed Shahabadi, A. Reyhani, Optimization of operating conditions in ultrafiltration process for produced water treatment via the full factorial design methodology, Sep. Purif. Technol. 132 (2014) 50–61, http://dx.doi.org/ 10.1016/j.seppur.2014.04.051. [6] S.H. Park, T. Yamaguchi, S.I. Nakao, Transport mechanism of deformable droplets in microfiltration of emulsions, Chem. Eng. Sci. 56 (11) (2001) 3539–3548, http://dx.doi.org/10.1016/S0009-2509(01)00047-1. [7] T. Bilstad, E. Espedal, A. Haaland, M. Madland, Ultrafiltration of oily wastewater, in: SPE Health, Safety and Environment in Oil and Gas Exploration and Production Conference, Society of Petroleum Engineers, 1994, pp. 1–5, http://dx.doi.org/10.2118/27136-MS. [8] A. Ezzati, E. Gorouhi, T. Mohammadi, Separation of water in oil emulsions using microfiltration, Desalination 185 (1–3) (2005) 371–382, http://dx.doi. org/10.1016/j.desal.2005.03.086. [9] J. Cakl, I. Bauer, P. Dole, P. Mikul, Effects of backflushing conditions on permeate flux in membrane crossflow microfiltration of oil emulsion, Desalination 127 (2000) 189–198. [10] S.H.D. Silalahi, T. Leiknes, Cleaning strategies in ceramic microfiltration membranes fouled by oil and particulate matter in produced water, Desalination 236 (1–3) (2009) 160–169, http://dx.doi.org/10.1016/j.desal. 2007.10.063. [11] S. Judd, H. Qiblawey, M. Al-Marri, C. Clarkin, S. Watson, A. Ahmed, S. Bach, The size and performance of offshore produced water oil-removal technologies for reinjection, Sep. Purif. Technol. 134 (2014) 241–246, http://dx.doi.org/10. 1016/j.seppur.2014.07.037. [12] K. Jepsen, M. Bram, S. Pedersen, Z. Yang, Membrane fouling for produced water treatment: a review study from a process control perspective, Water 10 (7) (2018) 847, http://dx.doi.org/10.3390/w10070847. [13] A. Abbas, Model predictive control of a reverse osmosis desalination unit, Desalination 194 (1–3) (2006) 268–280, http://dx.doi.org/10.1016/j.desal. 2005.10.033. [14] T. Darvishzadeh, N.V. Priezjev, Effects of crossflow velocity and transmembrane pressure on microfiltration of oil-in-water emulsions, J.

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

Membr. Sci. 423-424 (2012) 468–476, http://dx.doi.org/10.1016/j.memsci. 2012.08.043. K. Guerra, J. Pellegrino, J.E. Drewes, Impact of operating conditions on permeate flux and process economics for cross flow ceramic membrane ultrafiltration of surface water, Sep. Purif. Technol. 87 (2012) 47–53, http://dx. doi.org/10.1016/j.seppur.2011.11.019. S.R.H. Abadi, M.R. Sebzari, M. Hemati, F. Rekabdar, T. Mohammadi, Ceramic membrane performance in microfiltration of oily wastewater, Desalination 265 (1–3) (2011) 222–228, http://dx.doi.org/10.1016/j.desal.2010.07.055. H.J. Tanudjaja, V.V. Tarabara, A.G. Fane, J.W. Chew, Effect of cross-flow velocity, oil concentration and salinity on the critical flux of an oil-in-water emulsion in microfiltration, J. Membr. Sci. 530 (February) (2017) 11–19, http://dx.doi.org/10.1016/j.memsci.2017.02.011. M. Stoller, R.S. Mendes, Advanced control system for membrane processes based on the boundary flux model, Sep. Purif. Technol. 175 (2017) 527–535, http://dx.doi.org/10.1016/j.seppur.2016.09.049. J. Busch, W. Marquardt, Run-to-run control of membrane filtration in wastewater treatment – an experimental study, IFAC Proc. Vol. 40 (5) (2007) 195–200, http://dx.doi.org/10.3182/20070606-3-MX-2915.00080. B. Espinasse, P. Bacchin, P. Aimar, On an experimental method to measure critical flux in ultrafiltration, Desalination 146 (1–3) (2002) 91–96, http://dx. doi.org/10.1016/S0011-9164(02)00495-2. S.G. Lehman, L. Liu, Application of ceramic membranes with pre-ozonation for treatment of secondary wastewater effluent, Water Res. 43 (7) (2009) 2020–2028, http://dx.doi.org/10.1016/j.watres.2009.02.003. B. Espinasse, P. Bacchin, P. Aimar, Filtration method characterizing the reversibility of colloidal fouling layers at a membrane surface: analysis through critical flux and osmotic pressure, J. Colloid Interface Sci. 320 (2) (2008) 483–490, http://dx.doi.org/10.1016/j.jcis.2008.01.023. M. Hung, J. Liu, Microfiltration for separation of green algae from water, Colloids Surf. B: Biointerfaces 51 (2) (2006) 157–164, http://dx.doi.org/10. 1016/j.colsurfb.2006.07.003. S. Munirasu, M.A. Haija, F. Banat, Use of membrane technology for oil field and refinery produced water treatment – a review, Process Saf. Environ. Prot. 100 (2016) 183–202, http://dx.doi.org/10.1016/j.psep.2016.01.010. S. Nazirah Wan Ikhsan, N. Yusof, F. Aziz, N. Misdan, A review of oilfield wastewater treatment using membrane filtration over concentional technology, Malays. J. Anal. Sci. 21 (3) (2017) 643–658, http://dx.doi.org/10. 17576/mjas-2017-2103-14. A. Fakhru’l-Razi, A. Pendashteh, L.C. Abdullah, D.R.A. Biak, S.S. Madaeni, Z.Z. Abidin, Review of technologies for oil and gas produced water treatment, J. Hazard. Mater. 170 (2–3) (2009) 530–551, http://dx.doi.org/10.1016/j. jhazmat.2009.05.044. J.M. Dickhout, J. Moreno, P.M. Biesheuvel, L. Boels, R.G. Lammertink, W.M. de Vos, Produced water treatment by membranes: a review from a colloidal perspective, J. Colloid Interface Sci. 487 (2017) 523–534, http://dx.doi.org/10. 1016/j.jcis.2016.10.013. M. Padaki, R. Surya Murali, M.S. Abdullah, N. Misdan, A. Moslehyani, M.A. Kassim, N. Hilal, A.F. Ismail, Membrane technology enhancement in oil–water separation. A review, Desalination (2015) 197–207, http://dx.doi.org/10. 1016/j.desal.2014.11.023. G. Ferrero, H. Monclús, G. Buttiglieri, J. Comas, I. Rodriguez-Roda, Automatic control system for energy optimization in membrane bioreactors, Desalination 268 (1–3) (2011) 276–280, http://dx.doi.org/10.1016/j.desal. 2010.10.024. P. James, S. Vigneswaran, H. Hao, R. Ben-aim, H. Nguyen, Design of a generic control system for optimising back flush durations in a submerged membrane hybrid reactor, J. Membr. Sci. 255 (2005) 99–106, http://dx.doi.org/10.1016/j. memsci.2005.01.026. J. Busch, W. Marquardt, Run-to-run control of membrane filtration processes, IFAC Proc. Vol. 39 (2) (2006) 1003–1008, http://dx.doi.org/10.3182/ 20060402-4-BR-2902.01003. A.R. Bartman, C.W. McFall, P.D. Christofides, Y. Cohen, Model-predictive control of feed flow reversal in a reverse osmosis desalination process, J. Process Control 19 (3) (2009) 433–442, http://dx.doi.org/10.1016/j.jprocont. 2008.06.016. T. Mc Avoy, Y. Arkun, R. Chen, D. Robinson, P. Schnelle, A new approach to defining a dynamic relative gain, Control Eng. Pract. 11 (8) (2003) 907–914, http://dx.doi.org/10.1016/S0967-0661(02)00207-1. F. Ferranti, Y. Rolain, A local approach for the modeling of linear parameter-varying systems based on transfer function interpolation with scaling coefficients, IEEE Instrum. Meas. Technol. Conf. (2015), http://dx.doi. org/10.1109/I2MTC.2015.7151337. M.Y. Jaffrin, Dynamic shear-enhanced membrane filtration: a review of rotating disks, rotating membranes and vibrating systems, J. Membr. Sci. 324 (1–2) (2008) 7–25, http://dx.doi.org/10.1016/j.memsci.2008.06.050. D.C. Sioutopoulos, A.J. Karabelas, Correlation of organic fouling resistances in RO and UF membrane filtration under constant flux and constant pressure, J. Membr. Sci. 407–408 (2012) 34–46, http://dx.doi.org/10.1016/j.memsci.2012. 03.036. B. Hu, K. Scott, Microfiltration of water in oil emulsions and evaluation of fouling mechanism, Chem. Eng. J. 136 (2-3) (2008) 210–220, http://dx.doi.org/ 10.1016/j.cej.2007.04.003. L.G. Fernandez, C.O. Soria, C.A. Garcia Tourn, M.S. Izquierdo, The study of oil/water separation in emulsion by membrane technology, SPE Latin

K.L. Jepsen et al. / Journal of Process Control 81 (2019) 98–111

[39]

[40]

[41]

[42]

[43]

[44]

[45]

American and Caribbean Petroleum Engineering Conference, no. v. Society of Petroleum Engineers (2001), http://dx.doi.org/10.2118/69554-MS. K. Jepsen, L. Hansen, C. Mai, Z. Yang, Power consumption optimization for multiple parallel centrifugal pumps, 1st Annual IEEE Conference on Control Technology and Applications, CCTA 2017, vol. 2017-Janua (2017), http://dx. doi.org/10.1109/CCTA.2017.8062559, ISBN 9781509021826. Y. Zhu, Estimation of an N-L-N hammerstein-wiener model, IFAC Proc. Vol. (IFAC-PapersOnline) 15 (1) (2002) 247–252, http://dx.doi.org/10.1016/S00051098(02)00062-6. P. Bacchin, B. Espinasse, P. Aimar, Distributions of critical flux: modelling, experimental analysis and consequences for cross-flow membrane filtration, J. Membr. Sci. 250 (1–2) (2005) 223–234, http://dx.doi.org/10.1016/j.memsci. 2004.10.033. R.W. Field, D. Wu, J.A. Howell, B.B. Gupta, Critical flux concept for microfiltration fouling, J. Membr. Sci. 100 (3) (1995) 259–272, http://dx.doi. org/10.1016/0376-7388(94)00265-Z. F. Wicaksana, A.G. Fane, P. Pongpairoj, R. Field, Microfiltration of algae (Chlorella sorokiniana): critical flux, fouling and transmission, J. Membr. Sci. 387–388 (1) (2012) 83–92, http://dx.doi.org/10.1016/j.memsci.2011.10.013. M. Gander, B. Jefferson, S. Judd, Aerobic MBRs for domestic wastewater treatment: a review with cost considerations, Sep. Purif. Technol. 18 (2) (2000) 119–130, http://dx.doi.org/10.1016/S1383-5866(99)00056-8. P. Blanpain-Avet, N. Doubrovine, C. Lafforgue, M. Lalande, The effect of oscillatory flow on crossflow microfiltration of beer in a tubular mineral

[46]

[47]

[48] [49]

[50]

[51]

[52]

111

membrane system – membrane fouling resistance decrease and energetic considerations, J. Membr. Sci. 152 (2) (1999) 151–174, http://dx.doi.org/10. 1016/S0376-7388(98)00214-2. E.H. Bristol, On a new measure of interaction for multivariable process control, IEEE Trans. Autom. Control AC11 (1) (1966) 133–134, http://dx.doi. org/10.1109/TAC.1966.1098266. C.-C. Yu, M.K. Fan, Decentralized integral controllability and D-stability, Chem. Eng. Sci. 45 (11) (1990) 3299–3309, http://dx.doi.org/10.1016/00092509(90)80221-Y. S. Skogestad, I. Postlethwaite, Multivariable Feedback Control: Analysis and Design, 2nd ed., Wiley & Sons, 2005, ISBN 978-0-470-01167-6. M. Hovd, S. Skogestad, Simple frequency-dependent tools for control system analysis, structure selection and design, Automatica 28 (5) (1992) 989–996, http://dx.doi.org/10.1016/0005-1098(92)90152-6. J. Lee, T.F. Edgar, Conditions for decentralized integral controllability, J. Process Control 12 (7) (2002) 797–805, http://dx.doi.org/10.1016/S09591524(01)00056-7. ˜ Arranz, W. Birk, G. Nikolakopoulos, A survey on control M. Castano configuration selection and new challenges in relation to wireless sensor and actuator networks, IFAC-PapersOnLine 50 (1) (2017) 8810–8825, http://dx. doi.org/10.1016/j.ifacol.2017.08.1536. S. Glad, Extensions of the RGA concept to nonlinear systems, European Control Conference, ECC99, no. 2 (2000) F332.