1. business and industrial
2. energy
3. oil

Systematic approach in tuning PID controller for a
high purity control of Binary Distillation column
based on reduced order model: A Case study
Vetriselvi V1, N. Pappa2, Anuj Abraham3
Department of Instrumentation Engineering
MIT Campus, Anna University
Chennai, India
1
vetriselvi.mv@gmail.com, 2npappa@rediffmail.com, 3anuj1986aei@gmail.com
Abstract— This paper describes a systematic way for tuning
PID control structure for a Binary Distillation column process.
An approach to various Model Order Reduction (MOR)
techniques is compared with the original distillation column
process. The modeling and simulation of nonlinear model of
multivariable is implemented in MATLAB and a comparative
study is analyzed for full order linearised model using Taylor’s
Series expansion, Jacobian linearization method and reduced
order linear model namely Balanced Truncation, Singular
perturbation, Hankel Norm approximation methods. Also
controller performance evaluated for conventional PID
controller for simulated model. The result shows consistently
good servo and regulatory response.
distillate flow rate, bottom flow rate, cooling water flow rate
and control variables are concentration of distillate and bottom,
level of condenser and reboiler.
Index Terms—Model Order Reduction, Binary Distillation
Column, Balanced Truncation, Singular perturbation, Hankel
Norm.
I. INTRODUCTION
Distillation column plays a very important role in chemical
processes mainly petro-chemical and refinery industries, for an
effective way to separate mixtures of liquid or vapor mixture
of two or more substances into its component fractions of
desired purity[16].
Binary distillation column is used to separate two
components. The schematic diagram of a general distillation
column is shown in Fig. 1. Distillation column consists of a
vertical column, where plates or trays are used to increase
the component separations. Distillation column is separated
into two sections. They are stripping section and rectification
section. The trays above the feed tray is called stripping
section. The trays below the feed tray are called rectification
section. Reboiler and condenser are used as heat duties.
Condenser is used to condense distillate vapor and reboiler is
used to provide heat for the necessary vaporization from the
bottom of the column. Condensed vapor is collected in reflux
drum and require amount of it is used as a reflux.
The L-V (Liquid-Vapor) structure [9] is known as the
energy balance structure and can be considered as the
standard control structure for a dual composition control
distillation. The various manipulated variables considered for
distillation process in literature are reflux rate, boilup rate,
Fig. 1. Schematic diagram of distillation column
In this paper the three most important manipulated
variables taken are reflux rate, boilup rate. Also, the control
variables chosen are concentration of distillate and bottom. The
disturbances associated with distillation column are feed flow
rate and concentration of feed. The vertical column is designed
for 14 trays and the list of process parameters considered is
shown in Fig. 2.
Tray n (n=2 to 7):
Mx n = V(y n-1 - y n ) + (L + L F )(x n+1 - x n )
(6)
Re-boiler (n=1):
M B x1 = (L + L F )x 2 - Vy1 - Bx1
(7)
C. Process data for distillation column under normal
operating condition
The process data are based on a real petroleum project
Petro Vietnam Gas Company [1]. The input feed consist of
LPG and Naphtha. The relative volatility α under operating
conditions is 5.68. The properties and variations of the feed
includes molar weight, liquid density, feed composition of feed
under operating conditions are listed in Table I.
TABLE I. PROPERTIES OF FEED
Properties
Fig. 2. Process parameters of Distillation Column
II. MATHEMATICAL MODELING
A. Assumptions
The various assumptions considered for the distillation
column modeling are given below: [1]
a) The relative volatility α is constant throughout the
column. This means the vapor liquid equilibrium
relationship can be expressed by,
αx n
yn =
(1)
1 + (α -1)x n
th
b)
c)
d)
e)
f)
xn is liquid composition on n stage,
yn is vapor composition on nth stage,
α is the relative volatility.
The overhead vapor is totally condensed in a condenser.
The liquid holdups on each tray, condenser, and the reboiler are constant and perfectly mixed
The holdup of vapor is negligible throughout the system.
The molar flow rates of the vapor and liquid through the
stripping and rectifying Sections are constant.
The column is numbered from bottom (n=1 for the reboiler, n=2 for the first tray, n=f for the feed tray, n=N+1
for the top tray and n=N+2 for the condenser)
B. Dynamic model of distillation column process
Based on the assumptions described in section 2.1, the
dynamic models of distillation process are expressed by the
following component material balance equations:
Condenser (n=16):
M D x n = (V + VF ) y n-1 - Lx n - Dx n
(2)
Tray n (n=10 to 15):
Mx n = (V + VF )(y n-1 - y n ) + L(x n+1 - x n )
(3)
Tray above the feed flow (n=9):
Mx n = V(y n-1 - y n ) + L(x n+1 - x n ) + VF (y F - y n )
Tray
below
the
feed
flow
Mx n = V(y n-1 - y n ) + L(x n+1 - x n ) + L F (x F - x n )
(4)
(n=8):
(5)
LPG
Naphtha
Molar weight
54.4-55.6
84.1-86.3
Liquid density (kg/m3)
570-575
725-735
Feed composition (vol %)
38-42
58-62
TABLE II. OPERATING CONDITIONS OF DISTILLATION COLUMN PROCESS
Stream
Feed
LPG
Naphtha
118
46
144
Pressure (atm)
4.6
4
4.6
Density (kg/m3)
670
585
727
Volume flow rate (m3/h)
22.76
8.78
21.88
Mass flow rate (kg/h)
15480
5061
10405
Plant capacity (ton/year)
130000
43000
87000
Temperature (oC)
TABLE III. NOMINAL
VALUES
DISTILLATION COLUMN PROCESS
Variable
FOR
Stream
PROCESS
PARAMETERS
Molar flow
Unit
31.11
Kmole
5.8
Kmole
13.07
Kmole
VF
Liquid holdup in the
column base
Liquid holdup on a
tray
Liquid holdup in the
reflux drum
Vapor rate in feed
98.5152
kmole/h
LF
Liquid rate in feed
104.2491
kmole/h
V
Internal vapor rate
66.3407
kmole/h
MB
M
MD
L
Internal liquid rate
75.638
kmole/h
D
Distillate flow rate
92.7597
kmole/h
B
Bottoms flow rate
110.9235
kmole/h
OF
The operating conditions of distillation column process
and nominal values of process parameters are summarized in
Table II and Table III respectively.
D. Linearization of nonlinear differential equation
a) Taylor’s Series Expansion:
Multivariable binary distillation column is nonlinear
model, which are linearized to perform a simulation and
stability analysis. The nonlinear part of the model equations is
yn. It can be approximated using Taylor series approximation
around the steady-state operating point ( x n =x n )
f(x n )=f(x n )+f (x n )(x n -x n )+
f  (x n )
2
(x n -x n ) +... (8)
2!
For x sufficiently close to x , higher order terms will be very
close to zero, so we neglect the quadratic and higher order
terms to obtain the approximation[1].
y n  f(x n )+K n (x n -x n )
where, f(x n )=
αx n
1+(α-1)x n
, Kn =
(9)
A
( A, B, C )    11
  A21

