Development of on-line optimization-based control strategies for a starved-feed semibatch copolymerization reactor

ARTICLE IN PRESS Control Engineering Practice 18 (2010) 131–139

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Development of on-line optimization-based control strategies for a starved-feed semibatch copolymerization reactor L.F.T. Perea a, O.E. Okorafo b, R.A. Hutchinson b, M. Guay b, a b

Transmin Metallurgical Consultants, Lima, Peru Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada

a r t i c l e in fo

abstract

Article history: Received 22 January 2009 Accepted 13 October 2009 Available online 28 November 2009

A real-time nonlinear optimization scheme is implemented for the high temperature starved feed semibatch solution polymerization used to produce acrylic resins for automotive coatings. The system input and state trajectories are determined on-line to minimize a cost function defined in terms of deviation from target molecular weight and composition, and batch time. The performance of the proposed strategy is examined through numerical simulation for several study cases for semibatch butyl methacrylate and styrene (BMA/STY) free radical copolymerization. Results show that polymer quality can be assured in reduced batch time while achieving constraint satisfaction over the entire run. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Adaptive control Real-time optimization Model predictive control

1. Introduction Acrylic resins for automotive coatings are synthesized from a mixture of monomers selected from the methacrylate, acrylate and styrene families. In order to operate at reduced solvent levels as well as to control polymer composition and molecular weight (MW), the polymer is manufactured under higher temperature conditions. A ‘‘starved feed strategy’’, in which the monomer and initiator are fed at a fixed rate over a period of several hours, is utilized. Operating at these conditions of high instantaneous conversion (low monomer concentrations) ensures the production of polymer with relatively constant quality properties requiring minimal on-line measurement. In particular, the result of this feeding strategy is that copolymer composition and average molecular weight remain relatively uniform during the reaction. For batch and semibatch reactors the typical approach to improve performance has been divided into two parts. First, a recipe for the feed flowrate of one or more reaction components, typically based upon heuristics and experience, has been fixed. In a second step a controller is designed to implement the determined feed flowrate in the presence of disturbances (see e.g. Abel & Marquardt, 1998). Potential for optimization comes from the fact that both steps can be executed simultaneously through the use of novel control algorithms. Indeed, many applications in engineering can be formulated as optimization problems whose objective and constraints are often dictated by cost, safety, quality and environmental criteria. The optimization then proceeds through the manipulation of certain process  Corresponding author.

E-mail address: [email protected] (M. Guay). 0967-0661/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.conengprac.2009.10.003

variables such as flow rates, temperature, heat duty, etc. However, these process variables are not always adjustable and are often limited by physical bounds. Thus, the optimal operation of any process involves the determination of an extremum for the objective function by adjusting the process variables while satisfying the constraints on the process (Vemuri, 2004). The robust reactor operation comes at a price in terms of long batch time, and drift in both polymer molecular weight and composition still occurs in the early stages of the batch. Thus, there is an incentive for improving reactor operation and the efficiency of quality monitoring which can translate to higher production rates, improved product quality, safer operation, and, as the ultimate consequence, improved profits. The control strategies used for semibatch systems and a description of the particular approach used in this work is provided in Section 3 of this paper, with the specific copolymerization system tailored to the optimization strategy based on a dynamic model developed to capture the essence of the starved-feed free-radical copolymerization process (Section 2). The model of a butyl methacrylate and styrene (BMA/STY) system is used to test the on-line optimization algorithm, which is shown to reduce batch time while improving product uniformity. Several case studies are presented in Section 4 to illustrate the performance of the optimization approach.

2. System description The system under study in this work is a solution free-radical well-mixed polymerization reactor operating isothermally under starved-feed conditions. Initiator and monomers are continuously fed in the desired mass ratio to provide composition and molecular

ARTICLE IN PRESS 132

L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

weight control while maintaining low monomer concentration in the reactor. A representation of the process is given in Fig. 1. Solution free-radical BMA/STY copolymerizations have been performed at 138 3 C with monomer added at a constant rate over 6 h at different mass ratio (feed compositions); initiator was fed in at a constant mass ratio of 2.0 wt% relative to monomer feed. The experimental study of Li and Hutchinson (2006) showed that the free STY and BMA monomer levels in the semibatch reactor remain low throughout the course of these reactions such that the monomer and polymer compositions remain relatively constant. A detailed description of how polymer molecular weight and polymerization rate depend on reaction conditions and species concentrations from a defined set of kinetic mechanisms has been developed for semibatch high temperature BMA/STY copolymer-

