Calibration-free estimates of batch process yields and detection of process upsets using in situ spectroscopic measurements and nonisothermal kinetic models: 4-(Dimethylamino)pyridine-catalyzed esterification of butanol

Article abstract of DOI:10.1021/ac035356i

The use of an near-IR fiber-optic spectrometer with a high-speed diode array for calibration-free monitoring and modeling of the reaction of acetic anhydride with butanol using the catalyst 4-(dimethylamino)pyridine in a microscale batch reactor was presented. Nonlinear fitting of a first-principles model directly to the reaction spectra gave calibration-free estimates of time-dependent concentration profiles and pure component spectra. Real-time modeling of a batch reaction could be achieved by augmenting prerun batch data with recently acquired spectral data followed by fitting of a multiway kinetic model. In real-time applications, accurate initial guesses of the parameter estimates from prerun data can be used to facilitate rapid convergence of the iterative nonlinear fitting process.

Full text of DOI:10.1021/ac035356i

Anal. Chem. 2004, 76, 2575-2582
Ca lib ra t io n -Fre e Es t im a t e s o f Ba t c h P ro c e s s Yie ld s
a n d De t e c t io n o f P ro c e s s Up s e t s Us in g in S it u
S p e c t ro s c o p ic Me a s u re m e n t s a n d No n is o t h e rm a l
Kin e t ic Mo d e ls : 4 -(Dim e t h yla m in o )p yrid in e -
Ca t a lyze d Es t e rific a t io n o f Bu t a n o l
Pa ul Ge m pe rline ,*,† Gra e m e Puxty,‡ Ma rc e l Ma e de r,‡ Dw ight Wa lke r,§ Fra nk Ta rc zyns ki,§ a nd
Ma ry Bos s e rm a n†
Department of Chemistry, East Carolina University, Greenville, North Carolina 27858, Department of Chemistry, University of
Newcastle, Newcastle, NSW, Australia, and GlaxoSmithKline, 5 Moore Drive, Research Triangle Park, North Carolina 27709
streams and batches over a large spectral range with high signal-
to-noise ratios. The large volume of data produced by these
instruments has typically been used to develop multivariate
calibration models using techniques such as principal component
regression1 and partial least squares (PLS).2-4 The major goals
of these types of applications are to detect and diagnose process
upsets and monitor product quality and yield.5 In some cases, the
output of the calibration model is used for feedback control.5 One
drawback of the multivariate calibration approach is that a
significant amount of effort is required to develop and maintain
the calibration model. For example, instrument maintenance
events (lamp changes, fiber replacement, etc.) or process changes
(change in raw material suppliers or product formulation) may
degrade or even invalidate a calibration model, requiring a
significant amount of effort to update or revise it.
For batch chemical processes, fitting of a first-principles kinetic
model offers a calibration-free approach to monitoring and
modeling that does not require the significant effort of multivariate
calibration. In this approach, a reaction mechanism is postulated
giving a system of simultaneous ordinary differential equations
(ODEs). Using known initial conditions, the ODEs are numerically
integrated to estimate time-dependent concentration profiles. The
concentration profiles are fitted directly to the reaction spectra
by nonlinear optimization of the rate constants. Neither calibration
spectra nor off-line reference measurements are required in the
model-fitting process. For process analysis applications, ap-
proximate model parameters, in the form of rate constants can
be determined from prerun batches. During subsequent batch
runs, the prerun batch spectra can be augmented with the new
batch spectra as they are acquired. The model parameters can
then be updated in real time by nonlinear fitting methods and
used to monitor the progress of the batch, detect process upsets,
and forecast batch end points, product yield, and product quality.
In this paper, we report the use of an NIR fiber-optic
spectrometer with a high-speed diode array for calibration-
free monitoring and modeling of the reaction of acetic
anhydride with butanol using the catalyst 4-(dimethylami-
no)pyridine in a microscale batch reactor. Acquisition of
spectra at 5 ms/ scan gave information relevant for model-
ing these fast batch processes with a single multibatch
kinetic model. Nonlinear fitting of a first-principles model
directly to the reaction spectra gave calibration-free esti-
mates of time-dependent concentration profiles and pure
component spectra. The amount of catalyst was varied
between different batches to permit accurate estimation
of its effect in the multiway model. A wide range of
different models with increasing complexity could be fit
to each batch individually with low residuals and apparent
low lack of fit. However, only one model properly esti-
mated the concentration profiles when all five batches
were fitted simultaneously in a multiway kinetic model.
Inclusion of on-line temperature measurements and use
of an Arrhenius model for the estimated rate constant gave
significantly improved model fits compared to an isother-
mal kinetic model. Augmentation of prerun batches with
data from an additional batch permitted model-based
forecasts of reaction trajectories, reaction yield, reaction
end points, and process upsets. One batch with added
water to simulate a process upset was easily detected by
the calibration free process model.
During the last 20 years, a dramatic increase has been
observed in the use of process analytical chemistry (PAC) to
monitor and control modern chemical processes. PAC techniques
employ on-line sensors to measure physical properties of pro-
cesses in real time. There has been an increasing trend toward
the use of multivariate sensors, such as fiber-optic spectrometers,
capable of rapidly measuring the optical properties of process
(1) Wise, B. M.; Gallagher, N. B. J. Process Control 1 9 9 6 , 6, 329-348.
(2) Kourti, T.; MacGregor, J. F. Chemom. Intell. Lab. Syst. 1 9 9 5 , 28, 3-21.
(3) Nomikos, P.; MacGregor, J. F. Chemom. Intell. Lab. Syst. 1 9 9 5 , 30, 97-
108.
(4) Burnham, A. J.; Viveros, R.; MacGregor, J. F. J. Chemom. 1 9 9 6 , 10, 31-
45.
* To
whom
correspondence
should
be
addressed.
E-mail:
gemperlinep@mail.ecu.edu.
† East Carolina University.
‡ University of Newcastle.
§ GlaxoSmithKline.
(5) Wise, B. M.; Ricker, N. L. AIChE Symp. Ser. 1 9 8 9 , 85, 19-23.
10.1021/ac035356i CCC: $27.50 © 2004 American Chemical Society
Published on Web 03/27/2004
Analytical Chemistry, Vol. 76, No. 9, May 1, 2004 2575