By truncating the discard able states, the truncated reduced
system is then given by (A11 ,B1 ,C1 )
B. Singular perturbation (SP)
The full-order linear model which represents a two inputstwo outputs plant in equation can be expressed as a reduced
order linear model as,
1
x D 
L
=
G(0)
 x B  1+τ s
V
c
α
(1+(α-1)x n )
A12   B1 

,   , C1 C2  

A22   B2 
2
(14)
where, G(0) is the steady state gain defined by,
b) Jacobian Linearization Process:
-1
G(0)=-CA B
Vector form of nonlinear model
x=f(x,u)
(10)
z=g(x,u)
Nonlinear model of distillation process has 16 state
variables (x1 ,x 2 ,...,x16 ) , 2 input variables (L,V) , and 2 output
τc =
MI
+
Is lnS
M D (1-x D )x D
+
M B (1-x B )x B
Is
(15)
Is
where, τ c is the time constant,
M I is the total holdup of liquid inside the column defined by,
variables (x D ,x B ) .
14
M I =  Mi
x1 =f(x1 ,...,x16 ,L,V)
(16)
i=1
I S is the Impurity sum defined by,
x16 =f(x1 ,...,x16 ,L,V)
IS =D(1-x D )x D +B(1-x B )x B
(11)
z1 =x D
z 2 =x B
State space form of linearized model
x=Ax+Bu
z=Cx+Du
Elements of the linearization matrices
A ij 
f i
x j
g
Cij  i
x j
Bij 
x ,u
x ,u
f i
u j
g
Dij  i
u j
S is the Separation factor defined by,
(1-x B )x D
S=
(1-x D )x B
(17)
(18)
C. Hankel norm approximation
(12)
The hankel norm of a system is defined, as in Eq.19,
G=  A,B,C,D   H 

x ,u
(13)
G
x ,u
2
 y (t) dt
2
H
 sup 0
2
 u (t) dt
