The Rutgers Scholar
An Electronic Bulletin of Undergraduate Research

Reduction and characterization of cyclohexane oxidation

Rinku Parikh* and Marianthi Ierapetritou1
Department of Chemical and Biochemical Engineering,
Rutgers University, Piscataway New Jersey 08854

1Author to whom correspondence should be addressed. Tel: 732.445.2971
*Rutgers Undergraduate Research Fellow

Keywords: Kinetic model reduction, Feasibility and flexibility, Range of validity


A detailed reaction network for the liquid-phase oxidation of cyclohexane has been proposed in the literature. The aim of this work is to select from the detailed kinetic mechanism a reduced set of reactions that captures the essential features of the reactive process within a user-defined accuracy with the aid of mathematical programming techniques. The reduced model is characterized with respect to the range of initial conditions using the ideas of flexibility analysis.


The kinetic mechanisms that are required to describe reaction processes in the fields of combustion, polymerization, environmental science, materials technology, and microelectronics processing are characterized by a large numbers of species (NS) and elementary reactions (NR). Due to the increased complexity of the reaction mechanism, the computational time to solve these problems, especially when complex fluid dynamics are considered, becomes substantial. Therefore, methods involving mechanism reduction have become an important issue (Androulakis, 2000).

This paper is concerned with the reaction model of liquid-phase oxidation of cyclohexane. This reaction is part of a larger industrial process known as the CYCLOPOL process, which produces caprolactum, adipic acid, and polyamide fibers and plastics (i.e. nylon 6 nylon 66). The cyclohexane mechanism is a fairly complicated process given that the products are also intermediates in the series of reactions. The mechanism of cyclohexane oxidation is a multistage, free-radical chain reaction with degenerated chain, comprising different chain inhibition, chain propagation and chain termination steps (Ryszard, 2001).

One reaction mechanism containing 154 reactions was constructed by Tolman (Pohorecki et al. 2001). Another model containing 19 reactions and 10 species was developed by Kharkova et al. (1989). This work attempts to further reduce the mechanism with the aid of mathematical programming techniques. Once a reduced mechanism is found, flexibility analysis ideas are utilized to determine the range of validity in terms of nonisothermal conditions where the mechanism is valid.

The paper is organized as follows. In sections 2 and 3, brief reviews of the reduction methodology and of the flexibility approach are presented. The results obtained are shown in section 4 for the cyclohexane oxidation mechanism, while section 5 summarizes the work and outlines our continuing research direction in this project.

Kinetic model reduction approach

The methodology used for the reduction of the cyclohexane mechanism is the one proposed by Androulakis (2000) and extended by Sirdeshpande et al. (2000). The technique is based on eliminating reactions from the full set shown in Table 1, without allowing the predictions of the reduced mechanism to diverge from the full model. The measure of variance is based on a user-defined tolerance. The first stage involves selecting a reactor type and a discrepancy function. The discrepancy function is a measure of the error in eliminating reactions from the full mechanism. A user defined tolerance is selected for the system and the optimization problem of minimizing the number of reactions is solved to give the reduced mechanism. The tolerance is defined as the acceptable value of the discrepancy function.

Table 1: Sequence of reactions and rate constants, k, for oxidation of cyclohexane (Kharkova et at., 1989)

No. Reaction @ 423 K
1 RH + O2 --> R* + HO2* 1.80E-08
2 R* + O2 --> RO2* 1.00E+09
3 RH + RO2* --> ROOH + R* 12
4 RH + RO* --> ROH + R* 7.65E+07
5 ROOH --> RO* + OH* 5.07E-05
6 ROOH + ROH --> 2RO* + H2O 1.91E-03
7 ROOH + R=O --> 0.7R=O + 0.3P' + RO* + OH* 9.47E-04
8 ROOH + RH --> RO* + R* + H20 4.52E-06
9 ROOH + RO* --> ROH + RO2* 7.65E+07
10 ROOH + RO2* --> R=O + OH* + ROOH 7.44
11 ROH + RO2* --> R=O + RO2* 4.56
12 R=O + RO2* --> P' + RO2* 8.32
13 ROOH + R* --> ROH + RO* 4.59E+04
14 ROOH + R* --> RO2* + RH 1.64E+06
15 2RO2* --> R=O + ROH + O2 1.55E+06
16 ROOH + P' --> R=O + H2O + P' 2.18E-03
17 P' + ROH --> P 1.67E-03
18 2RO2* --> P' 1.55E+06
19 RH + OH* --> R* + H2O 7.65E+08

