# Model Predictive Control Strategies Performance Evaluation over a Pipeline Transportation System

Authors: Jorge Cardenas-Cabrera; Luis Diaz-Charris; Andrés Torres-Carvajal; Narciso Castro-Charris; Elena Romero-Fandiño; José David Ruiz Ariza; Javier Jiménez-Cabas
Journal: Journal of Control Science and Engineering
 (2019)

---

## Abstract

In several industries using pipelines to transport different products from one point to another is a common and indispensable process, especially at oil/hydrocarbon industries. Thus, optimizing the way this process is carried out must be an issue that cannot be stopped. Therefore, the performance of the control strategy implemented is one way of reaching such optimal operating zones. This study proposes using Model Predictive Control strategies for solving some issues related to the proper operation of pipelines. It is proposed a model based on physics and thermodynamic laws, using MATLAB® as the development environment. This model involves four pumping stations separated by three pipeline sections. Three MPC strategies are developed and implemented. Accordingly, the results indicate that a centralized controller with an antiwindup back-calculation method has the best results among the three configurations used.


---

## Body

## 1. Introduction

Transporting various fluids through pipelines is one of the most used transport mechanisms worldwide, mainly because it allows transporting large volumes in relatively short periods [1]. It is estimated that there are more than 3.500.000 (km) of pipelines in about 120 countries of the world and around 134.866 (km) of pipelines are planned and under construction worldwide [2]. Given the fatalities and severe environmental and economic consequences of potential failures or downtime, pipeline transportation systems have high requirements regarding safety, reliability, and profitability. Thus, it is an ongoing concern for pipeline system operators to guarantee increasingly better performance.Some works have been performed by researchers to enhance the performance of pipeline transportation systems [3, 4] by introducing a set of sufficient conditions to guarantee asymptotic stability for hydraulic pipeline systems described by the so-called water hammer equations. These conditions are then integrated into a receding horizon optimal control scheme. Although the proposed approach was validated in simulations, the dynamics of essential elements as sensors and actuators were not taken into account. Yuzhanin et al. [5] addressed the problem of maintaining pressure below the maximum rated load capacity in a linear section of an oil pipeline. The proposed algorithm can suggest the automatic control action that will allow dampening potentially dangerous pressure surges in the pipeline. The simulations results show an efficient performance of the algorithm. Osiadacz and Chaczykowski [6] described an algorithm for optimal control of a gas network. Since the compressors station, which forms a significant part of the operational plant of a gas transportation system, is driven by gas turbines which use natural gas as fuel, authors considered optimality concerning fuel consumption. Thus, minimizing this fuel usage is a primary objective in the optimal control problem considered. A PLC-based Fuzzy-PID controller is proposed in [7] to regulate the flow rate of the petroleum products at destination by controlling the multiple pressure signals range along the pipeline. Unlike the previously mentioned cases, whose results were tested through simulations, the performance of the Fuzzy-PID is assessed in a lab-scale experimental setup. Real-time experimental results showed that Fuzzy-PID controller provides better performance than a Cascade-PID controller.Motivated by the fact that model predictive control (MPC) is the industry chosen technology for advanced control in many processes even over other advanced process control techniques such as constraint control, nonlinear control and fuzzy logic control [8]; in this work, its performance is studied on a pipeline transportation systems. MPC control strategies are based on using a model to predict process behavior shortly (moving horizon window), and then the future manipulated variable trajectory is computed over the same horizon by solving a constrained optimization problem. Only the first sample of the calculated trajectory is applied to the actuator, and the horizon window is moved one sample forward, and the optimization process is repeated using new measurement data. Most popular MPC strategies found in literature are DMC (Dynamic Matrix Control) [9], GPC (Generalized Predictive Control) [10], IMC (Internal Model Control) [11], RHC (Receding Horizon Control) [4], NMPC (Nonlinear Model Predictive Control) [12, 13], among others [14].The essential contribution of this work is to explore and to implement in simulation three different MPC control configurations for a pipeline transportation process:(1) two decoupled control systems, an MPC based control system for pressure and an on-off control system for level, (2) two decoupled control systems for pressure and level based on MPC strategy, and (3) one MPC control system for pressure and level. Additionally, it has been proposed an alternative option for constraint handling which involves applying a back-calculation scheme. The addition of back-calculation as an antiwindup method allows reducing the control algorithms execution time significantly while maintaining the same performance as the achieved with the optimization based constraint handling method.The rest of the paper is organized in the following pattern. Section2 presents the proposed linear model that describes the fluid behavior in a pipeline and its validation by using data of a real pipeline used for transporting Liquefied Petroleum Gas. On the other hand, models for instrumentation and additional systems are presented too. Section 3 exposes the used Discrete Model Predictive Control strategy and the back-calculation scheme proposed as constrains handling method. Section 4, on the one hand, presents a comparison between the performances of the three control configurations proposed. Finally, Section 5 provides some conclusions.


