Using Taguchi's Method and Orthogonal Function Approximation to

2226

Ind. Eng. Chem. Res. 2002, 41, 2226-2237

Using Taguchi’s Method and Orthogonal Function Approximation to Design Optimal Manipulated Trajectory in Batch Processes Junghui Chen* and Rong-Guey Sheui Department of Chemical Engineering, Chung-Yuan Christian University, Chung-Li, Taiwan 320, Republic of China

A novel experimental scheme is proposed to integrate orthogonal function approximation with Taguchi’s method (orthogonal array) for designing the optimal manipulated trajectory of a batch process. The orthogonal function approximation finds a set of orthonormal functions as the basis to represent the batch trajectory. The optimal trajectory can be obtained if the location of the coefficients of the orthonormal functions is properly adjusted in the function space. The Taguchi approach is used to design and analyze the effect of each coefficient on reaching the optimal objective (quality) function. Because the coefficients are implicitly related to the objective function, they simply vary over two levels in a systematic way, and they would be moved into the optimal design condition. A search procedure for the optimal design coefficients is also proposed. As opposed to model-based design, the proposed method utilizes the simplicity of Taguchi methods to determine the potentially available knowledge of dynamic batch processes. For its potential applications, this combinational method is demonstrated through two simulation case studies. 1. Introduction As the current market is highly competitive, product life cycles are becoming shorter and shorter. Time-tomarket, quality, and individualization are the main competitive factors in the introduction of new products into the global marketplace. This is especially true for batch processes dealing with small volumes of specialized, value-added, and high-quality chemical products (e.g., polymer reactors, pharmaceutical manufacture, biochemical reactors, batch distillations, and etching processes). Batch processes are characterized by the prescribed processing of raw materials into products in finite duration. The operation of batch processes is quite different from that of continuous processes, because the former is always in unsteady-state operation. From an engineering point of view, there are two major objectives for operating batch processes: optimization of the yield rate and control of the product quality. Generally, the major control objective is to optimize the performance index at the end of the batch run. With proper design of the product quality set point, the desired performance index can be obtained. Therefore, obtaining an optimal reference trajectory for the manipulated variable is crucial to batch process operation. To solve the optimized set-point problem, Pontryagin’s maximum principle has often been used.1,2 However, this optimization solution is computationally expensive for most practical problems. Alternatively, the optimal control problem can be converted into a dynamic programming problem. The differential equation problem can be transformed into nonlinear algebraic equation problem by a discretization method (e.g., orthogonal collocation, finite difference).3 However, the convergence to solution can be slow, or it might even fail. Iterative dynamic programming was developed to solve limited dynamic programming problems.4,5 However, the fatal * To whom all correspondence should be addressed. Email: [email protected]. Tel.: 886-3-4563171 ext 4107. Fax: 886-3-4563171.

flaw of these methods lies in the fact that the above optimization can be obtained only with a reliable process model. For real processes (e.g., polymer manufacture, biochemical reactors, and etching processes), such explicit knowledge is hard to acquire. Take a batch reactor, for example. One of the greatest difficulties is determining the reaction rate of the system, including the reaction mechanism and the related parameters. To remedy the mismatch between the actual process and the model, several methods have been proposed. Luus (1974) derived the optimal control vector by employing a direct search procedure on the feedback gain matrix.6 Cuthrell and Biegler (1989) used a sequential quadratic programming approach to solve the optimal control problem.7 Other approaches have been proposed to reduce the sensitivity to model uncertainly via a state feedback law and obtain the optimal trajectory.8-10 Still, a mathematical description of this kind of system should be derived in a deductive manner using first-principles models. With shorter product life cycles, one can not afford to spend much time in determining detailed mathematical models. It might even be considered unrealistic or impossible to obtain a sufficiently accurate model in this way. The effectiveness of optimal trajectories in real applications would also be degraded. Some researchers have proposed on-line updating of the model using past and present process measurements and computation of the optimal trajectory using the updated model.11,12 Judging from the above discussion, a practical design should be inferred from a set of data collected from experiments on the system of interest. Even though this can be a useful short cut for deriving mathematic models, engineers need to have a free hand to conduct practical experiments that bring the system around its entire range of operation. With product life cycles becoming shorter and with the increasing emphasis on time-to-market, a systematic design strategy is needed to collect experimental data and explore possible quality-improvement designs. Randomly generating experi-

10.1021/ie010457x CCC: $22.00 © 2002 American Chemical Society Published on Web 04/05/2002

Ind. Eng. Chem. Res., Vol. 41, No. 9, 2002 2227

model-based design. The effectiveness of the proposed method is shown through two detailed simulation studies that illustrate the potential applications of the proposed method. 2. Problem Formulation The problem considered here is the end-point optimization of quality in batch processes. It can be mathematically formulated as

max J(y(tf)) u

Figure 1. Black-box design.

mental data that lack appropriate information for quality improvement wastes time and effort. Thus, proper procedures for experimental design are strongly needed. More importantly, the experiments should be designed to determine the optimal trajectory of the batch process and the effects of variations on the different batch design trajectories. These procedures can produce a set of guidelines that would be beneficial to practitioners. In this research, a new experimental scheme using Taguchi’s method (orthogonal array) and orthogonal function approximation is proposed. It is developed to design the optimal manipulated trajectory of the batch process. With the available experimental data, it can determine the possible optimal set-point trajectory. A systematic search technique is also established to diminish the size of the search region of each experimental design without abundant experimental time. Assume that the trajectory is a continuous, smooth curve. It is possible to find a set of orthonormal functions as the basis to represent the process trajectory. The trajectories of process measurements in batches are mapped onto the new feature coefficients in the function space. The optimal trajectory can be obtained if the location of the coefficients is properly adjusted in the function space. Once the coefficients for the approximated trajectory have been chosen, the next step is to define a fractional factorial experiment using the Taguchi approach to analyze the effect of each design coefficient on the optimization of an objective (quality) function and on the search for the optimal set of design coefficients. The design procedure can also minimize variations (such as external noise, batch-to-batch variations, and internal noise) to improve quality. By using the orthogonal array, the number of experiments can be reduced, and the effects of each design coefficients can be identified by performing one set of experiments. The method includes a set of tables that enables the coefficients to be investigated in a minimum number of trials. The coefficients with the greatest impact on the objective function are identified. These coefficients simply vary over two levels in a systematic way. Results are analyzed in terms of the design operation region. Control over achieving the objective is best obtained by changes in these primary coefficients. The true power of Taguchi methods comes from their simplicity of implementation. Rather than using model-based design, the proposed method determines the potentially available knowledge of the batch system to explore the optimal feasible region. This reduces the experimental time taken when the study is conducted. To our best knowledge, this is the first investigation for batch processes of the optimal trajectory based on experimental design, rather than

(1)

subject to

x3 ) f(x,u), 0 e t e tf y ) h(x,u)

(2)

and

L e y(t) e U where J is an objective function; x is the state vector; u is the manipulated variable; tf is the duration of each batch run; U and L are the lower and upper bounds of the manipulated trajectory, respectively; and f and h give the description of the dynamic model of the batch process. Numerous researchers have proposed methods to solve these problems using mathematical models based on the known physical principles (mass balances and energy balances) and engineering empirical relations. In fact, optimization based on a reliable process model is difficult. In contrast, engineers usually treat the mathematical model as a black-box problem. Traditionally, researchers have used statistical techniques to solve such optimization problems. The experimental design and Taguchi’s method,13 for examples, are two common approaches to the optimization of factor design in practice. In Figure 1, a black box represents either a product or a process, and various types of input and output variables are applied to different fields. It requires the determination of factor settings for manipulating inputs to yield an optimal quality performance via an orthogonal array and a data analysis procedure. Many applications of the Taguchi method have been presented, including manufacturing in the automobile, microelectronic, chemical process,14 and even service industries.15 However, the Taguchi method is not a panacea for all parameter design problems. As for batch design, the manipulated variable follows the prespecified trajectory within the course of a batch. It is not practical to add the design experiments to the nominal operating conditions at every time interval for a wide range of variation in the batch trajectories. It is hard to determine how often the process trajectory should be changed to obtain the desired quality. Prior knowledge about the duration is needed before a more appropriate trajectory experiment can be built. 3. Representation of the Trajectory Using an Orthogonal Function Curve-Fitting Technique In batch processes, the time evolution of a batch follows a prespecified trajectory. There are several wellknown mathematical methods for parametrization of the trajectory, such as Fourier series or Laguerre

2228

Ind. Eng. Chem. Res., Vol. 41, No. 9, 2002

functions. In this kind of approach, the coefficients of the function are estimated by dynamic optimization based on the process model. The desired trajectory can then be obtained.16,17 However, such approaches are useful for only limited applications as batch processes are highly nonlinear, time-varying, and seriously interconnected with uncertainties. It becomes a very difficult task to build the process model with first-principles equations that accurately describe the process behavior. Thus, in this research, a simple Taguchi scheme with orthogonal functions is developed for the optimization of dynamic batch processes. The conventional orthogonal method that uses powerseries expansions is efficient and reliable for function approximation. Let {φ0(t), φ1(t), ..., φN-1(t)} be a set of orthogonal functions over Ω ) [a, b]. From the definition of orthogonality, these functions has the following property

〈φi, φj〉 )

∫Ωφi(t) φ(t)j dt ) {0, 1,

i*j i)j

(3)

Hence, the linearly independent φi(t) functions span the function space and constitute a set of bases for the function space. By the approximation theorem, a function f(t) that is assumed to be continuous and squareintegrable in the range Ω ) [a, b] can be spanned in terms of these bases N-1

f(t) ≈ f(CN,t) )

∑ cnφn(t)

(4)

n)0

where the coefficients are CN ) {cn}n)0,1,...,N-1. This combination is called the best approximation to f(t) with respect to the linear combinations of {φn}n)0,1,...,N-1. CN ) {cn} is the projection of f(t) onto each basis function

cn )

∫Ωf(t) φn(t) dt

(5)

The coefficients, CN ) {cn}, are also called the Fourier coefficients of f(t) with respect to {φn}n)0,1,...,N-1. The advantage of orthonormal function approximation is that the old coefficients {cn}n)0,1,...,N-1 keep the same values even if a new term is added, and these coefficients are uncorrelated. Because the process variable is only measured at times t1, t2, ..., tK, by using the Gram-Schmidt method, the stepwise procedures can compute the coefficients ci

Table 1. Description of Two-Level Taguchi Experimental Design c1 trial

c2

c3

value level 1 level 2 level 1 level 2 level 1 level 2

1 2 3 4 avg effect

mapped onto the coefficient variables {cn} in the function space. This point in the function space spanned by N independent basis functions represents one batch run. Therefore, if the coefficients are properly adjusted in the function space, the optimal trajectory should be obtained. 4. Design of Orthogonal Coefficients Using the Taguchi Method To provide scientific discipline and conduct batch experimental design and statistical calculation, Taguchi’s method is used here. The goal of Taguchi’s method in this research is to determine which of the coefficients (i.e., orthogonal coefficients) is the most effective in the optimization of the set-point trajectory for batch product quality as the shape of the trajectory varies with the coefficients. 4.1. Taguchi’s Method of Experimental Design. A strategic design experiment can expose the coefficients to varying environmental conditions and determine the best value for a coefficient. Taguchi’s method, which uniquely harmonizes the statistical design of experiments and analysis of variation methods into industrial design and production, can be thought of very much as a prototyping methodology. The coefficients from the orthogonal function can be assigned to the columns of an orthogonal array. Suppose an experimental objective outcome J is a function of several coefficients [c1 c2 ... cn] whose values can be controlled. Write J ) f(c1,c2,...,cn). The adjusted coefficients [c1 c2 ... cn] are called factors. The goal is to find the factors that will have the most beneficial effect on J when the factors change. This can be done by varying each factor independently of the others and recording the corresponding change in J. The Taguchi methods are disciplined ways of varying two or more factors

e0(tk) ) f(tk) cn ) [ΦnTΦn]-1ΦnTEn

(6)

en+1(tk) ) en(tk) - cnφn(tk), k ) 1, 2, ..., K; n ) 0, 1, ..., N - 1 where

En ) [en(t1) en(t2) ‚‚‚ en(tK)]T and

Φn ) [φn(t1) φn(t2) ‚‚‚ φn(tK)]T By the orthogonal function approximation, the measurement variable {f(tk)} in each batch trajectory is

Figure 2. Hierarchical regions of the coefficient structure. The selection order of the coefficients is from top to bottom. The bold line within each horizontal line represents the maximum operating region, and the gray, the selected operating region.


Figure 3. Region elimination method that explores the possible optimal outcomes of the quality objective function.

simultaneously. In one of the simplest Taguchi methods, each individual factor varies between two levels (Table 1). Table 1 depicts the Taguchi experiment using three factors to explain the design strategy. The factors are denoted by c1, c2, and c3, each of which has two level values: level 1 and level 2. There are four trails of the experiment, each using level 1 and level 2 combinations in the columns c1, c2, and c3 (without shading) of the row containing the trial number. For example, in trail 3, the value of c1 is level 2, that of c2 is level 1, and that of c3 is level 2. The value of the experiment for each trail is recorded in the value column. Averages are taken of the white boxes in each column. The values in the effect row are found by subtracting the average outcome for level 1 values from the average outcome for level 2 values in the columns c1, c2, and c3. This provides a numerical value for the average effect on the outcome

of moving a factor from its level 1 value to its level 2 value. Note that Taguchi has tabled 18 basic orthogonal arrays,18,19 including two-level arrays, three-level arrays and mixed two- and three-level arrays. For easy application of the design strategy, only two-level arrays are used here. If additional factors (ci, i ) 4, 5, ...) are needed, the corresponding maximum number of factors should be properly selected in two-level arrays. Taguchi used a single measure, the signal-to-noise (S/ N) ratio, that combined both the predictable performance and the variance of unpredictable performance to evaluate the variation of the system’s performance. Because there are three different quality-loss functions available depending on the type of characteristic functions, the maximum-value method is used here to conduct the transformation of the quality-loss function and obtain a consistent objective function

2230


(i) cp ) /cp

(1) smaller is better J)-

1

n

∑yi2 ni)1

(7)

lower bound of cq

[ [

N-1

(2) larger is better J)

1

n

1

∑ ni)1

(8)

2

upper bound of cq

(3) nominal is better J)-

1

n

∑

ni)1

2

(yi - ys)

2

cU,1 ) min q

∑ cnφn(t) e U

0 e t e tf

(10)

where L and U are the lower and upper bounds, respectively, of the manipulated trajectory. If the values of ci, i ) 0, 1, ..., p - 1, p + 1, ..., N - 1, are fixed except for cp, the upper and lower bounds of cp can be obtained as follows:

lower bound of cp

[ ] [ ] N-1

L-

cL,1 p ) max upper bound of cp

(11)

N-1

cU,1 p ) min

ciφi(t) ∑ i)0 i*p

φq(t)

(14)

lower bound of cq

[ [

N-1

i*p

cL,2 q ) max upper bound of cq

] ]

ciφi(t) - /cpφp(t) ∑ i)0

L-

φq(t)

N-1

U-

(15)

ciφi(t) - /cpφp(t) ∑ i)0 i*p

φq(t)

(16)

U,2 Thus, cq ∈ [cL,2 q ,cq ] if cp ) /cp and the rest of the values of ci are fixed will satisfy L e f(t) e U. Therefore, the operating range [/cq,/cq] should be selected subject to U,1 L,2 U,2 [/cq,/cq] ⊂ ([cL,1 q ,cq ] ∩ [cq ,cq ]) (see Figure 2) to satisfy L e f(t) e U. Following the same procedure, the range of the coefficient cs, [/cs,/cs], can be defined as

2M-1

U,i ∩ [cL,i s ,cs ]

i)1

φp(t)

U-

i*p

[/cs,/cs] ⊂

ciφi(t) ∑ i)0 i*p


(ii) cp ) /cp

cU,2 ) min q

n)0

N-1

(13)

U,1 / Thus, cq ∈ [cL,1 q ,cq ] if cp ) cp and the rest of the values of ci are fixed will satisfy L e f(t) e U.

N-1

L e f(CN,t) )

φq(t)

U-

(9)

where n is the number of repetitions and ys is the target response. Now, we can obtain the maximum of J with minimum error variance for any combination of factor levels. Accordingly, we can analyze the data by using the J value as the objective property and making levelto-level comparisons for the factors assigned by an orthogonal array. 4.2. Upper and Lower Bounds of the Orthogonal Coefficients. Before considering the levels of Taguchi’s method, the upper and lower bounds of the coefficients should be computed to keep the approximated function within the feasibility region during each batch run. For operation and process safety, it is common to establish upper and lower bounds for the manipulated variable in many batch processes. The bounds can be obtained from the approximated trajectories, provided that the number of coefficients in the approximated function is sufficient.

i*p

cL,1 q ) max

yi

] ]


L-

(12)

φp(t)

Thus, if the operating range [/cp,/cp] is selected (Figure U,1 2), the condition [/cp,/cp] ⊂ [cL,1 p ,cp ] must be satisfied to ensure L e f(t) e U. Likewise, the upper and lower bounds of cq, q * p, can be obtained if the other coefficients are fixed for the following two cases:

Figure 2 shows the hierarchical plot for the range U,z coefficients. cL,z represent the upper and lower s and cs bounds, respectively, in the interval z for cs, and [/cs,/cs] represents the selected operating range for cs. Note that Figure 2 indicates that the selected operating ranges marked by the gray lines decrease in size on going from the upper selected coefficient to the lower one. As for the selection order of the coefficients, data collected from the orthogonal experimental array by effect analysis can be used. The effect value can separate the total variability of the data into the contributions from each of the coefficients. This can help determine which of the coefficients has a significant effect. The more significant effect the coefficient of the quality function has, the earlier the corresponding orthogonal function should be selected.


4.3. Optimization Search. A two-phase optimization procedure is developed to obtain the optimal coefficient combinations. First, maximize the J value to reduce the variation of the performance. Second, bring the performance onto target by adjusting factors. Phase 1 has been discussed previously in this paper. In this subsection, the region elimination method is proposed for phase 2. Region elimination methods mean that the linear search for a coefficient can delete a calculated portion of the range of CN ) {cn} from consideration on each successive stage of the search for an extremum of the optimal quality J. Two basic elements underlying the region elimination method are interval location and interval reduction. The purpose of interval location is to determine the initial interval that contains an optimal location. Interval reduction then is used to reduce the size of the initial interval until the desired quality has no significant change within the remaining interval. The following procedure illustrates how to locate the optimal point. Assume that the objective function is unimodal. This means that only a single extremum exists. First, evaluate the quality index by adjusting the coefficient ci at the initial two-level points, represented by /ci,j and /ci,j in Figure 3. Note that j represents the jth run. These two points are two corresponding levels of the coefficient ci in the experimental design. There are two possible outcomes for the quality of the two points: f(/ci,j) g f(/ci,j) and f(/ci,j) < f(/ci,j). In reality, the equal condition rarely occurs because of the round-off error inherent in digital computations. The search procedure for the two conditions can be the same except using opposite searching directions, so only the case f(/ci,j) > f(/ci,j) is discussed here. In Figure 3a, the region at the right of /ci,j contains values of f(ci) > f(/ci,j). The next step is to evaluate the quality at a third point located to the right of /ci,j, represented by ci,k, which is a certain distance from the point having the higher quality (Figure 3b1). The quality of the third location would have two possible outcomes: f(ci,k) > f(/ci,j) and f(ci,k) < f(/ci,j). In Figure 3c1, f(ci,k) > f(/ci,j). This indicates that the region to the right of ci,k must have much better values. Therefore, the previous procedure is repeated to obtain a new point in the search direction where the quality is enhanced. In Figure 3c2, f(ci,k) < f(/ci,j). At this point, we know that the maximum is bracketed by the two points /ci,j and ci,k. In these subplots, the gray background color shows the possible optimal region. The bold line with an arrow represents the level location of the current experimental design, the bold line indicates the optimal location of the experimental run up to now, and the dashed line(s) indicate(s) the location(s) of the previous experimental point(s). To reduce the size of the interval [/ci,j, ci,k], an internal point ci,l located either to the right or to the left of /ci,j is evaluated (Figure 3d1 or d2, respectively). Three possible conditions for each side can occur, as shown in Figure 3e1-e3 and Figure 3e4-e6, respectively. If f(/ci,j) > f(ci,l) > f(ci,k) (Figure 3e1) (or f(/ci,j) < f(ci,l) < f(*ci,j), Figure 3e6), then the optimum occurs in the interval [/ci,j,ci,l] (or[ci,l,ci,k]). Conversely, if f(/ci,j) < f(ci,l) and f(ci,l) > f(ci,k) (Figure 3e3) (or f(/ci,j) < f(ci,l) and f(ci,l) > f(/ci,j), Figure 3e4), then the optimum occurs in the interval [/ci,j,ci,k] (or[/ci,j,/ci,j]). In the above two cases, the same procedures are repeated to reduce the inter-

vals. The conditions depicted in Figure 3e2 and e5 do not arise under the assumption of unimodality. If the boundary of the coefficient is assigned, the optimal search rules for the boundary condition need to be included. That is, f(/ci,j) > f(/ci,j), and /ci,j reaches the upper boundary of its coefficient. This means that the optimal region falls in the interval [/ci,j, /ci,j]. The new point ci,k is located between these boundaries (Figure 3b2). Figure 3c3-c5 shows three possible outcomes for the new point. In Figure 3c3, the quality at the new point is better than f(/ci,j) and f(/ci,j). This indicates that the optimum still occurs in the interval [/ci,j, /ci,j]. The interval should be reduced. If the quality at the new point is worse than f(/ci,j) and f(/ci,j) (Figure 3c4), then the optimum can occur in two possible intervals: [/ci,j, ci,k] or [ci,k, /ci,j]. However, the conditions in Figure 3c4 do not arise under the assumption of unimodality. If f(/ci,j) > f(ci,k) > f(/ci,j) (Figure 3c5), then the optimum occurs in the interval [ci,k, /ci,j]. The same procedure is repeated to search the interval location. To determine the optimal set-point trajectory of the batch process, Taguchi’s method and orthogonal function approximation are integrated into an optimal search approach to optimize a single quality response. The above procedure of the proposed approach is summarized as follows: Step 1. Define the operating boundaries of the manipulated trajectory of the batch process (U and L). Step 2. Set a reference operating trajectory based on the prior knowledge. Select the maximum number of orthogonal function terms to be used to approximate the operating trajectory (Nmax) and the maximum number of runs to be conducted in the batch experiments (rmax). Set r ) 1. Note that the Legendre polynomial20 is employed here because its finite interval approximation is suitable in this batch study. Step 3. Approximate the operating trajectory using the proper number of orthogonal functions CN ) {cn}n)0,1,...,N-1 based on eq 4. Step 4. Compute the boundaries of each coefficient to satisfy the operating region (eq 10). Step 5. Construct the Taguchi experimental array and the levels of each coefficient. Step 6. Use the region elimination method and the effect analysis of the Taguchi table to search the optimal region of the quality response. Step 7. Select the current optimal point (Jopt)r in the current experimental design. Step 8. If (Jopt)r is converged to a certain threshold value, then stop the experimental design; otherwise, adjust the levels of each coefficient and go to step 9. Step 9. If r > rmax, then stop the experimental design because of the limitation on the maximum number of experimental runs; elseif the optimal quality response in the current run is significantly improved at a certain threshold ratio (R) compared with the previous run, that is, if

(Jopt)r - (Jopt)r-1 > R(Jopt)r-1

2232


Table 2. Feasibility Regions of the Coefficients for Example 1 run 1 U L

run 2

run 4

run 8

c0

c1

c0

c1

c0

c1

c0

c1

c2

-8.1650 8.1695

-4.9365 12.0075

-3.8532 10.9243

-5.0553 12.7334

-9.9313 9.9313

-4.2373 3.9896

-9.1789 9.1789

-0.9585 8.0295

-9.2828 9.2828

run 10 U L

run 12

c0

c1

c2

c3

c0

c1

c2

c3

0.4077 5.6151

-11.5890 5.8946

-1.6118 1.1859

-9.8434 3.1597

2.0289 5.0422

-12.7591 5.5489

-0.8463 0.3845

-9.775 3.2281

Table 3. First Run of Taguchi Experimental Design for Example 1 c0 trial

value

level 1 0.0000

c1 level 2 5.3149

level 1 0

level 2 2.5979

1 2 3 4

-9.00 × 104 -9.00 × 104 -9.00 × 104 -5.54 × 104 -5.54 × 104 -5.54 × 104 -6.35 × 104 -6.35 × 104 -6.35 × 104 -2.89 × 104 -2.89 × 104 -2.89 × 104

avg

-5.95 × 104 -7.27 × 104 -4.66 × 104 -7.68 × 104 -4.21 × 104 2.61 × 104

effect

3.46 × 104

Table 4. Second Run of Taguchi Experimental Design for Example 1 c0 trial

value

level 1 5.3149

c1 level 2 8.0296

level 1 2.5979

level 2 5.9403

1 2 3 4

-2.89 × 104 -2.89 × 104 -2.89 × 104 -2.01 × 104 -2.01 × 104 -2.01 × 104 -5.58 × 104 -5.58 × 104 -5.58 × 104 -4.70 × 104 -4.70 × 104 -4.70 × 104

then go to step 3 and set r ) r + 1;

avg

-3.80 × 104 -2.45 × 104 -5.14 × 104 -4.24 × 104 -3.36 × 104

elseif N g Nmax, then

effect

Figure 4. Desired trajectory for example 1.

stop the experimental design due to the limitation on the maximum number of orthogonal functions; else increase the number of terms of the orthogonal functions, N ) N + 1. Go to step 3.

-2.69 × 103

function approximation quickly converges without wasting the numbers of experimental designs. Assume that no prior information about the function is available; the quality function J is formulated to minimize the difference between the actual function f(t) and approximated function ˆf(t) during the operating time

∫03600(f(t) - ˆf(t))2 dt

J)-

End

5. Illustrative Examples Two examples were selected to illustrate the operation and performance of the proposed method. For a simple explanation, a mathematical function that partially consists of the orthogonal functions is used first because it helps illustrate the search path of the optimal value. The other example is concerned with an optimal batch trajectory problem in a nonlinear fermentation chemical process. 5.1. Example 1: Fitting Mathematical Curve Problem. An artificial function with three terms of orthogonal functions is represented by

f(t) ) 4φ0(t) + 5φ1(t) + 3φ3(t)

(17)

defined for 0 e t e 3600 and -10 e f(t) e 15 (Figure 4). The function contains the Legendre polynomial functions φi, i ) 0, 1, 3, without φ2. The aim of this problem is to show how Taguchi’s method with the orthogonal

8.80 × 103

(18)

where ˆf(t) denotes the orthogonal function approximation (eq 4). Thus, the larger the value of J, the better the approximation of the original function. First, only the two-term orthogonal function is used, and a straight line is set as the reference approximated trajectory. Using eqs 11-16, the boundaries of the coefficients can be determined (Table 2). The Taguchi table is reproduced at the first run (Table 3) to evaluate the effects of the different experimental coefficients. Only four combinations of experimental coefficients are needed for these two coefficients. The columns labeled c0 and c1 represent the two experimental coefficients. The J value of each experimental coefficient is shown in the column labeled “value”. The row labeled “effect” can be used to evaluate the results of the experiments. Table 3 indicates that c0 and c1 should increase to enhance the average J value and c1 should be the first selection to define the range of c1 because the effect of both is c1 > c0. In the second run, the new Taguchi table (Table 4) is conducted using the search direction in the previous run. Only three trials are performed because trial 1 in


Figure 5. Locations of the coefficients ci, i ) 0-3, and the search regions for different runs in example 1.

Table 4 is the same as trial 4 in Table 3. From the effect values in Tables 3 and 4, c0 lies between 0 and 8.0296, and c1 needs to be further increased. The same procedures are kept until the quality performance is not significantly improved (R ) 0.05) at run 4. Then, a new orthogonal term φ2(t) is added to increase the search space. Because the φ2(t) term makes no contribution, a new term φ3(t) should be added after three more runs. The optimal point search based on the proposed experimental design is depicted in Figure 5. Coefficients ci, i ) 1, 2, 3, 4, respectively, lay out the location of the experiments. The corresponding J value is also included.

The search region is based on the proposed region elimination method. For example, the optimal region of coefficient c0 in run 1 includes either 0 e c0 e 5.1349 or 5.1349 e c0 because J(c0)5.1349) > J(c0)0). Thus, the levels of c0 in run 2 range from 0 to a value larger than 5.1349. The same procedures are repeated . A total of 67 experiments are conducted to determine the optimum trajectory. The results of the 1st, 2nd, 4th, 8th, 10th, and 12th runs are shown in Figure 5. It is obvious that the search region with shaded area of each coefficient ci, i ) 0, 1, 2, is systematically reduced and the contribution of the fourth term c3 is not significant from

2234

Ind. Eng. Chem. Res., Vol. 41, No. 9, 2002 Table 5. First Run of Taguchi Experimental Design for Example 2 c0 level 2 35

trail

value

1 2 3 4

1.3233 1.3233 1.2850 1.285 1.4348 1.4657

avg

1.3772 1.30415 1.45025

effect

Figure 6. Optimal approximated trajectories for different runs along with the actual trajectory function.

level 1 31.64

c1 level 1 -2.0412

c2 level 2 0

level 1 0

1.3233

1.3233 1.285

1.4348 1.4657

1.285 1.4348

1.4348 1.4657

1.4657

1.37905 1.37535

1.3945 1.3599

-0.0037

0.1461

level 2 1.5359

-0.0346

around the optimal point. The optimal region elimination method helps us position the proper experimental location to fit the original function. The search directions and iterates for these four coefficients are separately shown in two subplots (Figure 7a and b), each of which represents the contour of the function with two coefficients fixed at the optimal values. 5.2. Example 2: Optimal Temperature Trajectory for Penicillin Fermentation. In this case study, the proposed method is applied to solve the two-point boundary-value problem of a batch penicillin fermentation process. The differential equations describing the state of the system in batch penicillin fermentation21 are

cell mass production dy1 b1 ) b1y1 - y12 dt b2

(19)

penicillin synthesis dy2 ) b3y1 dt

(20)

with

b1 ) w1

b2 ) w4

b3 ) w5

Figure 7. Execution of the coefficient search: (a) c0 vs c1, (b) c2 vs c3.

run 10 to run 12. The stopping criterion (J e 50) is satisfied after run 12. The optimal approximated functions at the various runs are plotted with the actual function (Figure 6). Note that there is not much experimental information in the first few runs. After eight runs, the locations of the coefficients are within the optimal region, and then the approximated function slowly approaches the original function with the new experimental data added (Figure 6). The coefficients locate more and more points

[ [ [

1.0 - w2(θ - w3)2

] ] ]

1.0 - w2(25 - w3)2 1.0 - w2(θ - w3)2

1.0 - w2(25 - w3)2 1.0 - w2(θ - w6)2

1.0 - w2(25 - w6)2

where w1 ) 13.1, w2 ) 0.005, w3 ) 30 °C, w4 ) 0.94, w5 ) 1.71, and w6 ) 20 °C. y1 and y2 represent the dimensionless concentrations of cell mass and penicillin, respectively; t is the dimensionless time, 0 e t e 1; and θ is the temperature in degrees Celsius. The parameters bi, i ) 1-3, are functions of temperature, bi g 0. The parameter-temperature functions have shapes typical of those encountered in microbial or enzyme-catalyzed reactions. The simulation conditions and relevant parameters are the same as those of Constantinides et al.22 The initial concentrations are specified as

y1(0) ) 0.03 and y2(0) ) 0.0 The performance index to be maximized is the concen-


Figure 8. Locations of the coefficients ci, i ) 0-4, and the search regions for different runs in example 2.

tration of penicillin at the end of the fermentation, i.e., t)1

J ) max y2(t)1) θ(t)

(21)

Because the Taguchi two-level array of L4 allows up to three factors, a three-term orthogonal function is used first to reduce the number of experiments and the number of runs. Initially, a straight-line shape is chosen for the temperature trajectory for the setting of the

2236


) 1.4780), the results with the actual solution (ymax 2 from the proposed method are satisfactory. 6. Conclusion

Figure 9. (a) Optimal trajectory and (b) maximum concentration of penicillin for different runs in example 2.

boundary values of each coefficient. The Taguchi orthogonal array and the corresponding data collection plane are illustrated in Table 5. Because the trajectory depends only on the coefficients, only the coefficients ci need to be optimized for J to be maximized. Judging from the main effects of each operating coefficient, we conclude the the setting c0(level 2), c1(level 1), and c2(level 1) would give the highest value of J. As in example 1, the proposed design procedure for this case is depicted in Figure 8. A new coefficient corresponding to the fourth term of the Legendre function is added after the third batch run. Then, the same procedure is repeated. The fifth term of the Legendre function is added after the sixth batch run, but the objective function is not significantly improved. The experiment is then stopped, and a near-optimum is found. A total of 50 experiments are required to reach an optimum. Figure 9 shows the resulting optimal trajectory and the optimal yield at each run. In Figure 9a, the optimal trajectories of runs 4-6 overlap; the optimal trajectories of runs 7 and 8 overlap as well. The shape of the optimal trajectory is changed after the new coefficient is added to increase the search space. As can be quite clearly seen, the optimal policy would produce remarkable improvements in the yield of the desired product in the fourth run. The maximum principle has been applied to this case to determine the optimal temperature trajectory.22 However, in that case, the given complicated models involving the chemical reaction are required. Compared

A wide variation of complicated behavior is often expected in batch process manufacturing. In practice, it is impossible to solve the optimization set-point problem corresponding to such a process using a reliable process model. In addition, the solution of the resulting optimization problem using existing techniques is computationally expensive. Instead, in this work, a statistical approach, Taguchi’s method, and a function approximation tool, the orthogonal function approximation, are applied to evaluate the effect of various trajectories and to discover the optimal operating batch trajectory in a batch process. Taguchi’s method of experimental design is beneficial in pinpointing the factors of the design process. It is most helpful in controlling the quality once the levels of the factors have been set. The Taguchi approach, however, is more appropriate for dealing with static process relationships rather than dynamic batch processes among the design factors. It cannot establish the appropriate levels of batch trajectories to determine the optimal operating conditions. In this research, the orthonormal function approximation is used to eliminate the problem resulting from the dynamic set-point trajectory. It can turn implicit information into key coefficients that contain the necessary part of the operating conditions for each variable in each batch. Taguchi’s method for designing the coefficients can then be used to develop the optimal trajectory in the batch process. Additional studies for searching for the optimal condition are also proposed. Finally, it is to be noted that all of the results presented in this paper are based on the simulation model. Although the numbers of experiments in the first and second examples here are 67 and 50, respectively, these experiments are designed without any prior knowledge of the system. The required number of experiments can be significantly reduced if the initial trajectory is explored based on historical operations. The aims of developing the proposed method are to construct the framework for batch process design systems and to finally perform the developed method on the actual batch plant data in the future. Further batch process design studies will be our topic of next research, including integration of the proposed method with the soft-computing model to overcome the limitations of multiple optima of batch trajectories, dynamic path constraint problems, and design of the optimal duration of a batch. Acknowledgment The authors thank the reviewers for their helpful suggestions for revising this paper. Also, this work is supported by the National Science Council, R.O.C. Nomenclature ci ) projection of the trajectory onto the basis function φi cU i ) upper bound of the coefficient ci cLi ) lower bound of the coefficient ci /cU ) selected upper bound of the coefficient c i i U /ci ) selected lower bound of the coefficient ci J ) quality objective function K ) number of measurement points in each trajectory

Ind. Eng. Chem. Res., Vol. 41, No. 9, 2002 2237 L ) lower bound of the manipulated trajectory n ) number of repeated experiments N ) number of orthogonal basis functions rmax ) maximum number of batch runs R ) threshold ratio for quality improvement U ) upper bound of the manipulated trajectory ys ) target response Greek Symbol φi ) orthogonal function

Literature Cited (1) Soroush, M.; Kravaris, C., Optimal Design and Operation of Batch Reactors. 2. A Case Study. Ind. Eng. Chem. Res. 1993, 32, 882. (2) Thomes, I. M.; Kiparissides, C. Computation of NearOptimal Temperature and Initiator Polices for Batch Polymerization Reactor. Can. J. Chem. Eng. 1984, 62, 284. (3) Tsang, T. T. C.; Clarke, D. W. Generalized Predictive Control with Input Constraints. IEEE Proc. D. 1988, 135, 451. (4) Luus, R.; Rosen, O. Application of Dynamic Programming to Final State Constrained Optimal Control Problems. Ind. Eng. Chem. Res. 1991, 30, 1525. (5) Luus, R. Optimal Control by Dynamic Programming Using Systematic Reduction in Grid Size. Int. J. Control 1990, 51, 995. (6) Luus, R. Optimal control by direct search on feedback gain matrix. Chem. Eng. Sci. 1974, 29, 1013. (7) Cuthrell, J. E.; Biegler, L. T. Simultaneous Optimization Methods for Batch Reactor Control Profiles. Compt. Chem. Eng. 1989, 33, 8, 1257. (8) Eaton, J. W.; Rawlings, J. B. Feedback Control of Chemical Processes Using On-Line Optimization Techniques. Compt. Chem. Eng. 1990, 14, 469. (9) Zafriou, E.; Zhu, J. M. Optimal Feed-Rate Profile Determination for Fed-Batch Fermentation in the Presence of ModelPlant Mismatch. In Proceedings of the American Control Conference; IEEE: Piscataway, NJ, 1989; p 2006. (10) Modak, J. M.; Lim, H. C. Feedback Optimization of FedBatch Fermentation. Biotechnol. Bioeng. 1987, 30, 528.

(11) Wilson, D. J.; Agarwal, M.; Rippin, D. W. T. Experiences Implementing the Extended Kalman Filter on an Industrial Batch Reactor. Compt. Chem. Eng. 1998, 22, 1653. (12) Duppen, D.; Bonvin, D.; Rippin, D. W. T. Implementation of Adaptive Optimal Operation for a Semi-batch Reaction System. Compt. Chem. Eng. 1997, 22, 185. (13) Robert, H. L.; Joseph, E. M. Designing for Quality, An Introduction to the Best of Taguchi and Western Methods of Statistical Experimental Design; ASQC Quality Press: Milwaukee, WI, 1990. (14) Fowlkes, W. Y.; Creveling, C. M. Engineering Methods for Robust Product Design; Addison-Wesley: Reading, MA, 1995. (15) Kumar, A.; Motwani, J.; Otero, L. An Application of Taguchi’s Robust Experimental Design Technique to Improve Service Performance. Int. J. Qual. Reliab. Manage. 1996, 13, 85. (16) Mujtaba, I. M.; Macchietto, S. Efficient Optimization of Batch Distillation with Chemical Reaction Using Polynomial Curve Fitting Techniques. Ind. Eng. Chem. Res. 1997, 26, 2287. (17) Ray, W. H. Advanced Process Control; McGraw-Hill: New York, 1982. (18) Phake, M. S. Quality Engineering Using Robust Design; Prentice Hall: Upper Saddle River, NJ, 1991. (19) Diamond, W. J. Practical Experiment Design for Engineers and Scientists; John Wiley & Sons: New York, 1981. (20) Sansone, G. Orthogonal Functions; Dover: New York, 1991. (21) Constantinides, A.; Spencer, J. L.; Gaden, E. L., Jr. Optimization of Batch Fermentation Processes. I. Development of Mathematical Models for Batch Penicillin Fermentations. Biotechnol. Bioeng. 1970, 12, 803. (22) Constantinides, A.; Spencer, J. L.; Gaden, E. L., Jr. Optimization of Batch Fermentation Processes, II. Optimal Temperature Profiles for Batch Penicillin Fermentations. Biotechnol. Bioeng. 1970, 12, 1081.

Received for review May 22, 2001 Revised manuscript received September 24, 2001 Accepted February 11, 2002 IE010457X

Using Taguchi's Method and Orthogonal Function Approximation to

Recommend Documents