Loop Interaction Analysis and Control Structure Selection: Application to a Fluid Catalytic Cracking Unit

A simulation approach for the assessment of variables interaction and consequent control structure selection of a fluid catalytic cracking unit (FCCU) is presented in this paper. The simulator which was implemented in Matlab draws from an earlier mathematical model of the FCCU, was used as a virtual FCCU for studying the dynamic response of the riser temperature (Trx), the regenerator temperature (Trg) and the regenerator flue gas oxygen concentration (Od) to step changes in air flow rate (Fa), regenerated catalyst flow rate (Frc), gas oil feed rate (Fgr). The results show strong interaction in FCCU variables, with Fa affecting Trg and Od; Frc affecting Trx, Trg and Od; Fgr affecting Trx, Trg and Od. A linearised state-space model based on the first-principle model was deduced and transformed to a 3x3 input-output model. Three channel interaction measures: Relative Gain Array (RGA), Effective Relative Gain Array (ERGA) and the Normalized Relative Gain Array (RNGA) were applied to the selection of FCCU control structure. All the measures point to a diagonal scheme with the following pairings: (Trx/Fgr), (Trg/Fa) and (Od/Frc) ,for the decentralized control of the riser temperature, the regenerator temperature and the flue gas oxygen concentration respectively. The suggested control structure offers a high promise of stability, with a Niederlinski index (NI) of 101.79. DOI: 10.7176/CTI/8-03