Changes in instrument response due to maintenance events such
as lamp changes are inconsequential as long as the instrument
response continues to provide sufficient signal to accurately and
precisely update parameter estimates.
decomposition does not represent Y exactly. We define a residuals
matrix R that is the difference between the measurements and
the model. The process of nonlinear least-squares fitting gives
estimates Cˆ and Aˆ that best represent Y by minimizing the sum
of squares of the residuals matrix Rˆ .
Traditional near-infrared (NIR) measurements have found
widespread use in process analytical application because of the
method’s good sensitivity, high information content, and low noise.
Commercially available scanning instruments or Fourier transform
(FT-NIR) instruments require long scan times. Some batch
processes take place on a much shorter time frame. For example,
4-(dimethylamino)pyridine (DMAP) is one of the best-known
catalysts available for acetylation with acetic anhydride, and
commercially available NIR scanning or FT instruments may be
too slow to give measurements relevant for estimation of rate
constants. Recent advances in semiconductor technology has
made possible indium gallium arsenide charge-coupled device
diode arrays with excellent sensitivity in the range from 1100 to
2200 nm and high scan speeds of 5 ms.6
In this paper, we report the use of an NIR fiber-optic
spectrometer with a high-speed diode array for monitoring and
modeling the reaction of acetic anhydride and butanol with DMAP
as the catalyst in a microscale batch reactor. The spectrometer is
capable of acquiring spectra at the rate of 5 ms/ scan, which gives
information relevant for modeling the batch process with a kinetic
model. A total of five different batches were modeled in a multiway
kinetic model. Augmentation of four prerun batches with data from
a fifth batch permitted updating of parameter estimates as well
as forecasting reaction trajectories, reaction yield and end points.
Data Analysis Method. The use of multivariate absorption
measurements for fitting multiway kinetic models has received
considerable attention in the past few years.7-10 In most of these
approaches, nonlinear least-squares estimation of model param-
eters is required. In this paper we chose the nonlinear least-
squares (NLLS) Levenberg-Marquardt algorithm11 to fit the
activation parameters of our proposed reaction mechanism using
a hard-modeling approach. The NLLS fitting of multivariate
absorption data has been used for some time,12 and the procedure
is briefly outlined here.
Rˆ ) Y - Cˆ Aˆ
(2)
Calibration-free modeling of batch reactions is summarized by the
following steps: (1) A reaction mechanism is postulated giving a
model composed of a system of simultaneous ODEs. (2) Numer-
ical integration of the ODEs produces an estimate of Cˆ . (3) Least-
squares fitting of Cˆ to Y produces calibration-free estimates of
pure component spectra, Aˆ . (4) The model parameters (rate
constants) are iteratively adjusted by a nonlinear estimation
method until no further reduction in Rˆ 2 is obtained.
Specialized software with a graphical user interface was written
in C++ to perform the kinetic modeling. After postulating a
chemical model, it is encoded into strings of text representing
the proposed reaction mechanism and input into the computer
program. For example, the reaction between acetic anhydride and
butanol might be represented by the following string: “AcOAc +
BuOH + DMAP > BuOAc + HOAc + DMAP”. The program uses
an intelligent model parser to extract the number of species,
species names, their corresponding stoichiometric coefficients,
and the number of reactions. The program produces a list of
parameters, their reaction coefficients and constructs a system
of ODEs that describes the change in concentration of each
species with time.13 This approach is completely general so that
a system of ODEs of arbitrary complexity can be automatically
generated for any number of coupled reactions. With knowledge
of the initial concentrations of the species in the chemical model,
the differential equations are integrated to yield the concentration
of each species at any desired time. Only the simplest system of
ODEs describing chemical models can be integrated explicitly;
consequently, the Bulirsch-Stoer numerical integration technique
was used,11 which is capable of integrating complex systems of
stiff ODEs to any desired level of accuracy.
The on-line spectroscopic measurements consist of s spectra
measured at w wavelengths arranged into a data matrix Y with
dimensions s × w. According to the Beer-Lambert law, this
matrix can be decomposed into the product of a concentration
matrix C (s × n) and a matrix of molar absorptivities A (n × w),
where n is the number of absorbing species.
The parameters to be fitted associated with a chemical model
are the rate constants of each step. Under isothermal conditions,
it is assumed that the rate constants of all reaction steps remain
constant for all experiments. The experiments reported in this
paper were performed under nonisothermal conditions, neces-
sitating the use of the Arrhenius model to describe the rate
constants, k, as a function of the absolute temperature, T:
Y ) CA + R
(1)
k ) Ae-E / RT
(3)
A
However, due to the noise inherent in any measurement and other
sources of measurement error such as concentration errors, this
where A is the preexponential factor, EA is the activation energy
of the reaction, and R is the universal gas constant. The two fitted
parameters, A and EA, are different by several orders of magnitude
in most typical applications, which causes significant problems
during nonlinear estimation. To alleviate this problem, the Ar-
rhenius equation was reparametrized to the form shown in eq 4,
(6) Hoogeveen, R. W. M.; van der A, R. J.; Goede, A. P. H. Infrared Phys. Technol.
2 0 0 1 , 42, 1-16.
(7) Dyson, R.; Maeder, M.; Neuhold, Y.-M.; Puxty, G. Anal. Chim. Acta 2 0 0 3 ,
490, 99-108.
(8) Neuhold, Y.-M.; Maeder, M. J. Chemom. 2 0 0 2 , 16, 218-227.
(9) Bijlsma, S.; Louwerse, D. J.; Smilde, A. K. J. Chemom. 1 9 9 9 , 13, 311-329.
(10) Bijlsma, S.; Louwerse, D. J.; Windig, W.; Smilde, A. K. Anal. Chim Acta
1 9 9 8 , 376, 339-355.
(11) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical
Recipes in C: The Art of Scientific Computing; Cambrige University Press:
Cambridge, U.K., 1992.
(12) Maeder, M.; Zuberbu¨ hler, A. D. Anal. Chem. 1 9 9 0 , 62, 2220-2224.
(13) Dyson, R.; Maeder, M.; Puxty, G.; Neuhold, Y.-M. Inorg, React. Mech. In
press.
2576 Analytical Chemistry, Vol. 76, No. 9, May 1, 2004