BMA

STY

TBPA

Fig. 1. Schematic of semibatch reactor system. The flowrates of butyl methacrylate (BMA) and styrene (STY) monomers and tert-butyl peroxyacetate (TBPA) initiator can be independently manipulated.

ization (Li & Hutchinson, 2006). This full model can be used for product development and off-line optimization of the current starved-feed operating strategy. In this work a dynamic mathematical model has been built based on a reduced set of kinetic mechanisms embedded into overall material balances. Our aim is the development of a simplified model that captures the fundamental process behavior to design an efficient and robust system to control polymer properties while fully exploiting online optimization-based control techniques. Therefore, predictions on the intended controlled variables should be consistent and represent the dominant dynamics of the system. To this regard, the kinetics represented by the terminal model, which assumes that the reactivity of the polymer radical depends solely on the nature of its terminal monomer unit, was employed. Other mechanisms not considered in the reduced model include styrene thermal initiation and methacrylate depropagation. Under standard assumptions such as perfect mixing, constant physical properties, quasi-steady-state assumption (radical stationarity) and long chain hypothesis, mass balances on the initiator, solvent and monomers, and moment balances on the radical (live) and dead polymer chains yield the isothermal mathematical model for the semibatch system summarized as follows: dx1 A B ¼ u1 ðtÞ  x1 ðkpAA lo ðtÞ þ kpBA lo ðtÞÞ dt

ð1aÞ

dx2 B A ¼ u2 ðtÞ  x2 ðkpBB lo ðtÞ þkpAB lo ðtÞÞ dt

ð1bÞ

dx3 ¼ u3 ðtÞ  kd x3 ðtÞ dt

ð1cÞ

dx4 A B ¼ x1 ðkpAA lo ðtÞ þkpBA lo ðtÞÞ dt

ð1dÞ

dx5 B A ¼ x2 ðkpBB lo ðtÞ þ kpAB lo ðtÞÞ dt

ð1eÞ

dx6 ¼ VðtÞRm0 ðtÞ dt

ð1fÞ

Table 1 Kinetic rate coefficients and model parameters for butyl methacrylate (A)/styrene (B) copolymerization (Li & Hutchinson, 2006). Rate expression

Value at 138 3 C

Initiation

kd ðs1 Þ ¼ 6:78  1015 expð17714=TÞ f ¼ 0:515

1:32  103

Propagation

kpAA ðL mol

1

s1 Þ ¼ 3:80  106 expð2754:2=TÞ

kpBB ðL mol s1 Þ ¼ 4:266  107 expð3910=TÞ kp kp rA ¼ AA ¼ 0:42; rB ¼ BB ¼ 0:61 kpAB kpBA 1

Termination

ðL mol kcop t

1

4:69  103 3:16  103

s1 Þ ¼ 4:5  108

a ¼ 0:65; b ¼ 0:01; g ¼ 0:33 2.345

Transfer to

4 kpAA ksol tr AA ¼ 5:0  10

solvent

ksol tr BB

Transfer to

kmon tr AA ðL mol

1

s1 Þ ¼ 1:56  102 expð2621=TÞ

0.27

monomer

kmon tr BB ðL mol

1

s1 Þ ¼ 2:31  106 expð6377=TÞ

0.427

kmon tr ij Density

4

¼ 8:0  10

¼

kpBB

2.528

kpij kmon trjj kpjj

rBMA ðkg L1 Þ ¼ 0:91545  9:64  104 T ð3 CÞ rST ðkg L1 Þ ¼ 0:9193  6:65  104 T ð3 CÞ rpol ðkg L1 Þ ¼ 1:19  8:07  104 T ð3 CÞ rsol ðkg L1 Þ ¼ 0:892  1:3  103 T ð3 CÞ

Molecular

M A ðkg mol

weights

M B ðkg mol

1

Þ ¼ 0:1422

1

Þ ¼ 0:10415

1

M I ðkg mol

M sol ðkg mol

Þ ¼ 0:13216

1

Þ ¼ 0:10617

0.782 0.827 1.079 0.713

ARTICLE IN PRESS L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

133

0.40

[STY] (mol/L)

0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00

0

3000

6000

0

3000

6000

9000

12000 15000 Time (s)

18000

21000

24000

0

3000

6000

9000

12000 Time (s)

18000

21000

24000

