INTRODUCTION
The industrial scale 2ethylhexenal hydrogenation process produces 2ethylhexanol oxoalcohol, used as raw material for dioctylphthalate (DOP) or 2ethylhexylphthalate (DEHP) production. DOP and DEHP represent around 50% of the global consumption in the production of poly(vinyl chloride) (PVC) (ICIS,2011;WRI,2009;SRI,2011).
Because the hydrogenation reaction is highly exothermic, the reactor output temperature is of great importance. Like most of chemical and petrochemical plants, the 2ethylhexenal hydrogenation process presents large time delays and is multivariable by nature. For multivariable time delay processes, the control solution will revolve around a multivariable Smith Predictor. A very important step for the multivariable Smith Predictors is the decoupling of the process.
The present paper is structured in 8 parts. After the introduction, the second part describes the hydrogenation process and also details the actual control and the new proposed control system. The third part presents the developed dynamic mathematical models of the hydrogenation process: the nonlinear distributed parameter mathematical model and the linear simplified mathematical model. Both models are validated using plant data, acquired courtesy of S.C. Oltchim S.A., Ramnicu Valcea, Romania. In the control algorithm design procedure for multiinput multioutput (MIMO) systems, it is necessary a measure of interactions within MIMO systems, which can be provided by the Relative Gain Array (RGA) analysis presented in section four. Sections five and six present a multivariable PID controller design, respectively IMC control design. In section seven the comparative simulation results between the two control strategies are presented. The conclusions represent the last part of the paper.
PROCESS DESCRIPTION
The 2ethylhexenal hydrogenation process is complex, which can take place both in liquid phase or gas phase in the presence of a catalyst and, besides the main product (2 ethylhexanol), can also form other sideproducts like: nbutanol or isobutanol from secondary sidereactions. Technological solutions present both advantages and disadvantages. The gas phase 2ethylhexenal hydrogenation reaction is characterized by high energy consumption and great dimensions of equipment, but has an important advantage: the possibility to regenerate the catalyst. In comparison to the gas phase hydrogenation reaction, the advantages of the liquid phase hydrogenation reaction of 2ethylhexenal are: low energy consumption, small dimension of equipment, and small volumes of catalyst. However, the main drawback is the impossibility of catalyst regeneration. Industrial scale liquid phase hydrogenation is predominant in chemical industries and represents the main focus of the present work.
An important role in the hydrogenation reaction is played by the catalyst. For the present case study, a nickel catalyst deposited on silica support was chosen based on the BASF license. This type of catalyst for the 2ethylhexenal hydrogenation reaction is also recommended in specific literature (^{Smelder, 1989}).
Several researchers like ^{Bozga (2001)} and Cocker (2011) studied the hydrogenation reaction mechanism. Hence, the reaction mechanism stages for the present case study are (^{Bozga et al., 2001}; ^{Cocker, 2011}): the transport of hydrogen from the gas phase to the gasliquid interface; the transport of hydrogen from the gasliquid interface to the liquid phase; the transport of reactants 2ethylhexenal and hydrogen from the liquid phase to the liquidsolid interface; reactant adsorption to the active sites of the catalyst; the chemical reaction between adsorbed reactant molecules; product desorption from the active sites of the catalyst; the product transport from the catalyst surface to the liquid phase. The reaction pathway is expressed as follows (^{Collins et al., 2009}):
where: A is 2ethylhexenal (reactant), B is 2ethylhexanal (intermediate product) and C is 2ethylhexanol (final product).
The main parameters that influence the two consecutive hydrogenation reactions (the intermediate product being 2ethylhexanal) are: the pressure, temperature (input/output), reactor loading, reflux rate and purity of input reactants (hydrogen, 2ethylhexenal). The temperature is a critical parameter due to the fact that the hydrogenation reaction is highly exothermic with a reaction heat of 25.27 x 10^{3} kcal/kmol (^{Smelder,1989}). High temperatures favor sidereactions, but also an input temperature of the reactants between 90°C and 110°C is necessary in order to start the reaction. The maximum technological recommended temperature inside the reactor is in the range of 160 °C  180 °C, depending on the catalyst degree of activity. The solution chosen to maintain the temperature below the critical value is the dilution of the 2ethylhexanal flow with the product itself (the reflux rate). A typical industrial process flow diagram of the 2ethylhexenal hydrogenation process is presented in Figure 1.
As presented before, a parameter of great importance is the reactor loading. This influences the hydrogenation reaction itself and also the temperature. It was determined experimentally that a specific flow ratio between the 2ethylhexenal flow and the recirculated 2ethylhexanol flow must be maintained to carry out the reaction with good parameters and also to maintain the temperature in the recommended range. The catalyst degree of activity acts as a disturbance because it influences the hydrogenation reaction and implicitly the output temperature of the product. The necessary steps to counteract the effect of a decrease in the catalyst activity are: to change the flow ratio (the recirculated 2ethylhexanol flow must be decreased) and also to increase the input temperature of the reactants to start the hydrogenation reaction.
At the present time, the hydrogenation reactor is operated using an open loop control system, so the main parameters like 2ethylhexenal input flow rate, recirculated 2ethylhexanol flow and input temperature of the reactants are maintained at the imposed values, without feedback control of output temperature or 2ethylhexanol output concentration. This control solution is characterized by a main disadvantage: it is not faulttolerant, being difficult to apply correction methods if one control loop fails. However, the hydrogenation reactor will operate in a satisfactory manner if the auxiliary control loops for flow control of the reactants and input temperature control of the reactants are working in good parameters.
At present an operator sets manually the reference values for the 2ethylhexenal input flow rate, recirculated 2ethylhexanol flow and input temperature of the reactants. The 2ethylhexanol output concentration is measured periodically and, if it is not satisfactory, the reference values for the main input variables are adjusted by the accumulated expertise of the operator. In order to optimize the production, the authors propose a feedback control system using the existing infrastructure. The implementation costs are reduced to the costs of a temperature transducer used to measure the output temperature of the product.
In order to ensure safe and normal conditions of operation for the hydrogenation process there are two main goals of the control strategy. The first goal is to maintain the temperature at an acceptable value. This goal can be achieved by assuring a specific volumetric flow ratio between the 2ethylhexenal and 2ethylhexanol flow rates. The recirculated 2ethylhexanol flow rate is considered to be the manipulated variable. The second goal is to ensure a high conversion of 2ethylhexenal to 2ethylhexanol. This goal can be reached by assuring the necessary input temperature of the reactants to start the reaction. The input temperature of the reactants must be increased as the catalyst degree of activity decreases in time. The input temperature of the reactants also has an important influence on the output temperature. Overall, the two considered control inputs are the 2 ethylhexanol flow rate and the input temperature of the reactants. The output temperature and product concentration are considered to be the measured outputs and the input hydrogen pressure (P), and 2ethylhexenal flow (Qenal) act as disturbances.
MATHEMATICAL MODELS
Distributed Parameter Nonlinear Mathematical Model
An accurate mathematical model is needed to evaluate the operational challenges and to understand the processes that occur inside the reactor and also to develop a more efficient control strategy. In the specific literature, a number of kinetic studies of the hydrogenation reaction are reported (^{Smelder, 1989}; ^{Niklasson, 1987}, ^{1988}). However, no mathematical model of the process is given.
Motivated by these reasons, in previous works (^{Both, 2013}), a first principle based mathematical model was developed starting only from some equations describing the hydrogenation reaction kinetics. The developed model (^{Both, 2013}), consisting of numerous partial differential equations and several algebraic ones, was implemented and validated using plant data acquired courtesy of S.C. Oltchim S.A. This model is useful in the design, optimization and operation of the chemical reactor and especially in control design and testing. The reaction rate equations, transport models, energy and mass balances were coupled and were solved with respect to time and space with algorithms suitable for partial differential equations (^{Attou et al., 1999}; ^{Burghardt et al., 1995}, ^{Iliuta et al., 1997}).
The equations describing the nonlinear mathematical model used are (^{Both, 2013}):