ref)EA/ R
k ) krefe-(1/ T-1/ T
(4)
mented this alternative merit function and found it to be very
sensitive to the initial parameter estimates. Even for a simple
simulated second-order reaction, where a single parameter was
fitted, it converged to the correct value only when the initial
parameter estimate was close to the optimal value, whereas using
function 6 resulted in convergence from a broader range of initial
estimates. Also, we found function 7 to be very sensitive to the
number of significant eigenvectors used. When applied to the
nonisothermal data, the fitted activation parameters did not agree
with those from the NLLS fit and produced unrealistic values.
Furthermore, the quality of the NLLS fits using function 6
indicated that the measured data were not suffering from any
significant systematic deviations from the model estimated mea-
surements.
where kref is the rate constant at temperature Tref. This form of
the Arrhenius model was used during nonlinear least-squares
fitting to reduce correlation between the parameters A and EA,
thereby improving nonlinear fitting performance.14 Once fitting
was complete, a straightforward calculation was used to recover
A from kref and Tref
.
Before nonlinear model fitting was performed, starting condi-
tions were specified, e.g., the initial concentrations of the species
involved in the reactions and initial estimates for the model
parameters kref and EA. Numerical integration of the system of
simultaneous ODEs to the spectral acquisition times was per-
formed with these initial estimates and the known initial concen-
trations to calculate the time-dependent matrix of concentration
profiles, Cˆ . After integration, calibration-free estimates of the molar
absorptivities, Aˆ , were obtained by a linear regression step,12
where Cˆ + is the pseudoinverse Cˆ .
EXPERIMENTAL SECTION
Equipment. All reactions were performed in an auto-MATE
reactor system (H.E.L. Inc., Lawrenceville, NJ), a miniature (50
mL) computer-controlled multiple reactor system. Process condi-
tions, including reaction temperature, jacket temperature, and
agitation, were controlled by WinISO software from H.E.L. running
on a 333-MHz Pentium II computer. The heads of the auto-MATE
reactors were modified to accept fiber-optic retroreflection Near-
IR probes from Equitech International (New Ellenton, SC).
Near-IR spectra were collected over the range of 1400-2100
nm using a CP-140 short focal length, 140 mm, near-IR fiber-optic
spectrophotometer fitted with a 120 grove/ mm grating (JY Horiba,
Edison NJ) coupled to a high-speed 256-element, InGaAS array
from Hamamatsu (Bridgewater NJ). Near-IR illumination was
provided by a tungsten filament lamp and adjustable power supply
(Thermo Oriel, Stratford CT). The retroreflection probe had a
1-mm path length housed in a stainless steel probe 12.5 mm in
diameter and 140 mm in length. The probe was fitted with two
400-µm core low-OH fused-silica fibers.
Spectrometer data acquisition was performed using in-house-
written software using LabVIEW version 6 software (National
Instruments, Austin TX) and an instrument control library supplied
by Hamamatsu. The power supplied to the lamp was adjusted to
ensure optimum use of the diode array and A/ D converter
dynamic range with a maximum digitized intensity in the range
of 6.5 × 104 counts.
Aˆ ) Cˆ +Y
(5)
The Levenberg-Marquardt method of nonlinear least squares was
used to minimize the sum of squares of the resulting residuals,
Rˆ , by iteratively adjusting the nonlinear parameters kref and EA.
Rˆ ) Y - Cˆ Aˆ
(6)
To achieve a more robust fit with reduced correlation between
the nonlinear parameters and also to take advantage of the fact
that a number of separate experiments were made under different
initial conditions, fitting was carried out in a second-order global
analysis mode.15,16 In this mode, a global chemical model was fitted
simultaneously to all the experiments. Ideally, a single matrix of
molar absorptivities would be calculated in such an analysis;
however, due to baseline shifts between measurements, it was
necessary to calculate individual Aˆ for each experiment. Furusjo¨
et al.14 have carried out similar fitting of nonisothermal in situ
infrared and Raman measurements. They argued for a different
function to be minimized:
Reaction Conditions. The reaction used to study the capabili-
ties of the combined calorimeter/ NIR spectrograph and batch
analysis software was the acetylation of 1-butanol with acetic
anhydride with DMAP as a catalyst. A 20-mL aliquot of anhydrous
1-butanol (Aldrich Milwaukee, WI) and 21 mL of acetic anhydride
(Aldrich) were added to a dry 50-mL reactor at room temperature.
(Warning, acetic anhydride is a corrosive and combustible liquid,
incompatible with strong oxidizing agents, water, strong bases,
and alcohols). All reagents came from freshly opened bottles. The
reactor jacket was cooled to -25 °C, and power was applied to
the internal reactor heater to warm the contents to 20 °C.
Accurately weighted amounts of DMAP in the range of 20-80
mg (Aldrich) were added directly into the glass reactor vessel to
initiate the reaction. Near-IR spectra were collected every second
using a dry N2 reference. The reactions were run until the solution
temperatures returned to 20 °C. An additional two reactions were
performed with 400 and 800 µL of water added to the reactor prior
to initiation of the batch reaction to simulate process upsets. An
Rˆ ) (I - UUT)Cˆ
(7)
where U contains the trimmed set of left eigenvectors from the
singular value decomposition of Y. For nonisothermal data, the
advantage of this function is that it is resistant to absorption peak
shifts with temperature that may occur for mid-IR and Raman data.
This function, however, relies on the assumption that Cˆ is
contained in the space spanned by the columns of U. This is only
true near the optimum solution, suggesting that if poor initial
parameter estimates are used, the routine may struggle to
converge to the correct C. To test this hypothesis, we imple-
(14) Furusjo¨ , E.; Svensson, O.; Danielsson, L.-G. Chemom. Intell. Lab. Syst. 2 0 0 3 ,
66, 1-14.
(15) Bugnon, P.; Chottard, J.; Jestin, J.; Jung, B.; Laurenczy, G.; Maeder, M.;
Merbach, A. E.; Zuberbu¨ hler, A. D. Anal. Chim Acta 1 9 9 4 , 298, 193-201.
(16) Dyson, R. M.; Kaderli, S.; Lawrance, G. A.; Maeder, M.; Zuberbu¨ hler, A. D.
Anal. Chim. 1 9 9 7 , 353, 381-393.
Analytical Chemistry, Vol. 76, No. 9, May 1, 2004 2577

