Nonlinear Black-box Modeling of a Reactive Distillation Process

Nonlinear Black-box Modeling of a Reactive Distillation Process

Abdulwahab GIWA and SüleymanKARACAN

Ankara University, Graduate School of Natural and Applied Sciences, Department of Chemical Engineering, Ankara, Turkey

Abstract

Two types of nonlinear black-box (treepartition and sigmoid network NARX) modelswere developed for thereactive distillation processused for the production of ethyl acetate, water being the by- product,from the esterification reaction between acetic acid and ethanol. The data used for the modeling were generated by operating the reactive distillation column setup. The model orders were determined using MDL and AIC criteriaand compared. The unit numbers of the nonlinear estimators were optimized and the optimized values were further used to develop the models. The good conformities between the experimental top segment temperatures and NARX models predicted onesrevealed that the models successfully represented the behavior of the ethyl acetate reactive distillation process. In addition, owing to the higher fit value and the lower loss function observed from the developed sigmoid network NARX model, it was found to be better than treepartition NARX model for this process.

Keywords: Reactive distillation process,NARX model,treepartition, sigmoid network, ethyl acetate.

1. INTRODUCTION

   Esters are of great importance to chemical process industries. Among them, acetate esters are important organic solvents widely used in the production of varnishes, ink, synthetic resins, and adhesive agents. They are produced from the reactions of acid and alcohols under an acidic condition. A key issue in the production of these esters is the low conversion from the reactions. As a result, heavy capital investments and high energy costs are inevitable. The reactive distillation is a very attractive way to reduce these investments and energy costs [8].

   Reactive distillation is a process that combines both separation and chemical reaction into one unit. It has a lot of advantages for those reactions occurring at temperatures and pressures suitable for the distillation of the resulting components. Its main advantages are derived from the elimination of equipment and the constant removal of products [12]because it is known that increased overall conversion can be achieved in equilibrium reactions if the products are continuously removed from the reaction zone.

   The combination of both chemical reaction and separation in a single unit results in the complexity of this process which makes its modeling a very challenging task to chemical engineers because the unavailability of process models capable of reliably describing the several complexes and interrelated phenomena including simultaneous reaction and separation has made the safe scale-up from laboratory experiments to industrial plants quite challenging[7], especially when a packed column is the type involved in the reactive distillation.

   There are researches already carried out on the modeling of reactive distillation column using the first principles approach that normally incorporates many assumptions; most of them focused on tray columns. Those concerning packed column are few. In fact, very few researches were found regarding using plant data to develop models for a reactive distillation plant. From the information obtained from the literature, Giwa and Karacan (2012c) developed ARX and ARMAX models for ethyl acetate

   reactive packed distillation process. Even though they were able to obtain good models that were capable of representing the process very well, the models they developed were linear. These linear models may not be able to account for some complexities happening in the reactive distillation process under certain conditions. Based on this, it was felt that it would be better to represent this process using nonlinear models because, naturally, the process is a nonlinear and a complex one.

   Therefore, this work is aimed at developing nonlinear black-box models for ethyl acetate reactive process using tree-partition and sigmoidnetwork as the nonlinearity estimatorswith the aid of System Identification Toolbox of MATLAB 2012a [10 and 11].

2. Process Description for Data Acquisition

   The data used for the development of the models in this work were acquired from the experiments that were carried out in a reactive packed distillation column (RPDC) set up as in Figure 1 and as also shown and described in the previous researches [1, 2, 3, 4, 5 and 6]. The column, excluding the condenser and the reboiler, had a height of 1.5 m and a diameter of 0.05 m. It consisted of a cylindrical condenser of a diameter and a height of 5 and 22.5 cm respectively. The main column section of the plant was divided into three subsections of 0.5 m each. The upper, middle and lower sections were the rectifying, the reaction and the stripping sections respectively. The rectifying and the stripping sections were packed with raschig rings while the reaction section was filled with Amberlyst 15 solid catalyst that had a surface area of 5300 m2/kg, a total pore volume of 0.4 cc/g and a density of 610 kg/m3. The reboiler was spherical in shape with a volume of 3 Litre. The column was fed with acetic acid at the top (between the rectifying section and the reaction section) while ethanol was fed at the bottom (between the reaction section and the stripping section) with the aid of peristaltic pumps which were operated with the aid of a computer via MATLAB/Simulink program. All the signal inputs (reflux ratio (R), feed ratio (F) and reboiler duty (Q)) to the column and the measured outputs (top segment temperature (Ttop), reaction segment temperature (Trxn) and bottom segment temperature (Tbot)) from the column were sent and recorded respectively on-line with the aid of MATLAB/Simulink computer program and electronic input-output (I/O) modules that were connected to the equipment and the computer system. The esterification reaction occurring in the column was an equilibrium type that is given as:

   3 2 5 3 2 5 2

   CH COOH C H OH Keq CH COOC H H O

   (1)

   Applying the input data shown in Table 1 below, the experimental reactive distillation process together with the column described above was run to generate the experimental data used for the development of the nonlinear black-box models for the process with the aid of System Identification Toolbox contained in MATLAB[10 and 11].

   Table 1. Experimental input values

   Parameter	Signal Type	Initial Level	Final Level
   Reflux ratio (R)	Step	4	2
   Feed ratio (Fa/Fe)	Step	1	1.25
   Reboiler duty (Q, kJ/s)	Step	0.56	0.49

   In developing the models, the experimental data were first converted to iddata type using iddata command before being fed into the System Identification Toolbox via MATLABmfile. The main command used for the development of the models in the MATLAB environment was nlarx. After developing the models, they were simulated and the simulated (model) top segment temperatures were compared with the experimental ones using the compare command of MATLAB.

   Figure 1. Reactive Packed Distillation Column

