Article pubs.acs.org/IECR
Nonlinear Model Predictive Control Based on Wave Model of HighPurity Internal Thermally Coupled Distillation Columns Xinggao Liu,* Lin Cong, and Yexiang Zhou State Key Laboratory of Industrial Control Technology, Control Department, Zhejiang University, Hangzhou 310027, China ABSTRACT: Although the internal thermally coupled distillation column (ITCDIC) is much more energy-efficient than the conventional distillation column, the complex dynamics stemming from thermal coupling between the rectifying column and the stripping column have restrained the achievement of control strategies and application in chemical industry. To reduce the difficulties, a novel nonlinear wave model of ITCDIC process is proposed based on the wave theory in the separation processes to simplify the mechanism model of ITCDIC. A nonlinear model predictive control based on the proposed wave model (NWMPC) of the ITCDIC process is further established to handle the multivariable interactions and improve the control performance. A dynamic optimization scheme characterized by an improved control parametrization method is introduced to solve the problem generated by the NWMPC model. An IMC controller is also explored as a comparison research. In the example of a benzene−toluene system, NWMPC presents great advantage over IMC in the servo and regulatory control.
1. INTRODUCTION The internal thermally coupled distillation column (ITCDIC) has attracted much attention in the world during the last decades. Many research works have revealed its large energysaving potential.1−4 A few test applications of ITCDICs in the plant have been successfully carried out, which further confirm the energy-saving capability.5−19 However, due to the complex nonlinear behaviors and high degree of interactions between the column sections, the control design of ITCDICs becomes a troublesome problem, which also becomes a barrier on the application of ITCDICs in a real plant.20−22 Although a proportional integral derivative (PID) algorithm could get stable control of an ITCDIC,23,24 it could not deal with the strong interactions between the two columns, and it could not work efficiently under various disturbances such as the feed flow rate, feed component fraction, and so on in the high-purity control. As a result, a model-based control algorithm is much more popular. Nakaiwa et al.25 proposed two model-based controllers, internal model control (IMC) and generic model control (GMC). Zhu et al.26,27 investigated an IMC scheme for an ITCDIC, which is based on a second-order model. However, the establishment of an efficient and simple model especially a nonlinear model of ITCDICs is still a hard assignment. In the control research of conventional distillation columns, the linear approximation model is widely used for its simplicity and convenience. However, under high-purity conditions, a linear approximation model could not capture the nonlinear dynamics any more.14,28 Up to now, the controller design of ITCDICs is mainly based on the data-driven model. Naikawa25 and Zhu et al.26 both established the IMC controller for ITCDICs. Compared with PID control algorithms, IMC attains better control performance when set-points transfer. However IMC is sensitive to the change of the operating conditions, and its regulatory control performance becomes invalid. Although an IMC controller is improved by choosing different transfer function models, the control performance could be improved hardly because the model mismatch always emerges when the operating conditions change.27 An advanced nonlinear model is © 2013 American Chemical Society
required to describe the system, considering the complex dynamics of ITCDIC. Nonlinear wave theory is first developed in the fixed-bed sorption and chromatography.29−40 Systems with distributed parameters tend to exhibit dynamic phenomena that resemble traveling waves. Luyben41,42 proposed profile position control of distillation columns with the propagation of temperature profiles. Marquardt et al.43 developed nonlinear model reduction for distillation using nonlinear wave theory and derived expressions for the wave propagation velocity and the trial function of the wave profile. Hwang et al.44−46 studied the wave propagation dynamics in binary distillation columns, such as the propagation, reflection, superposition, and selfsharpening behaviors of concentration waves. Later Hankins47 proposed nonlinear wave model in a conventional distillation column by considering the influence of variable molar flows for dynamics and disturbance propagation. On the basis of the nonlinear wave model, lots of advanced controllers were developed for the conventional distillation columns, such as reactive distillation columns, multiseparation columns, and so on.48−51 A nonlinear wave model and an advanced nonlinear model predictive control (NMPC) strategy was proposed by Henson and co-workers52,53 for cryogenic distillation columns. An accurate low-order dynamic model was established based on nonlinear wave phenomenon, and the combination of the wave model and the MPC strategy formed an outstanding control strategy. The composition profile expression (the trial function) used in their work is also applicable in ITCDIC and proves to be very effective. The difference is that their derivation of the wave velocity involved the assumption of constant molar flow rate, which is reasonable in the cryogenic distillation columns Received: Revised: Accepted: Published: 6470
January 4, 2013 April 10, 2013 April 16, 2013 April 16, 2013 dx.doi.org/10.1021/ie400033h | Ind. Eng. Chem. Res. 2013, 52, 6470−6479
Industrial & Engineering Chemistry Research
Article
but does not hold true in ITCDIC. So the wave velocity needs to be rederived for ITCDIC. The conventional nonlinear wave model cannot be applied to ITCDICs directly. Liu et al.54 proposed a nonlinear wave model of an ITCDIC and implemented a GMC strategy for an ITCDIC based on the proposed model. Although the performance of the controller is very good, the method has some limitations. Usually products of different purities or grades are required in the same column, and considering the significant nonlinear behaviors and sensitivity to the operating conditions, the conventional control schemes must be tuned at each operating condition to achieve good control performance. The values of the parameters of controllers are troublesome problems and when the state of the operating conditions changes largely, the controllers may become invalid. Considering the complex dynamics of ITCDIC, it is also required to employ an advanced nonlinear multivariable control scheme to satisfy the control performance along with an advanced nonlinear model. Nonlinear model predictive control possesses a strong potentiality in improving control and operation of nonlinear processes. In this work, a simpler nonlinear wave model of ITCDIC is established compared to that proposed by Liu et al.54 The wave velocity with variable molar flow rate is derived first, which plays the most important role in the dynamic wave modeling. The further modeling is carried out based on the profile trial function43 and thermally coupled relations. The wave model is much simpler than the model proposed by Liu et al.54 and maintains the accuracy, which is more convenient for solving dynamic optimization in MPC. On the basis of the proposed wave model, a nonlinear wave model predictive controller (NWMPC) is developed for the ITCDIC, and an improved control parametrization method is introduced to solve the dynamic optimization problems to improve the performance of the NWMPC algorithm, which is mainly extended from a twopoint gradient and variational method. The benzene−toluene system is studied as an illustration, where the proposed NWMPC is compared with an IMC controller previously established by Zhu et al.26 In the servo and regulatory control, NWMPC exhibits substantial improvement on eliminating model mismatch, and also, NWMPC presents a great advantage in the restraint of various disturbances. What is more, NWMPC still works properly when the set-point of the top product changes to 0.9999, a very high purity.
Figure 1. Schematic diagram of an ideal ITCDIC.
in the rectifying section and stage j + f − 1 in the stripping column is the following: Q j = UA(Tj − Tj + f − 1),
(j = 1, ..., f − 1)
(1)
where UA is heat transfer rate and Tj is the temperature of stage j. According to Antoine equation, Raoult’s law, and Dalton’s law, Tj = b/(a − ln Pvp, j) − c
(2)
Pvp, j = PYj/(γjXj)
(3)
where γj is the activity coefficient at the jth stage, a, b, and c are the coefficients of the Antoine equation, Pvp,j is the vapor saturated pressure of stage j, P is the pressure of the rectifying column (Pr) or the pressure of the stripping column (Ps), and Xj is the mole fraction of liquid at stage j. In the benzene− toluene system in this work, eq 3 can be written as Pvp, j = P /[Xj + (1 − Xj)/α]
(4)
where α is the relative volatility. Liquid flow rates in the rectifying section are j
Lj =
∑ Q k /λ ,
(j = 1, ..., f − 1) (5)
k=1
where λ is latent heat. And, vapor flow rates in the rectifying section are
2. MECHANISM MODEL OF THE ITCDIC An ITCDIC has neither a reboiler nor a condenser, and the manipulation of internal thermal coupling is accomplished through heat exchange between the rectifying section and the stripping section. In order to provide the necessary thermal driving force for the heat transference from the rectifying section to the stripping section, the former must be operated at a higher pressure than the latter. Figure 1 shows the schematic diagram of an ideal ITCDIC. Because of the high degree of thermal coupling relations, the interactions seem to become intensified.26 The traditional datadriven model and linear approximation model have a great difficulty in the description of the distinct nonlinear dynamics. Before the development of nonlinear wave model, the basic equations of the ideal ITCDIC are shown as follows (The definitions of the variables are listed in the Nomenclature): The configuration of an ITCDIC in this work is a full symmetrical structure, and the heat transferred between stage j
Vj + 1 = V1 + Lj ,
(j = 1, ..., f − 1)
(6)
V1 is the vapor flow rate of the top product leaving stage 1, which can be derived from the total mass balance equation: V1 = F(1 − q)
(7)
where F is the feed flow rate and q is the feed thermal condition. Liquid flow rates in the stripping section are j
Lf + j − 1 = Lf − 1 + Fq −
∑ Q k /λ ,
(j = 1, ..., f − 2)
k=1
(8)
Ln = F − V1
(9)
The vapor flow rates in the stripping section are 6471
dx.doi.org/10.1021/ie400033h | Ind. Eng. Chem. Res. 2013, 52, 6470−6479
Industrial & Engineering Chemistry Research
Article
j
Vf + j = Vf − F(1 − q) −
algebraic stat vectors. The model contains 40 differential equations and 220 algebraic equations.
∑ Q k /λ , k=1
(j = 1, ..., f − 2)
3. NONLINEAR WAVE MODELING OF THE ITCDIC 3.1. Nonlinear Wave Velocity of the ITCDIC. Dynamics of separation processes could be described by the propagation of concentration or temperature profiles. In conventional distillation columns, a nonlinear wave is commonly defined as a spatial structure moving with constant propagation velocity and constant shape along a spatial coordinate. Wave velocity and the describing function of the profile are two important aspects of the description of the moving front. The derivation of wave velocity for ITCDIC is given as follows: The approximate expressions of the component mass balances (eqs 11−14) are expressed as follows:
(10)
Component mass balances: H dX1/dt = V2Y2 − V1Y1 − L1X1
(11)
H dXj /dt = Vj + 1Yj + 1 − VjYj + Lj − 1Xj − 1 − LjXj , (j = 2, ..., n − 1 and j ≠ f )
(12)
H dX f /dt = Vf + 1Y f + 1 − Vf Y f + Lf − 1X f − 1 − Lf X f + FZf (13)
H dX n/dt = −VnYn + Ln − 1X n − 1 − LnX n
(14)
where H is the stage holdup, Yj is mole fraction of vapor at stage j, and Zf is the mole fraction of feed. Vapor−liquid equilibrium relationships: Yj = KjXj
(j = 1, ..., n)
h
(j = 1, ..., n)
(16)
∫0
The mole fraction of vapor at stage 1 (Y1) and the mole fraction of liquid at stage n (Xn) are the controlled variables (the outputs). Normally, the temperature is measured, and according to (2), (3), and (14), Y1 and Xn are derived. The pressure of rectifying column Pr and the feed thermal condition q are the operation variables in the control strategy. The feed flow rate F, the feed composition Zf, and the pressure of the stripping column Ps are the main disturbances in this work. So the values of the above variables need to be measured. Some of the states maintain fixed values in this work (n, f, a, b, c, λ, α, H) which are listed in Table 1, and other states can be derived from
40
0.501
pressure of rectifying section (Pr) [Mpa] pressure of stripping section (Ps) [Mpa] heat transfer rate (UA) [W/K] latent heat of vaporization (λ) [J/mol] liquid holdup (H) [kmol]
feed stage ( f)
21
feed flow rate (F) [kmol/h] feed composition (Zf) feed thermal condition (q) relative volatility (α) antoine constant (b)
100
1.5
2.317 2788.51
Antoine constant (a) antoine constant (c)
15.9008 −52.36
0.5
f −1
h
∂X dz = ∂t
n
∫0
f − 1 ⎛ ∂(VY )
⎜ ⎝ ∂z
⎛
n
∫f −1 h ∂∂Xt dz = ∫f −1 ⎜⎝ ∂(∂VYz ) −
−
∂(LX ) ⎞ ⎟ dz ∂z ⎠
∂(LX ) ⎞ ⎟ dz ∂z ⎠
(18)
(19)
Wave velocity describes how fast the wavefront travels. However it is difficult to describe the exact position of the wavefront; usually the sharpest point on the profile is chosen as a representation of the front. Denote wavefront positions in the rectifying and stripping sections by S1 and S2, respectively. Define ε1(t) = z − S1(t) and ε2(t) = z − S2(t), and transfer eqs 18 and 19 into wave traveling coordinate systems which are centered at the wave position (ε1 = 0 and ε2 = 0). After the algebraic transformations, the wave velocity of the ITCDIC is derived as follows:
Table 1. Operating Conditions of the ITCDIC Studied stage number (n)
(17)
where z = z/Δz, Δz is the height equivalent to a theoretical ̅ plate (HETP), and z̅ is a spatial coordinate. Then z becomes a dimensionless spatial coordinate belonging to [0 n], h = H/Δz. Integrate eq 17 from the top to the feed tray, and the feed tray to the bottom, respectively:
(15)
where Kj is the equilibrium vaporization ratio at the jth stage. In the benzene−toluene system, eq 15 can be written as Yj = αXj /[(α − 1)Xj + 1],
∂(VY ) ∂(LX ) ∂X = − ∂t ∂z ∂z
0.3387 0.101 9803
∫0
30001.1
f −1
h
∂X̃ dε1 ∂t 1 f − 1 − S1 ∂X̃ ∂ε1 = h dε1 −S1 ∂ε1 ∂t
∂X dz = ∂t
f − 1 − S1
∫−S
h
∫
(20)
where X̃ denotes the new liquid mole fraction in the moving wave coordinates. Note that ε1(t) = z − S1(t), so dε1/dt = −dS1/dt, and by substitution into eq 20, we can get the following expression:
the solution of the dynamic model and do not need to be measured. Normally, the value of the heat transfer rate UA is constant and does not need to be measured, but in this work, we also regard it as a disturbance like F, Zf, and Ps to test and verify the performance of the controller. All the system state vectors can be divided into two classes, the differential state vectors and the algebraic state vectors. The corresponding model contains both differential equations and algebraic equations. In this work, the stage number n is 40, and from the mechanism model above, we can see that there are 40 components of differential state vectors and 220 components of
∫0
f −1
h
f − 1 − S1
∂X dz = − ∂t
∫−S
=−
dS1 dt
=−
1
h
∂X̃ dS1 dε1 ∂ε1 dt
f − 1 − S1
∫−S
h
1
∂X̃ dε1 ∂ε1
dS1 h(X f − 1 − X1) dt
(21)
From eqs 17 and 11−13, we can obtain 6472
dx.doi.org/10.1021/ie400033h | Ind. Eng. Chem. Res. 2013, 52, 6470−6479
Industrial & Engineering Chemistry Research
∫0
f −1
h
∂X dz = ∂t
∫0
f − 1 ⎛ ∂(VY )
⎜ ⎝ ∂z
−
Article
∂(LX ) ⎞ ⎟ dz ∂z ⎠
= Vf Y f − V1Y1 − Lf − 1X f − 1
Combined with the material balances and the thermally coupled relations, the nonlinear wave model of the ITCDIC is established, including the thermally coupled relations (eqs 1−4), vapor and liquid molar flow rates (eqs 5 and 9 when j = f − 1, eqs 7 and 8), vapor−liquid equilibrium relationships (eq 15), the wave velocity (eqs 23 and 24), and the trial function of the profile (eqs 25−26). From the nonlinear wave model, we can see that there are only 2 components of differential state vectors and 184 components of algebraic state vectors, and the model contains only 2 differential equations and 184 algebraic equations, much less than the mechanism model, which contains 40 differential state vectors, 220 algebraic state vectors and 40 differential equations, 220 algebraic equations. The reduced wave model changes the large-scale problem to a smaller one, which makes the solution of the dynamic optimization problem used in the model predictive control much easier and more efficient. 3.3. Model Test. Dynamics of the ITCDIC are mainly characterized by the changes of the mole fractions of different stages when operation conditions change. Therefore, the mole fractions of different stages need to be tracked efficiently. Also, the movement of the sharpest point on the profile should be tracked efficiently; that is, the wave velocity by the wave model should be correspondent with the velocity derived by the observation of the real concentration profile observed from the mechanistic model. In this work, the benzene−toluene system is taken as an illustration. The system with the operating conditions presented in Table 1 is disturbed by a +10% step change of the feed flow F. Figure 2 shows the tracking of the wave velocity. From the figure, we can see that the wave velocities got by the wave
(22)
Combine eqs 21 and 22, the wave velocity of the rectifying section is derived as follows: Vf Y f − Lf − 1X f − 1 − V1Y1 dS1 = dt h(X1 − X f − 1)
(23)
In the same way, the wave velocity of the stripping section is the following: −Vf Y f − LnX n + Lf − 1X f − 1 + FZf dS 2 = dt h(X f − 1 − X n)
(24)
3.2. Nonlinear Wave Model of the ITCDIC. Due to the highly distorted concentration distributions, the profile of the ITCDIC presents two obvious characteristics: (1) asymptotic properties at the two ends; (2) only one inflection point exists both in the rectifying section and stripping section. Hence an elegant trial function derived by Marquardt43 is employed to describe the shape of the front as follows: X̂ (j) = X r _min +
X r _max − X r _min 1 + exp( −k r(j − S1))
j = 1, 2, ..., f −1 X̂ (j) = Xs _min +
(25)
Xs _max − Xs _min 1 + exp( −ks(j − S2))
j = f , f + 1, ..., n
(26)
where Xr_max, Xr_min, kr, Xs_max, Xs_min, and ks are the six profile parameters which are generated from the composition concentration data. Xr_min and Xr_max denote the asymptotic limits of the rectifying section when the profile extends to an infinite distance; Xs_min and Xs_max denote the asymptotic limits in the stripping section. Xr_max and Xs_min approximate to the liquid mole fraction at the top stage and the bottom stage respectively. The terms kr and ks characterize the tangent of the inflection points S1 and S2 (kr and ks do not equal the tangents). Nonlinear least-squares estimation is employed to update the six parameters online formulated as follows (take the rectifying section as an example): f −1
min
∑ (X̂i − Xi)2 i=1
Figure 2. Tracking of wave velocity (left is the rectifying section and right is the stripping section).
X1̂ − X1 = 0 X̂ f − 1 − X f − 1 = 0
model and by the profile observation are quite closer to each other. So we give the difference of the two velocities in the figure. At the very initial time when the feed flow changes, there is a small error existing for a very short duration, and after 1 h, the error comes close to zero, which means that the tracking of the wave velocity is very efficient. Figure 3 shows the comparison of liquid mole fraction at different stages. The average errors between the two tracking curves at different stages are 8.07 × 10−4, 1.35 × 10−3, 2.55 × 10−3, 3.13 × 10−3, 4.82 × 10−3, 2.48 × 10−3, 1.95 × 10−3, and 2.37 × 10−4, respectively. The error is smaller as the stage gets
X r_min − X r_max < 0 0 < X r_max
