Data reconciliation and optimal operation of a catalytic naphtha reformer

The naphtha reforming process converts low-octane gasoline blending components to high-octane components for use in high-performance gasoline fuels. The reformer also has an important function as the producer of hydrogen to the refinery hydrotreaters. A process model based on a unit model structure, is used for estimation of the process condition using data reconciliation. Measurements are classified as redundant or non redundant and the model variables are classified as observable, barely observable or unobservable. The computed uncertainty of the measured and unmeasured variables shows that even if a variable is observable it may have a very large uncertainty and may thereby be practically unobservable. The process condition at 21 data points, sampled from two years of operation, was reconciled and used to optimize the process operation. There are large seasonal variations in the reformer product price and two operational cases are studied. In case 1, the product price is high and throughput is maximized with respect to process and product quality constraints. In case 2, the product price is low and the throughput is minimized with respect to a low constraint on the hydrogen production. Based on the characteristics of the optimal operation, a ”self optimizing” control structure is suggested for each of the two operational cases.


Introduction
The naphtha reforming process converts low-octane gasoline blending components to high-octane components for use in high-performance gasoline fuels. "Octane" or, more precisely the octane number, is the measure or rating of the gasoline fuels antiknock properties. "Knocking" occurs in an engine when the fuel self detonates due to high pressure and temperature before it is ignited by the engine spark. Permanent damage of the engine cylinder and piston parts is a likely result of * This article was originally published as : T. Lid and S. Skogestad, "Data reconciliation and optimal operation of a catalytic naphtha reformer", Journal of Process Control, 18, 320-331 (2008). Reprinted with permission from Elsevier.
persistent "knocking". The most common measure of the octane number is the RON (Research Octane Number). By definition iso-octane (2,2,4 trimethyl pentane) is given an octane number (RON) of 100 and n-heptane an octane number of 0. A fuel with 95 RON has, by use of this measure, equal anti knock properties to a mixture of 95% of iso-octane and 5% n-heptane.
A simplified process model of a semiregenerative catalytic naphtha reformer, involving five pseudo components, was presented by Smith (1959) and validated against plant data. The same model was used in Bommannan et al. (1989), where reaction parameters were estimated from two sets of plant data, and in Lee et al. (1997) where a process with continuous catalyst regen-eration was modeled. In all three cases above, good agreement with plant data was reported. These models are used for simulation and design purposes except in Taskar and Riggs (1997) where optimal operation during a catalyst cycle, is considered. Taskar and Riggs (1997) developed a more detailed model of a semiregenerative catalytic naphtha reformer, involving 35 pseudo components. They claimed that the simplified model is an oversimplification of the process but no details of the practical consequences of the discrepancies where presented.
In this paper the simplified model of Smith (1959) is used for modeling a catalytic naphtha reformer with continuous catalyst regeneration. The model uses the unit model structure of Lid and Skogestad (2007). Scaling is applied to the process model variables and equations to improve its numerical properties. The process model is fitted to 21 data sets from the naphtha reformer at the Statoil Mongstad refinery. These data where collected in a two year period and include feed and product analysis and process measurements. The current state of the process is estimated using data reconciliation (Tjoa and Biegler, 1991), where redundancy of measurements, observability of variables and uncertainty of the estimate are examined. The same model is also used for computation of optimal operation and economical analysis of operational cases. Based on this analysis, a model predictive controller (MPC) for "optimal" operation of the process is suggested.

Data reconciliation
In this section, we summarize the equations used in this paper. For more details, it is referred to the references and the thesis of Lid (2007).
Data reconciliation is used to estimate the actual condition of the process and is here obtained as the solution of min z J(y m , z) where J(y m , z) is the objective function for data reconciliation, f (z) = 0 represents the process model, A r z = b r is used to specify known values and z r min ≤ z ≤ z r max physical constraints. The n y measured values are collected in the measurement vector y m .
If the measurement error is normally distributed N(µ, σ) and has a zero mean measurement error (µ = 0). The maximum likelihood estimate is achieved using a quadratic objective function where e = y m − y, and the measured variables represent the estimated values of the measurements y m . The measurement mapping matrix U has U (j, i) = 1 if variable j is measured and the measured value is located in y m (i), The weighting matrix Q is the inverse of the measurement error covariance matrix Σ m . If the measurement error is normally distributed N(µ, σ) with nonzero mean µ the quadratic objective function will result in a biased estimate. In data reconciliation, a mean measurement error µ = 0, is called a gross error.
In this work, the Combined Gaussian distribution of Tjoa and Biegler (1991) is used to handle data sets with gross errors, see figure 1.  The Combined Gaussian distribution is described by the following objective function which has two adjustable parameters, p and b. In sum-mary, data reconciliation is based on the Combined Gaussian objective (4), whereas the Gaussian objective (2) is used for analysis of the uncertainty in the estimate.