(19)
0
III. MODEL ORDER REDUCTION (MOR)
The requirement of model order reduction is that the
reduced order model so obtained should retain the important
and key qualitative and quantitative properties such as
stability, transient and steady state response etc. of the original
system[10].
A. Balanced Truncation (BT)
An arbitrary system (A,B,C,D) can be transformed into a
balanced system (A,B,C,D) via a state-space transformation.
The balanced systems states are ordered (descending) by how
controllable and observable they are, thus allowing a portion
of the form:
0
where, y(t)   C e

A(t-s)
B u(s)
Hankel norm gives how much energy can be transferred
from past inputs into future outputs through the system G. In
control theory, eigenvalues define system stability and Hankel
singular values define the "energy" of each state in the system.
Its characteristics in terms of stability, frequency, and time
responses are preserved by keeping larger energy states of a
system based on the Hankel singular values. They can achieve
a reduced-order model that preserves the majority of the
system characteristics. Mathematically, given a stable statespace system (A,B,C,D) its Hankel singular values are defined
as in Eq.20,
G
2
H
= λ max (PQ)=σ i
(20)
A. Decoupling control system
As shown in Fig.4, to compensate for process interactions
additional controllers are used and reduces control loop
interaction[3].
where, σi is the hankel singular values,
The controllability and observability grammians P and Q
respectively satisfies,
T
T
AP+PA =-BB
T
(21)
T
(22)
A Q+QA=-C C
One defines the Hankel operator ΓG of the system G(s) by,
0
ΓG :L2 (, 0] : (ΓG u(t))   C e
A(t-s)

B u(s) ds, t>0 (23)
This method also guarantees an error bound on the infinity
norm of the additive error G-G red  for well-conditioned
Fig. 4. Example of a figure caption. (figure caption)
model reduced problems as in balanced truncation method.
n
G-G red

 2  σi
(24)
K+1
where, σ i are singular values of a given system G(s).
IV. CONTROLLER DESIGN OF BINARY DISTILLATION COLUMN
Figure 3 shows the composition control diagram of binary
distillation column. In this control configuration, the vapour
flow rate V and the liquid flow rate L are the control inputs to
maintain the specification of the product concentration outputs
XB and XD (controlled variable) due to disturbance F (feed
flow) and XF (feed concentration).
B. Decoupler design equations
For multi-loop control, cross-controller T12 is used to
cancel the effect of U2 on Y1.
T12 GP11 U22 + GP12 U22 = 0
Because U22≠0 in general, then
G
T12 =- P12
(25)
G P11
4
T12 =
3
2
0.0053s +0.2816s +3.8374s +17.5437s+6.2479
4
3
4
3
2
0.0024s +0.1298s +1.558s +7.6499+4.152
Similarly, cross controller T 12 is used to cancel the effect of U1
on Y2.
T21 GP22 U11 + GP21 U11 = 0
G
T21 =- P21
(26)
G P22
T21 =
2
0.0735s +2.1553s +16.4563s +22.4s+4.932
4
3
2
0.0753s +2.2082s +17.0774s 25.1257s+7.4866
V. RESULTS AND DISCUSSION
A. Quality of top and bottom product
Fig. 5 and Fig. 6 shows the liquid and vapor concentration of
top product on each tray for nominal operating conditions of
the binary distillation column.
Fig. 3. Composition control of binary distillation column
Fig. 5 Response of Liquid concentration of LPG on each
tray
Fig.8 Product quality for nominal and ±10% change
Reflux rate
When reflux rate is decreased by 10% from its nominal value,
the quality of the distillate product increases by 2.5% and the
quality of the bottoms product increases by 2.1%. In contrast,
when reflux rate is increased by 10% from its nominal value,
the quality of the distillate product decreases by 4.3% the
quality of the bottoms product decreases by 3.6%.
8.4 Simulation with 10% decreasing and increasing Boilup
rate:
Fig. 6 Response of Vapor concentration of LPG on each
tray
From the simulation result the Purity of top product is
obtained as 96.45% and Impurity of Bottom product as 3.13%.
Simulation with the nominal values of stream, the purity of the
distillate product is 96.45% and the impurity of the bottoms
product is 3.13%. Fig. 9 shows the product quality for nominal
and ±10% change in boilup rate to the binary distillation
column.
8.2 Simulation with 10% decreasing and increasing feed flow
rate:
Fig. 7 shows the product quality for nominal and ±10%
change in feed flow rate to the binary distillation column.
Fig. 7 Product quality for nominal and ±10% change in
feed rate
When feed flow rate decreased by 10% from its nominal
value, the quality of the distillate product decreases by 6.9%
the quality of the bottoms product increases by 2.7%. In
contrast, when feed flow rate increased by 10% from its
nominal value, the quality of the distillate product increases by
9% and the quality of the bottoms product decreases by 8.6%.
8.3 Simulation with 10% decreasing and increasing Reflux
rate:
The product quality for nominal and ±10% change in reflux
rate to the binary distillation column is shown in Fig.8.
Fig. 9 Product quality for nominal and ±10% change
Boilup rate
When boilup rate is decreased by 10% from its nominal value,
the quality of the distillate product decreases by 4.8% and the
quality of the bottom product decreases by 4.01%. In contrast,
When boilup rate is increased by 10% from its nominal value,
the quality of the distillate product increases by 3.5% and the
quality of the bottoms product increases by 2.9%.
8.5 Response of nonlinear and linearized model
Based on the linearized model the simulation were carried out
for nominal values. Response of Top product liquid
concentrations of linearized model and simulated binary
distillation column process are shown in Fig. 10
8.8 Reduced order linearised model (Gh) based on Hankel
norm
The Hankel singular values in decreasing order are
3
3
4
5.2172  10 , 2.606  103 , 1.038  10 , 1.4281  10 ,
5
3.7081  10 and all remaining singular values are smaller
than 106. It is observed from the Fig. 11, that the Hankel
singular values that the 5th order of the system captures the
majority of the input–output behavior of the system.
Fig. 10 Response of nonlinear and linearised model
The linearized model response tracks the simulated binary
distillation, column process at steady state condition with
small deviation in transient condition is shown in Fig. 9.
Table 5 Steady state values of xn, yn, Kn on each tray
Trays
16
15
14
13
12
11
10
9
xn
0.9645
0.9216
0.8316
0.6824
0.5110
0.3804
0.3087
0.2763
yn
0.9936
0.9853
0.9656
0.9243
0.8558
0.7772
0.7173
0.6844
Kn
0.1868
0.2012
0.2373
0.3230
0.4938
0.7347
0.9502
1.0803
Trays
8
7
6
5
4
3
2
1
xn
0.2657
0.2617
0.2525
0.2325
0.1944
0.1373
0.0765
0.0313
yn
0.6727
0.6682
0.6574
0.6325
0.5781
0.4747
0.3199
0.1551
Kn
1.1284
1.1473
1.1931
1.3027
1.5576
2.1054
3.0806
4.3214
8.6 Reduced order linearised model (Gb) based on BT in state
space:
The linearised fifth order model (Gb) based on BT in state
space is obtained as,
5th order:
1.456
3.084
0.155 
8.006 3.227
 6.353 4.165