Introduction
Variable interaction refers to a situation in which a change in the numerical value of an input variable of a process produces a proportional or more than proportional change in more than one output variable, in a simulation study. One process with a very high promise of interaction is the Fluid Catalytic Cracking Unit (FCCU) of a petroleum refinery. Fluid catalytic cracking is the breaking down of gas oils and certain atmospheric residues, in a fluidized bed reactor, into more valuable components such as gasoline and light gases. The process is encumbered by the attendant deposition of coke on the catalyst, causing activity decline. In view of its impact on overall refinery economics, the FCCU is often one of the first units considered for the application of advanced control and optimization strategies. The FCC unit is difficult to control due to its strongly interacting and highly non consist of several loops of Single-Input-Single-Output (SISO) controllers that function independently. Although they enjoy the merit of well established control theories, control loop interaction posses the threat to -linear multivariate features which result in complex dynamics 1 . Multiple-Input-Multiple-Output (MIMO) controller design procedures are better suited for addressing FCCU control issues, however, due the intricacies associated with the design, tuning and maintenance of the MIMO structure, decentralized PID control systems are largely preferred in the process industry 2,3 . Such decentralized schemes desirable performance. In the face of a large number of manipulated inputs, some of which affect more than one measured output, the challenges posed by variables interaction has attracted keen research interests 4 . In order to eliminate or at least reduce the deleterious effects of interaction on controller performance, it is expedient to properly address control structure. Control structure selection is a sub-problem of control system design and is concerned with the appropriate choice of controlled variables, manipulated variables and the pairing of such variables. Although it is a twin problem in which the existence of loop interaction needs to be established before appropriately pairing the variables, most research efforts focus on the later. Bristol 5 proposed the Relative Gain Array, RGA, as a measure of interaction. It is based on the steady state properties of process transfer function and gives clues to solving pairing problems in Multi-Loop, Single-Input-Single-Output (SISO) control structures. The application of RGA was popularized by Shinskey 6 8 . This is an improvement on the work of Bristol 5 as it is based on the dynamic properties of process transfer functions rather than on their steady state properties. Several other measures of interaction, such as the Partial Relative Gain (PRGA) the µ interaction index, Performance Relative Gain Array (PRGA), the Participation Matrix, the Hankel Interaction Index Array (HIIA) and the Effective RGA , have also been reported in the literature (9,10,11,12,13,14,15,16) . In particular, Xiang et al. 17 applied the Effective Relative gain Array (ERGA) as a measure of Multiple Input, Multiple output(MIMO) systems interaction. The quantification of interaction was based on steady state gain and critical frequency variations as prescribed elsewhere 16 An extension of the ERGA was and proposed tested by Monshizadeh-Naimi et al. 18 .The interaction measure which is known as Effective Relative Energy Array (EREA) is, according to the authors, an energy-based compromise between steady-state gain and bandwidth of the system under investigation. The method was claimed to have out-performed the conventional ERGA. However, the desired closed-loop characteristics of the system must be assumed apriori. The implication is that the method does not guarantee global convergence to variable pairing since the desired closed-loop system characteristics vary with control objectives. Hamid and Komaraji 19 employed a Gramian-based approach for the selection of control configuration. With MIMO transfer functions that were drawn from the open literature, the method of Salgado and Cornley 13 was adopted in the evaluation of the controllability and observability gramians as well as the participation matrix A two-step variant of the ERGA was presented in Amit and Babu 20 in which the effective gain of the plant was obtained as scalar product of the steady state gain matrix and the bandwidth matrix. The ERGA was then obtained as a scalar product of the effective gain matrix and the transpose of its inverse, in the same manner as the RGA. Again, transfer functions drawn from literature 18 were used to buttress the superiority of the method over the conventional RGA.
Ajayi and Oboh 21 extended a ratio method for Two-Input, Two-Output (TITO) control configuration to higher dimension square systems by the introduction of multiple ratios. MIMO transfer functions for distillation case studies were drawn from the literature to exemplify the application of the method as well as serving as basis for comparison with previous works that were based on RGA. The traditional ratio method and its extension basically scale the static gain matrix as in the normalized gain array in Cai et al. 22 . In order to measure the loop and variable interactions that are associated with the fluid catalytic cracking unit, the development of a control-relevant dynamic model is needful. Early papers on the FCCU are mainly concerned with reaction kinetics. They basically focused on a simplified lumping approach due to the difficulty in the complexity of the chemical structure of the gas oil feed. Weekman and Nace 23 , presented the oldest and simplest lump model (the three-lump model) to describe catalytic cracking reactions. Several variations of this model are now available in the literature 24,25,26,27 . They include the four-lump model, the five-lump model, the ten-lump model and the 19-lump model. An elaborate review of FCCU kinetic models has been given elsewhere in the literature 28 .
Regarding FCCU process models that establish relationships between variables, there are basically three types that appear in the open literature. The first type of models consists of variants of one dimensional mass, energy and chemical species balances. The simplest kind of these models is the homogenous version where both catalyst and gas-phase species are assumed to move with the same velocity and the gas oil feed is considered to enter the riser totally vapourized 1,28,29,30,31,32,33 Elsewhere in the literature 34,35,36,37 , the slip between the vapour phase and the catalyst was considered. The second type of models 38, 39 , is semi-empirical and is basically described as core-annulus models. They are incapable of correctly capturing the features of cracking processes at conditions different from those of model parameters estimation. More detailed models that are based on phenomenological concepts 40,41,42,43,44 , constitute the third category. They often times require simultaneous solution of the mass, energy, momentum and species conservation equations. These models consider particle flow characteristics in the framework of the kinetic theory of granular flow, they are detailed and capture the dynamics of the FCCU without leaning on simplifying assumptions that are common to the other classes of models. However, the computational load in terms of time and computer resources is heavy. In the light of this, they are not suited for on-line control and optimization studies. Ahari et al. (45) presented a static mathematical model of the FCCU with attention to the riser reactor. The four-lump model (24) , was used as a descriptor of the riser reactions while correlations for riser hydrodynamics were taken from Patience et al. (38) . Due to the combinatorial nature of this selection problem, the number of candidate control structures may be large and favorable candidates are easily overlooked. A systematic approach that does not solely rely on linear systems theory is required to handle this problem. This paper aims at applying a simulation approach on the FCCU control loop interaction and control structure selection using a basic mathematical model of the unit that captures the control-relevant dynamics. The ultimate goal is to assess control system designs based on the selected structure. Model Development Three types of models namely non-linear dynamic, linearised dynamic and input-output models were employed in this work. While the details of model development are as presented elsewhere 46,47 , the highlights are presented in the sub-sections that follow. 2.1: Non-Linear Dynamic Models 2.1.1: Riser Temperature Model: The rate of coke combustion may be obtained, using equation (3), as given by (Hovd and Skogestad 12 ): Although the dynamics of the stripper is not relevant to the subject matter, equation (8) and (9) where solved at steady state to obtain expressions for Cst and Tst respectively as given below