## 2. Fluid Dynamics Modeling

### 2.1. Pipeline Model

By assuming negligible convective changes in velocity, as well as constant cross-sectional area, the momentum and continuity equations governing the dynamics of the fluid in a horizontal pipeline can be expressed as [15–18]:(1)dQt,zdt+APρt,zdPt,zdz+fQt,zQt,z2ϕAP=0(2)dPt,zdt+ρt,za2APdQt,zdz=0where z,t∈0,Lp×0,∞gathers the space m and time s coordinates, respectively, Lp is the length of the pipe, P is the pressure Pa, Q is the flow rate m3/s, f is the Darcy–Weisbach friction factor, AP is the cross-sectional area of the pipe m2, ϕ is the inside diameter of the pipe m, a is the wave speed in the fluid m/s, and ρ is the fluid density.Due complexity to solve this set of equations, several linearized versions of this set have been proposed [19–22]. At the current study, the linearized model (3) is used to design control strategies.(3)PLPQ0ρLPs=1cosh⁡nLP-ZKtanh⁡nLP-Q0,tTAPs1cosh⁡nLP-e-τs1ZKtanh⁡nLP1cosh⁡nLP-Q0,tTZKAPstanh⁡nLP00e-τs·P0QLPρ0swhere ZK=Ls+R/CsR=ρfQ/ϕAPC=1/ρa2T=λQQ/2ϕAP2L=ρ/AP, n=Ls+RCs, and the subindexes 0 and LP indicate input and output.In this work, the Liquefied Petroleum Gas (LPG) pipeline presented in [22] would be considered (see Figure 1).Figure 1
Liquefied Petroleum Gas pipeline.The physical parameters of the pipeline are shown in Table1.Table 1
Pipeline model parameters.
Variable Symbol Value Units Variable Symbol Value Units Density ρ 530 Kg/m3 Pipelines Diameter d 0.5 m  Friction Coefficient λ 0.00965  Tanks Height H m a x 15.7920 m  Fluid Speed of Sound a 2490 m/s Tanks Diameter D t 52.9032 m  Pipelines Lengths L p 1 10326  Pipelines Inclinations α 1 1.93×10−4 L p 2 18560 m α 2 3.77×10−4 rad L p 3 14980  α 3 5.47×10−4LPG is a flammable mixture of hydrocarbon (primarily propane and butane), gaseous at standard temperature (15°C) and pressure (0.1MPa), but easily liquefied by increasing pressure (see Figure2). Due to its ease of handling, cleanness of combustion, high purity, and high heating values; LPG is used as fuel in a large variety of industrial, commercial and domestic uses [23].Figure 2
Pressure-Temperature Diagram of (40% butane–60% propane) LPG.


### 2.2. Instrumentation and Additional Systems

In order to get closer to a realistic situation, the pipeline model was enriched with models of storage tanks, sensors/transmitters, and actuators. As it refers to fuel storage modeling, membrane tanks are commonly used. These tanks are designed to provide a constant inside pressure, which avoids phase changes in the system. Inlet and outlet taps in the tanks are positioned at the same height in order to have pressure changes only caused by the liquid column. The dynamic model for this kind of systems can be expressed as(4)dhtdt=1ATQint-Qouttwhere ht is the level in the tank in m, AT is the cross-sectional area of the tank m2, Qin and Qout are the flows of liquid into and from the tank, respectively m3/s.On the other hand, first-order models were considered for instruments (sensors/transmitters and actuators). There were added centrifugal pumps at the beginning of each pipeline section, control valves at each tank's inlet section, and level transmitters at each tank, as well as manual valves. Figure3 shows the piping and instrumentation diagram (P&ID) for a pipeline section. The total system is constituted by three of the sections shown in Figure 3.Figure 3
P&ID for a pipeline section.Table2 shows critical parameters of the instruments used to model its dynamical behavior. Taking into account the pipeline and instruments parameters, the obtained model corresponds to (5)PL1PL2PL3︸CPL=G11s00G21sG22s0G31sG32sG33s︸GPLsMp1Mp2Mp3︸Mp where PLi is the pressure of pipeline sectioni at the outlet and Mpi is pressure drop in pumpi. The form of the model representing the effect of valves position, Mvi on the pressure at outlet section of each pipeline section, is represented in (6)PL1PL2PL3︸CPL=G14s00G24sG25s0G34sG35sG36s︸GdsMv1Mv2Mv3︸MvTable 2
Sensors and final control elements features.
Instrument Symbol Value Units Pumps Δ p m a x 1 MPa τ p 30 s m - s 50 %  Pressure Transmitter span 6 MPa τ P T 1 s  Level Transmitter span 20 m τ L T ≈ 0 s  Valves K v 1 6×10-3 m 3/(s.%) K v 2 6×10-3 m 3/(s.%) K v 3 2.8×10-3 m 3/(s.%) τ p 45 s m - 50 %Moreover, the model that links valve position and tanks level takes the form shown in (7). Finally, it is important to mention manipulated (Mp,Mv), and controlled variables (CPL,CL) move in the range of 0 to 100%.(7)H1H2H3︸CL=G44sG45s00G55sG56s00G66s︸GLsMv1Mv2Mv3︸MvBy combining (5), (6), and (7) one obtains the perturbed system form(8)CPLs=GPLsMps+GdsMvsCLs=GLsMvsDiscrete-time state space models that approximate the continuous-time models (8) were computed, with sampling time 1s.(9)
G
P
L
z
=
A
12
×
12
P
L
B
12
×
3
P
L
1
2
1
2
C
3
×
12
P
L
0
3
×
3
G
d
z
=
A
12
×
12
d
B
12
×
3
d
1
2
1
2
C
3
×
12
d
0
3
×
3
G
L
z
=
A
10
×
10
L
B
10
×
3
L
1
2
1
2
C
3
×
10
L
0
3
×
3In (9) it shows the matrix dimensions for each model, where Gz=ABCD=CzI-A-1B+D.