Scaling of the variables and model
To improve the numerical properties, the process model f (z) = 0 and linear constraints Az = b are scaled according to the scaling procedure proposed in Lid and Skogestad (2007).
First, every equation is paired with one variable. The equation-variable pairing may be regarded as "equation i is used for computation of the value of variable j". It is written in a matrix P, where P (i, j) = 1 if variable j is paired with equation number i. All other values equal zero. This is done both for the nonlinear process model f (z) and the linear constraints A.
Second, all variables z are scaled z = S vz , such that the scaled variablez has a value close to one. S v is a n z × n z fixed diagonal scaling matrix.
Finally, the equation scaling matrices of the process model and the linear constraints, S f and S l , are computed as where × denotes element by element multiplication so that S f and S l are diagonal matrices. The scaled model is written If the model equations are properly scaled, the condition number of should be reasonable low (< 1 × 10 6 ).
It should be noted that the variable scaling has some pitfalls. A simple input-output mass balance of a two component process stream is used as an example.
(12) The condition number of H is in this case ≈ 5.3. If the feed composition specifications are changed to x 1 = [0.01 0.99] T the condition number of H is ≈ 6.7. This shows that small values of the variables x 1 (1) and x 2 (1) are not a problem. However, if variable scaling is added, such that the scaled variables have a value of ≈ 1 the condition number of H is ≈ 7.4 × 10 3 . That is, we have by improper variable scaling created an "ill conditioned" model.
On the other hand, if the molar flow F 1 is increased from 1 to 100 the condition number of H is ≈ 2.8×10 4 . If the flow variables are scaled such that the scaled variable has a value ≈ 1, and the equations are scaled according to the procedure above, the condition number of H reduces to ≈ 8.2. The "rule of thumb", which was applied to this model, is: be careful by assigning large variable scaling factors to variables with values close to zero. Typically, all molar fractions are in [0 1] and by definition close to one and are scaled by a factor equal to one. Scaling the reformer model according to the procedure above reduces the condition number of H from 2.3 × 10 12 to 3.6 × 10 4 . The maximum absolute value of the elements in H is reduced from 4.8 × 10 5 to 7.6 and all values of H corresponding to the equation-variable pairing has a value equal to one. 4 Case study: Naphtha reformer 4.1 Process description and model structure The feed to the naphtha reformer is a crude oil fraction from the refinery crude unit with a boiling range of ≈ 100 − 180 • C and a density of ≈ 763kg/m 3 . The products are high-octane naphtha, also called "reformate", "gas" (C 2 − C 4 ) and hydrogen. The increase in octane number is due to a conversion of paraffins and naphthenes to aromatics. The amount of catalyst in the four reactors is approximately in the ratio 1:1:2:3. The reactor inlet temperatures are in the range 770K-800K.
The overall reaction is endothermic and there is a significant temperature drop from the inlet to the outlet of the reactors. In order to compensate for this temperature drop, the reactor is separated into four sections with intermediate reheating, see figure 2. The fresh feed is mixed with hydrogen rich recycle gas and is preheated in the reactor effluent heat exchanger (E1). The feed is further heated in heater no. 1 (H1) before it enters reactor no. 1 (R1), and so on. The hot reactor product enters the feed pre-heater (E1) and is further cooled with cooling water before it enters the separator. Hydrogen rich gas is compressed, except for a small purge stream, and recycled. The liquid product from the separator (D1), a mixture of reformate and gas, is separated in a downstream distillation column.
The components in the process are lumped into five pseudo components. These are hydrogen (H), "Gas" C 2 − C 4 (G), paraffines (P), naphthenes (N) and aromatics (A). A description of the thermodynamic properties of these pseudo components can be found in Lid (2007). The justification for this simplification is that the carbon number of the molecules does not change in the two reactions (13) and (14). For example, a C 7 naphthene is converted to a C 7 aromatic and a C 7 paraffin is converted to a C 7 naphthene. This conversion is described by four main reactions (Smith, 1959 where the columns refer to the components H, G, P, N and A. The reaction rates are, where p x is the partial pressure of component x and p is the total reactor pressure.
For the forward and reverse rate constants, k f and k r , an Arrhenius type of rate expression is assumed where the activation energy E is dependent on the catalyst and k 0f is dependent of the molarity of the reaction (Bommannan et al., 1989). R is the universal gas constant. Reaction 1 is endothermic and reaction 2-4 are exothermic. Reaction 1 dominates such that the overall reaction is endothermic.
The structure of the reformer model is shown in figure  2. The liquid feed S 1 is mixed with recycle gas S 55 . The resulting vapor S 2 and liquid S 3 outlet stream are preheated in the reactor effluent heat exchanger E1 and then enter the first heater and reactor section. The heaters are modeled using direct heat input and each of the four reactors is modeled using ten CSTRs in series with even distribution of catalyst. Heat exchanger E2 and separator D1 is modeled using the same flash unit model .
In addition, variables and equations for the reformate octane number (RON), R1 feed hydrogen to hydrocarbon ratio, and some mass flows are added as internal variables in a "dummy" unit model. The mass flows are for the feed, reformate, gas and hydrogen products and recycle gas.

Process model
The model equations are organized in a unit model framework (Lid and Skogestad, 2007). For the CSTR elements the mass balances, energy balance, mole fraction summation and pressure drop relationship is written where the process stream variables x, T , P and F represents the molar composition, temperature, pressure and molar flow respectively. The CSTR inlet and outlet streams are in this case marked with subscript 1 and 2. In addition m c is the mass of catalyst, and A c is a catalyst activity parameter.
This gives N C + 3 equations for each reactor element. Similar model are formulated for the other units; heater, separator with cooling, compressor, heat exchanger, stream mix and stream split. For details, together with thermodynamics data (enthalpy, entropy, vapor-liquid equilibria), the reader is referred to the thesis of Lid (2007).
The resulting model and specifications are written As seen from Tables 1 and 2, the model f (z) = 0 contains n z = 501 variables z and n f = 442 equations. The first requirement for a unique solution is that n z − n f = 59 variables are specified. These specifications are added as n s = 59 rows in A s with the corresponding specification values in b s . Table 3 lists 23 of the specifications. The remaning 36 come from the catalyst efficiency factors for the CSTRs wich are assumed equal within one reactor. This is incorporated as 36 linear constraints in A s .   (Matlab, 2000). In order to reduce the computational load in solving the model, the first order derivatives where calculated analytically.
Unit model      Figure 3 shows for a typical case the molar flows of each component in the four reactors as a function of the normalized catalyst mass. There is a net production of hydrogen and gas. The largest amount of hydrogen is produced in reactor one and the largest amount of gas is produced in reactor four. The main reaction in reactor number one is conversion of naphthenes to aromatics. The main reaction in reactor number four is conversion of paraffines to naphthenes. The large temperature drop in reactor one is due to the large heat of reaction required for the conversion of naphthenes to aromatics.

Nominal operation
Other key variables like heater duties and product yields are listed in

Data reconciliation results
In the data reconciliation we want to estimate the 23 remaining degrees of freedom (rather than specifying them as we did in the simulation case in table 3).
The naphtha reformer process has n y = 26 measured values. These are from the feed, product and recycle gas analyzers, feed product and recycle gas mass flow measurements and various temperature measurements. All the measurements are listed in table 5. The values for the standard deviations are based on typical measurement uncertainties. For flow measurements the uncertainty are assumed to be 3% of the measurement range. For temperature measurements a fixed value of 3 • C is assumed. The standard deviation for the analyzers of 1% are based on instrument specifications except for the recycle gas H 2 analyzer which has a higher standard deviation (10%) due to a large modeling error in this section (see discussion section).
The feed hydrogen and gas content is known to be almost zero and specifications x 1 (1) = 0 and x 1 (2) = 0 are added in the linear constraints A r of the data reconciliation problem in equation (1). The remaining degrees of freedom then equal 21.
The observability of all variables, given the process model (f (z) = 0), linear constraints and specifications (A r ) and measurements (U ), is verified by the rank of When Γ has full column rank (equal to the number of variables n z = 501) the values of all variables are observable (Stanley and Mah, 1981). In this case the rank of Γ equals 498, which indicates that there are three unobservable variables.
One of these is the condenser liquid outlet pressure, which needs to be specified, as the liquid stream is not connected to any downstream units. In addition, there are no measurements of the cooling water inlet or outlet flow or temperature. Thus, in order to make all variables observable the values of P 52 , F CW and T CWi where specified by adding three linear constraints in A r . The degrees of freedom are now reduced from 21 to 18.
It is verified, using the definition of redundancy in Crowe (1989), that all measurements in the reformer process are redundant.
Data reconciliation using equation (1) and (4) was ap-plied to 21 data sets from the plant collected over a period of two years. The results are given in figures 4-7 and detailed results for data set no. 12 are shown in table 5.The uncertainty of the estimated values are computed using the method from Romagnoli and Stephanopoulos (1981) and are shown in table 5.
There is almost no reduction of uncertainty in the estimate of the reactor inlet or outlet temperatures, compared with the uncertainty of the measured values. This is probably because there is in practice little redundancy in the reactor section measurements (only inlet and outlet temperatures are measured). The feed (F 1 ) and product mass flow (F 52 ) uncertainty is reduced by approximately 30%. The compressor inlet temperature (T 54 ), separator outlet temperature (T 52 ) and in particular the recycle gas hydrogen content (x 54 (1)) has a large reduction of uncertainty. This is probably due to the oversimplification in the modeling of the separator and recycle gas system (i.e. model error).
The values and standard deviations of the heat exchanger heat transfer coefficients and reactor and compressor efficiency are shown in table 6. On average the uncertainties in these variables are 10-35% of the actual value except for the estimate of U 2 . The estimated uncertainty in U 2 shows that this variable is not practically observable and indeed the estimate of U 2 = 200W/m 2 /K is equal to its initial value.  Gross errors (non-zero bias) according to the criterion given in Tjoa and Biegler (1991) are detected for the measured values marked with * . For data set 12 we detect gross error for reactor 1 outlet (T 15 ), reactor 4 outlet (T 48 ) and E1 hot side outlet temperature (T 50 ).
The two latter (T 48 and T 50 ) have a gross error detected in all 21 data sets. The outlet temperatures of reactor 1 has gross errors detected in data sets 12 and 13 and the outlet temperature of reactor 4 has gross errors detected in data sets 14. The compressor mass flow has a gross error detected in three data sets and the feed temperature has a gross error detected in one data set. inlet temperatures for all 21 data sets. The adjustments of the catalyst efficiency factors contribute to an almost perfect fit to the measured data. We have the highest reaction rate, and thus the highest influence on the other measured values, at the inlet of the reactor and this may be one reason why the error in temperature drop over each reactor is assigned to the reactor outlet temperatures. There are large predicted measurement errors in the reactor outlet temperatures, as shown in figure 5. The outlet temperature of reactor one and two have gross errors in most data sets but some data points have almost zero measurement error. The outlet temperature of reactor number four has an almost fixed bias in all data sets. As a curiosity, the outlet temperature of reactor three is "accepted" as an untrustworthy measurement at the refinery. However, this is not supported by our results which show close to zero measurement error in all data points.
The estimated catalyst efficiencies for all data sets are    Ideally, the catalyst efficiency factors A c should be close to one in all data sets but due to variation in the catalyst circulation some changes in A c are expected.
In periods, where the catalyst regenerator is shut down, the unit may run for several days with no catalyst circulation . In these periods the catalyst efficiency will decrease due to coke build up on the catalyst.
The values of A c show large deviations in excess of 1 in data points 5, 10, 17 and 19. There is no clear reason for this and the data at these points does not differ significantly from the others. An observation is that the measurement error of reactor one outlet temperature is almost zero at these points but this is also true for data points 1, 2, 3 and 14.
From figure 7, we find the average deviation between the measured and reconciled values for the mass flows of feed, reformate and gas are 0.7t/h, -1.93t/h and 1.59t/h respectively. The average deviation for octane is -0.25. The reconciled gas mass flow is persistently lower than the measured value and even if no gross errors where detected in the measured value the presence of a systematic error is clear.
where J(z) = −p(z) T z. In our case p is a vector of fixed prices of feed, products and utilities, see table 7.
Fixed variables include feed data (composition and temperature), heat transfer coefficients and compressor efficiency and are set equal to their reconciled values using linear equality constraints A opt z = b opt in (31).
Operating constraints like maximum feed flow, maximum pressure, maximum temperature and minimum product octane are added as upper and lower bounds on the variables in z opt min and z opt max , see table 8.
The naphtha reformer is the main producer of hydrogen at the refinery and may not be shut down even if the product price is low and the unit profit is negative. Thus, to secure the availability of hydrogen a lower bound is added on the reformer unit hydrogen production.
The number of degrees of freedom for the optimization is n z − n f − n opt = 7. This follows because the number of variables is n z = 501, the number of equations is n f = 442 and the number of rows (fixed values) in A opt is n opt = 52. Specifically, the 52 specified (fixed) values added in A opt are 40 catalyst efficiency factors, 2 heat exchanger heat transfer coefficients, compressor efficiency, feed temperature and N C = 5 feed mole fractions, reformate outlet pressure, cooling water flow and cooling water inlet temperature. Note that the feed rate is not specified so its optimal value is obtained as part of the optimization.

Optimization results
Two operational cases, which both are common operational regimes for a naphtha reformer unit in a refinery, are analyzed.
• Case 1. The product (reformate) price is high and throughput (feedrate) is maximized, subject to satisfying constraints.
• Case 2. The product price is low and throughput is minimized subject to meeting the production demand on hydrogen.
The detailed results from the optimization for case 1 and 2 are shown in table 8. In both cases the minimum reformat RON of 103 is an active constraint. This is expected because reformate is the most valuable product of the three, and we want to avoid "give away". The maximum separator pressure of 10 bar and the minimum H2/HC ratio of 3 in reactor 1 are also active constraints in both cases.
In case 1, the operation is in addition constrained by the maximum heater duties. The improvement in profit, compared to the reconciled solution, is 245$/h (2.1 × 10 6 $/year). This comes as a result of an increased feed flow, and a reformate yield improvement of 0.43%. The yield improvement is mainly due to reduced temperatures in the reactors and reduced reformate RON.
In case 2, the operation is in addition constrained by the minimum hydrogen product mass flow of 3.5t/h.

Implementation of optimal operation
In order to operate the process optimally the seven degrees of freedom have to be specified or fixed. These specifications are implemented as controlled variables. The degrees of freedom can be thought of being related to the heat input to the four heaters, the feed, the compressor work (recycle flow) and the H 2 product flow (purge). The basic control layer includes heater duty control, feed flow control and pressure control.
In case 1, there are seven active constraints and implementation is obvious: the seven active constraints are selected as controlled variables.
In case 2, there are four active constraints and these are selected as controlled variables. It is less obvious what to select as controlled variables for the remaining three unconstrained degrees of freedom. The problem is that the optimal value of the unconstrained variables depend on the disturbance, and also that there is a implementation error associated with control of the unconstrained variables (Skogestad, 2000). The objective is to find "self-optimizing" control variables which are insensitive to disturbances and control errors, that is, which result in a small economic loss. A closer analysis shows that the optimal variation in the inlet temperatures to the for reactors (which are between 780.4K and 798.8K in case 2) are not important. In fact specifying that the four reactor inlet temperatures to be equal (which corresponds to adding three specifications) only marginally decreases the profit by 0.005$/h. This is shown by the column "Case 2 (same T)" in table 8.This is also consistent with the equal marginal values of the heater duties in case 1 shown in table 9.
In summary, "self optimizing control" is achieved by adding three reactor difference temperatures as controlled variables with a zero set point. The actual reactor inlet temperatures will be indirectly determined by the four active constraints.