Ta b le 1 . Re a c t io n Co n d it io n s a n d Mo d e l Es t im a t e d
In it ia l Ra t e s a
mass
DMAP
(mg)
concn
DMAP
initial rate
at 20 °C
exotherm
range (°C)
batch
(mol L-1
)
(mol L-1 min-1
)
1
2
3
4
5
80.6
51.7
53.1
55.0
21.2
0.0160
0.0102
0.0105
0.0109
0.00420
1.22
19.7-32.9
19.5-26.5
19.3-26.4
19.7-27.6
20.0-21.7
0.783
0.804
0.833
0.321
a
Initial concentration of acetic anhydride and butanol for all batches
was 5.41 and 5.33 mol/ L, respectively. The mass of DMAP added was
varied.
Fig u re 1 . Estimated pure component spectra for 80.6-mg batch.
The solid line is a combined spectrum for the reactants (acetic
anhydride and butanol), and the dash-dot line is the products (acetic
acid and butyl acetate).
unexpected absorption band was observed in the batch with 400
µL of water that suggested some additional source of contamina-
tion besides water was present; thus, it was not used in the
calculations described below.
batches 2-5 had smaller exotherms. Furthermore, a wide range
of different models with increasing complexity could be fit to each
batch individually with low residuals and apparent low lack of fit.
However, only the model described in eq 8 estimated the
concentration profiles with low residuals for all five batches
simultaneously in a multiway kinetic model. As expected, fitting
one model simultaneously to many different experiments permit-
ted elimination of incorrect models and unambiguous identification
of the correct model; this is due to the increased robustness on
globalized analyses.
RESULTS
A third-order reaction mechanism (8) was fitted to all five
batches simultaneously. The initial conditions for each batch are
listed in Table 1. There is some error in the determined initial
concentrations as no consideration was given to the change in
density of the reactants after mixing. The catalyst, DMAP, was
included in the mechanism as it had a different initial concentra-
tion for each batch. Its inclusion was necessary to fit the model
simultaneously to all batches. Even though the reaction is third
order overall, in each batch, the reaction is effectively pseudo
second order because DMAP is liberated at the same rate at which
it is consumed, thereby causing the concentration of DMAP to
remain fairly constant throughout the course of the reaction.
Ideally if there are no baseline shifts or other spectroscopic
changes between batches, it is possible to estimate a single matrix
of pure component spectra in the linear least-squares step
described in eq 5. However, it was apparent there were both
baseline shifts and other spectral inconsistencies between batches,
which was not unexpected since a single-fiber NIR probe was
used. Single-fiber probes tend to be sensitive to changes in fiber
position that can sometimes produce different amounts of baseline
offset. While we were able to remove the baseline shift between
batches,17 small spectral inconsistencies between batches were
still apparent. In this circumstance, it was necessary to use a local
spectral model for each batch instead of a global spectral model.
In a local spectral model, the estimated pure spectra are allowed
to be different from batch to batch. Also a reduced wavelength
range was used for fitting, 2016.2-2235.1 nm, to avoid regions
where obvious instrumental inconsistencies occurred such as no
signal or large noise spikes. The estimated pure component
spectra for the 80.6-mg batch are shown in Figure 1. The negative
regions are due to the presence of small, negative absorption
measurements in the measured reaction spectra. The estimated
pure component spectra for the other batches are very similar
but show some minor differences. Additionally, the first 30 s of
measurements from each batch was excluded from fitting. This
was to allow for the dissolution and mixing of DMAP since it was
added as a solid.
AcOAc + BuOH + DMAP f
BuOAc + AcOH + DMAP (8)
The amount of DMAP was varied over a wide range between
batches, giving significantly different initial rates and significantly
different reaction exotherms. For example, the initial reaction rate
estimated for batch 1 was 1.22 and 0.321 mol L-1 min-1 for batch
5. The starting temperature for each batch was ∼20 °C. The peak
temperature during batch 1 reached 32.9 °C, whereas the peak
temperature was only 21.7 °C for batch 5. This is remarkable
considering the reactor jacket was chilled to -25 °C and the
reactor set point was 20 °C. At steady-state conditions, the reaction
calorimeter heater controller supplied ∼18 W of power to the
auxiliary heater in the reactor to maintain this temperature
differential. During batch 1, the evolution of heat from the
chemical reaction was sufficiently rapid that the controller reduced
the auxiliary heater power to zero for nearly 3 min and still the
heat removed by the jacket at -25 °C was insufficient to maintain
constant temperature.
The fitted preexponential factor and activation energy and their
respective uncertainties are listed in Table 2. We found the
minimum to be well defined with the NLLS routine converging
Fitting the postulated model shown in eq 8 to each experiment
individually did not provide sufficient information to simulta-
neously estimate kref and EA with good reliability, except for batch
1. The relatively large exotherm of batch 1 ensured good
observability for EA and was defined by this batch alone, whereas,
(17) Maeder, M.; Neuhold, Y.-M.; Olsen, A.; Puxty, G.; Dyson, R.; Zilian, A. Anal.
Chim. 2 0 0 2 , 464, 249-259.
2578 Analytical Chemistry, Vol. 76, No. 9, May 1, 2004

