Semi-batch process optimization under uncertainty: Theory and experiments

Pergamon

Computers chem. Engng Vol.22, No. 1-2,pp. 201-213. 1998 Copyright© 1997ElsevierScienceLtd Printedin GreatBritain.All tightsreserved PII: S0098-1354(96)00359-6 0098-1354/97$17.00+0.00

Semi-batch process optimization under uncertainty: Theory and experiments Peter Terwiesch*, Dag Ravemark, Benedikt Schenker and David W. T. Rippinl" Systems Engineering Group, LTC, ETH Zurich, CH-8092 Zurich, Switzerland

(Received 21 September 1994; revised 11 September 1995) Abstract

Batch and semi-batch processes provide a flexible means of producing high value-added products in the chemical, biotechnical, and pharmaceutical industries. Unlike continuous processes, batch processes are inherently transient and typically also nonlinear, and a nonlinear dynamic end-point optimization problem needs to be solved in order to determine the optimal operating strategy. Inter- and intra-run variations and lack of better measurement information typically only allow us to build an imperfect model of the process, and the remaining uncertainties can be large. Nominal optimization techniques, requiring an exact process model, may thus not always be suitable for batch process optimization. This contribution suggests an alternative approach in which modeling and identification uncertainty are explicitly accounted for during the process optimization. As opposed to a deterministic quantity that is a function of only the nominal model, a probabilistic measure of success is optimized, leading to robustness of the desired objective to uncertainties and variations. Both simulation and experimental results are given to demonstrate the idea and the application of the proposed approach and to highlight the benefits that can be expected from its industrial application. © 1997 Elsevier Science Ltd Keywords: Semi-batch; End point; Optimization; Practical; Uncertainty

1. Introduction Due to their operational flexibility, batch and semibatch processes are frequently used in the production of high-value added products in the chemical, biotechnical, and pharmaceutical industries. Unlike continuous processes they are inherently transient with typically nonlinear dynamics, thus requiring solution of a dynamic end-point optimization problem for their optimal operation. This problem is traditionally stated and solved assuming an exact deterministic process model is available Denbigh (1958). However, because of the innovative nature and the comparatively small quantity of many (semi-)batch products, typically only few modeling runs are conducted. Inter- and intra-run variations are frequent and end-point reproducibility, even in the laboratory, is often not better than 5% Terwiesch et al. (1994). Optimization with only the nominal process model seems inappropriate under these circumstances. *Corresponding author. New address: ABB Corporate Research, Information Technology Dept., 69115 Heidelberg, Germany. t Deceased 26 June 1994.

The approach taken here is to optimize a probabilistic measure of success that accounts for possible uncertainties or variations found in the identification step. Such probabilistic measures are often either of the expected value type, which is adequate when a good average outcome is desired, or of the risk type, e.g., when safety issues need to be considered or when the risk of making unacceptable product is non-negligible. Since in many current industrial applications only insufficient online measurement information is available, this contribution focuses on open-loop policy optimization. Obviously, provided on-line measurement information and sufficient computing resources are available, the presented technique can also be used to continuously adjust the operational variables during the whole of batch processing time. A simple simulation example with single parameter uncertainty is used to demonstrate the basic ideas of the suggested approach. As a second example, a practical case of a semi-batch process optimization under uncertainty is studied experimentally. The authors have deliberately chosen a non-proprietary problem, allowing a transparent discussion of all necessary details from the modeling to the production runs. The aim of the investigated process is to produce a red dye, while

201

202

P. TERWIESCH e l

keeping both an undesired side-product (blue dye) and the amount of unreacted raw material below specified bounds in order to guarantee the desired properties of the final product. A third example examines the effects of structural mismatch between the assumed model and a simulated "true" process. A structurally simple model is fitted to the more complex "true" process dynamics and two parameters and their respective uncertainties are identified. Subsequent optimization considering the probabilistic uncertainty description is found to yield an operating policy that is more robust also with respect to deterministic variations in the true parameters, i.e., to structural mismatch. In all cases policy parameterization is kept simple, reflecting the fact that industrial implementation typically requires a minimum of policy sophistication. While more complex parameterizations, e.g., dynamic profiles in the manipulated variables, can also be treated in a similar manner Terwiesch and Agarwal (1995), Ruppen et al. (1994), this contribution assumes that the feeding policy is parameterized with only 3 degrees of freedom: a constant feed level, u~, the duration of the feed addition phase, tj, and the duration of a pure batch phase, t2. This class of feed addition profiles is probably most frequently used in industry and allows the optimization problem to be kept simple for our demonstration purposes.

1.1. Objectives Expected-value objectives express the desire to achieve the best average result, e.g., maximization of expected yield, or to reduce variations in final product composition, corresponding to, e.g., minimization of end-point variance. Possible risk objectives, on the other hand, include maximization of the probability for making acceptable product, maximization of that profit that one can be 90% sure of achieving, or determining the most rigorous product specifications that can be met in 95% of all batches. While such robustness objectives are commonly found in practice, they remain unaddressed by nominal optimization techniques. A more detailed discussion of potential objectives is given in Section 2.1. 1.2. Constraints In industrial batch process operations, three types of constraints can be distinguished: Product Specification: definition of final product quality, e.g., upper bounds for by-products or unreacted feed, lower bounds for concentration of desired product. Safety Regulations: definition of what is considered safe operation of the process, e.g., an upper bound for accumulation of hazardous intermediates that could lead to a thermal runaway. Manipulated Variable Bounds: range that the manipulated variables can be chosen from, such as upper and lower bounds for jacket temperature due to heating and cooling limitations or upper bound