3. Modeling Procedure 3.1.Model Structure

   The nonlinear black-box model structure (which is referred to as Nonlinear AutoRegressive with eXogenous inputs (NARX) model) used in thiswork for the prediction of top segment temperature of ethyl acetate reactive distillation process is an extension of a linear ARXmodel. A linear SISOARX model has the structurethat is given as:

   Ttop t a1Ttop t 1 a2Ttop t 2 anaTtop t na

   1

   2

   nb

   b Rt b

   Rt 1 b

   Rt nb 1 et

   (2)

   wherena is the number of past top segment temperatures(output) and nb is the number of past reflux ratios(input). The input delay nkhas been set to zero in the formulation in order to simplify the notation.This structure shown in Equation (2) implied that, with only the linear ARX model, the current top segment temperature (Ttop(t))would be predicted as a weighted sum of past top segment temperature values and current and past reflux ratio values. Rewriting the equation as a product:

   Ttopp t a1 ,a2 ,,ana , b1 , b2 ,, bnb

   Ttop

   t 1, Ttop

   t 2,, Ttop

   t na, Rt , ut 1,, Rt nb 1T

   (3)

   where Ttop(t1),Ttop(t2),…,Ttop(tna),R(t), R(t1),…,R(tnb1) are delayed output and input variables, called regressors. The linear ARX model thus predicted the current output Ttoppas a weighted sum of the regressors.

   This structure wasfurther extended to create a nonlinear form by utilizing a more flexible nonlinear mapping function instead of the weighted sum that represented a linear mapping. The form of the nonlinear ARX model is:

   Ttopp t f Ttop t 1,Ttop t 2,Ttop t 3,, Rt, Rt 1, Rt 2,

   (4)

   wheref is a nonlinear function. Inputs to f are the model regressors.The nonlinear functions used in this work were tree partition and sigmoid network.

   Shown in Figure 2 below is a block diagram representing the structure of theNonlinear ARX model used for the development of the nonlinear black-box models involving the reflux ratio (input) and top segment temperature (output).

   Regressors

   R(t)

   Estimator

   Nonlinear Function

   Ttop(t)

   Linear Function

   Figure 2. Structural block diagram of aNonlinear ARXModel

   As can be observed from the figure, the structure of the NARX model consisted of both a linear and a nonlinear function. In the calculation aspect of the nonlinear ARX model,the top segment temperature was estimatedby first computing the regressors from the current and past input values and past output data using both the linear and the nonlinear function blocks of the nonlinearity estimator. Further, the nonlinearity estimator block mapped the regressors to the model output using a combination of nonlinear and linear functions.

   The nonlinearity estimator block included linear and nonlinear blocks in parallel as in [11]:

   f x LT x r d gQx r

   (5)

   xis a vector of the regressors. LT(x)+d is the output of the linear function block and is affine when d0. dis a scalar offset. g(Q(xr))represents the output of the nonlinear function block. ris the mean of the regressors x. Q is a projection matrix that makes the calculations well-conditioned. The exactforms of f(x) were based on treepartition and sigmoid networks that were used in this work.

   Tree partition, coded using the treepartitioncommand in MATLAB R2012b, is a class representing binary-tree nonlinearity estimator for NARX models whereas sigmoid network which was coded in the same MATLAB version using the sigmoidnetcommand has a structure which is given as:

   n

   gx k k x k

   k 1

   k k k

   where s es 11 is the sigmoid function. is a row vector such that x is a scalar.

   (6)

   Estimating a NARX model is the computation of the model parameter values, L, r, d, Q, and other parameters specifying g. The resulting models developed using the nlarx command are idnlarx objects that stored all the model data, including model regressors and parameters of the nonlinearity estimator.

   1. Model Order Selection

      The model orders (na, nb and nk) are the model parameters used to predict the current output (top segment temperature). The determinations of the optimum values of the orders were found to be very important in order to avoid any under-fitting or over-fitting of the models. As such, the model orders used were determined by minimizing Akaike Information Criterion (AIC) and Rissanens Minimum Description Length (MDL) criterion.

      The mathematical expression for the Akaike Information Criterion (AIC) is given as[9]:

      AIC logV 2d

      N

      (7)

      whereV is the loss function, d is the number of estimated parameters, and N is the number of values in the estimation data set.The loss function V is defined by the following equation:

      1 N

      T 2d

      V det t,N t,N logV

      (8)

      N

      1

      N

      where

      N represents the estimated parameters.

      For d<<N:

      AIC

      2d

      (9)

      logV 1

      N

      while that of the Rissanens Minimum Description Length is as well described as:

      MDL

      d logN

      (10)

      V 1

      N

      Similarly, in Equation 10 above,V is the loss function, d is the total number of parameters in the structure, and N is the number of data points used for the estimation.

   2. Optimization of Nonlinear Model Estimator Units

      At present, as there is not any basic rule for choosing the number of units of the nonlinear estimators, the unit numbers of the estimators were determines by varying their numbers from 1 to 170 for the tree-partition and from 1 to 30 for the sigmoid network. The unit number with the highest fit value for each of the models was then chosen and used to develop the models.

   3. Model Parameter Estimation

      The estimation of the model parameters was carried out in MATLAB Environment by minimizing the error between the measured temperature and the model output temperature. That is,

      min et min Ttope t Ttopp t