2: Linearised Dynamic Models
The riser model was reduced to an ordinary differential equation in time using the approximation that is given in equation (9).
Where rx T  represent the temperature drop over a displacement of L  along the riser. The reduced form of equation (1), equation (2) and equation (6) were then transformed to obtain the dynamic linearised model of the form given in equation (13). Where: The state transition matrix, A, the input matrix, B and the disturbance matrix,  were obtained by partially differentiating the reduced form of equation (1) , along with (2) and (6) with respect to each element in X,U and D respectively. Thus, The elements that constitute equation (17) were implemented in the symbolic math toolbox of MatLab ® 2008. Accordingly, the linear state-space model in compact form becomes:

2.3: Input-Output Model
The state-space-to-transfer-function feature of Matlab was exploited in the transformation of the model in the preceding section to the input-output form. Table 1 shows the elements of the resultant 3x3 FCCU.

3:
Dynamic Simulation of uncontrolled FCCU Two sets of simulation were conducted in this study and the details and outcomes are as presented in the sub-sections below

3.1: Nominal Yields and Control variables Profiles
The riser-regenerator modeling equations were solved for product yields (coke on catalyst, gasoline, light -gases), unconverted gas-oil, regenerator oxygen concentration, regenerator temperature and riser temperature variation with time. A dynamic solver-simulator 47 that was implemented in MatLab was employed as a solution tool. With the operating data and unit sizes taken from literature 11,46,47 as input to the simulator, the dynamic behaviour of the FCCU in the absence as well of presence of hick-ups was studied in simulation mode. Figure 1 and figure2 respectively show the profiles of product yields and control variables, du ring nominal operations. It could be readily deduced from the figures that these variables assume and maintain the steady state values in the absence of hick-ups in input (manipulated) variables.

3.2: Response to step change in Input (manipulated) Variables
Step changes of +5 percent of the nominal values followed by -5 percent of the nominal values of the manipulated inputs were applied to the virtual FCCU at t=1000 minutes and t=1500 minutes respectively. Figure 3, Figure 4 and Figure 5 each show the riser temperature variation, regenerator temperature variation, variation of oxygen concentration in the regenerator and the variation of coke concentration on the regenerated catalyst with time. A 5% increase in the nominal value of air flow rate (Fa) at t= 1000 minutes causes a shift to a new regenerator steady state temperature that is higher than the nominal. However, a 5% reduction in the nominal Fa value at t-15000 minutes caused the regenerator temperature to attain a new steady state value that is lower than the nominal. Without operator intervention, the regenerator maintains this new steady state temperature (Figure 3).While coke concentration and oxygen concentration profiles are the direct opposite of regenerator temperature, there are no notable changes in riser temperature profile following step changes in air flow rate, in the light of Figure 3.
Step changes in gas oil flow rate (Fgr) and regenerated catalyst flow rate (Frc) each produced responses ( Figure 4 and Figure 5) in regenerator temperature, oxygen concentration and coke concentration that are similar to the observed trends in Figure 3. However, the riser temperature assumed a new steady state value as opposed to the observed trend in Figure 3 where the effect of changes in Fa as not noticeable . Based on the observed profiles following step changes in the respective input variables, strong variable interaction can easily be inferred. The observed profiles during nominal steady operations were compared with the step response plots and the results are as summarized and tabulated in table 2. The cause-effect relationships given in table 2 are consistent with the functional relationships between the state variables and the input variable in equation (1), equation (2) and equation (5) respectively.

Control Structure Selection
Three channel interaction measures: The Relative Gain Array (RGA), Effective Relative Gain Array (ERGA) and Relative Normalized Gain Array (RNGA) were applied to control structure selection in this work. RGA was selected due to its wide acceptability as a first choice while the other two were employed to check the structure that is suggested by the RGA. Only highlights of interaction measures are presented here as details are available elsewhere.

4.1
Relative Gain Array (Bristo l5 , Hamid 19 ) Consider a multivariable (n xn ) plant that is represented by the transfer function given in equation (19)      (19) The RGA () is obtained as Where G(0) is the static gain matrix: nn n n n n n g g g g g g g g and  denotes element-by-element product.