9000

12000 15000 Time (s)

18000

21000

24000

0.35 0.30 [BMA] (mol/L)

0.25 0.20 0.15 0.10 0.05 0.00

9 8

Mn (kg/mol)

7 6 5 4 3 2 1 0 15000

Fig. 2. Comparison of reduced model predictions (lines) and the experimental data (points) for STY (top) and BMA (middle) monomer concentrations and polymer number-average molecular weight Mn (bottom) from Li and Hutchinson (2006).

where x1 and x2 represent the total mass of unreacted BMA and STY in the reactor, while x3 is the mass of unreacted initiator. States x4 and x5 represent the mass of BMA and STY contained in the polymer chains, and the number of moles of polymer is represented by x6 ; u1 , u2 and u3 are the input mass flowrates of the three species. The total moles of radicals in the system is calculated assuming radical stationarity:

ltot o



2kd f ¼ MI

1=2 

rtot ðtÞ kcop t

1=2 

1=2 x3 Mr ðtÞ

where system density rtot ðtÞ is calculated assuming volume additivity of the reactor components and the total mass in the reactor is calculated as Mr ðtÞ ¼ x1 þ x2 þ x3 þx4 þ x5 þ msol where msol is the initial charge of solvent in the reactor at the start is the termination rate of polymerization. The variable kcop t coefficient in the copolymerization system and was adjusted to fit reduced model to the experimental data. The complete set of

model parameters is summarized in Table 1 and the model fit is shown in Fig. 2. The concentration of radicals of BMA (superscript A) and STY (superscript B) is A lAo ðtÞ ¼ ltot o ðtÞfR ðtÞ;

B lBo ðtÞ ¼ ltot o ðtÞfR ðtÞ

where the radical fraction of each type is calculated according to the long-chain hypothesis: fRA ðtÞ ¼

kpBA M B x1 kpBA M B x1 þ kpAB M A x2

fRB ðtÞ ¼ 1  fRA ðtÞ The rate of production of new polymer chains, required to calculate polymer molecular weight, is defined by sol 1 Rm0 ðtÞ ¼ Rinit ðtÞ þ Rmon tr ðtÞ þRtr ðtÞ  2Rtc ðtÞ

ARTICLE IN PRESS 134

L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

where the rate of radical generation from initiation, transfer to monomer, transfer to solvent, and termination by combination is   1 2kd f rtot ðtÞx3 Rinit ðtÞ ¼ Mr ðtÞ M I   x1 x2 Rmon þ m ðtÞ ¼ m ðtÞ ðtÞ ltot r ðtÞ 1 2 tf Mr ðtÞ Mr ðtÞ o tot Rsol tr ðtÞ ¼

sðtÞ tot l ðtÞrtot ðtÞ Mr ðtÞ o A

B

A

B

2 2 Rtc ðtÞ ¼ kcop t ðð1  aÞðlo ðtÞÞ þ ð1  bÞðlo ðtÞÞ þ 2ð1  gÞlo ðtÞlo ðtÞÞ

with m1 ðtÞ ¼

m2 ðtÞ ¼ sðtÞ ¼

1 MA 1 MB

msol M sol

A mon B ðkmon tr AA fR ðtÞ þ ktr BA fR ðtÞÞ

B mon A ðkmon tr BB fR ðtÞ þ ktrAB fR ðtÞÞ

A sol B ðksol trAA fR ðtÞ þ ktr BA fR ðtÞÞ

and the fraction of termination events that occur by disproportionation for AA, BB, and AB termination given by a, b and g, respectively. The reduced model derived in this work capture the essential dynamics of the full model and the experimental data (Li & Hutchinson, 2006). The small mismatch in BMA level early in the reaction is caused by the neglect of depropagation in the reduced model. The primary references for the rate coefficients and parameters summarized in Table 1 are contained in the work of Li and Hutchinson (2006). Various quantities of interest can be calculated from the states. For the semibatch reactor, monomer conversion is defined as the mass of polymer in the reactor divided by the mass of monomer and polymer:

Xm ðtÞ ¼

x4 þ x5 x1 þ x2 þ x4 þ x5