represented in the training batches was used for batch 5; thus,
the forecast of batch 5 reaction profiles represents an extrapolation
of the model to new conditions. Slightly different parameter
estimates for the preexponential factor, A, and the activation
energy, EA, were obtained when only batches 1-4 were included
in the model-fitting procedure. However, accurate predictions of
batch 5 concentration profiles were obtained (see Table 4). At
15.67 min, the model forecasted concentration of reactants and
product deviated from the values obtained from the full model by
less than 2%.
Ta b le 2 . Fit t e d Va lu e s fo r t h e P re e x p o n e n t ia l Fa c t o r
a n d t h e Ac t iva t io n En e rg ya
A (L2mol-2min-1
1.94 × 1 0 7 (1)
)
EA (kJ)
38.5 (2)
a
Uncertainties are given in parentheses and represent two standard
deviations of the last significant digit.
Detection of P rocess Upsets. To illustrate the use of kinetic
modeling for detecting process upsets, a kinetic model was fitted
to batches 1-5 plus one batch with 55.1 mg of DMAP catalyst
and 800 µL of adder water to simulate a process upset. The
addition of the 800 µL of water resulted in a significant change to
the model-estimated parameters, the overall standard deviation,
and the estimated pure component spectra (see Table 5).
Additionally, the standard deviation of the residuals for the
adulterated batch was significantly larger than for the other
batches, indicating the adulterated batch deviates from the batches
1-5. Examination of the estimated pure component spectra
revealed that the spectra of the adulterated batch differ signifi-
cantly from those of the five good batches. The most obvious
difference is that the adulterated batch contains an additional peak
due to water. Part of the peak can be seen at the lower
wavelengths in Figure 5.
Speed of P rocess Upset Detection. To study the speed at
which process upsets can be detected by kinetic modeling, the
model parameters were fixed at the values determined for the
five good batches. Spectra were then added from the adulterated
batch one at a time until enough spectra were present such that
detection of a process upset was evident. The first 30 s of data
from the adulterated batch was omitted to allow for DMAP mixing
times and to remain consistent with the data analysis of the five
good batches. Following the addition of two spectra, the standard
deviation of the residuals for the spectra from the adulterated
batch was 6.25 × 10-3, ∼6 times that of the overall standard
deviation prior to the addition of the spectra from the adulterated
batch. The standard deviation of the spectra from the adulterated
batch remained high, about 4-5.5 times greater than the overall
standard deviation prior to the addition of the adulterated batch.
The presence of the process upset is detected after the inclusion
of two spectra from the adulterated batch.
Fig u re 2 . Nonisothermal fits at 2202.9 (upper) and 2212.6 nm
(lower) for 80.6 mg of DMAP. The dotted lines are the measured data
and the solid lines the calculated data.
to the listed values from a number of different initial starting
points. This fact, coupled with the excellent quality of the fits (0.5%
RSD in the fitted parameters, kref and EA, 0.5% RSD in fitted
absorbance values), indicates the observability of the estimated
parameters was excellent. The listed uncertainties are estimated
by statistical means, assuming a local linearization of the nonlinear
model near the minimum gives a good approximation to the actual
error covariance matrix.11 They are most certainly an underesti-
mate of the true uncertainties as no account has been taken of
departures from the above assumptions or the presence of
experimental errors such as those in mass and volume. The quality
of the fit is good, and representative fits at two different
wavelengths are shown in Figure 2. The overall standard deviation
in the residuals from the fit was 1.13 × 10-3 absorbance units,
which is well within the instrumental noise level. A plot of the
estimated reaction profiles is shown in Figure 3.
To illustrate the benefit of fitting the activation parameters with
an Arrhenius model, we also carried out the fit by assuming
isothermal conditions. In this situation, only a rate constant is fitted
rather then the activation parameters of the Arrhenius equation.
The fitted rate constant is assumed to be invariant throughout
the reaction. Figure 4 shows the temperature profile for the 80.6-
mg batch and fits at the same wavelengths as Figure 2, assuming
isothermal conditions. The isothermal model underestimates the
initial reaction rate early in the batch when the temperature is
higher and overestimates the rate later in the batch when the
temperature is lower. This illustrates the benefit of incorporating
temperature into the model.
CONCLUSIONS
In this paper, we have demonstrated that fitting multivariate
kinetic models to in situ spectroscopic measurements can be a
powerful calibration-free tool for batch process monitoring and
control. First-principles models coupled with in situ spectroscopic
measurements can be used to determine reaction mechanisms,
time-dependent batch concentration profiles, reaction yields, and
reaction end points. Costly, time-consuming off-line reference
measurements with multivariate calibration methods such as PLS)
are not required to obtain accurate real-time estimates of concen-
tration profiles. Real-time modeling of a batch reaction can be
achieved by augmenting prerun batch data with recently acquired
spectral data followed by fitting of a multiway kinetic model. In
real-time applications, accurate initial guesses of the parameter
P rediction of Batch 5 from Batches 1 -4 . To illustrate the
use of kinetic modeling to forecast batch profiles and detect batch
end points, a kinetic analysis of batches 1-4 was performed (see
Table 3). The resulting model was used to forecast the reaction
profiles for batch 5. A significantly lower level of catalyst not
Analytical Chemistry, Vol. 76, No. 9, May 1, 2004 2579