## 2.1. Pipeline Model

By assuming negligible convective changes in velocity, as well as constant cross-sectional area, the momentum and continuity equations governing the dynamics of the fluid in a horizontal pipeline can be expressed as [15–18]:(1)dQt,zdt+APρt,zdPt,zdz+fQt,zQt,z2ϕAP=0(2)dPt,zdt+ρt,za2APdQt,zdz=0where z,t∈0,Lp×0,∞gathers the space m and time s coordinates, respectively, Lp is the length of the pipe, P is the pressure Pa, Q is the flow rate m3/s, f is the Darcy–Weisbach friction factor, AP is the cross-sectional area of the pipe m2, ϕ is the inside diameter of the pipe m, a is the wave speed in the fluid m/s, and ρ is the fluid density.Due complexity to solve this set of equations, several linearized versions of this set have been proposed [19–22]. At the current study, the linearized model (3) is used to design control strategies.(3)PLPQ0ρLPs=1cosh⁡nLP-ZKtanh⁡nLP-Q0,tTAPs1cosh⁡nLP-e-τs1ZKtanh⁡nLP1cosh⁡nLP-Q0,tTZKAPstanh⁡nLP00e-τs·P0QLPρ0swhere ZK=Ls+R/CsR=ρfQ/ϕAPC=1/ρa2T=λQQ/2ϕAP2L=ρ/AP, n=Ls+RCs, and the subindexes 0 and LP indicate input and output.In this work, the Liquefied Petroleum Gas (LPG) pipeline presented in [22] would be considered (see Figure 1).Figure 1
Liquefied Petroleum Gas pipeline.The physical parameters of the pipeline are shown in Table1.Table 1
Pipeline model parameters.
Variable Symbol Value Units Variable Symbol Value Units Density ρ 530 Kg/m3 Pipelines Diameter d 0.5 m  Friction Coefficient λ 0.00965  Tanks Height H m a x 15.7920 m  Fluid Speed of Sound a 2490 m/s Tanks Diameter D t 52.9032 m  Pipelines Lengths L p 1 10326  Pipelines Inclinations α 1 1.93×10−4 L p 2 18560 m α 2 3.77×10−4 rad L p 3 14980  α 3 5.47×10−4LPG is a flammable mixture of hydrocarbon (primarily propane and butane), gaseous at standard temperature (15°C) and pressure (0.1MPa), but easily liquefied by increasing pressure (see Figure2). Due to its ease of handling, cleanness of combustion, high purity, and high heating values; LPG is used as fuel in a large variety of industrial, commercial and domestic uses [23].Figure 2
Pressure-Temperature Diagram of (40% butane–60% propane) LPG.


## 2.2. Instrumentation and Additional Systems