al.

for feed addition rate limited by pump throughput. The definition of product specifications itself is often part of the objective. Their violation on parts of the parameter space can be tolerated as long as the overall result is satisfactory. Safety and manipulated variable constraints may not be violated anywhere.

1.3. Structure of this contribution The following Section 2 discusses a number of objectives with uncertainty and lists an algorithm that can be used for the numerical optimization with probabilistic uncertainty and constraints. Section 3 characterizes the investigated process examples, their models, uncertainties, and optimization objectives. Numerical optimization and simulation results for all three processes are given in Section 4, while Section 5 is reserved to experimental results. Section 6 summarizes and concludes this paper.

2. Mathematical approach 2.1. Objectives with uncertainty In a probabilistic uncertainty notation, typical objective functions can be written in the form max u ~oJ(u, O)p(O)dO

(1)

where u, 0, (9, J, and p represent operating policy, parameter vector, parameter space, profit associated with a particular outcome, and probability, respectively. Several different types of objectives are common:

Threshold Risk: The profit function J is of the binary type, i.e., J= 1 for admissible product and J = 0 for off-spec product. The objective then is to maximize the probability of making acceptable product. a,%-Risk Threshold: The probability or relative frequency a of acceptable batches is fixed while some function of the operating policy is optimized. This can be formulated as max Jl(u,s,O) u.s s't.~0~o J z(u,s, O)p(O)d0>----ot where J2 is of the binary type, J~ is some function of the operating policy, and s defines the product specifications. For example, the product specification itself can be subject to optimization and the question could be "How tight can we make the specifications while still making acceptable product with a probability of at least 90%?". Or, Jj can be chosen as a quantitative measure, e.g., the amount of product to be maximized, while the J~ constraint guarantees quality. Expected Value: Find an operating policy that maximizes the expected profit of operation. Exampies may include penalties for violation of soft constraints, for example the cost of reworking or disposing of unacceptable batches, or simply J can

203

Semi-batch process optimization under uncertainty be equal to concentration or total amount of produced product. Minimum Variance: Find an operating policy that reproduces a target product specification as tightly as possible despite variations and uncertainty in the system parameters. 2.2. Numerical algorithm We use an algorithm similar to the one suggested in Straub and Grossmann (1993) in the context of plant design optimization. It is illustrated here for the case of two uncertain parameters 0, and 0,: 1. Determine the maximum interval 8, = [ Of,13:]where a 0, exists so that the constraints g, are met. 2. Collocate 8, with n, points f?, and collocation weights wi. 3. For each of these points t?,, find the largest interval @i= [ 0?, e,“i] where the constraints g, are met. 4. Collocate each @i with n, points f@ and collocation weights w(. 5. Compute the probability of a successful batch run given a strategy u as

J(e’,&) wj2p(@&).

(2)

This integration of the probabilistic measure of success is performed within an outer loop for the optimization of the operating policy which is a NelderMead simplex algorithm Nelder and Mead (1965). J(e’,,#) represents a weighting of the probability, e.g., profit. 3. Investigated processes

1 (3) where the rate constants have the nominal value k, =OS, k,=OS, and where the symbols [A], [B], V, and u denote concentration of A and B, current relative volume, and feed rate, respectively. Initially, the reactor contains [A](O)=O.2 moles per liter of A, no B ([B](O)=O) and is filled to 50% (V(O)=O.5). We assume here that the rate constant of the side reaction, k2, is not as well determined as that of the desired reaction. Its uncertainty can be described by a normal distribution around the nominal value pk2=0.5 with standard deviation a,,=O.OS Such an extent of uncertainty also represents a typical magnitude for many practical problems. 3.2. Nominal optimum Using traditional nominal optimization over the described class of staircase feed profiles, one obtains u, =O.Ol, t, =50 and tz= 120, i.e., maximum feed addition until the reactor is filled, then wait as long as possible to let the remaining B react. Even though the example was originally used to demonstrate the benefits of more sophisticated feeding profiles, it reaches more than 97% of what could be reached with the theoretically optimal singular feeding policy (Young (1980), Terwiesch (1994)). The final product is acceptable when the following specifications are met: l

This section introduces three semi-batch optimization problems that are used to demonstrate and study the proposed approach. 3.1. A simple simulation problem Consider a chemical reaction system kI

A+B+C B+B2D

conducted in an isothermal semi-batch reactor. While simple, this example characterizes an important class of industrial problems. Since it is also a common benchmark for optimal control approaches, we will assume the same numerical values as given by Young (1980). The aim is, through addition of reactant B, to convert as much as possible of reactant A to the desired product, C, in a specified time t,= 120. It would be suboptimal to add all B initially as the second-order side-reaction yielding the undesired species D will be favoured at high concentrations of B. To keep the effect of this undesired reaction small, the system is operated in semi-batch mode and B is added in a feed stream with concentration br,=0.2. The rate equations for this system can be expressed as

l

l

final concentration of unreacted feed B no greater than 0.0035 moles per liter; final concentration of undesired side product D no greater than 0.018 moles per liter; final concentration of desired product C at least 0.058 moles per liter.

From a practical point of view, since the nominally optimal policy satisfies these constraints even with a certain margin, one could be tempted to run the reaction with this policy even with the given model uncertainty. However, as Monte Carlo simulations show, the nominal policy yields acceptable product in only 35% of all runs. This value is undesirably small, so that normally more modeling runs would need to be conducted to reduce the model uncertainty, which would be expensive with respect to various resources (manpower, laboratory cost, time to market, material). 3.3. A real process As mentioned in the introduction, the aim of this practical process is to make red dye. However, since the raw material is produced in a previous reaction step, comparatively large variations in its composition may occur. These changes cannot be detected with the existing equipment and make it difficult to determine an optimal open-loop feeding policy. Chemical knowledge

P. TERWlESCHet al.

204

n~ "o~

so~"

m''"~ + I-I+

"o~

~m~

__f3÷

+ I.-I+ N2*

-o~s

-o~s

sol

suggests that the reactions of interest occur in the sequence (see top of page) To simplify the subsequent discussion we introduce the following abbreviations. A represents the raw material, 4,5-Dihydroxynaphthalene-2,7-disulfonic acid (chromotropic acid), that acts as a coupling reagent. B stands for the diazotized aniline (diazonium salt), C for the mono- and D for the di-substituted (bisazo) product. The solution is buffered to a pH of 10, so the concentration of H ÷ remains constant during the reaction and can be omitted in the further discussion. The mechanism can then be expressed through the short notation kl

A + B "* C C+B"

k2

D.

These reactions take place in a single well-mixed homogeneous liquid phase. If the initial amount of A was constant and/or known, a one-parametric trial-and-error search could be conducted and an operating strategy could be found without a kinetic model. As this is not the case, heuristic determination of a suitable operating strategy would require a large number of trial experiments. The following sections of this paper are concerned with the use of mathematical modeling and optimization methods in order to find a suitable operating strategy with reduced experimental effort. Based on the chemical knowledge about the investigated system (Szele and Zollinger (1983) and Zollinger (1991)), the following kinetic model structure is hypothesized:

d[A]ldt - k, [A] [B] (d[Blldt~ [ - k,[AltBI - k=[C][B]~ d[C]ldt)=t k,[A][B]-k,[C][B] )+( d[D]ldt /

dVIdt

~

k=[C][B]

0

"

- [A]/V

b f ~ - [B]

_[C]Iv )u, -

[D]/V

1

where u, bf~, k~, and k2 denote feed flow rate, feed concentration and two rate constants, respectively. [A], [B], [C], and [D] denote molar concentrations of the species A, B, C, and D, The hypothesis that the rates for

sth'

both reactions depend linearly on each of the involved reactants is based on Szele and Zollinger (1983) and can be confirmed experimentally. In the next step, an attempt is made to identify the values of the rate constants k~ and ks using UV-VIS spectroscopy. No reference spectra are available to the authors for C and D, and time considerations forbid the synthesis of pure components for calibration purposes. Consequently, both spectra and rate constants need to be estimated using nonlinear regression on four experimental modeling runs. bf~d=0.39 mol/1 is a known constant. Based on the modeling runs the spectra can be estimated reasonably well. Figure 1 shows the obtained component spectra for the diluted samples and an optical length of 1 cm. Estimates for the two rate constants are overly correlated to be useful, as the available data does not provide sufficient information to estimate k~ and k2 separately. The reason is that although a numerical minimum for the estimation objective function is found at kc5.25 mol/1 s, k2=0.645mol/l s, only a lower bound but no upper bound on both rated constants can be determined with sufficient statistical significance. Consequently, the model is reduced by one parameter and only the ratio O=k2/kj is estimated. The lower bound on both rate constants is directly anticipated in the planning of the operating policy. Comparison of modeling runs with different stirrer speeds and feed rates yields the insight that the reaction is not limited by mixing. The remaining uncertainty of 0 is found to be adequately represented by a normal distribution with mean /~0=0.1228 and standard deviation tr0=0.01. The raw material A is desired to be transformed, through addition of diazonium salt B, into the red dye (monoazo) product C. The extent of the second reaction to D must be kept small in order not to spoil the desired red monoazo product C with blue bisazo product D. To avoid ageing and separating problems, the amount of final A allowed to be contained in the product is limited to [A]r-~0.01 mol/1. Traces of blue dye must stay below [D]r.~0.0055 mol/I in order not to spoil the visual impression of the red dye. Undetectable variations in initial raw material concentration [A]0 are modeled by a normal distribution with mean/hA~0--0.05 mol/l, as is nominally required by the recipe. The standard deviation is o-iA~--0.005mol/l and anticipates various possible sources of variation in

Semi-batch process optimization under uncertainty

205

70 60 ID U

50 4O 30 20

XO

5~o

5~o

660

wavdength ,X/n,n

Fig. I. Absorptionspectra for C (solid) and D (dashed). upstream operations. Since no opposite information is available, no correlation between 0 and [A]o is assumed. 3.4. Mismatch problem

The purpose of the third example is to study the influence of mismatch between the assumed model structure and the true process dynamics. From a theoretical point of view, it is not clear to what extent structural model-plant mismatch can adequately be described as parametric uncertainty. In practice, on the other hand, even in the most simple batch kinetics studies, verification of residuals between model prediction and observation typically shows clear trends rather than white noise, indicating the presence of structural mismatch rather than ideal measurement noise. Consider the following reaction mechanism: A +/~ - - ~ . -1_ ~ C

(4)

k~.

k,u

B + B ~ D (order 2.5).

This mechanism is subsequently simulated for the true kinetics. The reaction is carried out in an isothermal semi-batch reactor and the simulated true dynamic model is - k~,[A][B] + k~[/]

d[A ]/dt

d[/]/dt /

k,,[A][B]-k2,[/] 0

\

dV/dt

,]

- [A]/V

-

analyzed for the four components A, B, C, and D. The short-lived intermediate I cannot be observed. A mechanism is postulated that agrees relatively well with the observations, i.e., produces a good fit in A and C concentration errors. kl

A+B-* C kl

B+B--~ D.

This mechanism is the same as in the first example and is one that is frequently identified in fed-batch chemical reactions; in our case, it fits the observations well. As in most practical examples, only part of the residuals appear to be due to noise and the remaining part contains a slight trend (not shown). Using this nominal model rather than the true kinetic model, standard regression techniques (Bard (1974) and Bates and Watts 0988)) yield the following estimates and uncertainties*: k, =0.5051±0.035, k2=0.5088±0.094, and correlation coefficient p=0.95. Since this yields the same nominal model as in the first example section, we keep the quality constraints as suggested there. From a practical point of view, the obtained model structure and parameters are simple and explain the observations. There is no need for further refinement or increased model complexity. Even the violation of the mass balance, where the quantities of A and C do not exactly add up to a constant over time, is very small, since the intermediate I is only short-lived and thus exists in concentrations so small that the difference is easily masked by measurement noise.

[II/V

1

At time t=0, only A is present with concentration [A]o=0.2. The true values of the rate constants are kl,= 1.1355, k~=5.0, ka,=4.2539, and k4,=3.5880. Two modeling runs are simulated. During the run, A and C are measured at 10 points. The final product is

4. Numerical results Numerical optimization results for the three examples are given in this section. * Given as expectedvalue ± one standarddeviation.

P. TERWIF_~CH et al.

206

4.1. Simple problem Three different optimization attitudes are demonstrated for the first example.

4.1.1. Risk minimization First, a minimization of the risk of making unacceptable product is conducted. This optimization problem can be formulated as:

max p([B]r~[B]/.moxA[C]/~[C]/.m~,A[D]:[D]:,,,~)

(ul,tl.t2)

s.t.

(6)

O--
u~tl~0.5 0~t~<--t2-< 120, where [A] denotes concentration of A and indexfdenotes value at final time. For the specifications defined above, [B]/,.~=0.0035, [C]/.m~,=0.058, and [D]/..~=0.018. The three other constraints express a limitation of the feed stream, the maximum relative volume and the timing requirements. Note that with the given objective there is

no sense in reducing the batch time to less than its upper bound, so that t2= 120 is always optimal. Figure 2 compares the obtained risk-minimizing feed profile (solid) with the nominal feed profile (dashed). Figure 3 shows the final concentration of C as a function of the uncertain parameter ke for both policies together with the constraint [C]:,ml,. Figure 4 does the same for the other two constraints. The acceptable region of/q, i.e., the interval for which a given policy produces a final product that meets the specifications, is shown in Fig. 5 together with the probability density (PDF) for p(k2). One can clearly see (Figs 3 and 4) that the conservative risk-minimizing strategy (solid) produces a slightly smaller concentration [C] for the nominal k2, but extends the acceptable product region substantially (98.21%, solid, compared to 34.70%, dashed). Using a Nelder-Mead simplex algorithm, a fast integrator (with relative error <--10-s), and seven collocation points for each scalar normal distribution, the problem is solved within 2 minutes on a DEC /~VAX.

0.012

0.01

-

0.008

0.006

0.004

0.002

0

20

0

40

60

80

100

t

Fig. 2. Comparison of feed addition profiles: nominal (. . . . ) and risk-minimizing(

).

0.072, 0.07

0.068 .,., ~

O.OGG 0.064

0.062

0.~ 0.058 0.05~

=

o~s

Fig. 3. [CJ/as a function of k2 for risk-minimizing ( ~

.

o

o:4

o:s k~

o~6

0.7

) and nominal (.... ) feed addition profile.Constraint [C]/,=~,also

Semi-batch process optimization under uncertainty

207

0.025

0.02

f. f

/

/

/

/

/

/ ..°...°/ .... / °°--/ ..<-......~-:~"~"~ /

/

0.01~

O.O1

0 . ~

o'.3

0'.4

Fig. 4. [B]:and [D]: as a function of k2 for risk-minimizing (

4.1.2. Risk threshold optimization The optimization problem is now reformulated to illustrate the concept of a,%-risk threshold optimization. In the previous section, only product quality issues were addressed, while of course in realistic problems product quantity is also an important issue. One could take that into account either by solving a series of problems with variable lower bounds on product quantity or, as is demonstrated here, by optimizing the product quantity that can be reached given a specified risk. Mathematically, the optimization problem then becomes max V/

(7)

(/~ldl,/2)

s.t.pf[B]/~ [B]/,m~a[C]~- [C]f,minA [D]:-~[D]/,m.)~ a

o:s

o'.6

ks

o'.~

) and nominal (. . . . ) feed profile. Constraints [B]/.. . . .

[D]f.max"

tion optimum from the previous section, which is the upper limit on a for which a solution exists. It fills the reactor only to VI=0.8680. The nominal solution (filling the reactor with h=50, u~=0.01) results as a limiting case on the other side, as for 0t=34.70% the reactor is filled to VI= 1. The insight that can be obtained from this trade-off curve is that the risk of making unacceptable product can be reduced to about 5% requiring only a small reduction in final volume. When the objective is to maximize the actual quantity of C produced, while taking into account that part of the batches need to be rejected, the optimum is located at tr=0.93 on the trade-off curve and the reactor is filled to VI=0.9613.

4.1.3. Quality variance minimization

0--
u:~<--0.5 0--
When reproducibility of final product is the main goal, the minimum variance objective is used. Let us assume that reproducibility of [C]iis the desired goal. Obviously, if the solution space is not constrained by any further specifications, ut =0 is a trivial solution that produces no C with no variance. One obtains more meaningful results if the following, more realistic specifications, are made:

5 o..5 4

3.5

3 2.5 2 1.5 1 0.5

~,'.2

J

0.3

o14

ols

ole

i

0.7

k2 Fig. 5. Comparison of risk-minimizing (

) and nominal (. . . . ) region of acceptable product: PDF p(k2), also

o.a

P. TERWIESCHet al.

208 1 0.98 0.ge O.94

O.Q2 O.G

0.158 ! i !

0,89

0.8~

i

I

i

i

0.82

0.84

0.SQ

0.tt8

0 g

|

|

/

!

0.92

0.94

0.96

0.98

Fig. 6. Final volume V/( ) as a function of risk threshold a. (. . . . ) limiting case of risk minimization, O: maximum quantity of C considering rejected batches.

min ([C]/-E{ [C]/})2

(8)

(Ul,tl.t2)

s.t. p([B]:[B]/.mRxA[C]~[C]/.m,,A[D]~-[D]y.,~x)>-90%

0--
uncertain, the overall effect of product quality constraints is visualized through constant-risk contour lines (outer 10%, inner 95%). To the lower left, the amount of fed B is insufficient to produce enough C, so that the C:,m~.constraint is violated. Near the top right, the feed is stopped too late, so that too much B remains in the final product--a problem that could be dealt with by increasing the batch time (if not forbidden due to plant cycle time constraint). Also to the right, near where the hard feed rate constraint becomes active, the amount of final D violates the D/,m~, constraints as feed is added too quickly, thus temporarily increasing the concentration of B and favoring the polymerization to D.

4.2. Practical problem This section states numerical optimization results for the experimental example and suggestions for alternative operating policies for experimental comparison in Section 5. The requirements and assumptions translate into the following optimization problem:

350

23O

200

150 IO0

~0

8.0~1

0.032

0.034

0.056

0.038

0.06

0.062

[Clj

0.064

0.066

0.068

0.09

Fig. 7. Probability distribution for final product concentration [C]j, Left: minimum variance result, right: nominal solution.

Semi-batch process optimization under uncertainty

209

80 "70 60

,,..7

50 4O 30

20 10

0.005

0.01

0.015

0.02.

0.025

'U,

Fig. 8. Decision plane (u,fl) with lines of constant risk, feed and volume constraints. max p([A]/-~0.01 moi/IA[D]~0.0055 mol/l)

(u,J,,t2)

s.t.

(9)

0-
Although this optimization problem is already of minor complexity, it can be reduced further by heuristically apparent choices of u, and t2, as is discussed in the sequel. In order to avoid the risk of exothermic runaway, feed-limited semi-batch operation is chosen (Hoppe (1992)). An upper bound on the feed rate has to be obeyed such that the generated heat of reaction can always be removed isothermally. The information on the lower rate constant bound is taken into account by always leaving a pure batch phase t2=60 s between the last feed addition and reactor shutdown. This time is sufficiently long to guarantee vanishing of practically all B and is sufficiently short for not being relevant to the policy optimization problem. Given the kinetic mechanism, the fast reaction rates, the length of the final phase t2, and the upper feed constraint Um~= 10-4 I/S, potential operating strategies can differ only in the feed addition time h. This is convenient for the purpose of this comparative study, as it keeps the example simple and limits the necessary experimental verification effort. Note that similar to Meadows and Rawlings (1991) only the integral amount of added feed influences the final product quality, whereas the exact profile is not important. As adding feed at the upper feed rate bound is efficient in terms of reaction time, the feed rate is set to u=um,~ during

O
algorithm as ~=1230 s, the third strategy is based on equimolar addition of B, i.e., d = 1280 s, and the fourth policy waits for the maximum of [C], as would be optimal if unconstrained product quantity was maximized, i.e., ~ = 1330 s. The algorithm is used to evaluate the cumulative probabilities of making acceptable product. It yields the probabilities pj=55.2%, p2=63.6%, p3=50.8%, and P4=33.5% with n,=5, n2=7 collocation points per parameter.

4.3. Mismatch problem Based on the estimated parameter values, nominal optimization, risk minimization, and risk-threshold optimization are compared for the mismatch problem. The nominal policy is the same as in the previous section. For the other two, uncertainties in parameters k~ and k2 are anticipated in the policy optimization by obtaining the probability of making feasible product by integrating over their joint probability distribution. Figures 9-11 show the obtained regions of acceptable product in the (kj, ks) plane for each of the three strategies. The axes are chosen corresponding to the nominal value :1:3 standard deviations. A full ellipse, as is almost reached in the risk minimizing case (Fig. 11), would correspond to the case in which acceptable product is made for all parameter values within a region bounded by AX2<-9.0. This is a region that corresponds roughly to ±3tr in the scalar case, and is used as a collocation bound here. It is clearly advantageous over ±3tr collocation in each parameter as long as the true distribution can be approximated as normal since parameter estimates in batch processes are often highly correlated. In all three figures, the constraint that bounds the acceptable region to the lower left is the upper limit on [B]~ which is reached because ka is too small to consume B. It slightly declines to the right as an increase in kmcan compensate by consuming some of the B to C. In Fig. 9, an additional constraint appears at the upper left. This means that the nominal feeding policy, since it involves a high feedrate, favours the second-order side

210

P. TERWIESCHet al.

0.7

0.6

0.5

0.4 /,,, ,-****'**o...,o..., , , , . - * ' ' ' ' * * " t.,° . . . . . --0.3 0.45

0.$

kl

0.55

0.6

Fig. 9. Region of acceptable product for nominal operating strategy in the (k.k2) plane for ±3 standard deviations. constraints, ©: collocation points.

0.7

0.6

0.5

0.4.

0.3 0.45

0.5

kz

0.55

0.6

Fig. 10. Region of acceptable product for risk threshold operating strategy in the (k.k2) plane for ±3 standard deviations. constraints, C): collocation points.

0.7

0.6 t~ ,.kl 0.5

0.4

0.3

o.~s

ols

kz

o.~s

ole

Fig. 1I. Region of acceptable product for risk minimizing operating strategy in the (k,,k2)plane for ±3 standard deviations. constraints, O: collocation points.

Semi-batch process optimization under uncertainty Table 1. Summary of results

u t, v/ Qt F~ktlkl 8k21k2

~k~/k3 8kJk4

Nominal

a-Risk

Robust

0.01 50 1.000 41.2% +~ - 0.0497 +0.4236 -1 +oo - 0.3319 +0.3006 -0.2118

0.01 46.5 0.965 90% +oo - 0.1433 + 1.3349 -1 +~ - 0.894 + 1.2190 - 1

0.0083 45.5 0.878 99.2% +~ - 0.1393 + 1.2925 -1 +e¢ - 0.8715 + 1.2537 - 1

reaction to an extent at which the by-product D overshoots its upper quality limit. The larger the value of k~, the more B is also used to make 6', such that larger values of k2 can be tolerated. As it is sometimes held against probabilistic approaches that robustness to probabilistic uncertainty might only in pa~ be significant for performance on the true process, deterministic robustness margins are also determined. This is achieved by keeping three of the four parameters of the true model constant while finding a lower and an upper bound for the fourth. These bounds signify to what extent a parameter can vary without leading to unacceptable product. Table 1 summarizes the results obtained with the three different operating policies. First, a characterization in terms of the operating variables u and fl is given, leading to final volume V:. Then, the probabilistic measure of performance a, i.e., the probability of making acceptable product is given. In the final part of the table, a deterministic robustness measure, the relative changes in each parameter that are allowed while still making acceptable product, is stated. For example k~ variations are not bounded above since an increase in kl with any strategy produces more and improved product. The lower bound on k~ in the nominal case is given by a 5% decrease, while the other two operating policies are robust up to a 14% decrease. Similar results are found on all other parameters, where the policies obtained through probabilistic optimization based on a model with structural mismatch are considerably more robust to parameter changes than the nominal policy. The value of - 1 as a lower bound on any ~ / k i is not a consequence of a quality constraint but simply the bound of physical meaning for this parameter, i.e., ki=0. It is clear that these results do not have the quality of a mathematical proof. An operating strategy optimized for robustness to changes in two parameters with probabilistic uncertainty must not always yield robustness to deterministic variations in any four parameters of any structurally different model. However, one can adapt the point of view that in most batch processes the true process structure and its full complexity remain unknown, and that the identified model is simply a good local approximation to a more complex physical process. CACE 22-1/2-)1

211

Assuming conceptually that the model parameters depend on the true parameters through a smooth mapping, incorporation of probabilistic uncertainty information at least permits the avoidance of extreme sensitivity in a neighborhood of what appears to be the best approximate parameter values.

5. E x p e r i m e n t a l results

This section describes the procedure and equipment used in the experimental verification. Reaction data is given and analyzed and an interpretation is given. 5.1. Experimental setup

All experiments are conducted in a 1.5 1 jacketed stirred glass reactor at constant temperature T=5°C. The raw material, A, is purchased in the disodium dihydrate form from Aldrich Chemie, Buchs, Switzerland. The feed material, B, is prepared according to Fierz-David and Blangey (1952) using analytical-grade reagents from Fluka AG, Buchs, Switzerland. It is added continuously by means of a peristaltic pump. Feed consumption is monitored as a weight loss of the feed container on standard laboratory scales. UV-VIS spectroscopy together with the reference spectra generated during the modeling runs is used to analyze sample composition. Samples are prepared by diluting the original 2 ml sample in two steps. First the sample is diluted by a factor of 12.5 with a 2-normal buffer solution, which prevents further reaction. Then, in order to reach an area of satisfactory linearity and sensitivity of the UV-VIS instrument, a further dilution by a factor of 50 is performed with de-ionized water. 5.2. Experimental procedure

The aim of this experimental study is to evaluate the probability of making acceptable product with each of the investigated strategies. Consequently, experiments are designed where the usually uncertain and varying values of [A]0 are carefully controlled. The other uncertainty, O--/q/k~, is inaccessible to the experimenter, but the hypothesis of a constant but unknown O over all runs seems justified in view of the kinetic mechanism. In order to obtain accurate estimates of the probability of making acceptable product, 11 experimental runs are conducted, covering an interval of [A]oE[/hA~--20]A. ~0,/hAlo+2Oh~0]. The individual experiments are spaced such that the interval between any two experiments contains the same probability, so as to obtain reliable estimates. During each of the runs samples are taken at the time indicated by the appropriate strategy and are given one further minute before buffering and dilution, corresponding to the final batch phase t2=60 s. 5.3. Experimental data

Figure 12 summarizes the relevant experimental data. The vertical axis indicates the initial concentration [A]o

E TERWIESCHet al.

212 0.06

! !

'

ii

I

t t. . . . . . . . . . . . . . . . . . . . . . . . . . .

0.055

T

t /

(]i)®

@@

..

....... (:~!).~~.® .......................

0.05 ................... ~

\ \ \ \

®

::-.i.-i--i--i.. ::.. i.

,;) @ ®

@i)®®

0.045

X \ \ , 02

\

,

0.015

,

0.01

t~t~'

0.005~=''~'s

0.005

~--[A],

"~' ~aO'~'

'~'

[D], --~

Fig. 12. Experimental data, quality constraints, and theoretical probability distribution p([A]0). that was used in the 11 experiments. The horizontal axis from the center to the left shows the final concentration [A]/of unreacted raw material, intersected at [A]/=0.01 mol/i by the first product constraint. The horizontal axis from the center to the right displays the final concentration of blue dye, [D]/. The outcome of each of the I 1 experiments with each of the 4 operating policies is indicated by the circled policy numbers. For example, at the nominal initial concentration [A]0=0.05 mold, strategy no. 1 slightly violates the [A]/constraint, strategy no. 2 makes acceptable product, and strategies no. 3 and no. 4 violate the [D]/constraint. The experimental data is displayed as recorded. Apparently in the [A]0=0.0441 mol/l experiment samples no. 3 and no. 4 have been swapped accidentally during analysis. However, this is without consequence for the numerical integration of the probability of making acceptable product. The run with A]0=0.1M75 mol/l (fourth from below) appears to be an outliner that disadvantages strategy no. 2 with an early [D]/ constraint violation. However, since repetition of single runs with

undesirable results would be conceptually unsound, it is considered here as originally measured.

5.4. Interpretation Given the experimental data, the practical probability of making acceptable product with each of the investigated strategies is computed. This is done numerically by integrating the probability density p([A] 0) over those intervals for which the different strategies produce acceptable product. The bounds of these intervals are found using linear interpolation between the last experimental results within the acceptable region and the first one outside of it. The resulting cumulative experimental probabilities, the interval bounds, and the number of successful batches for each strategy are summarized in Table 2. Both in the theoretical policy optimization and the practical verification, policy no. 2 is most successful. It is no surprise that the predicted and the experimental probability differ considerably. The reason seems to be that the value of O=k2/k~ is apparently larger than the

Table 2. Comparison of the four operating strategies Strategy

tl/s

[A]o..,,/molll [A]o..Jmol/l

No. 1 No. 2 No. 3 No. 4

1180 1230 1280 1330

4.547.10-2 4.807.10-2 5.042.10-2 5.279.10-2

4.970.10 -2 5.200.10-2 5.369.10-2 5.634.10-2

p,~ 55.2% 63.6% 50.8% 33.5%

Pc\p,

P,xp2

Success

29.2% 31.1% 23.1% 18.8%

41.6% 57.6% 35.1% 12.6%

3 3 2 2

Semi-batch process optimization under uncertainty mean anticipated in the optimization step, as is suggested from the number of [D]f violations being larger than that of the [A]/ violations. Evidently the quality requirements in this contribution are hard and lead to a relative frequency of unacceptable batches that is significantly higher than in most industrial applications. However, a statement of the experimental optimization problem with the given specifications has the advantage that an experimental evaluation can yield a significant number of both acceptable and unacceptable batches for a small number of experiments. In realistic industrial situations only about one out of a hundred batches might turn out to be unacceptable, requiring a prohibitively large number of experimental runs for a verification study. Given a high frequency of unacceptable batches and fixed quality specifications, a next realistic step would be either to attempt to reduce the uncertainty in 0 through more modeling runs or regression analysis of past production runs or to reduce the standard deviation in raw material composition o'tAp. For example, assuming that one is able to reduce the uncertainty in raw material concentration o%l0 by 50%, the probability of making acceptable product with each of the investigated strategies would be as indicated in column Po~p2.

6. Conclusions Admittedly, the presented approach is limited to optimization for uncertainty in only a few model parameters, since computational complexity grows exponentially with the number of uncertain parameters. Also, with the given algorithm only comparatively simple policy parameterizations can be optimized. However, it is the authors' opinion that in most practical batch optimization problems parametric uncertainty can be lumped into a small number of parameters. Also, looking at current industrial operating strategies and equipment (Terwiesch et al. (1994), Ruppen (1994)), the limitation to simple strategies appears to be no severe drawback and the benefits that may result from application of the proposed technique to industrial problems can be significant. The first example has illustrated the proposed technique and has highlighted the importance of choosing the appropriate objective function and of using trade-off curves whenever no clear weighting of possibly conflicting objectives can be selected prior to optimization. The second example has demonstrated the sources and incorporation of two different types of uncertainty into the selection of an appropriate operating policy. While the given experimental results are insufficient to support optimality claims, they clearly indicate that the proposed optimized policy performs better than any of the three alternative heuristic strategies. The third examples supports the hope that statistical optimization approaches can also be useful when confronted with at least slight structural model-plant mismatch. A promising field for further research could be to find out whether or how infeasible path optimization methods

213

for simultaneous integration of the model, collocation of the probability density, and optimization of the operating profile could be used to improve computational efficiency.

Acknowledgements The authors want to thank Dr E Skrabal and Dr A. Klaus for their advice and practical support concerning the azo reaction. The first-named author acknowledges financial support from the Swiss National Science Foundation.

References Bard Y.(1974) Nonlinear Parameter Estimation. Academic Press, New York. Bates D. M. and Watts D. G. (1988)nonlinear Regression Analysis and its Applications. John Wiley and Sons. Denbigh, K. G. (1958) Optimum temperature sequences in reactors. Chem. Eng. Sci. 8, 125-132, New York. Fierz-David H.E. and Blangey L. (1952) Grundlegende Operationen der Farbenchemie. Springer-Verlag, Wien, 8th edition. Hoppe T. (1992) Use reaction calorimetry for safer process design. Chemical Engineering Progress, pp. 70-74, A1ChE, New York. Meadows E. S. and J. B. Rawlings (1991) Model identification and control of a semibatch chemical reactor. In Amer. Cont. Conf., pp. 249-255, ACC, San Francisco. Nelder, J. A. and Mead, R. (1965) A simplex method for function minimization. Computer Journal 7, 308-313. Ruppen D., Benthack C. and Bonvin D. (1994) Efficient computation of batch reactor control profiles under parametric uncertainties. In IFAC Adchem'94, pp. 79-84. Kyoto, Japan. Ruppen D. (1994) A theoretical and experimental contribution to the implementation of adaptive optimal operation for discontinuous chemical reactors. PhD thesis, Eidgen0ssische Technische Hochschule, Ztlrich. Dissertation ETH No. 10634. Stranb, D. A. and Grossmann, I. E. (1993) Design optimization of stochastic flexibility. Computers Chem. Engng. 17, 339-354. Szele, I. and Zollinger, H. (1983) Azo coupling reactions. Topics in Current Chemistry 112, 1-66. Terwiesch, P. and Agarwai, M. (1995) Robust input policies for batch reactors under parametric uncertainty. Chem. Eng. Comm. 131, 33-52. Terwiesch, P., Agarwal, M. and Rippin, D. W. T. (1994) Batch unit optimization with imperfect modeling: A survey. J. Proc. Control 4, 238-258. Terwiesch P. (1994) Dynamic Optimization of Batch Process Operations with Imperfect Modeling. PhD thesis, Eidgentssische Technische Hochschule, Ztlrich. Dissertation ETH No. 10857. Young M. J. (1980) Semi-batch reactions offer optimal yields. Processing, pp. 27-33, March. Zollinger H. (1991) Color Chemistry. VCH Verlagsgesellschaft mbH, Weinheim (FRG) and VCH Publishers, Inc., New York (U.S.A.), 2nd edition.