The total mass balances for the gas and liquid phases are (^{Both, 2013}):

The component mass balances for liquid and gas phases are (^{Both, 2013}):

The heat balances for the liquid and gas phases are (^{Both, 2013}):
where F_{L} and F_{G} are the liquid and gas phase flow, respectively, S is the crosssectional area, a_{v} is the specific gasliquid contact area, N_{HG} is the flux of hydrogen transferred from gaseous phase to liquid phase, C_{enal} is the concentration of 2ethylhexenal, C_{anal} is the concentration of 2ethylhexanal, C_{oct} is the concentration of 2ethylhexanol, is the concentration of hydrogen in the liquid phase and is the concentration of hydrogen in the gaseous phase. vph is the transferred hydrogen flow through the liquid film, adjacent to the catalyst pellet, onto its surface (^{Bozga, 2001}) and r_{i} is the rate of surface reaction i and Δ_{R}H_{i} is the reaction heat of reaction i.

The reaction rates are described by (^{Smelder, 1989}):
where K_{i} is the adsorption equilibrium constant for component i (i: enal, anal, H), C_{i} concentration of component i and k_{i} rate constant of surface reaction i,(i:1,2).
The rate of reaction incorporating the catalyst deactivation can be obtained as follows:
where: a is a fraction of active catalyst and r_{i} is the reaction rate of species i.
For validation of the developed model a considerable number of trials were chosen to compare simulation results with plant data (^{NIST, 2011}; ^{Gaspar et al., 2010}; ^{Perry et al., 1999}; ^{Tobiensen et al., 2008}, ^{Silva et al., 2003}) and the two most representative cases are presented in Table 1, which emphasizes the small differences between the plant data and simulated values.
Run  H_{2} concentration (kmol/m^{3})  2ethylhexenal concentration (kmol/m^{3})  2ethylhexanol concentration (kmol/m^{3})  Outlet liquid temperature (K)  