In order to get closer to a realistic situation, the pipeline model was enriched with models of storage tanks, sensors/transmitters, and actuators. As it refers to fuel storage modeling, membrane tanks are commonly used. These tanks are designed to provide a constant inside pressure, which avoids phase changes in the system. Inlet and outlet taps in the tanks are positioned at the same height in order to have pressure changes only caused by the liquid column. The dynamic model for this kind of systems can be expressed as(4)dhtdt=1ATQint-Qouttwhere ht is the level in the tank in m, AT is the cross-sectional area of the tank m2, Qin and Qout are the flows of liquid into and from the tank, respectively m3/s.On the other hand, first-order models were considered for instruments (sensors/transmitters and actuators). There were added centrifugal pumps at the beginning of each pipeline section, control valves at each tank's inlet section, and level transmitters at each tank, as well as manual valves. Figure3 shows the piping and instrumentation diagram (P&ID) for a pipeline section. The total system is constituted by three of the sections shown in Figure 3.Figure 3
P&ID for a pipeline section.Table2 shows critical parameters of the instruments used to model its dynamical behavior. Taking into account the pipeline and instruments parameters, the obtained model corresponds to (5)PL1PL2PL3︸CPL=G11s00G21sG22s0G31sG32sG33s︸GPLsMp1Mp2Mp3︸Mp where PLi is the pressure of pipeline sectioni at the outlet and Mpi is pressure drop in pumpi. The form of the model representing the effect of valves position, Mvi on the pressure at outlet section of each pipeline section, is represented in (6)PL1PL2PL3︸CPL=G14s00G24sG25s0G34sG35sG36s︸GdsMv1Mv2Mv3︸MvTable 2
Sensors and final control elements features.
Instrument Symbol Value Units Pumps Δ p m a x 1 MPa τ p 30 s m - s 50 %  Pressure Transmitter span 6 MPa τ P T 1 s  Level Transmitter span 20 m τ L T ≈ 0 s  Valves K v 1 6×10-3 m 3/(s.%) K v 2 6×10-3 m 3/(s.%) K v 3 2.8×10-3 m 3/(s.%) τ p 45 s m - 50 %Moreover, the model that links valve position and tanks level takes the form shown in (7). Finally, it is important to mention manipulated (Mp,Mv), and controlled variables (CPL,CL) move in the range of 0 to 100%.(7)H1H2H3︸CL=G44sG45s00G55sG56s00G66s︸GLsMv1Mv2Mv3︸MvBy combining (5), (6), and (7) one obtains the perturbed system form(8)CPLs=GPLsMps+GdsMvsCLs=GLsMvsDiscrete-time state space models that approximate the continuous-time models (8) were computed, with sampling time 1s.(9)
G
P
L
z
=
A
12
×
12
P
L
B
12
×
3
P
L
1
2
1
2
C
3
×
12
P
L
0
3
×
3
G
d
z
=
A
12
×
12
d
B
12
×
3
d
1
2
1
2
C
3
×
12
d
0
3
×
3
G
L
z
=
A
10
×
10
L
B
10
×
3
L
1
2
1
2
C
3
×
10
L
0
3
×
3In (9) it shows the matrix dimensions for each model, where Gz=ABCD=CzI-A-1B+D.


## 3. Implemented Control Strategy

By increasing pressures, above the bubble curve for a given temperature or density (see Figure2), LPG can be held in the liquid phase, which simplifies transporting, measuring, and storage. Once pressure is released, LPG vaporizes and becomes in an ideal energy source for fuel. In order to maintain LPG the liquid phase, tight pressure control is necessary to avoid possible phase changes during the transportation process. Additionally, it is well known that oscillations and rapid changes of pressure within the pipeline can cause undesirable phenomena as cavitation, slugging and hydraulic shock (water hammer). Then, these aspects must be considered in the control strategies performance analysis.


### 3.1. Discrete Model Predictive Control Strategy

The function-based discrete time model predictive control (DMPC) formulation proposed by Wang [24] was applied to the LPG pipeline. The DMPC strategy considered is based on state space models and was implemented to control the downstream pressure of each pipeline section (CPL) and the liquid level of each storage tank (CL), by manipulating the pressure drop in pumps (MP) and control valves positions (Mv), respectively, while respecting the physical limitations of the actuators. A set of discrete time Laguerre function is used as an orthonormal basis to describe the optimum future control trajectory (that is, the difference of the manipulated variable) within a prediction horizon. The use of the Laguerre function allows reducing the computational burden associated with MPC. Below the DMPC strategy is briefly presented.


#### 3.1.1. Augmented State Space Model

As mentioned above, at the current study the MPC strategy used is based on state space models. The general extended state space model in (12) is employed. Recall that we embedded integrators in the state space model for obtaining integral action.(10)Δxpk+1ek+1︸xk+1=Ap0CpAp1︸AeΔxpkek︸xk+BpCpBp︸BeΔukek=01︸CeΔxpkykwhere subindex p indicates process, Δxpk+1=xpk+1-xpkΔuk=uk-uk-1 and ek=csetk-ck. Moreover, the poles of the new matrix are the poles of the process and one pole located at 1. So, firstly it was confirmed the model was controllable and observable. Then, on the other hand, it was necessary to estimate the states, which are not measured; in this sense a Kalman filter was used [25], basically because, in the presence of noise and uncertainties, it produces state estimates that tend to be more precise. Kalman filter is commonly used for MPC vendor products to provide state estimation functionality [26]. Additionally, it is worth mentioning that the performance of the control strategies was evaluated by using the nonlinear pipeline model described by (1) and (2). By using the finite difference method, a spatial-discrete version of the nonlinear model [15, 27] was obtained, which was used in closed loop simulations. This way, Figure 4 shows the block diagram of the implemented strategy, G the process to be controlled, and GC the DMPC controller.Figure 4
A general block diagram of the implemented strategy.