4. Results and Discussions

   T (oC)

   top

   The acquired data set from the experimentcarried out is as shown in Figure 3 below.

   (11)

   75

   74

   73

   72

   0

   500

   1000

   1500

   2000

   2500

   3000

   6

   4

   500

   1000

   1500

   t (s)

   2000

   2500

   3000

   0

   R

   Figure 3. Dynamic response of top segment temperature to a negative step (from 7 to 4) in reflux ratio It can be observed from the experimental results shown in Figure 3 above that the top segment temperature was able to respond to the change in the reflux ratio. The existence of the change in the

   top segment temperature as a result of a change in the reflux ratio is one of the major reasons for using the reflux ratio as the input variable of the top segment temperature.

   The model orders obtained using the two criteria (MDL and AIC) are as shown in Table 2 below. From the results shown in Table 2, it can be noticed that the value of na given by the MDL criterion was higher than that of the AIC. This observation is in accordance with the information available in the literature because it has been discovered from the literature that MDL normally allows the shortest possible description of the observed data (Ljung, 1999). That is to say to that MDL can be chosen as the model orders selection criterion where a low-order model is required.

   Table 2.Model orders

   Parameters	Range used for optimization	Value given by MDL	Value given byAIC
   na	1:10	4	9
   nb	1:10	1	1
   nk	1:10	1	1

   After obtaining the model orders, another imprtant parameter considered were the unit numbers of the nonlinearity estimators (tree partition and sigmoid network) used. In order to determine the appropriate number of units for the estimators, the model orders obtained from the MDL criterion was used and the two nonlinearity estimators were used to develop the models. The fit values (describing the performances) of the models with the different unit numbers were calculated using Equation 12 which is given as:

   1 normTtopp Ttope

   tope

   fit normT

   • meanTtope

   100

   (12)

   The terms Ttopp and Ttopein the Equation (12) equation stand for experimental and predicted top segment temperaturerespectively.

   Fit (%)

   The calculated fit values were recorded and further plotted as shown in Figures 4 and 5 respectively for the treepartition and sigmoid network estimators.

   86

   84

   82

   80

   78

   76

   74

   72

   20

   40

   60

   80

   Unit number

   100

   120

   140

   160

   Figure 4. Performance of tree partition unit number

   From Figure 4, the maximum fit value of 87.03% was found to occur in the treepartition nonlinear estimator when the unit numberwas63. This value (87.03%) was found to remain constant throughthe

   unit number of 126 for this model. As such, 63 was chosen as the unit number used in developing the model with the treepartition nonlinear estimator.

   Fit (%)

   As can be seen from Figure 5, the performance of the sigmoid network was found to be zigzag with the highest performance value of 89.01% occurring with 3 units. This low value of units giving this good performance was found to be very encouraging in using the sigmoid network for the development of NARX model for the reactive distillation process being studied.

   88

   86

   84

   82

   80

   78

   76

   74

   5

   10

   15

   Unit number

   20

   25

   30

   Figure 5. Performance of sigmoid network unit number

   Based on the results obtained from Figures 4 and 5 above, 63 and 3 units were respectively used for the development of the tree partition and sigmoid network NARX model for the reactive distillation process. The results obtained from the models developed indicated that the two models had 5 regressors (4 past outputs and 1 past input).

   				Experimental Tree partition NARX model Sigmoid network NARX model

   75.5

   75

   74.5

   T (oC)

   top

   74

   73.5

   73

   72.5

   72

   71.5

   0 500 1000 1500 2000 2500 3000

   t (s)

   Figure 6.Experimental andpredicted top segment temperatures

   However, the loss function of the sigmoid network NARX model which was found to be approximately 0.00031 was discovered to be lower and better than that of the treepartition NARX

   model that had a loss function of approximately 0.00032, even though the difference between the loss functions of the two models was found to be very small.

   Simulating the developed models, comparisons were made between the experimental (measured) top segment temperature and the predicted (model) ones. The comparisons are as shown in Figure 6 above.According to the results shown in Figure 6, the top segment temperatures predicted using the two methods (treepartition and sigmoid network) were found to be in good agreement with each other and with the experimentally measured ones.