Number-average molecular weight and copolymer composition (BMA mole fraction), denoted by Mn and F, are calculated as follows:    x4 þ x5 x4 x4 x5 1 Mn ðtÞ ¼ ; FðtÞ ¼ þ x6 MA MA MB An objective of the starved-feed strategy is to keep these latter two quantities uniform throughout the course of the semibatch reaction, while maintaining high conversion in the reactor. Thus, the objective function is defined so as to minimize the batch time and maintain the copolymer molecular weight ðMn Þ and composition (F) at their targets values throughout the batch while also satisfying path and terminal constrains on the states and/or the controls.

3. Real-time optimization Industry and academia have shown an increasing interest in using measurement-based optimization methods such as model predictive control (MPC) for trajectory tracking (Jarupintusophon, Lelann, Cabassud, & Casamatta, 1995; M’hamdi, Helbig, Abel, & Marquardt, 1996; Nagy & Agachi, 1997; Peterson, Hernandez, Arkun, & Schork, 1992) and on-line re-optimization strategies (Fiber, Gesthuisen, & Engell, 2002; Ruppen, 1994; Terwiesch, 1995). This tendency goes hand in hand with recent developments in sensor technology that are becoming prevalent in many industrial settings. In this work, we consider an on-line optimization technique (Peters, 2005) to design a control structure that

meets the specific objectives of the fed-batch polymerization reaction system presented in the previous section. This technique incorporates the feedback measurements directly into the optimization procedure rather than into a low level tracking controller so it truly becomes an online optimization rather than an online re-optimization, thus avoiding the need for controller design and the related complexities. The driving idea behind the method is that as measurements come in from the batch, new profiles are determined on-line and implemented simultaneously without the need of tracking the calculated profiles. The proposed approach differs from existing online optimization techniques (Abel & Marquardt, 2003; Chen & Allgower, 1998; Fiber et al., 2002; Mayne, Rawlings, Rao, & Scokaert, 2000; Richards & Congalidis, 2006; Ruppen, 1994; Terwiesch, 1995) since it involves minimal computational requirements, does not require knowledge of a nominal optimal solution and provides stable closed-loop performance. The basic control problem addressed here is to find the input trajectory that solves the dynamic optimization problem for some user defined cost functional. A dynamic optimization problem is formulated describing the feed flow rates parameters and the time interval lengths as decision or optimizing variables. For tracking purposes the objectives for Mn and F have been aggregated into a single scalar function while defining the batch time as a terminal objective function:  2  2 ! Z tr Mn ðtÞ FðtÞ dt ð2Þ o1  1 þ o2 1 J¼ Msp Fsp 0 where J is the objective function to be minimized, and Msp and Fsp are the target values for molecular weight and composition, respectively. The set of weights to scale the objectives are defined by o1 and o2 . The optimization is subject to process dynamics (1) and the constraints described below. In this study, the monomer and initiator levels in the batch must not exceed certain maximum limits to meet safety constraints related to heat release. These constraints are described by ; 0 rxk rxmax k

xk þ 3 Z0;

k A f1; 2; 3g

ð3Þ

Placing the constraints on unreacted monomer in the reactor is a convenient way to limit potential heat release, as this quantity is directly related to the maximum temperature rise in the reactor that could occur in the case of a complete loss of heat removal capacity in the system. Control constraints are often dictated by actuator limitations, with non-negativity of flow rates, a common input constraint: ; 0 ruk ðtÞ r umax k

k A f1; 2; 3g

ð4Þ

Terminal constraints are commonly related to safety or productivity considerations. With regard to the latter, the desired mass of polymer at the batch end is specified as an endpoint constraint: x4 ðtf Þ þx5 ðtf Þ ¼ mpol ðtf Þ

ð5Þ

The control functions are parameterized in three distinct strategies. The first (Case 1) and simplest strategy considers using a single constant: ðuk ðtÞ ¼ uk , k ¼ 1; 2; 3Þ through the entire batch time horizon subject to a starved-feed operation. The corresponding control parameterization is given by Case1 u1 ðtÞ ¼ y1 ;

u2 ðtÞ ¼ y2 y1 ;

u3 ðtÞ ¼ y3 y1

ð6Þ

The second parameterization (Case 2) considered in this work allows for exponential feeding in a starved-feed operation. This parameterization is written as follows: Case2 u1 ðtÞ ¼ y1 ey2 t ;

u2 ðtÞ ¼ y3 u1 ðtÞ;

u3 ðtÞ ¼ y4 u1 ðtÞ

ð7Þ

The final parameterization (Case 3) allows for independent exponential feeding strategy. The objective is to verify the cost of

ARTICLE IN PRESS L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

utilizing a starved-feed strategy in practice. This slightly more complex parameterization is given by