Reactor Model

For our case, the reactor was assumed to be a perfectly-mixed stirred tank reactor under isothermal conditions. Therefore, the key variables are the concentration of the species. The initial conditions were chosen to be the same as those presented by Kharkova (1989). The set of differential equations is shown in Table 2.

Table 2: System of kinetic equations
(Kharkova et at., 1989)

r1 = k1[RH][O2] r11 = k11[ROH][RO2*]
r2 = k2[R*][O2] r12 = k12[R=O][RO2*]
r3 = k3[RH][RO2*] r13 = k13[ROOH][R*]
r4 = k4[RH][RO*] r14 = k14[ROOH][R*]
r5 = k5[ROOH] r15 = k15[RO2*]2
r6 = k6[ROOH][ROH] r16 = k16[ROOH][P']
r7 = k7[ROOH][R=O] r17 = k17[ROH][P']
r8 = k8[ROOH][RH] r18 = k18[RO2*]2
r9 = k9[ROOH][RO*] r19 = k19[RH][OH]
r10 = k10[ROOH][RO2*]  

d[RH]/dt = -r1 - r3 - r4 - r8 - r19 + r14

d[ROOH]/dt = -r5 - r6 -r7 - r8 - r9 - r13 - r14 - r16 + r3

d[ROH]/dt = -r6 - r18 - r17 + r4 + r9 + r13 + r15

d[R=O]/dt = -r12 - 0.3r7 + r10 + r18 + r15 + r16

d[P']/dt = r12 + 0.3r7 + r18 - r17

d[RO2*]/dt = r2 + r9 + r14 - r3 - r10 - 2(r15+r18)

d[RO*]/dt = r5 + 2r6 + r7 + r8 + r13 - r4 - r9

d[R*]/dt = r1 + r3 + r4 + r8 - r2 - r13 - r14

d[OH]/dt = r5 + r7 + r10 - r19

d[P]/dt = r17

Discrepancy function

The discrepancy function used to measure the error in eliminating reactions is the integral of the square error (ISE). The following is the representation of the ISE:


where M is the number of matched species, wj the weight assigned to the jth quantity, and Ej is the time dependent error between full and reduced models. The common form of Ej is a scaled residual as the following:


The u represents the value of the matched quantity, where the superscript f and r, stand for full and reduced respectively, and s is a scale factor defined as:


corresponding to no scaling, scaling the residuals of an observable with the maximum value of that observable over the entire time horizon for the full mechanism, and scaling residuals of an observable with the corresponding value of the full mechanism at each time step. k is a lower limit (typically 10-15) on s and avoids scaling by extremely small quantities (Sirdeshpande et al., 2000).

User-defined tolerance

The user-defined tolerance, which is represented by e, characterizes the degree to which the reduced mechanism should match the full mechanism for the quantities of interest. The quantities of interest are the molar concentration profiles for the matched species. The tolerance is chosen carefully in order to capture the essential features of the reactive process. If the tolerance is too small, which means the reduced model must match the full model very closely, reduction may not occur and all the reactions must be included in the reduced set. However, if the tolerance is too large, then the reduced set found may not adequately express the reactive process of the full model. Therefore, a tolerance value somewhere in between must be selected. Practically, the tolerance is selected by starting with a large value and sequentially decreasing it until the required reduction and error is achieved. Note that the matched quantities are dependent on the system under consideration. For example, for combustion systems, the ignition time and the pollutant concentration are of greatest interest.

Optimization problem

The reduced mechanism is obtained by minimizing the objective function, which is represented by the sum of the binary variables, li, for each reaction. If the reaction is present in the reduced mechanism, l is denoted as one. However, if the reaction is absent from the reduced set, the value for l corresponds to zero. Therefore, the lower bound of the objective function is equal to zero meaning no reactions are present, and has an upper bound equal to NR, implying the full mechanism is present. The objective function:


is subject to the following constraint:


The constraint implies that the discrepancy function, c(l), must be less than or equal to the user-defined tolerance, e.

The overall problem of mechanism reduction is an integer nonlinear programming (INLP) problem. The INLP was solved using a branch-and-bound (B&B) strategy (Floudas, 1995). The relaxed INLP at each node of the branch-and-bound tree was solved using a sequential quadratic programming (SQP) code. Both NPSOL (Gill et al. 1990) and FFSQP (Zhou et al. 1998) have been used.

Flexibility analysis

One of the main limitations of most of the existing reduction techniques is the fact that reduction is performed at specific initial conditions. Thus, there is no guarantee that the reduced model will perform well when different conditions are utilized. In this paper, the ideas of flexibility analysis will be utilized to determine the range of conditions where the reduced model is valid away from the nominal point where the reduction is performed (Balakrishnan et. al. 2002).

Let the initial conditions for ROOH and temperature be denoted by q1 and q2, respectively. Next, assume that Dqi+ and Dqi- are the maximum expected deviation from there nominal values, qiN . These choices lead to the expected range for each parameter as:


Now, in order to check the percentage of the expected range where the mechanism is valid, the delta factor must be incorporated into the inequality. This incorporation gives rise to the following range:


Therefore, the solution to the problem lies in finding the value of d. The flexibility index F for a given reduced mechanism set {lr} is the value of d that defines the critical parameter range responsible for the error function c to be equal to its desired value e. In this paper, F is determined by the vertex enumeration strategy (VES).

Vertex enumeration strategy

The vertex enumeration strategy solves four problems, one at each vertex of the parameter range. The coordinates of the vertices Vk for a two-parameter problem are


The value of delta can be defined as a scaled step length along each vertex direction. Moreover, it can be considered the fraction that can be traveled along each vertex direction until the error constraint is violated. Therefore at each individual vertex, the following problem is solved:


subject to:


Finally, the flexibility index is the defined as:


Therefore, the flexibility index, F, is the minimum of the four delta values. This is the value that will define the flexibility region. Note that the minimum over the vertices is selected to guarantee that the scaled hyper-rectangular box will be inscribed within the feasible region (Swaney and Grossmann, 1985).

Inscribed convex hull

Instead of using the minimum value of delta to construct the flexibility area, a different approach can be used to determine the inscribed convex hull region (Ierapetritou, 2001). This approach solves a number of additional problems equal to the number of uncertainty parameters along the directions of the axes. Then, these points are used to determine the inscribed convex region as shown in Figures 11 and 12 for the cyclohexane case study.

Feasibility analysis

In order to evaluate the results of the flexibility analysis, feasible region plots are generated by a grid search method. With the cyclohexane mechanism, the variables of interest are the initial amount of ROOH (cyclohexyl hydroperoxide), and the initial temperature. The reduced mechanism was generated with no initial amount of ROOH and constant temperature. Checking the feasibility at each point, the region over which the reduced mechanism is valid in terms of initial concentration and temperature is determined.

2500 equally-spaced points were generated in each parameter range. Then, the error constraint was evaluated at each of these individual points. Every point was classified as either pass or fail depending on whether the error constraint was satisfied or not.

Results for Cyclohexane Oxidation

Results from reducing the number of reactions included in the mechanism

The initial conditions used for the reduction are:

The watched species are the species that are considered in the discrepancy function and are the quantities that are to be matched. For the cyclohexane oxidation these are: [RH] (cyclohexane), [ROH] (cyclohexanol), [R=O] (cyclohexanone), and [P] by-product.

Table 3: Effect of user-defined tolerance to reduced set

e (user defined tolerance) Reaction Set NR Ns
1 1-12; 14-19 18 10
10 1,2,3,4,6,7,9,11,12,15,16,17,18,19 14 10
28 1,2,3,4,6,9,11,12,14,15,17,18 12 9
40 1,2,3,4,6,11,12,15,16,17,18 11 9
71 1,2,3,4,6,12,15,16,17,18 10 9
128 1,2,3,4,6,15,17,18 8 9

The results of the reduction optimization are summarized in Table 3 for different values of the tolerance parameter (e). As expected, as the tolerance increases, the number of reactions found in the reduced set decreases. However, a species reduction only occurs when the tolerance is greater than or equal to 28. At this value of e, it is noticed that the OH* disappears from the reduced set. The reduced set obtained from a tolerance value of 71 is the reduced model that is used in this paper. As shown in figures 4 and 5, further increase in the tolerance value resulted in unacceptable agreement with the full model. Figures 1, 2, and 3 show the concentration profiles for [RH], [ROH], [R=O], [P'], [ROOH], and [P].

Although, [ROOH] and [P] are not among the species involved in the constraint of the discrepancy function, the reduced model predicts the behavior of these species, too (Figure 3). This is of extreme importance since it illustrates the predictive capability of the reduced model.

Profile comparison for watched species

Figure 1: [RH] and [ROH] profiles for both reduced vs. full model for e = 71.

Figure 1

Figure 2: [R=O] and [P'] profiles for both reduced vs. full model for e = 71.

Figure 2

Profile comparison for unwatched species

Figure 3: [ROOH] and [P] profiles for both reduced vs. full model for e = 71.

Figure 3

Profiles for reduced set of larger tolerance

Figure 4: [RH] and [P] profiles for both reduced vs. full model for e = 128; T0=403 K. Reduced set consists of 9 reactions and 9 species.

Figure 4

Figure 5: [ROOH] and [R=O] profiles for both reduced vs. full model for e = 128; T0=403 K. Reduced set consists of 9 reactions and 9 species.

Figure 5

Figures 4 and 5 above demonstrate the effect in the profiles for the reduced model versus the full model when a large tolerance is allowed in the system. For example, from figure 4 the species profile for P in the reduced model decreases while the full model shows an increase in concentration. This observation reinforces the significance of the user defined tolerance.

Flexibility results

Figures 6 shows the feasible region (i.e. the region in which the reduced model produces acceptable results) corresponding to the reduced set of 10 reactions and 9 species for the cyclohexane problem. The grid search was performed for initial [ROOH] molar concentration in the range {0.0, 1.8} mole/L and initial temperatures in the range {373, 474} K.

Figure 6: Feasible region for reduced set of 10 reactions and 9 species corresponding to e = 71.

Figure 6

Figure 7: Feasible region for the 10 reaction and 9 species set when the tolerance, e, is increased from 71 (figure 6) to 85.

Figure 7

Flexibility Region

The following results correspond to the values of the d obtained in solving the problem at each of the four vertices for the expected range.

Figure 8: Feasible region for reduced set of 10 reactions and 9 species with different nominal point.

Figure 8

V1 = (391.75, 0.0) --> d = 0.625

V2 = (392, 0.992) --> d = 0.620

V3 = (457.3, 0.0) --> d = 0.686

V4 = (449.6, 0.8512) --> d = 0.532

The first point, labeled V1, corresponds to the vertex point in figure 9, with a d value of 0.625. Therefore, the coordinates of the point are obtained from the following relations:

V1x = 423 - (50) * d = 391.75 K

V1y = 0.0 - (0.0) * d = 0.0 mol/L

This is the same procedure used to determine the coordinates of the other three vertices. However, the flexibility region is defined by the minimum of these four d values, in this case 0.532. This factor is used to determine the coordinates of the flexibility box in Figure 9.

Note that the size of the flexibility region as defined in section 3, depends on the choice of the nominal point, which represents the point at which reduction is performed. Figures 9 and 10 show different flexibility regions for different nominal points of the same reduced model consisting of 10 reactions and 9 species. This change in the shape of the flexibility region occurs because the flexibility index, which defines the general shape of the flexibility box, is different for the two cases. In particular for Figure 10, the smallest value for d occurs along the vertex direction towards the point (0.0, 373) not along the vertex direction towards the point (1.6, 473) as in Figure 9.

Figure 9: Flexibility region for reduced set of 10 reactions and 9 species 0.0<= [ROOH]
[ROOH]<= 1.6 mol/L; 373 <= T0 <= 473 K, qN1 = 0.0 mol/L; qN2 = 423 K.

Figure 9

Figure 10: Flexibility region for reduced set of 10 reactions and 9 species with different nominal point 0.0<= [ROOH] <= 1.6 mol/L; 373 <= T0 <= 473 K, qN1 = 0.0 mol/L; qN2 = 403 K.

Figure 10

Figure 11: Inscribed Convex Hull for reduced set of 10 reactions and 9 species 0.0<= [ROOH] <= 1.6 mol/L; 373 <= T0 <= 473 K, qN1 = 0.0 mol/L; qN2 = 423 K.

Figure 11

Inscribed convex hull

Using the approach in section 3, a convex hull is inscribed within the range of validity of a given reduced model as shown is figure 11 and 12 for the two nominal points used in the flexibility analysis in the previous subsection. Note that the area predicted following this approach corresponds to a better representation of the actual region compared with the flexibility region shown in figure 9 and 10.

Figure 12: Inscribed Convex Hull for reduced set of 10 reactions and 9 species with different nominal point 0.0<= [ROOH] <= 1.6 mol/L; 373 <= T0 <= 473 K, qN1 = 0.0 mol/L; qN2 = 403 K.

Figure 12

Summary and future directions

A new model for liquid-phase oxidation of cyclohexane was developed in this paper. A mathematical programming technique approach was utilized to derive a more compact system involving only the key reactions and species. The reduced model was then characterized using the flexibility analysis approach that resulted in the determination of the range of validity for the reduced model in terms of initial conditions.

The feasibility plots showed that the system displays non-convex characteristics. Therefore, flexibility analysis and the inscribed convex hull approach will not always compute a region that lies within the calculated feasible solution (see figure 12). For example, different nominal conditions may produce flexibility regions that lie outside this desired feasible boundary. Thus, research is underway to investigate efficient approaches that can lead to the determination and characterization of the nonconvex regions (Goyal and Ierapetritou, 2002)


Financial support from the National Science Foundation is gratefully acknowledged (Grant # 9983406).


Androulakis, I. P., Kinetic Mechanism Reduction Based on an Integer Programming Approach. AIChE J., 46, 361, 2000.

Balakrishnan S., I. Banarjee and M.G. Ierapetritou, Uncertainty Considerations in the Complex Kinetic Mechanisms Reduction. Accepted for publication AIChE J, 2002.

Floudas, C. A., Nonlinear and Mixed-Integer Optimization, Oxford University Press, 1995.

Gill, P.E., W. Murray, M. A. Saunders, and M. H. Wright, "User's Guide for NPSOL (Version 4.0): A Fortran Package for Nonlinear Programming," Technical Reports SOL 86-2, 1986, Systems Optimization Laboratory, Dept. of Operations Research, Stanford University, Stanford, CA, 1986.

Goyal, V. and M.G. Ierapetritou, Determination of Operability Limits Using Simplicial Approximation. Accepted for publication. AIChE J, 2002.

Ierapetritou, M.G, A New Approach for Quantifying Process Feasibility: Convex and one Dimensional Quasi-Convex Regions. AIChE J., 47, 1407, 2001.

Ierapetritou, M. G., and I. P. Androulakis, Uncertainty Considerations in the Reduction of Chemical Reaction Mechanisms. International Symposium on Foundations of Computer-Aided Process Design, Brickenridge, CO, 1999.

Kharkova, T.V., Arest-Yakubovich, I.L., and Lipes, V. V., Kinetic model of the liquid-phase oxidation of cyclohexane. I. Homogeneous proceedings of the process. Kinetics and Catalysis, 30, 954-958, 1989.

Ryszard P., Jerzy B., Wladyslaw M., Wioletta P., Artur Z., and Piotr W., Kinetic model of cyclohexane oxidation. Chem. Eng. Sci., 56, 1285-1391, 2001.

Sirdeshpande, A.R., M.G. Ierapetritou, and I.P. Androulakis, Design of Flexible Reduced Kinetic Mechanisms. AIChE J., 46 (2), 361, 2000.

Swaney, R. E., and I. E. Grossmann, An Index for Operational Flexibility in Chemical Process Design. Part 1-Formulation and Theory. AIChE J., 31, 621, 1985a.

Swaney, R.E., and I.E. Grossmann, An Index for Operational Flexibility in Chemical Process Design. Part 2-Computational Algorithms. AIChE J., 31, 631, 1985b.

Pohorecki, R., J. Baldyga, W. Moniuk, W. Podgorska, A. Zdrojkowski, P.T. Andrzchowski, Kinetic Model of Cyclohexane Oxidation. Chem. Eng. Sci., 56,1285-1291, 2001.

Zhou, J. L., A. L. Tits, C. T. Lawrence, "User's Guide for FFSQP Version 3.7: A Fortran Code for Solving Constrained Nonlinear (Minimax) Optimization Problems, Generating Iterates Satisfying All Inequality and Linear Constraints," Electrical Engineering Dept. and Institute for Systems Research, University of Maryland, College Park, MD, 1998.

Copyright 2002 by Marianthi Ierapetritou
Current URL: