Dynamic Response Surface Models: A Data-Driven Approach for the

Article pubs.acs.org/IECR

Dynamic Response Surface Models: A Data-Driven Approach for the Analysis of Time-Varying Process Outputs Nikolai Klebanov and Christos Georgakis* Department of Chemical and Biological Engineering and Systems Research Institute for Chemical and Biological Processes, Tufts University, Medford, Massachusetts 02155, United States ABSTRACT: In a recent publication (Ind. Eng. Chem. Res. 2013, 52 (35), 12369) we generalized the classic design of experiments (DoE) methodology by introducing the Design of Dynamic Experiments (DoDE), allowing for the systematic design of experiments involving time-varying inputs. Here, we expand the response surface model (RSM) methodology, used in DoE and DoDE problems, so that it describes the time-evolution of the process, not just the results at the end of the experiment. We apply this generalized type of RSM model, to be denoted by DRSM, to three example processes; a nonisothermal batch reactor with a simple reaction, an isothermal semibatch reactor with several reactions, and a semibatch penicillin fermentation process. Using a limited number of online measurements at prespecified equidistant time instants, we are able to quickly and accurately represent the time evolution of the process output through these simple interpolative data-driven models. The ever-increasing availability of time-resolved measurements is expected to make the proposed approach widely useful.

1. INTRODUCTION Data-driven modeling is a crucial tool for the optimization and control of many manufacturing processes for which it is not possible to develop a detailed and accurate knowledge-driven model describing the inner workings of the process. Continuous processes, which are widely used in the chemical and petrochemical industries, have high production rates and have justified the development and use of detailed knowledge-driven models for process design, optimization, and control. On the other hand, batch and semibatch processes, such as batch fermentations employed by the pharmaceutical industry, have smaller production rates and their detailed modeling is a major challenge. A small fraction of such batch processes have been studied by means of a knowledge-driven model,2 while for their majority no quantitative model has been developed. This is because such a model requires a thorough knowledge of the process’s inner workings for the modeler to write the necessary material and energy balances. Such details include the reaction stoichiometry, the details of the thermodynamic equilibria, and the related kinetic or transport rate relationships. While such an approach is usually capable of providing significant insight, the development of the comprehensive knowledge-driven models requires a significant investment of time and money. Data-driven models are thus an attractive alternative as they utilize only experimental data from an existing process, and have reduced developmental costs compared to those inherent in knowledgedriven models. The design of experiments (DoE) methodology is a classical data-driven approach to guide systematic experimentation, as well as the development of data-driven models. DoE allows the experimenter to run a minimal number of experiments, properly © XXXX American Chemical Society

designed, to characterize the overall behavior of the process. The related modeling approach necessary to statistically estimate the relationship between the input variables, often called factors, and the output variables, is called response surface methodology (RSM);3,4 a robust and simple technique for the estimation of time-invariant process outputs, utilizing linear parameter estimation tools. This has found wide applicability in a variety of processes5−10 In a previous publication,1 the design of dynamic experiments (DoDE) methodology was introduced as a way to incorporate input variables that are time-dependent. This provides a reliable and efficient way to design experiments around input prof iles, such as a time-varying reactant feed rate or the reaction temperature. However, the related RSM model relates these inputs, regardless of their nature, to a single output, usually the end-point of the process. Although separate RSM models could be developed from experimental results at different time instants during each batch, this approach offers a fragmented view of the time evolution of a batch operation. In the present publication we propose a novel approach that generalizes the RSM for enabling it to handle time-resolved output data. As this methodology is equipped to model the dynamic behavior of a process, we call this the dynamic response surface methodology (DRSM). The DRSM is a data-driven model structure that captures the relationship between process inputs (time-variant or time-invariant) and time-resolved output Received: October 1, 2015 Revised: March 7, 2016 Accepted: March 13, 2016

A

DOI: 10.1021/acs.iecr.5b03572 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

whether one has time-resolved measurements in the preliminary set of experiments. One can then switch to the quadratic DRSM model once the experimental region has shifted near the point of optimality. Here, we use τ, 0 ≤ τ ≤ 1, as the dimensionless time in the (t0, tf) interval of interest, defined by τ = (t − t0)/(tf − t0). Quite often we will have t0 = 0. The unknown parametric functions for this model are β0(τ), βi(τ), βij(τ), and βii(τ) with i, j = 1, 2, ..., R; i < j. They need to be estimated by a finite set of measurements at instants τ1, τ2, ..., τK. Consequently, we aim for a finite approximation for the identification of the β0(τ), βi(τ), βij(τ), and βii(τ) functions of time. 2.1. Parameterization with Shifted Legendre Polynomials. A finite approximation of β0(τ), βi(τ), βij(τ), and βii(τ) is achieved by a polynomial expansion, using Shifted Legendre polynomials. We select them as the basis for the parametrization of each{βq(t)|q = 0 or q = i, ij, ii}function because of their orthogonality in the (0, 1) interval. The first five shifted Legendre polynomials are

variables. We describe how the DRSM is constructed; focusing on how the set of sampled measurements from each experiment is incorporated into the linear regression used for the estimation of the model parameters. The use of the statistical tool of the analysis of variance (ANOVA)11 ensures that the computed model has no insignificant parameters. We also define and apply an F-test to determine the minimum number of polynomials necessary to accurately represent the time evolving profiles of the measured process outputs during each batch. We demonstrate the performance of the DRSM approach by applying it in silico to three example processes: (i) a nonisothermal batch reactor with a simple reaction, (ii) an isothermal semibatch reactor with several reactions, and (iii) a semibatch penicillin fermentation process.

2. THE DYNAMIC RESPONSE SURFACE MODEL (DRSM) METHODOLOGY A quadratic RSM in both DoE and DoDE problems has the following form, where y is the response of interest and the variables (x1, x2, ..., xn) represent the process inputs, called factors. n

y = β0 +

n

n

∑ βi xi + ∑ ∑ i=1

P0(τ ) ≜ 1,

n

βijxixj +

i=1 j=i+1

∑ βiixi

2

P1(τ ) ≜ −1 + 2τ , (1)

i=1

P2(τ ) ≜ 1 − 6τ + 6τ 2

In the case of a DoDE, some of the xi variables are classical factors and some are dynamic subfactors corresponding to one or more dynamic factors.1 The above quadratic form of the RSM is not the only one of use in either the DoE or DoDE set of experiments. Depending on whether the experimental program is in the preliminary exploratory stage or not, at least two other RSM forms might be considered. One is the purely linear one and the other is often called the “two-factor interaction” (2FI) form. They are described in the following two equations, respectively. n

y = β0 +

n

∑ βi xi ;

y = β0 +

i=1

n

P4(τ ) ≜ 1 − 20τ + 90τ 2 − 140τ 3 + 70τ 4

R−1

β0(τ ) =

R−1

i=1 j=i+1

βi (τ ) =

y(τ ) = β0(τ ) +

i=1

j=1 i Fa , n1, n2

(24)

T (τ ) = 32.5 + 17.5z(τ )

Here n1 and n2, defined earlier, are the two degrees of freedom of the F distribution. For the usual 95% confidence level, the value of the parameter α is equal to 0.05. Instead of making the above inequality test, we can calculate the p = p(R, K) value corresponding to the F0((R, K)) statistic. With such a p value at hand, we can decide as follows:

(28)

We consider here the simplest DoDE design with two dynamic subfactors x1 and x2, representing constant or linearly varying temperature profiles, parametrizing z(τ) by a polynomial expansion with two shifted Legendre polynomials, z(τ ) = x1P0(τ ) + x 2P1(τ ) = x1 + x 2( −1 + 2τ )

(29)

resulting in the following temperature dependence in time T (τ ) = 32.5 + 17.5{x1 − x 2 + 2x 2τ }

If the conclusion of the above test is that the Null Hypothesis is rejected, this implies that, with a certainty of 95%, the part of the data not explained by the model and represented by SŜ un(R , K ), is larger than the normal variability of the process, represented by SŜ err . To certify that the model has represented all of the nonrandom variability in the data we wish to fail to reject the null hypothesis and thus need to have p(R, K) < 0.95. If the Null Hypothesis is rejected, then the time-resolved measurements used are not sufficient and the present value of K should be increased. An excellent source of all the statistical concepts utilized here can be found in the following textbook.12

The function z(τ) needs be constrained between −1 and +1, so we place the following constraint on x1 and x2 as detailed in our previous publication1 −1 ≤ x1 ± x 2 ≤ +1 Because the Shifted Legendre polynomials are bound between −1 and +1 when τ ∈ (0,1), the above set of inequalities on the xi variables guarantees that z(τ) is constrained between −1 and +1. A set of experiments specified by the Center Composite Design (CCD)11 with α = 2 will be performed, for the purpose of estimating a following quadratic DRSM model.

4. BATCH REACTOR WITH SINGLE REVERSIBLE REACTION 4.1. Definition of the Problem. Here, we apply the DRSM methodology to develop a data-driven model for a batch reactor with a reversible reaction of reactant A to product B with the following characteristics: k1

A⇄B k2

r = k1[A] − k 2[B]

y(τ ) = β0(τ ) + β1(τ )x1 + β2(τ )x 2 + β12(τ )x1x 2 + β11(τ )x12 + β22(τ )x 2 2

(31)

The values of the dynamic subfactors, x1 and x2, defining the temperature profiles linear in time, along with the characteristics of the temperature profiles, are given in Table 1. Table 1. Center Composite Design (CCD) DoDE Design with 13 Runs (Linear T Profiles)

(26)

where

coded variable coefficient

⎛ E ⎞ k1 = k10 exp⎜ − 1 ⎟ , ⎝ RT ⎠ k10 = 1.32 × 108 h−1, E1 = 10000 kcal, ⎛ E ⎞ k 2 = k 20 exp⎜ − 2 ⎟ , ⎝ RT ⎠ k 20 = 5.25 × 1013 h−1, E2 = 20000 kcal

(30)

(27)

This kinetic model is used in the simulation of the process, and the developer of the DRSM model does not have any knowledge of either the reaction stoichiometry or of the kinetic rates. We will consider the effects of temperature on the conversion of compound A. The DoDE methodology will be used to define several in silico experiments with different time-varying temperature profiles in the reactor. Although a time-varying input was selected for this simulation, traditional time-invariant inputs could also be used for the deployment of a DRSM model. We collect measurements of the reactant’s concentration at a series of equidistant time intervals K. To properly represent true experiments, a constant fractional measurement error will be

run #

x1

x2

temperature profile

1 2 3 4 5 6 7 8 9 10 11 12 13

−1.0 −0.5 0.0 −0.5 0.0 0.5 0.0 0.5 1.0 0.0 0.0 0.0 0.0

0.0 0.5 1.0 −0.5 0.0 0.5 −1.0 −0.5 0.0 0.0 0.0 0.0 0.0

at 15 °C 15 → 32.5 °C 15 → 50 °C 32.5 → 15 °C at 32.5 °C 32.5 → 50 °C 50 → 15 °C 50 → 32.5 °C at 50 °C at 32.5 °C at 32.5 °C at 32.5 °C at 32.5 °C

This set of 13 experiments consists of 9 unique experiments, 6 for the same number of model parameters and 3 for the estimation of the GoF statistic. Four additional replicate experiments are perforemd at the center point (x1 = 0; x2 = 0) and used to estimate the normal process variability. E


Article


more reasonable maximal value for R. In the present simple case the time dependency of the DRSM model is a second order polynomial in time.

The concentration measurements of reactant A during the batch are taken at K equidistant time instants, the first at τ = 0 and the last one at τ = 1, and are used to calculate the fractional conversion of compound A. The 9 unique temperature profiles versus time are presented in Figure 1. The resulting fractional conversion profiles (before the error is added) are given in Figure 2.

y(τ ) = {γ0,1P0(τ ) + γ0,2P1(τ ) + γ0,3P2(τ )} + {γ1,1P0(τ ) + γ1,2P1(τ ) + γ1,3P2(τ )}x1 + {γ2,1P0(τ ) + γ2,2P1(τ ) + γ2,3P2(τ )}x 2 + {γ12,1P0(τ ) + γ12,2P1(τ ) + γ12,3P2(τ )}x1x 2 + {γ11,1P0(τ ) + γ11,2P1(τ ) + γ11,3P2(τ )}x12 + {γ22,1P0(τ ) + γ22,2P1(τ ) + γ223P2(τ )}x 2 2

(32)

We collect 27 data points from the 9 unique experiments to estimate 18 parameters. The model coefficients estimated along with their 95% confidence intervals are

Figure 1. Three-level DoDE Center Composite Design (CCD) in Table 1 with two dynamic subfactors.

⎡ γ = + 0.621 ± 0.0187 ⎢ 0,0 β0 ⇒ ⎢ γ0,1 = + 0.379 ± 0.0085 ⎢ ⎢⎣ γ0,2 = − 0.225 ± 0.0139

⎡ γ1,0 = − 0.094 ± 0.015 ⎢ β1 ⇒ ⎢ γ1,1 = − 0.060 ± 0.015 ⎢ ⎢⎣ γ1,2 = + 0.035 ± 0.017 (33)

⎡ γ2,0 = − 0.020 ± 0.013 ⎢ β2 ⇒ ⎢ γ2,1 = − 0.064 ± 0.013 ⎢ ⎢⎣ γ2,2 = − 0.044 ± 0.015

⎡ γ , not significant ⎢ 12,0 β12 ⇒ ⎢ γ12,1 = − 0.112 ± 0.058 ⎢ ⎢⎣ γ12,2 = − 0.103 ± 0.055 (34)

⎡ γ11,0 = − 0.025 ± 0.024 ⎢ β11 ⇒ ⎢ γ11,1, not significant ⎢ ⎢ γ , not significant ⎣ 11,2

⎡ γ = + 0.026 ± 0.025 ⎢ 22,0 β22 ⇒ ⎢ γ22,1, not significant ⎢ ⎢ γ = − 0.048 ± 0.023 ⎣ 22,2

(35)

We note that several of the gamma parameters are not significant. For each of the significant ones, a confidence interval is estimated along with the nominal parameter value. Some of these confidence intervals in the above equations [33, 34, and 35] are unacceptably large, because we limited ourselves to only three measurements (K = 3). From the analytical solution of the model under constant temperature profiles, we know that the time dependency of the conversion of A is an exponential function, which is not accurately approximated by a second order polynomial. 4.3. Assessing the Accuracy of the DRSM Model. At this point, we would like to assess the accuracy of the DRSM model for different values of K and R. For this reason we examine a range of values for K, between K = 3 and K = 10. We are considering the case of limited measurements, which is characteristic of a situation when these composition values are obtained by GC, LHPLC, or any other time-laborious technique. In the case where spectroscopic means of measuring compositions are used, the value of K will be considerably larger and the modeling task easier. We expect that an increase in the number of time-resolved measurements K, which allows for a larger value of R polynomials, will improve the accuracy of the DRSM model. For each investigated pair of R and K values, we calculate the fractional sum of squares in eq 21 and we report the values in Table 2. We observe that for K = 10 there are several values of R (= 5, 6, 7, 8, and 9) that produce almost identical values for SSun and thus

Figure 2. Nine simulated profiles of fractional conversion of A to B, corresponding to the DoDE design in Table 1 and Figure 1. Solid, dashed, and dotted lines represent groups of profiles corresponding to three possible initial temperatures.

4.2. Illustrative Example of Estimating a Simple DRSM. To illustrate the character of DRSM results, we consider the simple case when only three measurements are taken (K = 3) at the beginning, midway, and at the end of the batch. In a practical situation this is a very small number of measurements to provide a reasonable model for the time evolution of the process. It is presented here merely for demonstrational purposes, since the number of model parameters is minimal and their estimated values can be used to show the results of the ANOVA. In this case R can be at most equal to 3. In later examples, where there will be more time-resolved measurements, we advise that K − 1 be a F


Article

Industrial & Engineering Chemistry Research Table 2. SSun Values for DRSM of Conversion A⇄ B in a Batch Reactor with Single Reversible Reactiona K

a

R

3

4

5

6

7

8

9

10

1 2 3 4 5 6 7 8 9

0.344 0.095

0.262 0.099 0.014

0.214 0.094 0.020 0.004

0.183 0.083 0.023 0.005 0.001

0.162 0.081 0.023 0.007 0.003

0.147 0.072 0.023 0.008 0.003 0.002 0.002

0.135 0.065 0.022 0.008 0.003 0.002 0.002 0.002

0.126 0.061 0.022 0.008 0.003 0.002 0.001 0.002 0.002

K = number of time-resolved measurements; R = number of polynomials.

are candidates for appropriate DRSM models. We can offer similar observations for K = 8 and 9. These are denoted with bold face numbers. For smaller values of K the evidence that an increase in the value of R does not much affect the value of SSun is not as evident. In such a case the goodness of fit test should be our guidance on the accuracy of the respective DSRM model. We now turn our attention to the GoF test defined in eqs 23, 24, and 25. The related p values are given in Table 3. Table 3. Corresponding F-Test p Values for DRSM of Conversion A ⇄ B in a Batch Reactor with Single Reversible Reactiona K R

3

4

5

6

7

8

9

10

1 2 3 4 5 6 7 8 9

1.00 1.00

1.00 1.00 1.00

1.00 1.00 1.00 0.97

1.00 1.00 1.00 1.00 0.51

1.00 1.00 1.00 1.00 0.99 0.82

1.00 1.00 1.00 1.00 0.90 0.34 0.47