Effective Relative Gain Array (Xiang et al. 19 )
Given the plant in equation (8), the ERGA is obtained as  = E  E -T (22) Where E = G(0)   (23)  is the matrix of the critical frequencies of the elements in equation (19) and as follows:

Relative Normalized Gain Array (He et al 15 )
The RNGA is computed according to equation (14) as = KN  KN -T (25) Where the normalized gain KNij for each element in equation (19) is computed as In (26), ij is the average residence time of the multivariable transfer function element, gij and it is a measure of how fast a controlled variable responds to a change in manipulated variable. ij is the sum of two components, the delay time, ij and the time constant, iij according to equation (16).
In this work, a spline function was implemented in MatLab to fit the transfer functions and the time taken to achieve a step response of 63.2% of the statistic gain was calculated as consistent with equation (27) However, the critical frequencies were obtained by finding the bandwidth of the close-loop transfer function.

4:
Application to FCCU Control Structure Selection Following section 3.1 to 3.3, the interaction measures RGA, ERGA and RNGA were implemented in MatLab R2008. The transfer functions given in table 1 were used as the only input while the other arrays such as the critical frequency () the static gain G(0), average residence time (Tar), RGA, ERGA and RNGA were returned as output from the computer programme (interact.m) that was developed in this work. Accordingly, the outputs are as given below The three interaction measures point to a diagonal pairing of input-output variables in which the riser temperature is controlled by manipulating the gas oil feed rate(Trx/Fgr),air flow rate is manipulated to control the regenerator temperature(Trg/Fa) and the regenerator oxygen concentration is controlled by manipulating the flow rate of regenerated catalyst(Od/Frc).

Conclusion
A basic mathematical model of the fluid catalytic cracking unit (FCCU) that captures the control-relevant dynamics of a fluid catalytic cracking has been used as a virtual plant to study the interaction of variables and control structure selection. Non-linear dynamic analysis in simulation mode for capturing variables interaction has been presented along with linearised input-output model to which the RGA,ERGA and RNGA were applied. A diagonal control structure with very reasonable row dominance that suggests stability has been selected and presented.

6.
Notation Ar = Area of riser, m 2 Cpa = Heat Capacity of air, kJ/kg-K Cpc = Catalyst heat capacity kJ/kg-K Cpf = Heat capacity of feed, kJ/kg-K Crc=Coke on regenerated catalyst, wt % ΔHrcb = Heat of coke combustion, kJ/kmol ΔHij = Heat of reaction , kJ/kmol EAcb= Activation energy for coke combustion, kJ/mol ɛgr = Hydrocarbon volume fraction in the riser Fa = Flow rate of regenerator air, kg/s Fgr = Gas oil feed rate, kg//s Frc = Flow rate of regenerated catalyst, kg/s Fsc = Flow rate of spent catalyst, kg/s Fsr = Catalyst flow rate, kg/s kcb = coke combustion rate constant,s -1 kij= Rate constant for riser reaction,1/s Lr = Length of riser, m Ma= Molecular mass of air, kmol/kg Mc= Molecular mass of coke, kmol/kg Mgr =Molecular mass of hydrocarbon in the riser, kmol/kg Mor= Air hold-up in regenerator, kg Mrg = Catalyst hold up in regenerator, kg Mst= Catalyst hold up in separator, kg Oin= Mole fraction of oxygen in air Od= Mole fraction of oxygen in regenerator flue gas ri =Rate of gas oil cracking reaction.s -1 rcb = Rate of coke combustion Rg = Gas constant, kJ/ kmol -K ῥc= catalyst density, kg/.m 3 ῥgR= Hydrocarbon vapour density, kg/.m 3 Ta= temperature of air entering regenerator, K Tc= temperature of regenerated catalyst, K Trx = riser outlet temperature,K Trg = Regenerator Dense Bed temperature, K Tst = Stripper temperature, K Tvapf= feed vaporization temperature, K yij =Initial and boundary values for riser product yields, wt fraction y1= Steady state gas oil Concentration in the riser y2 = Steady state gasoline concentration in the riser y4 =Concentration of coke on catalyst entering stripper, wt fraction z= Dimensionless Length 7.