good control of polymer quality (molecular weight and composition) using starved-feed operation (for Cases 1 and 2). To handle constraint, we transform the constrained problem to the minimization of the unconstrained cost:

Case3 u1 ðtÞ ¼ y1 ey2 t ;

u2 ðtÞ ¼ y3 ey4 t ;

u3 ðtÞ ¼ y5 ey6 t

ð8Þ

Jip ðyÞ ¼

Note that (8) can be reduced to either (6) or (7) for specific choices of the parameters ðy1 , y2 , y3 , y4 , y5 , y6 Þ. Cases 1 and 2, however, maintain the current plant convention of ratioing the comonomer and initiator feeds. The objective of the optimization technique is to steer the system to the local optimum of the user-defined cost (2) subject to a specific parameterization (i.e. from (6), (7) and (8)) while responding in real-time to changes in plant operation. Clearly the technique cannot be used to recover the actual optimal control. The quest for an optimal control strategy cannot be performed in realtime, or, at least, not in sufficient time to provide some responsiveness to disturbances and sudden process changes. The technique proposed provides suboptimal control strategies but it allows the adjustment of the parameters of the input parameterization in a way that provides continuous improvement with respect to the user-defined cost function. It also ensures constraint satisfaction. It results in a manageable tuning process for algorithm parameters and gain factors, and takes advantage of the generally

50:50 BMA:STY (wt ratio)

mfed TBPA (kg)

9:7  103 2:46  101

mfed STY (kg) tf (min)

2:46  101 360

m X

 2  2 Mn ðtÞ FðtÞ  1 þ o2 1 Msp Fsp 1 logðwj þ eÞA dt þMðwf Þ2

wk ¼  xk þxmax þ e; k

wk þ 3 ¼ xk þ e;

wk þ 6 ¼  uk þ umax þ e; k

k ¼ 1; 2; 3

wf ¼ x4 þ x5  mpol ðtf Þ Numerical values for the bounds of the path constraints are given with the cases studies. The integral part of the modified cost is considered as an additional state such that 0 1  2  2 m X   t t M ð Þ Fð Þ n  1 þ o2  1  m1 log wj þ e A x_ 7 ¼ @o1 Msp Fsp j¼1 Finding a locally optimal solution can often be achieved through appropriate initial guesses to the optimizer and the use of realistic constraints. However, to find a feasible set of initial guesses is sometimes not a trivial task due to the interior feasible point method (barrier function) used. The barrier method requires that the initial guess of the solution is strictly feasible. In addition, ill conditioning inherent to the nonlinearity of the reactor system may play a role in the sensitivity of the solution to initial conditions. The constrained parameter set

Case 2 Base Case

0.2

Case 2 Base Case

0.1

Setpoint 1

wk þ 9 ¼ uk þ e;

The end point constraints are as follows:

0.3

0

k ¼ 1; 2; 3

The control constraints are given by

6

0

ð9Þ

where the m path constraints are given by

0.4

2

o1

j¼1

8

4

tf

m1

F

Mn (kg/mol)

mfed BMA (kg)

Z 0

Table 2 Conventional starved feed strategy used in Li and Hutchinson (2006). Results

135

Setpoint 0

2 x 104

0

1

2 x 104

0.8

0.6

0.3

Mr (kg)

mpol (kg)

0.4

0.2 Case 2

0.1 0

0.4

Case 2 Base Case

Base Case 0

1

2 Time (s)

3 x 104

0.2

0

1 Time (s)

2 x 104

Fig. 3. Closed-loop reactor trajectories for Mn (kg/mol), F (BMA mole fraction), polymer mass (kg), and total reactor mass (kg), as calculated using the basecase model (no manipulation of inputs) and Case 2 dynamic optimization.

ARTICLE IN PRESS 136

L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

x1 (t)

0.04 Case 2

0.02

Base Case 0

0

0.5

1

1.5

2

2.5 x 104

x2 (t)

0.02

Case 2

0.01 0

Base Case 0

0.5

1

1.5

2

2.5 x 104

x 10−3

x3 (t)

1

Case 2

0.5 0

Base Case 0

0.5

1

1.5

2

2.5 x 104

Time (s)

Fig. 4. Closed-loop reactor trajectories for x1 (kg BMA), x2 (kg STY), and x3 (kg TBPA) for Case 2.

2.5

x 10−5

u1

2 1.5 1

0

2000

4000

6000

8000

10000

12000

14000

2000

4000

6000

8000

10000

12000

14000

2000

4000

6000 8000 Time (s)

10000

12000

14000

−5

2.5

x 10

u2

2 1.5 1

0 x

10−7

u3

10

5 0

Fig. 5. Closed-loop reactor trajectories for u1 (kg/s BMA feed), u2 (kg/s STY feed) and u3 (kg/s TBPA feed) for Case 2.

O ¼ fyn A RN JwðxðtÞ; fðy; tÞÞ Z 0; 8t A ½0; tf g describes a convex subset of RN . It is assumed that the parameters belong to a compact subset Op of RN , and that the cost functional J : RN -R is convex and continuously differentiable on Op . The assumptions made up to this point guarantees a local optimization of the constrained problem exists and that the gradient can be used to achieve that minimization (Peters, 2005). An interior point method with penalty function is used to include the constraint costs. The interior point method incorporating a log barrier function

enforces the state and input constraints (essentially converting the constrained optimization problem into an unconstrained one) while the end-point constraints are incorporated through a terminal penalty function (Peters, 2005). In the remaining equations obvious notation has been omitted. Thus, let the path cost with barrier function be stated as follows: Lðx; uÞ ¼ qðx; uÞ  m1

nc X i¼1

logðwi ðx; uÞ þ eÞ

ARTICLE IN PRESS

9

9

8

8

7

7

6

6 Mn (kg/mol)

Mn (kg/mol)

L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

5 4

5 4 3

3 2

0

0.5

1 Time (s)

1.5

1 0

2 x 104

Fig. 6. Comparison of Mn (kg/mol) based on the online dynamic optimizations calculated with Cases 1–3 input parameterizations.

Case 1 Case 2 Case 3 Base Case Setpoint

2

On−line Off−line Base Case Setpoint

1 0

137

0

1 Time (s)

1.5

2 x 104

Fig. 8. Comparison of Mn trajectories calculated using the basecase model, off-line optimization, and the online dynamic optimization for Case 2.

0.45

0.45

0.4

0.4

0.35

0.35

0.3

0.3

0.25

0.25 F

F

0.5

0.2

0.2

0.15

0.15

0.1

Off−line

0.05

Base Case

0

0.5

1 Time (s)

1.5

Case 2 Case 3

0.05

Setpoint

0

Case 1

0.1

On−line

Base Case Setpoint

0

0

2 x 104

0.5

1 Time (s)

1.5

2 x 104

Fig. 7. Comparison of F (copolymer composition) based on the online dynamic optimizations calculated with Cases 1–3 input parameterizations.

Fig. 9. Comparison of F (copolymer composition) trajectories calculated using the basecase model, off-line optimization, and the online dynamic optimization for Case 2.

To emphasize that the optimization is based on the current conditions, the new cost functional with interior point method and penalty function is stated as Z tt Jip ðt; yðtÞ; xðtÞÞ ¼ Lðxp ðtÞ; fðyðtÞ; tÞÞ dt þ MJxp ðtf Þ  xf J2