1.923
2.887 1.115 


A  0.988
0.915
0.792
1.903
0.882


0.944
12.23 0.5977 
 3.185 3.061
 0.9387 0.5377 0.6563 7.55 16.94 
 0.1918 0.2162 
 0.0567 0.2834 
 0.1192 0.0864 
 0.0669 0.1311 




B  0.0403
0.0049

 C  0.0211 0.0346 

0.0367
0.0463
0.00315
0.059




 0.0121 0.0333 
 0.0352 0.0043 
D
T
Fig.11 Hankel Singular Values
5th order:
20.69 4.391 4.548 1.85 1.562
 0
12.62 1.194 1.375 1.621


A
0
6.285 1.983 1.667
 0

0
0
1.378 0.397 
 0
 0
0
0
0
0.397 
 0.0568 0.9297 
 0.1871 0.1725
0 0 
D


0 0
B  0.113
0.1035


0.03752 0.09185 
 0.0309 0.00242 
C
0.03884 0.009042 0.007691 0.08733 0.006092 
0.09442 0.2501 0.1421
0.0944 0.03215
8.9 Comparative study of reduced order model
Based on the reduced order model the simulation is carried out
for nominal values.Response of linearized model and reduced
order model for top and bottom product is shown in Fig.12,
Fig.13
0 0 
0 0
8.7 Reduced order linearised model (Gs) based on SP in state
space:
The linearised second order linearised model (Gs) based on SP
in state space:is obtained as,
2nd order:
A
0 
0 
0.4639
0.0625
B
 0
 0
0.4639
0.0625
C
 0.0289 0.0445