Fig u re 3 . Estimated reaction profiles of acetic anhydride (s), butanol (- -), and acetic acid and butyl acetate (‚‚) for five batches. DMAP is
not included in the plots as its concentration is much lower than the other species and it remains constant throughout.
estimates from prerun data can be used to facilitate rapid
convergence of the iterative nonlinear fitting process.
overfitting. In this paper, we used fitting of a kinetic model to
several batches at significantly different conditions (different
amount of catalyst, different reaction exotherms, and temperature
profiles) to validate our model and guard against overfitting. While
it is possible to find several reaction mechanisms (models) that
give acceptable fits to individual batches, the likelihood is very
low that the wrong mechanism will correctly fit other batches at
significantly different conditions.
There are several difficulties with the multivariate kinetic
modeling approach to process monitoring and control; however,
the benefits (discussed later) outweigh the difficulties in many
circumstances. First, as in any model-fitting procedure, increasing
the model complexity will improve the quality of the fit. In
multivariate kinetic fitting, model complexity is increased by
increasing the complexity of the reaction mechanism. This usually
results in lower residuals for individual batches or experiments.
Some method of model validation is necessary to guard against
Another drawback of the multivariate kinetic modeling ap-
proach to process monitoring and control is that a significant
amount of effort is required to elucidate the correct mechanism.
2580 Analytical Chemistry, Vol. 76, No. 9, May 1, 2004