#### 3.1.2. DMPC Formulation

The procedure established at [24] is centered on an orthonormal-base set of discrete functions to approximate the future differences of the manipulated variable Δmki+i=ΔmkiΔmki+1⋯Δmki+Nc-1T, where NC indicates the control horizon. These functions are known as discrete Laguerre functions, and its z-transform is shown in (11).(11)Γkz=z-1-a1-az-1Γk-1z,Γ1z=1-a21-az-1The pole of the function a determines the velocity at which lk converges to zero. Then, Δm is approximated by N Laguerre functions, such as in (12), where ki represents the current sampling period, k are the future sampling periods, N is the number of terms used at the expansion, and cj are the coefficients of the process which are the function of future sampling periods.(12)Δmki+k=∑j=1NcjkiljkLaguerre functions can be expressed in the state space form in (13), which reduce computational burden. Al is a N×N matrix, which is a function of a.(13)Lk+1=AlLk,L0T=1-a21-aa2-a3⋯-1N-1aN-1TTherefore (12) can also be represented as a function of the Laguerre vector defined at (13) such as(14)Δmki+k=LkTη,η=c1c2⋯cNTThus, the prediction of future states and outputs in an instantm is defined by (15).(15)xki+m∣ki=Amxki+∑i=0m-1Am-i-1BLiTηyki+m∣ki=Cxki+m∣kiThen, the cost function to be minimized can be written as (16), where Q=CTC and RL>0. It is necessary to note J is a function of η instead of ΔM, so minimization must be done in function of η.(16)J=∑m=1Npxki+m∣kiTQxki+m∣ki+ηTRLηResults of the minimization operation are shown in (17)η=-Ω-1ΨxkiΩ=∑m=1NpϕmQϕmT+RLΨ=∑m=1NpϕmQAmϕmT=∑i=0m-1Am-i-1BLiTTherefore, the control law is written at (18).(18)Δmk=kDMPCx^k


#### 3.1.3. Constraint Handling

To take into account physical actuators limitations, a constraint is imposed on the control signalΔmk, which is mathematically defined in (19)mmin≤mk≤mmax By using the Laguerre functions, it is possible to rewrite expressions in (19) as shown in (20)mmin≤∑i=0k-1LiTη+mk-1≤mmaxFurthermore, it is important to note that constrained control strategies require real-time optimization using quadratic programming [24], where the minimization function would be (21). This restricted optimization must be done at the entire prediction horizon, which increases the computational burden.(21)J=ηTΩη+2ηTΨxkiThus, at this research, a second option for constraining handling has been proposed which involves applying a back-calculation scheme, Figure 5.Figure 5
Constraint handling scheme using the antiwindup back-calculation method.Both methods were implemented, and simulations were run considering one change in the pressure set point of station 1. Results showed that both methods back-calculation and constraint handling based on quadratic optimization are very similar, as can be seen in Figure6. As it refers to the IAE value, it is observed that the back-calculation method is slightly better than the constrain-handling method for the setpoint change in station 1. However, even though the effect over this change over the other stations is low, the response of the optimization based constrain-handling method is somewhat better than the back-calculation method. So, provided both algorithms are equivalent, it was decided to use the method with the best performance from a computational point of view. In this case, the back-calculation method is much better than the optimization based constrain-handling method, having that it is approximately three times faster. Therefore, for this research, the back-calculation scheme was used. Some performance indices obtained with back-calculation scheme were overshoot of 0.1526%, the rise time of 39.056s, and settling time of 67.805s.Figure 6
Performance comparison between the Constrain and the Back-Calculation methods.


## 3.1. Discrete Model Predictive Control Strategy

The function-based discrete time model predictive control (DMPC) formulation proposed by Wang [24] was applied to the LPG pipeline. The DMPC strategy considered is based on state space models and was implemented to control the downstream pressure of each pipeline section (CPL) and the liquid level of each storage tank (CL), by manipulating the pressure drop in pumps (MP) and control valves positions (Mv), respectively, while respecting the physical limitations of the actuators. A set of discrete time Laguerre function is used as an orthonormal basis to describe the optimum future control trajectory (that is, the difference of the manipulated variable) within a prediction horizon. The use of the Laguerre function allows reducing the computational burden associated with MPC. Below the DMPC strategy is briefly presented.


### 3.1.1. Augmented State Space Model