1.00 1.00 1.00 1.00 0.95 0.48 0.47 0.58

1.00 1.00 1.00 1.00 0.70 0.03 0.01 0.01 0.14

Figure 3. Experimental data with 10% noise compared to the DRSM predictions of fractional conversion profiles at two different sets of conditions within the experimental design: x1 = −0.5, x2 = −0.5 on the left, x1 = 1.0, x2 = 0.0 on the right (corresponding to runs 4 and 9, respectively).

due to the substantial (10%) measurement error added to the simulated values. A reduction of this oscillatory behavior can be easily achieved by a further decrease in the value of R to 7, 6, or even 5, something that is supported by the corresponding values in Tables 2 and 3. The oscillatory behavior in Figure 3 points to a future improvement of the DRSM methodology, incorporating knowledge of the examined process to limit these small but unwanted oscillatory elements in the model’s predictions. Overall, the estimated profiles give valuable insight into the shape and magnitude of the response. Now we proceed to examine the efficacy of the DRSM methodology in analyzing processes with more complex stoichiometry and kinetics.

a K = number of time-resolved measurements; R = number of polynomials.

Those cases with a p value smaller than 0.95 are in bold numbers. We observe that we can have an acceptable DRSM model with as few measurements as 6 utilizing a fourth order polynomial (R = 5 total polynomials). The smaller is the value of p, the more accurate the DRSM models are, as the unmodeled sum of squares is the smallest. For any other value of K, several DRSMs can be easily estimated for different values of R and they can be quantitatively assessed from the corresponding p values. Furthermore, a qualitative assessment of the model can be performed by plotting the DRSM predictions against the experimental data used for its development (Figure 3). In all four subplots, we assume that we measure at 10 equidistant points. Each column relates to a separate run, #4 and #9, respectively, for the left and right column. The bottom row utilizes nine polynomials (R = 9) and the top row utilizes eight polynomials. We see that the DRSM represents the data quite accurately, despite the substantial measurement error we introduced. A slight drawback of the method is that the arrived model has a bit of a wavy nature that is not expected in this simple reaction. This has been caused by the wavy character of the data points

5. SEMIBATCH REACTOR WITH A COMPLEX REACTION NETWORK 5.1. Problem Definition. Here, we examine the time evolution of a reaction network with three reactions in a semibatch isothermal reactor with constant volume. The coreactant B is fed in semibatch mode to minimize the amount of byproduct produced. The reaction network is as follows: Rxn1: A + B → C , Rxn2: 2B → D , Rxn3: C → E , G

r1 = k1[A][B]

with k1 = 2 L g mol h−1

r2 = k 2[B]2

with k 2 = 1 L g mol h−1

r3 = k 3[C]

with k 3 = 1 h−1

(36)


Article


met. Since the experimental region shrinks to zero at the end of the batch, we will also impose the constraint w(1) = 0, which yields x0 + x1 + x2 = 0. A constraint must also be imposed on u(τ) to ensure that the amount of reactant B fed is the desired one: ∫ 10{u0(τ) + Δu(τ) w(τ)} dτ = 15 + 5a. This simplifies to

We assume a reactor volume of 10 L, and an initial concentration for A of 1.0 g mol/L. The batch time, tb, will be set at 1.0 h. We assume that we do not know the stoichiometry or the kinetics of the reaction network besides the need to feed B in semibatch mode. It is desirable to maximize the production of the desirable product C, which was examined in prior publications.1,13 The decision variables are the total amount of B to be fed, and the dependence of the feeding profile in time. We set 15 g mol as the reference value for the total amount of B fed, and we will consider a range of interest between 10 and 20 g mol. The material balances defining the model are given below:

( 43x

a = −2

−0.5 ≤ x1 ≤ +0.5,

−0.5 ≤ x 2 ≤ +0.5,

−0.5 ≤ x1 + x 2 ≤ +0.5, −1.5 ≤ 4x1 + 3x 2 ≤ +1.5

d[B] = u(t )/V − k1[A][B] − 2k 2[B]2 , dt

(43)

To illustrate that the DRSM methodology is applicable to any of the various choices for a DoE or a DoDE design, a D-optimal design will be used here rather than a CCD one used in the first example. A D-optimal design with the above constraints consists of a minimum of six experiments to estimate the six parameters of a quadratic model, three additional experiments to assess the GoF statistic, and four replicates (runs #9−13) very near the center point of the region to assess the inherent variability of the process (Table 3). As before, measurement error will be added to the results of the simulated experiments: ye = y(1 + σN(0,1)) with N(0,1) a random normally distributed number with zero mean and standard deviation equal to 1. This time we select σ equal to 0.02 (defined as σ = 0.04/1.96), yielding a simulated measurement error within 4% at 95% confidence. Figure 4 shows the nine unique feeding profiles of reactant B in the DoDE set of experiments. We are interested in developing a DRSM model for product composition (C). The nine distinct product concentration profiles are shown in Figure 5. 5.2. DRSM Estimation. Using the profiles generated above, parameters will be fitted to the following quadratic DRSM:

d[C] = k1[A][B] − k 3[C], dt d[D] = k 2[B]2 , dt d[E] = k 3[C] dt

(37)

with the following initial conditions [B]0 = [C]0 = [D]0 = [E]0 = 0 g mol/L (38)

The total amount of B fed, BT, will be parametrized by introducing decision variable a as follows: BT = 15 + 5a with the constraint −1 ≤ a ≤ +1. We define a nominal feeding profile u0(t), which uses the nominal amount of B fed (15 g mol) and a batch time set to 1.0 h. As the batch time is fixed at 1.0 h, the actual and the dimensionless times will have the same values, so u0(t) = u0(τ). To satisfy these conditions, u0(t) will be constrained by the following equation:

∫0

)

+ x 2 . We observe that of the four decision

variables, a, x0, x1, and x2, and the two variables a and x0 depend on x1 and x2. With the only independent factors being x1 and x2, the previously defined constraints are reformulated to the following four:

d[A] = −k1[A][B] dt

[A]0 = 1.0 g mol/L,

1

1

u0(t ) dt = 15

y(τ ) = β0(τ ) + β1(τ )x1 + β2(τ )x 2 + β12(τ )x1x 2

(39)

+ β11(τ )x12 + β22(τ )x 2 2

There are many possible choices for u0(t), given its timedependent nature, that satisfy the above constraint. Without any other prior process information, it is reasonable to choose the simplest meaningful profile, a linear one. As B is a reactant, we choose a linearly decreasing profile with a value of zero at the end of the batch: u0(τ ) = 30(1 − τ )

As before, K is the number of sampled measurements collected in each of the 13 experiments and R is the number of polynomials Table 4. D-Optimal DoDE Design for Semibatch Reactor with 3-Reaction Network

(40)

The flow rate of B for the experiments will be defined in terms of the nominal profile u0(τ): u(τ ) = u0(τ ) + Δu(τ ) w(τ )

(41)

We select Δu(τ) = 20(1 − τ) and we thus define the experimental region. This implies that −1 ≤ w(τ) ≤ +1. The dynamic factor w(τ) will be parametrized using the first three Shifted Legendre polynomials, following the DoDE methodology. We introduce three new decision variables, x0, x1, and x2, called the dynamic subfactors, to define w(τ ) = x0P0(τ ) + x1P1(τ ) + x 2P2(τ )

(44)

(42)

The constraint −1 ≤ w(τ) ≤ +1 on w(τ) is satisfied if we require that −1 ≤ x0 ± x1 ± x2 ≤ +1. Here again the constraints on the xi variables are selected so that the constraint in w(τ) is

a

H

run

aa

x0a

x1

x2

1 2 3 4 5 6 7 8 9 10 11 12 13

−1.000 0.993 0.333 1.000 −0.027 −0.333 −0.993 0.993 −0.020 −0.020 −0.020 −0.020 −0.020

−0.50 0.33 0.00 0.50 −0.01 0.00 −0.33 0.38 −0.01 −0.01 −0.01 −0.01 −0.01

0.00 −0.50 −0.50 0.00 0.01 0.50 0.50 −0.35 0.00 0.00 0.00 0.00 0.00

0.50 0.17 0.50 −0.50 0.00 −0.50 −0.17 −0.03 0.01 0.01 0.01 0.01 0.01

Dependent on x1 and x2. DOI: 10.1021/acs.iecr.5b03572 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


sum of squares, SSum, as calculated by eq 20. These values are given in Table 5. Table 6 gives the p values of the F-test described in eq 25. The values of Table 5 are smaller than the values of Table 2 because the added error in the experimental data was quite smaller using σ = 0.02, versus the initial value of σ = 0.05. Similarly the acceptable values for p in Table 6 are quite small. It is observed that a fourth order polynomial (R = 5) will do an excellent job when the number of measurements ranges from K = 6 to K = 10. In Figure 6 we present six subplots in which we compare the simulated profiles of product C and the estimated ones from the developed DRSM. Each of the columns represents different experimental conditions defined by the values of the a0, x0, x1, and x2 factors. In the first row of figures we have 6 measurements and we use 5 polynomials in the DRSM. In the second row we measure at 10 time instants and we again use 5 polynomials. The two figures of the third row also utilize 10 measurements but now the DRSM is allowed to try terms for 9 total polynomials. The approximation is excellent.

Figure 4. Nine unique feeding profiles for reactant B, generated by D-optimal DoDE design, as seen in Table 4 (note that the two linear profiles are very close and thus may be visually hard to discern).

6. PENICILLIN FERMENTATION 6.1. Definition of the Process. Here, we consider the penicillin fermentation model stated by Bajpai and Reuss,14 which has been successfully modeled with traditional RSM methodology using the DoDE approach.1 In contrast to the previously discussed batch processes, the penicillin fermentation model is characterized by a more complex dynamic behavior. The model used to simulate the process consists of balances for the volume (V), biomass (x), substrate (s), and product (p) as follows: dV = u(t ) dt xu(t ) dx = μx − dt V (s f − s)u(t ) mssx μx ρx ds =− − − + dt YX / S YP / S km + s V

Figure 5. Nine concentration profiles for product C corresponding to the nine input profiles in Table 4 and Figure 4. Solid, dashed, and dotted lines represent groups of profiles corresponding to initial feed rates of B near 10, 30, and 50 g mol h−1.

dp pu(t ) = ρx − kdp − dt V

(45)

with

used to parametrize each of the βq(τ) functions. A DRSM will be constructed for the time evolution of the product composition C. For each pair of R and K values, the assessment of the accuracy of the developed DRSMs is done as previously by first constructing a table of the fractional values of the unmodeled

⎛ s ⎞ μ = μmax ⎜ ⎟ ⎝ kxx + s ⎠

and

⎡ ⎤ s ⎥ ρ = ρmax ⎢ ⎢⎣ kp + s + (s 2 /kin) ⎥⎦ (46)

Table 5. SS Values for DRSM of Product [C] in Semibatch 3-Reaction Networka K R 1 2 3 4 5 6 7 8 9 a

3

4 −01

3.63 × 10 3.72 × 10−02

5 −01

3.25 × 10 3.02 × 10−02 2.94 × 10−03

6 −01

3.08 × 10 2.54 × 10−02 4.20 × 10−03 7.78 × 10−04

7 −01

2.83 × 10 2.46 × 10−02 4.60 × 10−03 1.12 × 10−03 5.45 × 10−05

8 −01

2.75 × 10 2.37 × 10−02 4.63 × 10−03 1.24 × 10−03 8.18 × 10−05 3.20 × 10−06

9 −01

2.70 × 10 2.37 × 10−02 4.23 × 10−03 1.31 × 10−03 1.07 × 10−04 4.55 × 10−06 1.19 × 10−06

10 −01

2.64 × 10 2.28 × 10−02 4.15 × 10−03 1.30 × 10−03 1.08 × 10−04 5.20 × 10−06 1.23 × 10−06 1.23 × 10−06

2.65 × 10−01 2.24 × 10−02 4.00 × 10−03 1.29 × 10−03 1.16 × 10−04 5.75 × 10−06 1.28 × 10−06 1.28 × 10−06 1.28 × 10−06

K = number of time-resolved measurements; R = number of polynomials. I


Article

Industrial & Engineering Chemistry Research Table 6. Corresponding F-Test p Values for DRSM of Product [C] in Semibatch 3-Reaction Networka K

a

R

3

4

5

6

7

8

9

10

1 2 3 4 5 6 7 8 9

1.00 1.00

1.00 1.00 1.00

1.00 1.00 1.00 1.00

1.00 1.00 1.00 1.00 2.33 × 10−06

1.00 1.00 1.00 1.00 8.06 × 10−06 3.18 × 10−26

1.00 1.00 1.00 1.00 2.89 × 10−05 1.22 × 10−29 6.89 × 10−36

1.00 1.00 1.00 1.00 1.45 × 10−02 5.18 × 10−27 9.60 × 10−38 9.60 × 10−38

1.00 1.00 1.00 1.00 2.22 × 10−05 6.14 × 10−37 1.07 × 10−50 1.07 × 10−50 1.07 × 10−50


We vary the initial biomass concentration between 1 and 2 gx/L, and express it through the decision variable w1 as follows: x(0) = 1.5 + 0.5w1 with −1 ≤ w1 ≤ +1. The initial concentration of the substrate is set to 500 g/L, and similarly we will fix the substrate concentration in the semibatch feed at the same value. We will use a dynamic factor to characterize the time-variant substrate feed, increasing the reactor volume to a final value of 10 L. As in the previous example, we will define a nominal feeding profile u0(t). The final reactor volume of 10 L will impose t a constraint on u0(t) as follows: V(0) + ∫ 0bu0(t) dt = V(tb). This constraint simplifies to

∫0

1

u 0( τ ) d τ =

10 − 7 3 = ≜ θ(L /h), tb 130

τ = t /tb (47)

As before, we assume the simplest meaningful time-dependent profile for u0(τ), a linear dependence on time, of the form u0(τ) = A + Bτ. For such batch fermentations, the initial amount of substrate fed should be kept high in order to promote rapid growth of the biomass. For this model, the feeding profiles will tend toward zero at the conclusion of the batch, u0(1) = 0. We define the nominal feeding profile as 6 u0(τ ) = A(1 − τ ), A = 2γ = 130 L/h . We select δu(τ) = A(1 − τ) such that all feeding profiles examined are between zero and 2u0(t). We define the feeding profile u(τ) in terms of the nominal feeding profile δu(τ) and the dynamic factor z(τ):

Figure 6. Simulated (experimental) data compared to the DRSM predictions of product C concentration profiles at two different sets of conditions: x1 = 0.03, x2 = 0.33 on the left, x1 = −0.09, x2 = −0.26 on the right.

Table 7. Model Parameters for Semibatch Penicillin Fermentation definition μmax ρmax kx kp kin kd km ms YX/S YP/S sf

maximum specific biomass growth rate maximum specific production rate saturation parameter for biomass growth saturation parameter for penicillin production inhibition parameter for penicillin production penicillin degradation saturation parameter for maintenance consumption maintenance consumption rate yield factor, substrate (S) to biomass (X) yield factor, substrate (S) to product (P) feed concentration of substrate (S)

value 0.1 h−1 5.5 × 10−3 gp/(gx h) 6.0 × 10−3 gs/gx 0.1 × 10−3 gs/L 0.1 gs/L 0.01 h−1 0.1 × 10−3 gs/L

u(τ ) = u0(τ ) + δu(τ ) z(τ )

where

− 1 ≤ z (τ ) ≤ + 1 (48)

z(τ) is parametrized here in terms of the first three Legendre polynomials, with decision variables x1, x2, and x3 as weights. z(τ ) = x1P0(τ ) + x 2P1(τ ) + x3P2(τ )

2.9 × 10−2 gs/(gx h) 0.47 gx/gs 1.2 gp/gs 50 g S/L

(49)

As before, the dynamic subfactors in eq 49 are constrained by the following inequalities. −1 ≤ x1 ± x 2 ± x3 ≤ +1

The substrate is fed to the bioreactor in a semibatch mode, by the incoming flow rate u(t). The model parameters used here are those used by Riascos and Pinto,15 given in Table 7. In this process, it is desirable to maximize the product yield. We will design a set of experiments around two input factors: the initial biomass concentration and the substrate feed profile. 6.2. DoDE Design. The batch time tb will be fixed at 130 h, and the initial volume of the biomass mixture will be set to 7.0 L.

(50)

The initial four degrees of freedom, associated with three dynamic subfactors x1, x2, and x3 and the classical factor w1 are reduced to two by the following two constraints. Because δu(1) = 0, we can impose the constraint z(τ) = 0 as there is no reason for it to be any other value. We thus have x1 + x2 + x3 = 0 and x3 = −(x1 + x2). The second constraint is imposed by the reactor volume that has a maximum value of 10 L, resulting in the J


Article

Industrial & Engineering Chemistry Research following constraint that must be satisfied by u(τ): V (0) +

∫0

We observe that they pose a significant challenge for our polynomial approximation. 6.3. DRSM Development. The discontinuous-like dynamic behavior that occurs at different time instances for each run is due to the switch from biomass growth to penicillin production. Efforts to develop a classical RSM to estimate the dependence of this transition time on the experimental conditions have not been very successful so far. This could have helped approximate the penicillin concentration profile with two different sets of polynomials, one per segment. Nevertheless, we will estimate a quadratic DRSM with the following form to represent the penicillin profile over the whole time interval of 130 h as a function of the two input factors A and B.

1

u(τ ) dτ

= V (0) +

∫0

1

[u0(τ ) + Δu(τ )z(τ )] dτ = V (τ )

(51)

∫ 10Δu(τ)

This simplifies to z(τ) dτ = 0, which implies x2 = 3x1. Thus, there are now two independent variables: w1 and x1, and for simplicity, we will denote them with w1 = A and x1 = B. A DoDE design is defined for the two factors and should satisfy the following constraints: −1 ≤ w1( = A) ≤ +1

and

− 0.125 ≤ x1( = B) ≤ +0.125

y(τ ) = β0(τ ) + β1(τ )A + β2(τ )B

(52)

+ β12(τ )AB + β11(τ )A2 + β22(τ )B2

The experimental design consists of a minimum of six experiments to estimate the parameters of a quadratic model, three additional experiments to assess the GoF statistic, and four replicates at the center point to assess the normal variability of the process. The values defining these experiments are listed in Table 6. A 4% error is added to the simulated profiles (σ = 0.02). The product profiles of different runs are presented in Figure 7.

We will now investigate a suitable combination of R and K values which results in an accurate model. We will assume that we are measuring the penicillin concentration at the K instants, where K ranges between 13 and 20. For each such value of K, several values of R will be used and the respective DRSM will be evaluated as in the prior examples. In Table 9 we present the dependence of the SSum values on the values of K and R. For each value of K considered in the range of 13 to 20 there are several values of R that provide the smallest values for SSum. They are denoted with bold face numbers. Table 10 provides the p values of the GoF test defined in eq 22. We note that there is a multitude of K and R value pairs with p values less than 0.95 indicating a reasonable accuracy of the respective DSRM. In Figure 8 a set of six subplots is given in which the penicillin product profiles estimated by the DRSM are compared to the simulated ones, illustrating their substantial accuracy of the DRSM model. Each column of subplots relates to two different sets of experiments corresponding to the indicated two sets of values for the A and B factors. In the first row we examine the case of a DRSM with R = 12 over 15 measurements, approximating the profile with good accuracy. In the second row of subplots we have the best DRSM approximation when R = 12 and the number of measurements is 20. With the same number of measurements and an increase in the values of R the accuracy deteriorated at initial and final times as shown in the third row of subplots. Interestingly enough, this could not have been predicted from the corresponding values in Tables 9 and 10. Of further interest is that there is no loss of predication accuracy around dimensionless time equal to 0.3 when the derivative of the penicillin profiles appears to have a step change. This indicates that the desirable results of Tables 9 and 10 cannot by themselves provide the final judgment on the quality of the DRSM and plots like the one in Figure 8 can help us substantially in selecting a suitable value of R for the most desirable DRSM model.

Table 8. DoDE Design for the Semibatch Penicillin Fermenter

a

run

w1 = A

x1 = B

x2a

x3a

1 2 3 4 5 6 7 8 9 10 11 12 13

1.00 1.00 1.00 −1.00 −1.00 −1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

−0.125 0.125 0.000 −0.125 0.125 0.000 −0.125 0.125 0.000 0.000 0.000 0.000 0.000

−0.375 0.375 0.000 −0.375 0.375 0.000 −0.375 0.375 0.000 0.000 0.000 0.000 0.000

0.500 −0.500 0.000 0.500 −0.500 0.000 0.500 −0.500 0.000 0.000 0.000 0.000 0.000

(53)

Dependent on factors w1 and x1.

7. POTENTIAL USAGES OF THE DRSM The usages of the DRSM model are expected to be as widespread as the usages of the classical RSM methodology. In a preliminary investigation one might develop DRSM models that include only two factor interaction terms, as in eq 2 but with time dependent βq(t). This helps one investigate the change(s) in the ranges of the DoE or DoDE factors that will bring us closer to the optimal condition of the process. Because of the dynamic character of the DRSM, as compared to the RSM models, time needs not be a

Figure 7. Nine unique profiles of penicillin concentration, simulation based on the DoDE design shown in Table 8. K


Article

Industrial & Engineering Chemistry Research Table 9. SS Values for DRSM of Penicillin Producta K

a

R

13

14

15

16

17

18

19

20

8 9 10 11 12 13 14 15 16 17 18 19

0.0009 0.0008 0.0006 0.0005 0.0004

0.0008 0.0006 0.0005 0.0005 0.0004 0.0004

0.0008 0.0006 0.0005 0.0004 0.0004 0.0004 0.0004

0.0008 0.0006 0.0005 0.0005 0.0004 0.0004 0.0003 0.0004

0.0008 0.0006 0.0004 0.0004 0.0004 0.0003 0.0003 0.0003 0.0003

0.0008 0.0007 0.0005 0.0005 0.0005 0.0004 0.0004 0.0004 0.0004 0.0004

0.0007 0.0006 0.0005 0.0005 0.0004 0.0004 0.0004 0.0004 0.0004 0.0003 0.0003

0.0008 0.0006 0.0005 0.0004 0.0004 0.0004 0.0004 0.0003 0.0003 0.0003 0.0004 0.0003


Table 10. Corresponding p Values for the DRSM of Penicillin Producta


factor in the DoE or the DoDE experimental design, substantially reducing the number of experiments. Deciding on whether the batch duration should be increased or decreased from the value initial considered can be easily deduced. After the preliminary experiments guide us to the region where the optimal operation is expected more detailed DRSM models, like in eq 3, can be estimated and the process optimized both with respect the values of the factors xi as well as time. One can thus claim that the DRSM model enables dynamic optimization in comparison to the static optimization possible with the traditional RSMs. Because the DRSM approach models the overall dynamic behavior of the process examined, one can envision a multitude of other usages besides optimization; for example, in the development of recursive dynamic models and their usage in modelbased control. For reacting systems, batch or semibatch, the development of a DRSM for each of the reacting species can be instrumental in the discovery of the stoichiometry of the active reactions and the estimation of the respective kinetic models. The possibilities are indeed numerous.

Figure 8. Simulated (experimental) data compared to the DRSM predictions of penicillin concentration profiles at two different sets of conditions: A = −0.626 and B = −0.003 on the left, A = 0.323 and B = −0.075 on the right.

8. CONCLUSIONS We have introduced a new data-driven methodology, named dynamic response surface methodology (DRSM), for the modeling of batch processes. It generalizes the traditional response surface methodology by making it possible to incorporate time-varying outputs. This methodology was demonstrated in conjunction with another generalization of the design of experiments framework, the newly proposed design of dynamic experiments (DoDE) methodology.1 The DRSM is also applicable in classical DoE problems when the output is measured at several time instants. We have introduced two statistical criteria to help identify the appropriate number of interpolating polynomials, the value of R, which is expected to produce the most accurate model for a given number of measurements K per batch. A visual inspection of the estimated profiles which meet the postulated statistical tests can avoid any that exhibit unnecessary small oscillatory overtones on the estimated profiles. The present numerical algorithm does not warrant against this happening. To illustrate the power and versatility of this new technique, the DRSM methodology was applied to a set of three sets of in silico experiments, corresponding to three different processes. The presented results clearly demonstrate that the DRSM provides an excellent representation of the dynamic behavior of the processes examined. With such a detailed model at hand, one can readily proceed in optimizing the operations of the process

K R

13

14

15

16

17

18

19

20

8 9 10 11 12 13 14 15 16 17 18 19

1.00 1.00 1.00 1.00 0.98

1.00 1.00 1.00 1.00 0.98 0.98

1.00 1.00 0.96 0.95 0.79 0.81 0.88

1.00 1.00 0.94 0.95 0.89 0.62 0.56 0.61

1.00 1.00 0.94 0.90 0.85 0.66 0.46 0.48 0.62

1.00 1.00 1.00 1.00 0.99 0.98 0.98 0.94 0.90 0.96

1.00 1.00 1.00 1.00 1.00 0.99 1.00 0.99 0.99 0.96 0.91

1.00 1.00 0.99 0.98 0.92 0.87 0.88 0.78 0.69 0.72 0.80 0.70

a

L


Article

Industrial & Engineering Chemistry Research with respect to a performance characteristic. Furthermore, the novel dynamic character of the DRSM also enables the selection of the optimal duration of the batch, as it will be described in a forthcoming publication. It was amply obvious in a subplot of Figure 8 that the method needs to be constrained from estimating oscillatory profiles when none are justified from the nature of the process. The current DRSM methodology assumes a uniform level of uncertainty applicable to all process measurements. This will not be true if a given process output is measured using two different measuring tools, one slow and accurate, like HPLC, and another fast and not as accurate. Furthermore, one needs to more systematically account for missing data as well as for measurements, which are not performed in equidistant time instants and for batch experiments which have varying durations.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: 617-627-2573. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Georgakis, C. Design of Dynamic Experiments: A Data-Driven Methodology for the Optimization of Time-Varying Processes. Ind. Eng. Chem. Res. 2013, 52 (35), 12369−12382. (2) Nolan, R. P.; Lee, K. Dynamic model for CHO cell engineering. J. Biotechnol. 2012, 158 (1−2), 24−33. (3) Myers, R. H.; Montgomery, D. C.; Anderson-Cook, C. M. Response surface methodology: process and product optimization using designed experiments, 3rd ed.; Wiley: Hoboken, N.J., 2009; p xiii, 680. (4) Box, G., Draper, N. R. Empirical Model-Building and Response Surfaces; John Wiley & Sons, Incorporated: New York, 1987; p 688. (5) Andersson, M.; Adlercreutz, P. Evaluation of simple enzyme kinetics by response surface modelling. Biotechnol. Tech. 1999, 13 (12), 903−907. (6) Bas, D.; Boyaci, I. H. Modeling and optimization I: Usability of response surface methodology. J. Food Eng. 2007, 78 (3), 836−845. (7) Beg, Q. K.; Saxena, R. K.; Gupta, R. Kinetic constants determination for an alkaline protease from Bacillus mojavensis using response surface methodology. Biotechnol. Bioeng. 2002, 78 (3), 289− 295. (8) Gholami-Seyedkolaei, S. J.; Mirvaghefi, A.; Farahmand, H.; Kosari, A. A. Optimization of recovery patterns in common carp exposed to roundup using response surface methodology: Evaluation of neurotoxicity and genotoxicity effects and biochemical parameters. Ecotoxicol. Environ. Saf. 2013, 98, 152−161. (9) Naik, L.; Mann, B.; Bajaj, R.; Sangwan, R. B.; Sharma, R. Process Optimization for the Production of Bio-functional Whey Protein Hydrolysates: Adopting Response Surface Methodology. Int. J. Pept. Res. Ther. 2013, 19 (3), 231−237. (10) Tuan, N.; Hazimah, A. H.; Wan, M. Optimization of reaction conditions for enzymatic synthesis of palm fatty hydrazides using response surface methodology. J. Oleo Sci. 2012, 61 (5), 297−302. (11) Montgomery, D. C. Design and Analysis of Experiments, 8th ed.; John Wiley & Sons: 2013. (12) Montgomery, D. C.; Runger, G. C.; Hubele, N. F. Engineering Statistics, 5th ed.; Wiley: Hoboken, N.J., 2011; p xxii, 515. (13) Troup, G. M.; Georgakis, C. Process systems engineering tools in the pharmaceutical industry. Comput. Chem. Eng. 2013, 51, 157−171. (14) Bajpai, R. K.; Reuss, M. A Mechanistic Model for Penicillin Production. J. Chem. Technol. Biotechnol. 1980, 30 (6), 332−344. (15) Riascos, C. A. M.; Pinto, J. M. Optimal control of bioreactors: a simultaneous approach for complex systems. Chem. Eng. J. 2004, 99 (1), 23−34.

M


Dynamic Response Surface Models: A Data-Driven Approach for the

Recommend Documents