Fig u re 5 . Spectra for the 80.6-mg batch. The solid lines show the
estimated pure component spectra for the five good batches. The
dash-dot lines show estimated pure component spectra of the
adulterated batch. Spectra were normalized to 1 to facilitate visual
comparison.
Fig u re 4 . Temperature profile of the reaction and isothermal fits at
2202.9 (upper) and 2212.6 nm (lower) for 80.6 mg of DMAP. The
dotted lines are the measured data and the solid lines the calculated
data.
models that accurately estimate reaction profiles and reaction
yields over a wide range of operating conditions. The time-
dependent concentration profiles of species not directly observable
in the spectroscopic measurements can be inferred from the
model. For example, a catalyst may be present at concentrations
too low to be analytically detectable, or there may be a lack of
useful absorption bands in the range studied. In this paper, we
have shown that the effect of a catalyst, however, can still be
elucidated if a sufficient number of properly designed experiments
are included in a multiway kinetic model.
Ta b le 3 . Fit t e d Mo d e l P a ra m e t e rs fo r Ba t c h e s 1 -4 a
A (L2mol-2min-1
3.34 × 107 (1)
)
EA (kJ)
39.9 (2)
a
Uncertainties are given in parentheses and represent two standard
deviations of the last significant digit.
The multivariate kinetic modeling approach is sensitive to the
initial concentrations of reagents and reagent purity. Accurate
values of the initial concentrations of reagents must be obtained
from the amounts of reagents used to charge the reactor. This
sensitivity, however, allows for very rapid and sensitive detection
of process upsets caused by impurities. When process upsets
occur, inspection of model results (estimated rate constants,
spectral residuals, and estimated pure component spectra) can
be compared to results from prior batches to give useful diagnostic
information. For example, changes in estimated rate constants
may indicate a significant change in catalyst activity, catalyst purity,
or the presence of a different reaction mechanism from contami-
nation, such as water as shown in this paper. A significant increase
in spectral residuals can also indicate the presence of contamina-
tion. The appearance of new unexpected absorption bands in
estimated pure component spectra can give diagnostic information
and help identify the source of contamination.
Ta b le 4 . P re d ic t e d a n d Fit t e d Co n c e n t ra t io n s fo r Ba t c h
5 a ft e r 1 5 .6 7 m in
predicted
(mol/ L)
fitted
(mol/ L)
error
(%)
acetic anhydride
butanol
acetic acid and butyl acetate
2.81
2.73
2.60
2.76
2.68
2.65
1.81
1.87
1.89
Ta b le 5 . Co m p a ris o n o f Re s u lt s Fo llo w in g Ad d it io n o f
t h e Ad u lt e ra t e d Ba t c h
A (L2 mol-2
min-1
EA
(kJ)
SD of
residuals
)
5 batches
1.94 × 107 (1)
38.5 (2)
68.5 (4)
1.13 × 10-3
2.30 × 10-3
5 batches + adulterated 1.96 × 1012 (2)
batch
Last, in multiway kinetic fitting, the use of one global kinetic
model and different local spectral models for individual batches
offers significant advantages (discussed below). One drawback
of this approach is that the estimated concentration matrices of
the batches, Ci, are usually rank deficient. For example, in the
reaction A + B f C + D, both reactants A and B disappear at the
same rate and both products are formed at the same rate.
Although the matrix C has four columns, its rank is only 2. This
means the solution to eqs 5 and 6 must be solved at lower rank,
and only two “pure component” spectra are estimated. The
estimated spectra thus represent “pseudospecies”, where the
spectra of the two pseudospecies represent the sum of the pure
component spectra, A + B and C + D, respectively.
Multiple batches or experiments are required. Good experimental
designs must be used with sufficiently different levels of controlled
variables to unambiguously identify the correct mechanism and
provide sufficient information to unambiguously identify and
estimate all important model parameters. Inadequate designs will
result in confounded model parameters. Sometimes many different
reaction mechanisms must be tried until the best mechanism is
found.
The benefits of multivariate kinetic modeling approach to
process monitoring and control however are significant. Fitting
multivariate kinetic models gives fundamental process insights.
We have shown that including important information in the model
such as the amount of catalyst and temperature gives robust
Analytical Chemistry, Vol. 76, No. 9, May 1, 2004 2581

The advantage gained by using different local spectral models
for individual batches is significant, however. We have shown in
this paper that different baseline offsets or even small shifts in
spectral response from batch to batch do not seriously impact
the quality of the calibration-free kinetic models. These batch-to-
batch differences appear as small differences in the estimated pure
component spectra. As long as the time variant spectral response
gives adequate information to accurately and precisely estimate
the model parameters, excellent calibration-free estimates of time-
dependent concentration profiles are obtained. This represents a
significant advantage compared to traditional multivariate calibra-
tion approaches such as PLS, where small changes in a spec-
trometer’s response caused by changing a burned-out lamp or
replacing a damaged fiber-optic probe or cable can invalidate a
calibration model.
ACKNOWLEDGMENT
P.G. acknowledges that this research was supported in part
by a grant from the National Science Foundation (CHE-0201014)
and a grant from the Measurement and Control Engineering
Center (MCEC) at the University of Tennessee, Knoxville, TN.
Received for review November 15, 2003. Accepted
February 23, 2004.
AC035356I
2582 Analytical Chemistry, Vol. 76, No. 9, May 1, 2004