As mentioned above, at the current study the MPC strategy used is based on state space models. The general extended state space model in (12) is employed. Recall that we embedded integrators in the state space model for obtaining integral action.(10)Δxpk+1ek+1︸xk+1=Ap0CpAp1︸AeΔxpkek︸xk+BpCpBp︸BeΔukek=01︸CeΔxpkykwhere subindex p indicates process, Δxpk+1=xpk+1-xpkΔuk=uk-uk-1 and ek=csetk-ck. Moreover, the poles of the new matrix are the poles of the process and one pole located at 1. So, firstly it was confirmed the model was controllable and observable. Then, on the other hand, it was necessary to estimate the states, which are not measured; in this sense a Kalman filter was used [25], basically because, in the presence of noise and uncertainties, it produces state estimates that tend to be more precise. Kalman filter is commonly used for MPC vendor products to provide state estimation functionality [26]. Additionally, it is worth mentioning that the performance of the control strategies was evaluated by using the nonlinear pipeline model described by (1) and (2). By using the finite difference method, a spatial-discrete version of the nonlinear model [15, 27] was obtained, which was used in closed loop simulations. This way, Figure 4 shows the block diagram of the implemented strategy, G the process to be controlled, and GC the DMPC controller.Figure 4
A general block diagram of the implemented strategy.


### 3.1.2. DMPC Formulation

The procedure established at [24] is centered on an orthonormal-base set of discrete functions to approximate the future differences of the manipulated variable Δmki+i=ΔmkiΔmki+1⋯Δmki+Nc-1T, where NC indicates the control horizon. These functions are known as discrete Laguerre functions, and its z-transform is shown in (11).(11)Γkz=z-1-a1-az-1Γk-1z,Γ1z=1-a21-az-1The pole of the function a determines the velocity at which lk converges to zero. Then, Δm is approximated by N Laguerre functions, such as in (12), where ki represents the current sampling period, k are the future sampling periods, N is the number of terms used at the expansion, and cj are the coefficients of the process which are the function of future sampling periods.(12)Δmki+k=∑j=1NcjkiljkLaguerre functions can be expressed in the state space form in (13), which reduce computational burden. Al is a N×N matrix, which is a function of a.(13)Lk+1=AlLk,L0T=1-a21-aa2-a3⋯-1N-1aN-1TTherefore (12) can also be represented as a function of the Laguerre vector defined at (13) such as(14)Δmki+k=LkTη,η=c1c2⋯cNTThus, the prediction of future states and outputs in an instantm is defined by (15).(15)xki+m∣ki=Amxki+∑i=0m-1Am-i-1BLiTηyki+m∣ki=Cxki+m∣kiThen, the cost function to be minimized can be written as (16), where Q=CTC and RL>0. It is necessary to note J is a function of η instead of ΔM, so minimization must be done in function of η.(16)J=∑m=1Npxki+m∣kiTQxki+m∣ki+ηTRLηResults of the minimization operation are shown in (17)η=-Ω-1ΨxkiΩ=∑m=1NpϕmQϕmT+RLΨ=∑m=1NpϕmQAmϕmT=∑i=0m-1Am-i-1BLiTTherefore, the control law is written at (18).(18)Δmk=kDMPCx^k


### 3.1.3. Constraint Handling

To take into account physical actuators limitations, a constraint is imposed on the control signalΔmk, which is mathematically defined in (19)mmin≤mk≤mmax By using the Laguerre functions, it is possible to rewrite expressions in (19) as shown in (20)mmin≤∑i=0k-1LiTη+mk-1≤mmaxFurthermore, it is important to note that constrained control strategies require real-time optimization using quadratic programming [24], where the minimization function would be (21). This restricted optimization must be done at the entire prediction horizon, which increases the computational burden.(21)J=ηTΩη+2ηTΨxkiThus, at this research, a second option for constraining handling has been proposed which involves applying a back-calculation scheme, Figure 5.Figure 5
Constraint handling scheme using the antiwindup back-calculation method.Both methods were implemented, and simulations were run considering one change in the pressure set point of station 1. Results showed that both methods back-calculation and constraint handling based on quadratic optimization are very similar, as can be seen in Figure6. As it refers to the IAE value, it is observed that the back-calculation method is slightly better than the constrain-handling method for the setpoint change in station 1. However, even though the effect over this change over the other stations is low, the response of the optimization based constrain-handling method is somewhat better than the back-calculation method. So, provided both algorithms are equivalent, it was decided to use the method with the best performance from a computational point of view. In this case, the back-calculation method is much better than the optimization based constrain-handling method, having that it is approximately three times faster. Therefore, for this research, the back-calculation scheme was used. Some performance indices obtained with back-calculation scheme were overshoot of 0.1526%, the rise time of 39.056s, and settling time of 67.805s.Figure 6
Performance comparison between the Constrain and the Back-Calculation methods.


## 3.1.1. Augmented State Space Model