0 0 
0.0349 0.0527  D  0 0
Fig.12 Comparison of linearized model and reduced order
models for top product quality
Fig.13
Comparison of linearized model and reduced order models
for bottom product quality
Based on the response of the reduced order model based on
Hankel norm approximation captures the majority of the input
output behavior of the system.
Fig. 14
Servo response of the simulated model based PID
controller for top product quality
From the controller responses, it is observed that the PID
controller is able to track the set point changes with smaller
overshoot and lesser settling time.
Table 6 ISE, IAE and MSE for original and reduced order
systems
Top
Bottom
Product
Order
5th order HN
2nd order SP
ISE
2.5*10-3
2.54
IAE
2.02
41.45
MSE
8.3*10-8
8.4*10-5
5th order BT
1.31
48.02
4.3*10-5
5th order HN
2nd order SP
5th order BT
2.02*10-4
5.11*10-1
8.10*10-1
0.61
18.13
383.18
6.7*10-9
1.7*10-5
2.6*10-3
From the Table 6 it is observed that the reduced order model
obtained using Hankel norm has minimum ISE, IAE and MSE
values for both top and bottom product quality.The transfer
function of the reduced order model is given below.
G11 
0.0024 s 4  0.1298s 3  1.558s 2  7.6499  4.152
s  42.1s 4  579.8s 3  2853.3s 2  3697.9s  1049.3
G21 
0.0735s 4  2.1553s 3  16.4563s 2  22.4 s  4.932
s 5  42.1s 4  579.8s 3  2853.3s 2  3697.9s  1049.3
G12 
0.0053s 4  0.2816 s 3  3.8374s 2  17.5437 s  6.2479
s 5  42.1s 4  579.8s 3  2853.3s 2  3697.9s  1049.3
5
Fig. 15
Servo response of the simulated model based PID
controller for bottom product quality
The simulation experiment was carried by increasing feed
flow rate by 3% from its nominal value. The responsed
obtained for the decoupled system with simulated model are
shown in Fig.16 and Fig.17 respectively.
0.0753s 4  2.2082 s 3  17.0774 s 2 25.1257 s  7.4866
8. G22  s 5  42.1s 4  579.8s3  2853.3s 2  3697.9s  1049.3
10 Analysis of controller performance
Fig.16 Regulatory response of the simulated model PID
controller for top product quality.
The simulation was carried out for the setpoint change of
Purity of top product from 0.964 to 0.99 and Impurity of
bottom product from 0.031 to 0.01. Both the controllers are
tuned based on ZN closed loop method. Servo response of the
simulated model for purity of top product and impurity of
bottom product based PID controller for binary distillation
column process are shown in Fig.14 and Fig.15
respectively[3].
Fig.17 Regulatory response of the simulated model based
PID controller for bottom product quality.
It is observed that the PID controller is able to reject the
disturbance with smaller overshoot and lesser settling time.
9. Conclusions
The first principle model of binary distillation column is
developed using governing equations and parameter values.
The simulated distillation column is validated under nominal
and steady state operating conditions. A linearized model of
order 16 is obtained using Taylor’s Series expansion, Jacobian
linearization Process. Three different model order reduction
techniques
namely
Balanced
Truncation,
Singular
Perturbation, Hankel Norm approximation are obtained and it
is observed that a 5th order reduced model obtained using
Hankel Norm captures the majority of behavior of the system.
PID controller was designed using simulated model. Designed
controller is able to provide good servo and regulatory
response.
Acknowledgement
The authors gratefully acknowledge Anna University, Chennai
for providing financial support to carry out this research work
under Anna Centenary Research Fellowship (ACRF) scheme.
References
[1] Vu TrieuMinhand John Pumwa “Modeling and Adaptive
Control Simulation For A Distillation Column’’ 14th
International Conference On Modelling And Simulation,
2012.
[2] Pradeep B. Deshpande, Charles A. Plank, “Distillation
Dynamics and Control”, Instrument Society of America,
1985.
[3] Francisco Vázquez Fernando Morilla, “Tuning
Decentralized Pid Controllers For Mimo Systems With
Decouplers” 15th Triennial World Congress, Barcelona,
Spain, 2002.
[4] Franks, R. G.E., Modeling and Simulation In Chemical
Engineering. Wileyinterscience, N.Y., (1972).
[5] Nelson, W. L., Petroleum Refinery Engineering.
Auckland Mcgraw-Hill, (1982).
[6] Joshi, M. V., Process Equipment Design. New Delhi,
Macmillan Company of India, (1979).
[7] Mccabe, W. L. , and Smith J. C., Unit Operations of
Chemical Engineering. N.Y Mcgraw-Hill, (1976).
[8] Wuithier, P., Le PetroleRaffinageEt Genie Chimique.
Paris Publications De L’institutFrancaise Du Petrole,
(1972).
[9] Stephanopoulos G., Chemical Process Control. Prentice
Hall International, (1984).
[10] Juergen Hahn, Thomas F. Edgar “An Improved Method
For Nonlinear Model Reduction Using Balancing of
Empirical Gramians” Computers and Chemical
Engineering 26 (2002) 1379– 1398.