5. Conclusion

   The good conformities between the experimental top segment temperatures and those predicted using the developed nonlinear black-box (NARX) models have shown that the models can be used to successfully represent the behavior of the ethyl acetate reactive distillation process. However, owing to the higher fit value and lower loss function observed in the case of the NARX model developed using sigmoid network, it (sigmoid NARX model) has been found to be better than treepartition NARX model for this process.

6. Acknowledgements

   Abdulwahab Giwa wishes to acknowledge the support received from the Scientific and Technological Research Council of Turkey (TürkiyeBilimselveTeknolojikAratrmaKurumu – TÃœBTAK) for his Ph.D. Programme. In addition, this research was supported by the Scientific Research Project Office of Ankara University (Ankara ÃœniversitesiBilimselAratrmaProjeOfisi AÃœ BAP) under Project No. 09B4343007.

7. Nomenclatures

   AIC Akaike Information Criterion

   F Feed ratio (mL s-1Fa / mL s-1 Fe)

   Fa Acetic acid feed rate (mL/s)

   Fe Ethanol feed rate (mL/s)

   MDL Rissanens Minimum Description Length

   NARX Nonlinear AutoRegressive with eXogenous Inputs Q Reboiler duty (kJ/s)

   R Reflux ratio

   RPDC Reactive Packed Distillation Column t Time (s)

   Tbot Bottom segment temperature (oC)

   Trxn Reaction segment temperature (oC)

   Ttop Top segment temperature (oC)

   Ttope Experimental top segment temperature (oC) Ttopp Predicted top segment temperature (oC)

8. References

1. Giwa A. and Karacan, S., Black-Box Modelling of Ethyl Acetate Reactive Packed Distillation Column,

   AU Journal of Technology, 2012c, 15: 172-178.

2. Giwa A. and Karacan, S., Decoupling model Predictive Control of a Reactive Packed Distillation Column, International Journal of Advances in Science and Technology, 2012e, 4: 39-51.

3. Giwa A. and Karacan, S., Decoupling PID Control of a Reactive Packed Distillation Column,

   International Journal of Engineering Research & Technology, 2012f, Volume 1, Issue 6.

4. Giwa A. and Karacan, S., Development of Dynamic Models for a Reactive Packed Distillation Column,

   International Journal of Engineering, 2012d, 6: 118-128.

5. Giwa A. and Karacan, S., Modeling and Simulation of a Reactive Packed Distillation Column Using Delayed Neural Networks, Chaotic Modeling and Simulation, 2012a, 1: 101-108.

6. Giwa A. and Karacan, S., Simulation and Optimization of Ethyl Acetate Reactive Packed Distillation Process Using Aspen Hysys, The Online Journal of Science and Technology, 2012b, 2: 57-63.

7. Klöker M., KenigE.Y., Górak A., MarkusseA.P., Kwant G. and Moritz P., Investigation of Different Column Configurations for the Ethyl Acetate Synthesis via Reactive Distillation, Chemical Engineering and Processing, 2004: 43, 791801.

8. Lai, I-K., Hung, S.-B., Hung, W.-J., Yu, C.-C., Lee, M.-J., Huang, H.-P., Design and Control of Reactive Distillation for Ethyl and Isopropyl Acetates Production withAzeotropic Feeds, Chemical Engineering Science, 2007, 62: 878-898.

9. Ljung L., System identification: Theory for the User, 2nd Edition, Prentice Hall, Upper Saddle River, New Jersey, 1999.

10. Ljung, L., System Identification Toolbox: Getting Started Guide, The MathWorks, Inc., Apple Hill Drive, Natick, 2012a.

11. Ljung, L., System Identification Toolbox: Users Guide, The MathWorks, Inc., Apple Hill Drive, Natick, 2012b.

12. SneesbyM.G., TadeM.O. Datta R. and Smith T.N., ETBE Synthesis via Reactive Distillation. 2. Dynamic Simulation and Cotrol Aspects, Industrial & Engineering Chemistry Research, 1997, 36:1870-1881.