As mentioned above, at the current study the MPC strategy used is based on state space models. The general extended state space model in (12) is employed. Recall that we embedded integrators in the state space model for obtaining integral action.(10)Δxpk+1ek+1︸xk+1=Ap0CpAp1︸AeΔxpkek︸xk+BpCpBp︸BeΔukek=01︸CeΔxpkykwhere subindex p indicates process, Δxpk+1=xpk+1-xpkΔuk=uk-uk-1 and ek=csetk-ck. Moreover, the poles of the new matrix are the poles of the process and one pole located at 1. So, firstly it was confirmed the model was controllable and observable. Then, on the other hand, it was necessary to estimate the states, which are not measured; in this sense a Kalman filter was used [25], basically because, in the presence of noise and uncertainties, it produces state estimates that tend to be more precise. Kalman filter is commonly used for MPC vendor products to provide state estimation functionality [26]. Additionally, it is worth mentioning that the performance of the control strategies was evaluated by using the nonlinear pipeline model described by (1) and (2). By using the finite difference method, a spatial-discrete version of the nonlinear model [15, 27] was obtained, which was used in closed loop simulations. This way, Figure 4 shows the block diagram of the implemented strategy, G the process to be controlled, and GC the DMPC controller.Figure 4
A general block diagram of the implemented strategy.


## 3.1.2. DMPC Formulation

The procedure established at [24] is centered on an orthonormal-base set of discrete functions to approximate the future differences of the manipulated variable Δmki+i=ΔmkiΔmki+1⋯Δmki+Nc-1T, where NC indicates the control horizon. These functions are known as discrete Laguerre functions, and its z-transform is shown in (11).(11)Γkz=z-1-a1-az-1Γk-1z,Γ1z=1-a21-az-1The pole of the function a determines the velocity at which lk converges to zero. Then, Δm is approximated by N Laguerre functions, such as in (12), where ki represents the current sampling period, k are the future sampling periods, N is the number of terms used at the expansion, and cj are the coefficients of the process which are the function of future sampling periods.(12)Δmki+k=∑j=1NcjkiljkLaguerre functions can be expressed in the state space form in (13), which reduce computational burden. Al is a N×N matrix, which is a function of a.(13)Lk+1=AlLk,L0T=1-a21-aa2-a3⋯-1N-1aN-1TTherefore (12) can also be represented as a function of the Laguerre vector defined at (13) such as(14)Δmki+k=LkTη,η=c1c2⋯cNTThus, the prediction of future states and outputs in an instantm is defined by (15).(15)xki+m∣ki=Amxki+∑i=0m-1Am-i-1BLiTηyki+m∣ki=Cxki+m∣kiThen, the cost function to be minimized can be written as (16), where Q=CTC and RL>0. It is necessary to note J is a function of η instead of ΔM, so minimization must be done in function of η.(16)J=∑m=1Npxki+m∣kiTQxki+m∣ki+ηTRLηResults of the minimization operation are shown in (17)η=-Ω-1ΨxkiΩ=∑m=1NpϕmQϕmT+RLΨ=∑m=1NpϕmQAmϕmT=∑i=0m-1Am-i-1BLiTTherefore, the control law is written at (18).(18)Δmk=kDMPCx^k


## 3.1.3. Constraint Handling

To take into account physical actuators limitations, a constraint is imposed on the control signalΔmk, which is mathematically defined in (19)mmin≤mk≤mmax By using the Laguerre functions, it is possible to rewrite expressions in (19) as shown in (20)mmin≤∑i=0k-1LiTη+mk-1≤mmaxFurthermore, it is important to note that constrained control strategies require real-time optimization using quadratic programming [24], where the minimization function would be (21). This restricted optimization must be done at the entire prediction horizon, which increases the computational burden.(21)J=ηTΩη+2ηTΨxkiThus, at this research, a second option for constraining handling has been proposed which involves applying a back-calculation scheme, Figure 5.Figure 5
Constraint handling scheme using the antiwindup back-calculation method.Both methods were implemented, and simulations were run considering one change in the pressure set point of station 1. Results showed that both methods back-calculation and constraint handling based on quadratic optimization are very similar, as can be seen in Figure6. As it refers to the IAE value, it is observed that the back-calculation method is slightly better than the constrain-handling method for the setpoint change in station 1. However, even though the effect over this change over the other stations is low, the response of the optimization based constrain-handling method is somewhat better than the back-calculation method. So, provided both algorithms are equivalent, it was decided to use the method with the best performance from a computational point of view. In this case, the back-calculation method is much better than the optimization based constrain-handling method, having that it is approximately three times faster. Therefore, for this research, the back-calculation scheme was used. Some performance indices obtained with back-calculation scheme were overshoot of 0.1526%, the rise time of 39.056s, and settling time of 67.805s.Figure 6
Performance comparison between the Constrain and the Back-Calculation methods.


## 4. Results and Discussions