sim  plant  sim  plant  sim  plant  sim  plant  
1  0.078  0.0549  0.0021  0.0020  5.346  5.32  441.6  441.3 
2  0.088  0.055  0.0022  0.0023  5.343  5.31  442.1  441.5 
The accuracy of the distributed parameter nonlinear mathematical model of the process is highlighted by the evolutions of 2ethylhexenal, 2ethylhexanol concentrations and the temperature evolution along the reactor height, presented in Figure 2 and Figure 3.
Operational Mathematical Model
Due to the fact that nonlinear mathematical models are more accurate, but also too complex for efficient use in controller design, the authors propose another approach, to use a simple operational model of the process which describes its most important properties (^{Iancu et al., 2010}; ^{Simson et al., 2001}).
The linear operational mathematical model was determined using experimental identification methods based on plant data, as a transfer function matrix of second order elements between the main input and output variables (^{Both, 2012}): the 2ethylhexanol recirculated flow (Qoct), the input temperature of the reactants (Tin), the output temperature of the product (Tout) and the output concentration of the product (Cout).
where
The input hydrogen pressure (P) and 2ethylhexenal flow (Qenal) act as disturbances. The transfer matrix describing the dependence between the output variables and the disturbances is:
where
Following the previous equations, the transfer matrix of the simplified hydrogenation reactor is:
Using the operational mathematical model, the plant data and the simulation results of the nonlinear validated mathematical model, analysis of the dynamical behavior of the simplified model is done. Figure 4 and Figure 5 present the evolution of the output temperature (Tout) for a +14% step variation of the 2ethylhexenal input flow and for a step variation in the 2ethylhexanol recirculated flow.
The validation of the simplified linear mathematical model is based on the comparison between equation results and plant data. The presence of an acceptable error is emphasized in the previous figures.
RELATIVE GAIN ARRAY (RGA) ANALYSIS
Considerable research effort was dedicated to the development of control techniques for MIMO processes, but it still represents a difficult task in practical engineering. At the present time, there are two main directions of approach in control design for multivariable processes. The first approach supports the idea to design a centralized MIMO controller, which implies the use of algebraic decoupling methods or optimal control theory. However, the result is a complex controller, difficult to implement. Another approach is to design a decentralized controller, which attempts to control a fairly multivariable process by separating the control problem into several SISO (single  input, single  output) systems and then uses conventional control for each loop. The advantages of this approach are very clear: simple to design, easy to implement or tune (fewer parameters); hence it is more often used in practical applications. However, the main disadvantage is represented by the substandard closed loop performance due to interactions between loops. A mandatory step is to determine the appropriate loop configuration in order to obtain minimal interactions between loops. This is achieved by proper pairing of the manipulated and controlled variables. In the chemical industry the process complexity reaches up to several hundred control loops, making input  output pairing a difficult task. This task is also essential because it can lead to an unstable overall system, even if the individual loops are stable.
In 1996, ^{Bristol (1966)} introduced the Relative Gain Array technique to determine a measure of the process interactions, giving advice in solving the pairing problem. This technique considers only the steadystate gains of the process to determine the best pairing solution. Due to its simplicity, the RGA technique was widely used in different applications and still has the widest application in industry. The main disadvantage of this technique is represented by the lack of dynamic information of the process, which may sometimes lead to an incorrect loop pairing decision. A dynamic version of the RGA technique was developed by ^{McAvoy (1983)}. The dynamic RGA (DRGA), in comparison with the RGA, uses the transfer function model instead of the steadystate gain matrix in order to consider the process dynamics, giving a more accurate interaction assessment. ^{Xiong (2005)} combined in his work the advantages of both RGA and DRGA by introducing the concept of an effective relative gain array (ERGA), which provides an accurate, detailed and easy to understand description of the dynamic interaction among individual loops without the need to specify the controller type.
In this study the RGA technique was used to determine the best control loop pairs and, at the end, the performances of the selected control loop pairs are studied. RGA is represented by a matrix with one column for each input variable and one row for each output variable. From this matrix one can easily compare the relative gains associated with each inputoutput variable pair and match the input and output variables that have the biggest effect on each other, while also minimizing undesired side effects. Standard analysis suggests that RGA elements corresponding to inputoutput pairings close to 1 should be preferred. Negative or large RGA elements are disadvantageous as they correspond to loops which may be nominally stable, but which become unstable if saturation occurs (e.g. (^{Bristol, 1966})).
The simplest approach was found to be the use of the developed linear mathematical model and and to compute RGA from the steadystate gain matrix. The steadystate gain matrix is thus obtained with its elements:
The RGA, for a nonsingular square matrix G, is a square complex matrix defined as RGA(G)=G×(G^{1})^{T}, where x denotes elementby element multiplication (Hadamard or Schur product).
For the obtained steadystate gain matrix H_{p0}
and
From the RGA matrix one can observe that there is a strong coupling between Qoct and Tout and Tin and Coct.
MULTIVARIABLE SMITH PREDICTOR PID CONTROLLER DESIGN
The previously presented linear simplified operational mathematical model of the hydrogenation process consisting of a transfer matrix Hp, which describes the influence of the control inputs u on the controlled outputs y (Tout, Cout), and a transfer matrix H_{dist}, describing the effect of disturbances in the feed z (the input 2ethylhexenal flow Qenal, hydrogen pressure P) on the outputs y is used in order to develop a multivariable PID controller.
Both matrices Hp and D are (2,2)matrices.
The hydrogenation process is characterized by the presence of time delays. A typical approach to deal with time delay is the nondelayed output prediction (^{Melo, 2008}; ^{Wang, 2008}). The nondelayed output may be estimated and the controller can be computed as for a process without delay. The most popular output predictor is the Smith Predictor, Figure 6.
For the controller design the transfer matrix of the process without disturbance is used:
The desired multivariable controller matrix has the following form [34]:
where HR_{11}, HR_{22} are designed for the direct control of the outputs and HR_{12}, HR_{21} are designed to compensate the coupling effect.
Considering the large time delays of the plant the Smith Predictor structure is used, presented in Figure 7.
All the controllers are computed using the imposed phase margin design method, for γ_{k}*=60º.
The obtained controllers are:
INTERNAL MODEL CONTROL DESIGN
As presented before, the most popular structure for time delay compensation is the Smith Predictor structure. A very important step in the multivariable Smith Predictor structure is the decoupling of the process. A modified multivariable Smith Predictor structure was proposed by Wang, Zou and Zhang (^{Wang et al., 2000}). In this case the decoupling matrix is computed in frequency domain. More recently the Internal Model Control (IMC) method was used in order to compute the decoupling matrix (^{Wang et al., 2002}). ^{Chen et al. (2011)} used the IMC method only as an intermediary step in order to compute the final controller  a matrix of PI controllers. The method used in this paper is based on the approach proposed by ^{Pop et al. (2011)} in which the IMC controllers are used as final controllers.
The model of the process is assumed to be equal to the process transfer function matrix Hp, presented before:
The next steps in order to design the IMC controllers are: a) Determine the pseudoinverse matrix; b) Determine the decoupled process; c) Approximate the elements on the first diagonal of the steadystate decoupled process matrix; d) Design the IMC controllers.
The method used to decouple the multivariable system is the pseudoinverse (Skogestad, 2011) of the steady state gain matrix:
The pseudoinverse is computed based on the following equation (Skogestad, 2011):
where is the pseudoinverse which will act as a precompensator matrix.
The decoupled multivariable process is done by:
The matrix (36) represents the steadystate decoupled process, in which all elements are weighted sums of the original transfer functions Hf_{ij}·e^{sτ} Due to the static decoupling, the steadystate matrix H_{D}(s) = 0 will be equal to the unit matrix. In this way all the elements which are not on the first diagonal will be equal to zero.
The next step is to approximate the elements on the first diagonal of the matrix H_{D}(s) with simple transfer functions:
The elements on the first diagonal, approximated by the transfer function given in (37), will be used for controller design.
The approximation of each element Hf_{iid}*(s) can be obtained by using genetic algorithms or graphical identification methods (^{Chen et al., 2011}).
The next step is to design the IMC controller as follows (^{Chen et al., 2011}):
where f(s) is the IMC filter and Hf_{iidinv}*(s) represents the invertible part of Hf_{iid}*(s). The IMC filter is computed as:
where λ_{i} is the time constant of the filter associated with each output and n is chosen such that the final IMC controller is proper. The parameters of the IMC controller will be tuned in accordance with the imposed performance values of overshoot and settling time. Figure 8 and Figure 9 show the closed loop response of Hf_{11d} and Hf_{22d} considering different values of the IMC filter, λ= 5, 10, 15, 20. Analyzing the simulation results in terms of settling time and overshoot, the final IMC controllers for the case study results in the following transfer function:
Figure 10 presents the closed loop control scheme using delay compensator and IMC controller. The MoorePenrose pseudoinverse is used as a precompensator for the process H_{p}(s), the model H_{m}(s) and the fast model H͂_{m}(s) (the process model without the time delays). H_{D}*(s) is the transfer function matrix of all elements Hf_{iid}*(s) and is the transfer function matrix H_{D}*(s) without the time delay.
SIMULATION RESULTS
For analysis purposes both control structures were implemented in Matlab/SIMULINK (Matlab User Guide, 2008). The simulation scenarios are presented in comparison for the two control strategies in order to conclude the results. The two main objectives of the scenarios are the setpoint tracking analysis and disturbance rejection analysis. Due to the fact that, for both control strategies, the controllers were designed based on the linear model of the process, it is necessary to test the robustness of the control strategies considering the variation of the model parameters. This analysis is performed in the third simulation scenario.
The first simulation scenario is focused on testing the setpoint tracking performances of the designed control strategies.
To this end a reference step variation of 10 K for the output temperature of the product at t=10s is considered. Figure 11 presents the evolution of the output temperature while Figure 12 shows the influence of the reference step variation for the output temperature on the 2ethylhexanol concentration.
For the same analysis a reference step variation of 0.25 [kmol/m^{3}] is considered for the 2ethylhexanol output concentration. The evolution of the 2ethylhexanol concentration is presented in Figure 13. and the influence of the reference step variation on the output temperature is detailed in Figure 14.
To reach the second objective of the section, the second simulation scenario is focused on the disturbance rejection analysis of the designed control strategies. Figure 15 and Figure 16 present the effects of a 0.25 [kmol/m^{3}] disturbance in 2ethylhexenal flow on the two considered outputs.
From all the above figures it can be concluded that the IMC control strategy presents better performances in comparison with the MIMO PID control. Another advantage of the method is that it is more straightforward and easy to design. The minor disadvantage is the LagLead form of the controllers opposite to the classical PID forms.
The last objective of this section is to test the robustness of the designed control strategies. To this end, the following figures present the evolution of the output temperature and 2ethylhexanol concentration considering a reference step variation for both variables and also the variation of the uncertain parameters of the linear plant.
The closed loop simulations using the multivariable Smith Predictor PI controller under nominal parameter values (solid black line), as well as the simulations considering the variation of the uncertain parameters of the linear plant (solid grey lines), for a reference step variation of the output temperature are presented in Figure 17 and Figure 18. The uncertain parameters considered are the time constants of the process as well as the gain parameter.
The closed loop simulations using the multivariable Smith Predictor PI controller under nominal parameter values (solid black line), as well as the simulations considering the variation of the uncertain parameters of the linear plant (solid grey lines), for a reference step variation of the 2ethylhexanol concentration are presented in Figure 19 and Figure 20.
In order to choose the control algorithm with the best performances, the same simulation scenario is considered for the IMC controller. To this end, the closed loop simulations using the IMC controller in Smith Predictor structure under nominal parameter values (solid black line), as well as the simulations considering the variation of the uncertain parameters of the linear plant (solid grey lines), for a reference step variation of the output temperature are presented in Figure 21 and Figure 22.
The closed loop simulations using the IMC controller in Smith Predictor structure under nominal parameter values (solid black line), as well as the simulations considering the variation of the uncertain parameters of the linear plant (solid grey lines), for a reference step variation of the 2ethylhexanol concentration are presented in Figure 23 and Figure 24.
Based on the simulation results presented, it can be concluded that the most appropriate control strategy with satisfactory results is the IMC control strategy in Smith Predictor structure.
It can be observed that the differences between the performances of the nominal and perturbed systems are smaller than in the case of classical control, Figure 17 to Figure 24. It must be mentioned that the simulation scenarios are considered for four extreme points of operation, implying a very large range of parameter variation. If the variation range of the uncertain parameters is smaller, in normal operation conditions, the differences between the performances of the nominal and perturbed systems are very small.
CONCLUSIONS
This paper describes the two control strategies of an important chemical process, the 2ethylhexenal hydrogenation process, in a catalytic fixedbed reactor. This process is a nonlinear system with distributed parameters. However, it is possible to describe the processes that occur inside the reactor by a linear nominal transfer matrix.
Several simulation scenarios are considered in order to choose the control algorithm with the appropriate performances. The simulation results presented in the sixth section of the paper show that the IMC control presents better performances for setpoint tracking and also for disturbance rejection tests, and better robustness, so it is recommended for implementation. Additionally, this method reduces the design burden of the controllers.
NOMENCLATURE
C_{i}  molar concentration (kmol/m^{3)} 
c_{pi}  specific heat capacity (kJ/kmol K) 
F  volumetric flow rate (m^{3}/s) 
T  Temperature (K) 
Δ_{R}H  reaction heat (kJ/kmol) 
K  adsorption equilibrium constant, component i (m^{3}/mol) 
k  rate constant of surface reaction j (mol/s kg) 
H_{H2}  Henry's law constant 
v  velocity (m/s) 
S  reactor crosssectional area (m^{2)} 
M  molecular mass (kg/kmol) 
a_{v}  specific gasliquid contact area (m^{1}) 
r  rate of surface reaction (mol/s kg) 
t  time (s) 
z  axial reactor coordinate 
Greek Letters  
ρ  density (kg/m^{3)} 
Subscripts/Superscripts  
G  gas phase 
L  liquid phase 
i  H_{2}, 2ethylhexenal, 2ethylhexanal, 2ethylhexanol 
j  2ethylhexenal hydrogenation reaction, 2ethylhexanal hydrogenation reactor 