guarantees that the parameters remain in a convex set Ow ¼ fy A Op jJyJ r Zg for some Z 40 while ensuring that the cost decreases until the optimal profile parameters are determined. A Lyapunov-based method was used to show convergence to the local minimizers of a user-defined cost functional (Guay & Zhang, 2003; Peters, 2005). The projection algorithm is given by ( g if JyJ r Z or ðJyJ ¼ Z and ry PðyÞ r 0Þ _y ¼ ð10Þ c otherwise

t

where t A ½t; tf  is the integration variable and nc is the number of inequality constraints. The parameter m1 4 0 is the barrier function parameter for the logarithmic term, M 40 is the penalty term and e 4 0 is a constraint relaxation factor that prevents the barrier term from singularity. The formulation uses a gradientbased method for the solution of the dynamic optimization problems in real time, with a straightforward diagonally scaled steepest descent parameter update law of the form:

y_ ¼ ProjðGry Jip ; Ow Þ where G 40 is a scaling matrix to be defined later. To avoid divergence of the update law a projection algorithm is used. This

T

where c ¼ g  gry PðyÞry PðyÞ=Jry PðyÞJ2 , PðyÞ ¼ y y  Z. The projection algorithm as defined guarantees that y A Ow and hence in Op . In addition, let us consider the total cost: Jtot ðtÞ ¼

Z

t

LðxðtÞ; fðyðtÞ; tÞÞ dt þ

0

Z t

tf

Lðxp ðtÞ; fðyðtÞ; tÞÞ dt þ MJxp ðtf Þ  xf J2

ð11Þ p

where x ðtÞ are state predictions of the model equations (1) evaluated over the interval t A ½t; tf  with initial conditions

ARTICLE IN PRESS 138

L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

xp ðtÞ ¼ xðtÞ where xðtÞ is the current measurement of the state variables from the fedbatch reactor. The first term of the RHS of (11) is the elapsed cost of the fedbatch reactor subject to the evolving input trajectories given by fðyðtÞ; tÞ for t A ½0; t. It can be shown that the projection algorithm yields a parameter update law, y_ , that guarantees that J_ tot r0:

ð12Þ

Thus choosing the parameter update formula in this way guarantees that the total cost is improved continuously. Of course, this special property relies on the assumption that the process model provides a precise description of the process dynamics. In fact, it is possible to show that the property (12) is robust to relatively small perturbations and process uncertainties. This is an inherent property of most optimization based control such as nonlinear model predictive control. In the context of model predictive control, the problem of robustness has been a central research theme over the last ten years. Using the results in DeHaan and Guay (2007), one can easily deduce the nominal robustness of any gradient based dynamic optimization routine such as the technique proposed in this study. In many aspects, the proposed technique is very similar to model predictive control. However, there are some important distinctions with existing approaches. The proposed technique allows one to use a low dimensional parameterization of the input trajectories. This contrasts with standard approaches that rely on a piecewise constant discretization of the entire input trajectories. Under the assumption of initial feasibility, the proposed technique guarantees continuous improvement of the cost function using even the simplest possible parameterization. Another important property of the proposed approach is that it does not require the solution of a dynamic optimization problem at every sampling instant as in other techniques. It relies on a single gradient update at each step which requires a minimal number of function evaluations. Overall the proposed approach provides a lot of the features of MPC such as constraints satisfaction and nominal robustness but implements the dynamic optimization as a standard adaptive controller using the minimal amount of computations required at each sampling instant. The effectiveness of the technique will be demonstrate in the following section. 4. Case studies Cases 1–3 illustrate the optimization of the BMA/STY system relative to the conventional starved feed strategy for the 50:50 BMA:STY mass ratio. Cases 1–3 follow experimental initial conditions (monomer and initiator feed rates) used in the experimental study of Li and Hutchinson (2006) and summarized in Table 2. In fed fed Table 2, mfed TBPA , mBMA and mSTY represent the total mass of initiator, BMA and STY fed to the reactor to the final batch time tF . The mass of polymer produced at the final batch time is mpol ðtf Þ, taken as the terminal (end-point) constraint in the optimization routine. For all three cases, the target values of F and Mn are 0.43 and 8, respectively. The constraints for the states and input variables are set to x1 ðtÞ A ½0; 0:03;

x2 ðtÞ A ½0; 0:03;

x3 A ½0; 5:0  103 

and ui ðtÞ A ½0; 0:01, i= 1,2,3. The algorithm parameters for all three cases are as follows. The adaptive gain (k) in the steepest descent law used by the algorithm is 0.0001. The end constraint parameter M is set to 108 . This large value ensures that endpoint values are strictly enforced. The barrier function parameters, m and e, are set to 0.1 and 107 , respectively. The initial conditions for the system are xi ð0Þ ¼ 0;

i ¼ 1; . . . ; 6

Figs. 3–5 show the simulation results for Case 2. In Fig. 3, the outputs of the system are shown including, Mn , F, the mass of polymer produced and the total reactor mass. The simulation results based on the nominal parameter values corresponding to the experimental conditions listed in Table 2 are labeled as basecase on the plots and correspond to the model fit of the experimental data with no optimization of input feeds. It is important to note that these trajectories correspond to the initial trajectory. The real-time optimization routine provides an effective and efficient mechanism to update the profile online using the process measurements. In this specific case, it results in a reduction of total batch time from 21 600 to 13 500 s. Fig. 4 shows the corresponding state measurements for x1 , x2 and x3 . The state trajectories for x1 and x2 are shown to approach their upper limit. The input feed trajectories are shown in Fig. 5 (see Figs. 6 and 7). In Figs. 8 and 9, we compare the simulation results for Mn and F for the nominal basecase simulation result (dashed line), the realtime dynamic optimization result (dotted line) and the optimum trajectory computed off-line (dash-dotted line). The off-line optimal trajectories are based on the input parameterization of Case 2. The composition trajectories are identical, but the off-line optimum is more aggressive and yields a slight overshoot in Mn . The final batch time of 10 500 s with the off-line optimization indicates that the real-time optimization technique, with a batch time of 13 500 s, results in a suboptimal result. However, it is possible to use the current input trajectories as a starting point for a subsequent batch, which will then yield an improved trajectory. In Figs. 6 and 7, the results for the three case studies are compared. All three cases show improvement with Case 3 providing the greatest improvement followed by Cases 2 and 1. This is as expected since Case 3 provides the most flexibility and degrees of freedom in the choice of input trajectories. Case 3 provides a final batch time of 12 900 s, while Case 1 yields batch time of 13 500 s.

5. Conclusions This paper presented the application of a real-time dynamic optimization method for the optimization of a starved-feed semibatch copolymerization reactor. Using a validated process model, the study demonstrates that the real-time optimization technique yields significant improvement in process productivity (reduction of batch-time by as much as 40%) while improving Mn control and maintaining starved-feed operation. The results show that the real-time optimization technique can be used in a continuous quality improvement program for a batch process control system. It guarantees continuous improvement in the process performance both within batch and from batch to batch.

Acknowledgments The authors would like to acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada. References Abel, O., & Marquardt, W. (1998). A model predictive control scheme for safe and optimal operation of exothermic semibatch reactors. In Proceedings of IFAC symposium Dycops5 (pp. 761–766). Abel, O., & Marquardt, W. (2003). Scenario-integrated on-line optimization of batch reactors. Journal of Process Control, 13, 703–715. Chen, H., & Allgower, F. (1998). Nonlinear model predictive control schemes with guaranteed stability. In Nonlinear model based process control (pp. 465–494). Dordrecht, The Netherlands: Kluwer Academic Publishers. DeHaan, D., & Guay, M. (2007). Nonconvex optimization and robustness in realtime model predictive control. International Journal of Robust and Nonlinear Control, 17, 1634–1650.

ARTICLE IN PRESS L.F.T. Perea et al. / Control Engineering Practice 18 (2010) 131–139

Fiber, S., Gesthuisen, R., & Engell, S. (2002). Re-optimization based control of copolymer quality in semi-continuous emulsion polymerization process. In Proceedings of the American control conference anchorage, Anchorage, Alaska. Guay, M., & Zhang, T. (2003). Adaptive extremum-seeking control of nonlinear systems with parametric uncertainties. Automatica, 30. Jarupintusophon, P., Lelann, M. V., Cabassud, M., & Casamatta, G. (1995). Realistic model-based predictive and adaptive control of batch reactor. Computers and Chemical Engineering, 18, s445–s449. Li, D., & Hutchinson, R.A (2006). High temperature semibatch free radical copolymerization of butyl methacrylate and styrene. Macromolecular Symposium, 243, 24–34. Mayne, D. Q., Rawlings, J. B., Rao, C. V., & Scokaert, P. (2000). Constrained model predictive control: Stability and optimality. Automatica, 36(6), 789–814. M’hamdi, A., Helbig, A., Abel, O., & Marquardt, W. (1996). Newton-type receding horizon control and state estimation—a case study. In Proceedings of the 13th World congress of IFAC (pp. 121–126).

139

Nagy, Z., & Agachi, S. (1997). Model predictive control of a pvc batch reactor. Computers and Chemical Engineering, 21, 571–591. Peters, N. J. (2005). Real-time dynamic optimization of nonlinear batch systems. MSc thesis, Queens University. Peterson, T., Hernandez, E., Arkun, Y., & Schork, F. J. (1992). A nonlinear dmc algorithm and its applications to a semibatch polymerization reactor. Chemical Engineering Science, 47(4), 737–753. Richards, J. R., & Congalidis, J. P. (2006). Measurement and control of polymerization reactors. Computers and Chemical Engineering, 30, 1447–1463. Ruppen, D. (1994). A contribution to the implementation of adaptive optimal operation for discontinuous chemical reactors. PhD thesis, ETH Zrich. Terwiesch, P. (1995). Cautious on-line correction of batch process operation. AIChE Journal, 41(5), 1337–1340. Vemuri, Y. J. (2004). Real-time optimization of semi-batch reactor. PhD thesis, Florida State University.