Three different control configurations (see Figure7) were implemented: (1) two decoupled control systems, a DMPC based control system for pressure, and an on-off control system for level, (2) two decoupled control systems for pressure and level based on DMPC strategy, and (3) one control system for pressure and level, in this configuration, the design model used is shown in (22). Table 3 shows the values of the controller's parameters used. The setpoint for level control was fixed in 7.89m for the three tanks. Operating process conditions evaluated were pressure set points increments of 10%, 10%, and 20% (at times 5, 200 and 500 seconds) forPL1,PL2, andPL3, respectively, an increment of 50% (at 1000 seconds) in the output flow rate in tank 3. The obtained performances are shown next.(22)PL1PL2PL3H1H2H3s=G11G12G13G14G15G16G21G22G23G24G25G26G31G32G33G34G35G36G41G42G43G44G45G46G51G52G53G54G55G56G61G62G63G64G65G66sMp1Mp2Mp3Mv1Mv2Mv3sTable 3
Pressure and level control parameters.
Configuration Pressure Control Parameters Level Control Parameters 1 a = 0.98 0.4 0.8     N=229  Np=70Q=CeTCe  RL=diag0.0010.10.1 A S-R flip-flops were used, and lower and higher levels were set to 7.88m and 7.9m, respectively. 2 LQR Control,Setpoint=7.89m  Q=0.01CeTCeRL=diag100100300  3 a = 0.4 0.4 0.4 0.4 0.4 0.98 N = 25 25 9 15 15 2   N p = 600   Q = C e T C e RL=diag0.10.10.1111Figure 7
Pipeline Section P&ID for (a) Configuration 1, (b) Configuration 2, and (c) Configuration 3.
(a) (b) (c)First response to be seen is the performance of the control strategy implemented under setpoint changes, so, in Figure8 one step change to each controlled pressure is shown. Thus, as it is appreciated that a setpoint change onPL1, let us call it PL1set, has a low impact on controlled pressuredPL2 andPL3. Referring to a step change at PL2set, it is shown that the effect onPL3 is larger, and its disturbance toPL1 is appreciable but smaller than change seen atPL3. On the other hand, a change inPL3 also has a low effect onPL2 and almost an imperceptible effect onPL1. Thus, even though the system is interacting, this interaction seems to be restricted to locally near the unit process. In other words, changes at station 1 have a significant effect on station 2 rather than 3, changes at station 2 impact stations 1 and 3 in a similar way and changes at station 3 affect mainly over station 2.Figure 8
Pressure Response, under setpoint changes, for Configurations 1, 2, and 3.In Figure8 it is seen that control strategies can maintain the process at its desired setpoint. The performance of configurations 1 and 2 is similar, as it was expected since the pressure controller used in both configurations is the same. Moreover, the interaction between stations is more prominent as they are closer. In what concerns configuration 3, its performance is close to that showed by the other configurations; however, it looks its response velocity is higher. As it refers to tanks levels, pressure changes do not have an impact on them; basically level in tanks remained constant under pressure setpoints changes.The strategies were also tested by producing a disturbance at the outlet flow rate at station 3. Thus, Figure9 shows how, for configuration 1, this disturbance impacts the pressure of the other stations. As it was expected response at station 3 is more intense and, as it is seen, previous stations variations are still evident but lower. For configurations 1 and 2, variations in tanks levels due to disturbances cannot be seen, and constant levels are obtained.Figure 9
Level and Pressure Response, under disturbances, for Configurations 1, 2, and 3.Table4 shows the IAE performance index obtained for each configuration, where it is definitively seen; configuration 3 gives the best results as it refers to IAE.Table 4
IAE performance index obtained.
Variable Configuration 1 Configuration 2 Configuration 3 PL 1 1634 203.1 154.7  PL 2 1589 182.9 154.3  PL 3 3181 344.3 288.7  H 1 5.527 × 105 4.579 41.9  H 2 5.528 × 105 10.54 68.2  H 3 5.528 × 105 393.6 106.7


## 5. Conclusions

Three control configurations were implemented for a liquid pipeline transportation system:(1) two decoupled control systems, an MPC based control system for pressure and an on-off control system for level, (2) two decoupled control systems for pressure and level based on MPC strategy, and (3) one MPC system for pressure and level. A linearized pipeline model was used as a design model; however, the performance of the control configurations was assessed with a nonlinear pipeline model. The addition of back-calculation as the antiwindup method allows reducing the control algorithms execution time significantly while maintaining the same performance as the achieved with the optimization based constraint handling method. In comparison with previous works [7], the strategies implemented show an enhanced performance. Lower IAE (0.1526% vs. 14.63%) and overshoot (177.3 vs. 8.4×105) were obtained. Higher rise time and settling time was obtained with the proposed control strategies, but it is due to dimension, length (10326m vs. 7m) and diameter (0.5m vs. 0.0254m), of the pipeline considered. Among the three control configurations presented at the current study, the last one (Centralized Pressure and Level Control) showed the best performance among the different test carried out in the face of setpoints changes and disturbances. The IAE performance index confirmed this result.



---

*Source: 4538632-2019-05-20.xml*