Multiobjective Dynamic Optimization of the Cell ... - ACS Publications

Article pubs.acs.org/IECR

Multiobjective Dynamic Optimization of the Cell-Cast Process for Poly(methyl methacrylate) Martín Rivera-Toledo,*,† Ehecatl Antonio Del Río-Chanona,‡ and Antonio Flores-Tlacuahuac† †

Departamento de Ingeniería y Ciencias Químicas, Universidad Iberoamericana Prolongación, Paseo de la Reforma 880, México D.F., 01219 México ‡ Facultad de Química, Departamento de Ingeniería Química, Universidad Nacional Autónoma de México (UNAM), México D.F., 04510 México ABSTRACT: During the normal operation of plastic sheet polymerization reactors, large monomer conversion and average molecular weight should be achieved. However, improving any one of these two operating objectives would result in worsening the remaining one, and vise versa. Because of physical constraints, it would not be possible to meet both operating objectives in similar extent. Therefore, trade-off solutions must be obtained from which the designer can pick up the one that best suits her/his design objectives. In this work we compute the whole set of trade-off solutions by solving a set of multiobjective optimization problems. In this work we select as the best trade-off solution, the point on the trade-off curve closest to the utopia region. We demonstrate that dynamic optimal warming temperature profiles improve the control on the weight-average molecular weight and monomer conversion simultaneously, and lead to a polymer sheet with better homogeneity characteristics, when trade-off solutions based on multiobjective optimization are used.

■

INTRODUCTION The world economy and high commercial pressure demand manufacturing better and cheaper products, reducing capital and operating costs. Hence, productivity maximization, without compromising product quality, is the top priority for process industries. In polymerization systems product quality is specified within a narrow range, and productivity is enhanced by keeping product specification within such a narrow range. To achieve these processing objectives, optimal operating policies can be used.1 In polymerization reactors the average molecular weight determine several physical properties, e.g., strength and impact resistance, and they can be treated as an indirect measurement of product quality, whereas the polymerization time is a direct measurement of productivity. Therefore, the simultaneous minimization of polymerization time and maximization of average molecular weight distribution (MWD) can lead to the goal of productivity maximization without affecting product quality. However, frequently target operating objectives are in conflict. Examples of industrial systems where conflicting objectives emerge are as follows: the structure and manufacturing processes of large composite plastic parts,2 design of bioprocesses,3 biochemical reaction network,4 for drawbead design in sheet metal forming5 and polymerization processes.6,7 In all of these practical applications mentioned above, two or more objectives are in conflict, and no objective can be improved without sacrificing the remaining ones. Multiobjective optimization (MO) refers to finding values of decision variables which correspond to and provide the optimum of more than one objective. MO can be applied to handle conflicting performance objectives. The goal is to obtain a set of equally good solutions, known as Pareto optimal solutions. In a Pareto set, no solution can be considered to be better than any other solution with respect to all objective functions. When one moves from one Pareto solution to © 2014 American Chemical Society

another, at least one objective function improves while at least one other one gets worse. Hence, MO involves special methods for considering more than one objective and analyzing the results obtained. In this work, we formulate a multiobjective optimization problem using conflicting performance objectives in polymerization systems. Specifically, we compute optimal operating policies such as to minimize the distribution of molecular weight (DMW) and to maximize monomer conversion. We have to emphasize that computing multiobjective optimal operating policies of the addressed polymerization system turns out to be a complex problem, because of the two-dimensional distributed and dynamic nature of the polymerization system. In this work we demonstrate that improved homogeneous properties in the polymerization process of plastic sheets can be obtained using dynamic optimal control policies based on a specific trade-off solution located on the corresponding Pareto front. This turns out to be particularly true when the operating policies are compared against single objective dynamic optimal solutions, or against empirical operating policies. Moreover, the dynamic optimal results were obtained by using an industrially validated dynamic mathematical model of the polymerization plastic sheet process, which give us confidence in that similar results can also be observed by applying the off-line optimal policies to the real industrial system.

■

LITERATURE REVIEW The most often exploited approaches to generate the Pareto set are the weighting method and the ε-constraint method,8 Received: Revised: Accepted: Published: 14351

April 4, 2014 July 9, 2014 August 27, 2014 September 5, 2014 dx.doi.org/10.1021/ie5014162 | Ind. Eng. Chem. Res. 2014, 53, 14351−14365

Industrial & Engineering Chemistry Research

Article

although more recent approaches are available.9 Despite surveys of MO methods being available, they are often incomplete in terms of comprehensive coverage, algorithm presentation, and general applicability to engineering design. A good review of this topic has been published by Marler and Arora.10 In the present work, we used a well-known strategy called normal boundary intersection (NBI)9 for getting the Pareto set for the MO of the cell-cast process for poly(methyl methacrylate) to be described later. Following, we discuss some of the references about MO in polymerization systems. BonillaPetriociolet et al.11 presented a study using genetic algorithms for the acrylonitrile and vinyl acetate copolymerization system with the objective of finding the characteristic correlations parameters that are used to determine the autoacceleration reaction phenomenon12 called “gel effect”. This phenomenon arises as a consequence of the viscosity increase into the reaction medium, due to the formation of polymer molecules when small solvent concentrations are presented. The result is a rate constant termination decrease and a conversion increase. These authors observed that temperature improves monomer conversion and its mole fraction. However, an increment in operation temperature increases polydispersity and reduces average molecular weight. Benyahia et al.,13 presented a multicriteria optimization formulation for styrene (St) and butyl acrylate to determine the dynamic optimal feed rate profile necessary to produce core−shell latex particles with specific end-use properties depending on the application, e.g., paints or adhesives. The model proposed by the authors predicts the global monomer conversion, the average molecular weight, the particle size distribution, and the residual monomers mass fraction. An evolutionary algorithm to determine the Pareto set was applied, and a unique solution was selected and experimentally implemented. Gomes et al.14 reported optimal control policies for emulsion terpolymerization of St, methyl methacrylate (MMA), and methyl acrylate (MA), which were determined in a semibatch reactor using the multiobjective dynamic optimization method. Six variables were used as manipulated variables: three monomers feed rates (St, MMA, and MA), surfactant and initiator feed rates and the reactor temperature as well. Their results showed that the optimization procedure could simultaneously minimize the reaction time and monomer flows to obtain the desired composition of the polymers, particle size distribution, and molecular weight. The work showed high temperature influence on molecular weight, due to the effect of the diffusion of the radical ends. In those cases the objectives were often conflicting, which means achieving the optimum for one objective requires some compromise on one or more other objectives.

Figure 1. Cell-cast process for PMMA plastic sheet manufacture.

observed that even operating the new reactor configuration under empirical isothermal conditions, improved homogeneous polymer properties were achieved. This observation led us to conclude that to maximize homogeneous polymer properties, dynamic optimal warm air temperature profiles should be used, to ensure product quality (through MWD) and to minimize the time required for polymerization product with uniform polymerization properties. As conflicting objectives, we have selected the maximization of monomer conversion and the minimization of MWD. Control of MWD is important in industrial polymerization processes because end-use properties of polymers are strongly dependent on its MWD. In concentrated systems a strong acceleration in the polymerization rate occurs when the reactions become controlled by diffusion phenomena called gel effect. In bulk conditions, the glass effect produces a rapid decrease in the polymerization rate and a limitation in the final conversion is reached. The chain length of the polymer produced during polymerization is directly related to the reaction kinetics, as a consequence also the number-average molecular weight (Mn), the weight-average molecular weight (Mw) and the molecular weight distribution are strongly affected by the diffusion phenomena. During the gel effect a strong increase in the values of the Mw is observed; as a consequence, the MWD becomes broader and a bimodality in the distribution can be reached. The glass effect produces a decrease in the average molecular weights in the final phase of the reaction. The values of Mn and Mw and the MWD curve depend on the choice of operating conditions, monomer and initiator concentrations, and the heating policy. The temperature profile used during the operation plays an important role. In the present study, our goal is to compute the air temperature profiles so to drive the plastic sheet temperature to ensure that the average molecular weight and monomer conversion of the final polymerization process have the desired target values as soon as possible (i.e., minimum polymerization time), taking into account that the two objectives are in conflict. In the open literature, there have been relatively few research papers dealing with the optimization of convective ovens for PMMA plastic sheet manufacturing purposes. Zhou et al.16 have used a multiobjective optimization approach for the continuous casting of thin films of PMMA using an adaptive genetic algorithm. In a previous work,17 we have presented an optimization formulation based on the dynamic representation of the two-directional heating, but we just used a single objective framework. In this work we claim that better homogeneous polymerization properties can be achieved if

■

PROCESS DESCRIPTION The traditional poly(methyl methacrylate) (PMMA) manufacturing process is based on placing plastic sheets inside warm water baths to allow polymerization reactions to proceed until monomer and initiators are exhausted, as can be seen in Figure 1. However, industrial experience indicates that the quality of the PMMA plastic sheets tend to display wide variations of their MWD due to nonuniform heating patterns. These distributions are observed along the thickness and length of the plastic sheet. To achieve uniform polymerization properties, it has been suggested to carry out the set of PMMA polymerization reactions using an oven where heating proceeds by circulating warm air. In pilot plant experiments,15 we have 14352

dx.doi.org/10.1021/ie5014162 | Ind. Eng. Chem. Res. 2014, 53, 14351−14365


Article

cally. Consequently, the objective functions should be transformed such that they are dimensionless. The most robust approach to transforming objective functions, regardless of their original range, is given as follows: ψ̂ k(x,u) = (ψk(x,u) − ψmin k )/ (ψmax − ψmin k k ), this approach is consistently referred to as normalization, where ψ̂ k(x,u) stands for the normalized objective function and it has values between zero and one. During the past decades, many methods have been proposed to deal with multiobjective optimization problems. Full reviews can be found in the books by Liu19 and Yann20 and references therein. Following, we provide a brief description of three wellknown methods to compute the Pareto front: weighted sum (WS), epsilon-constraint (EC), and normal boundary instersection (NBI) methods. Weighted Sum Method. A standard technique for generating the Pareto set in multicriteria optimization problems is to minimize (convex) weighted sums of the different objectives for different settings of the weights. In other words, k weights wk are chosen such that wk ≥ 0, k = 1, 2, ..., n, and n ∑k=1 wk = 1 and the following problem is solved: min ∑nk=1wkψk(x,u). However, it is well-known that this method succeeds in getting points from all parts of the Pareto set only when the Pareto curve is convex. Das and Dennis21 provide a graphical interpretation of the weighted sum method for twoobjective problems to address some of its deficiencies. Furthermore, these authors show that varying the weights consistently and continuously may not necessarily result in an even distribution of Pareto optimal points and an accurate, complete representation of the Pareto optimal set. Finally, they discuss in detail the necessary conditions to yield an even spread of points on the Pareto curve in the objective space. ε-Constraint Method. Haimes et al.22 introduced the εconstraint approach (also called the ε-constraint or trade-off approach). This technique minimizes the single most important objective function ψi, whereas the remaining n − 1 objective functions are added as inequality constraints of the form ψk(x,u) ≤ εk for all k = 1, 2, ..., n, k ≠ i where i ∈ 1, 2, ..., n. The vector of upper bounds, ε = (ε1, ε2, ..., εn), defines the maximum value that each objective can take. These inequalities can be interpreted as hyperplanes reducing the feasible criterion space. To obtain a subset of the Pareto optimal set, one must vary the vector of upper bounds along the Pareto front for each objective, and run a new optimization problem for each new vector. Normal Boundary Intersection Method. To remove the deficiencies of the weighted sum approach to compute the Pareto front, Das and Dennis9 introduced the normal boundary intersection (NBI) method. This method provides a means for obtaining an even distribution of Pareto optimal points for a consistent variation in the user-supplied parameter vector w, even with a nonconvex Pareto optimal set. This method essentially works by solving a set of NLPs of the following form:

compromise solutions, based on considering multiobjective optimization problems, are used. Basic Concepts and Multiobjective Optimization Methods. The general multiobjective dynamic optimization problem is posed as follows: min Ψ(x ,u) = [ψ1(x ,u), ψ2(x ,u), ..., ψk(x ,u)]T x ,u

(1a)

subject to dx = f (x ,u), initial conditions: t = t0 , x = x0 dt

(1b)

h(x ,u) = 0

(1c)

g (x ,u) ≤ 0

(1d)

xlb ≤ x ≤ x ub

(1e)

ulb ≤ u ≤ u ub

(1f)

where x ∈ 9 is the state vector, u ∈ 9 is the vector of manipulated variables, Ψ = 9 n × m → 9 k is the vector of objective functions, f = 9 n × m → 9 n is the map representing the dynamic process behavior, and h = 9 n × m → 9 ne + nu is the set of system constraints. In this formulation n is the number of states, m is the number of manipulated variables, k is the number of objective functions, ne is the number of equality constraints, nu is the number of inequality constraints, lb stands for lower bound, ub is an upper bound, and t is the processing time. Ψ(x,u) is a vector of objective functions ψk(x,u), also called objectives, criteria, payoff functions, cost functions, or value functions. The objective functions are assumed to be conflicting so that no one can be minimized without increasing the other, this situation gives rise to the concept of a Pareto solution. Definition 1. Pareto Solution. A feasible point x* for the multiobjective problem (1) is said to be Pareto optimal if and only if there exists no other feasible point (x) such that ψk(x) ≤ ψk(x*) ∀k ∈ Rn and ψk(x) < ψk(x*) for at least one index k ∈ Rn. The family of Pareto solutions form the Pareto front, which represents a limiting curve of performance in the objective space. Definition 2. Utopia Point. The utopia point is a point given by the solution of min ψk(x) ∀k ∈ Rn, and subject to the following constraints: dx/dt = f(x,u), h(x) = 0, g(x) ≤ 0, xlb ≤ x ≤ xub, and ulb ≤ u ≤ uub with respect to the vector of decision variables x. It is important to emphasize that the minimization is given for each individual objective function. The utopia point is unattainable because the objectives are conflicting. However, it can still be used as a reference point. For instance, it is possible to compute the closest point along the Pareto front to the utopia point, also known as the compromise solution. Definition 3. Compromise Solution. The compromise solution is a point (xcs) located on the Pareto front with objective ψ(xcs) given by the solution of the following minimum Euclidian distance problem: n

min{ ∑ [ψk(x ,u) − ψk(x up , u)]2 }1/2 x ,u

k ∈ Rn

m

min γ

(NBIF)

x ,u,γ

subject to dx = f (x ,u), initial conditions: t = t0 , x = x0 dt

(2)

However, it is not necessary to restrict closeness to the case of a Euclidean norm.18 In addition, if different objective functions have different units, the Euclidean norm or a norm of any degree becomes insufficient to represent closeness mathemati-

h(x ,u) = 0

g (x ,u) ≤ 0 14353



Article

xlb ≤ x ≤ x ub

effects. In the past 35 years, several models have been published dealing with the mathematical description of diffusioncontrolled kinetic rate constants in free-radical polymerization.26 The reaction mechanism adopted here consists of a simple approximation of the well-known bulk free-radical polymerization kinetics featuring straightforward initiation, propagation, chain transfer to monomer, and termination reactions. Although there are more complete kinetics polymerization models (VHW) proposed in the literature,27 we found that the model (CCS) proposed in ref 28 was enough to reproduce the experimental pilot plant data for the PMMA plastic sheet15 polymerization system. In Figure 2, the VHW

ulb ≤ u ≤ u ub

Φw + γv = Ψ(x ,u) − ψ up

which will be called the “normal boundary intersection formulation (NBIF)”. In this formulation, Φ is a n × n payoff matrix in which the ith column is composed of the vector Ψ(x*i ,u) − ψup, Ψ(x*i ,u) is the vector of objective functions evaluated at the minimum of the ith objective function. The diagonal elements of Φ are zeros; w is a vector of scalars such that ∑ni=1wi = 1 and w ≥ 0, v = −Φe, and e ∈ Rk is a column vector of ones in the criterion space. v is called a quasi-normal vector. Because each component of Φ is positive, the negative sign ensures that v points toward the origin of the criterion space. Accordingly, for any w, a solution point is independent of how the objective functions are scaled. As w is systematically modified, the solution of the optimization problem (NBIF) yields an even distribution of Pareto optimal points representing the complete Pareto set. In this work, we compute the whole Pareto frontier of the MMA polymerization reactor using the NBI method. In the following section, we present the problem statement, previously to the multiobjective dynamic optimization problem formulation. Problem Statement. Cell-Cast Process Description. In the typical casting of acrilic sheet, molds formed by two glass plates separated by a peripheral gasket sealer and clamped together are filled with casting syrup through a gap left in the gasket.23,24 The casting syrup is made up of partially polymerized monomer (20%) which, once placed in the mold, is inserted into a furnace that is heated by circulating warm air (Figure 1). It is important to control the progress of the polymerization reactions and to create suitable mild thermic conditions which, in turn, requires quick and effective dissipation of excess heat. However, due to the low heat capacity of air, effective control of the thermal conditions during the operation turns out to be difficult. In fact, the heating pattern is affected by the temperature of the circulating air as has been shown in ref 15. Therefore, through the proper control of the air temperature profile, acceptable molecular weight distribution and monomer conversion can be achieved. As previously stated, commonly these operating objectives are in conflict. Hence, a trade-off solution must be obtained which meets both operating objectives in the best possible manner. Hence, the problem to be solved in this work can be formulated as follows: “Given is a dynamic, two dimensional model of a MMA plastic sheet reactor. We would like to compute the set of optimal trade-off solutions (Pareto front), such that molecular weight distribution and monomer conversion targets can be picked up from this set of optimal solutions. Having computed the target values of the operating objectives, then dynamic open-loop optimal control policies will be calculated to reach the target values as soon as possible using air temperature as the control variable”. Kinetic Model. Besides chemical kinetics, physical phenomena related to the diffusion of various chemical reactive species are important in free-radical polymerization reactions. In fact, at high monomer conversions, almost all elementary reactions can become diffusion-controlled. Reactions influenced by diffusion phenomena include termination of live macroradicals,25 propagation of a growing chain, and chemical initiation reactions. Diffusion-controlled termination, propagation, and initiation reactions have been related to the gel, glass, and cage

Figure 2. Monomer conversion dynamic time evolution when a “series approach” (···) (VHW) and a “parallel approach” (−−) (CCS) are used for computing the diffusion controlled kinetic rate constants: kt, kp, and kf at 50, 70, and 90 °C and initiator loading of 0.0258 mol/dm3. The circles represent experimental data.29

and CCS kinetic models and experimental data for monomer conversion in a batch reactor at 50, 70, and 90 °C are shown. As can be observed, the CCS model correctly predicts the variation in monomer conversion until 70 °C, over the conversion range from 0 to 95% conversion. However, around 90 °C, the precision of the CCS model decreases. The kinetic mechanism of a chemically initiated free-radical polymerization process can be described in terms of the following elementary reactions.30 Initiation: k

I → 2R• ki

R• + M → R1• Propagation: kp

R •n + M → R •n + 1

Chain transfer to monomer: kf

R •n + M → P•n + R1• Termination by disproportionation: k td

R •n + R •m → Pn + Pm Termination by combination: 14354



Article

k td

of reaction. The heat of reaction is given by (−ΔHr)(−rp) where ΔHr is the molal specific heat of reaction and rp is the molal rate for propagation reaction step. ΔHr changes along the x and z axes and the time coordinates. The solution domain will be divided into a series of nodes (i.e., the discretization points). At each one of these nodes the mass balance equations read as follows:

R •n + R •m → Pn + m

The initiation reaction involves the decomposition of the initiator, I, to produce radicals R•, reacting with the monomer molecules, M, to initiate new live (radical) polymer chains R•1 . During the propagation step, monomer molecules are added to the live-polymer chains, R•n and R•m. The growth of the chains terminates when the propagating radicals lose their activity through either of the termination reactions, resulting in deadpolymer chains, Pn, Pm, and Pn+m. Plastic Sheet PMMA Mathematical Model. The plastic sheet polymerization mathematical model was derived assuming that the heating source resulting from the set of polymerization reactions is a function of the local temperature and reactant/product concentration. It was also assumed that polymer properties changed along the thickness (x-axis) and length (z-axis) of the plastic sheet, as depicted in Figure 3.

∂I εI = −kdI − λ 0(k p + k f )(1 − X ) ∂t 1 + εX

(4a)

∂X = λ 0(k p + k f )(1 − X ) ∂t

(4b)

∂λ 0 ελ 0 2 =− (k p + k f )(1 − X ) + 2fkdI − k tλ 0 2 ∂t 1 + εX (4c)

ελ λ ∂λ1 = − 1 0 (k p + k f )(1 − X ) + 2fkdI − k tλ 0λ1 ∂t 1 + εX 1−X + (k p + k f )λ 0 M 0 (4d) 1 + εX

ελ 2λ 0 ∂λ 2 =− (k p + k f )(1 − X ) + 2fkdI − k tλ 0λ 2 ∂t 1 + εX 1−X (2λ1 + λ 0) + k pM 0 (4e) 1 + εX ∂μ0 ∂t

=−

εμ0 λ 0 1 + εX

(k p + k f )(1 − X ) + k tdλ 0 2 +

1 k tcλ 0 2 2 (4f)

Figure 3. Control volume for energy balance. Ta is the warming air temperature, 2H is the thickness, L is the length, W is the width, q is the amount of heat transferred along each direction.

∂μ1 ∂t

=−

εμ1λ 0 1 + εX

(k p + k f )(1 − X ) + k tdλ 0λ1 + k tcλ 0λ1 (4g)

Hence, polymer properties were taken as constants along the yaxis. For simplicity the glass plates are considered homogeneous and isotropic (i.e., the properties are the same in all directions (thermal conductivity, heat capacity, density)), average values were taken for MMA and PMMA. To obtain the two-dimensional dynamic energy balance, the total heat entering and leaving at the x and z coordinates was modeled by the Fourier law, and the rate of change of energy within the control volume was obtained by applying the shell energy balance method.31 Because of symmetry considerations, only half of the sheet thickness (from the center to the external surface) was taken into account. Dynamic mass and energy balances coupled through polymerization kinetics30 describe polymer conversion and molecular weights time evolution. Within the reactor, forced convection warm air provides the required energy to raise the plastic sheet temperature until a point where significant polymerization rates take place. Within the polymer bulk the dominant heat transfer mechanism is conduction. Therefore, heat is transferred along the x and z axes giving rise to the following two-dimensional dynamic heat transfer equation: ⎛ ∂ 2T ∂T ∂ 2T ⎞ Q rxn = α⎜ 2 + ⎟+ ∂t ρCp ⎝ ∂x ∂z 2 ⎠

∂μ2 ∂t

=−

εμ2 λ 0 1 + εX

(k p + k f )(1 − X ) + k tdλ 0λ 2

+ k tc(λ 2λ 0 + λ12)

(4h)

The initial (IC) and boundary (BC) conditions of the model are given by IC: t = 0,

∀x ∈ [0, H ], z ∈ [0, L], T = T0 , I = I0 , X = X 0 , λ 0 = λ 00 , λ1 = λ10 , λ 2 = λ 20 , μ0 = μ0 , μ1 = μ1 , μ2 = μ2 0

0

0

(5a)

∂T =0 ∂x

BC1: t > 0,

x = 0, ∀z ∈ [0, L],

BC2: t > 0,

x = H , ∀z ∈ [0, L], − k

(5b)

∂T = ha(T − Ta) ∂x (5c)

BC3: t > 0,

z = 0, ∀x ∈ [0, H ], − k

∂T = ha(T − Ta) ∂z (5d)

(3)

BC4: t > 0,

where T is the polymer temperature, t is the polymerization time, x is the axis of the sheet thickness, z is the axis for the sheet length, ρ is the density, Cp is the heat capacity, α is the thermal diffusivity, and Qrxn stands for the polymerization heat

z = L , ∀x ∈ [0, H ],

∂T =0 ∂z

(5e)

where To and Mo are initial temperature and monomer concentration, respectively, Ta is the warming air temperature, X is the monomer conversion, λ0, λ1, and λ2 are the zeroth-, 14355



Article

first-, and second-order moments of growing radical distribution, respectively, whereas μ0, μ1, and μ2 are the zeroth-, first-, and second-order moments of dead polymer, respectively. The subscript “o” on λi and μi, i = 0, 2, 3 stands for initial values of these variables. I is the initiator concentration, k is the monomer average thermal conductivity, f is the initiator efficiency, ε is the volume expansion factor, ha is the heat transfer coefficient, kd, kp, kf, and kt are the kinetic coefficients for initiation, propagation, chain transfer to monomer, and termination steps, respectively. Information regarding the PMMA kinetic rate constants is displayed in Table 1. Following

To determine the full chain-length distribution of the polymer, a large number of mass balances must be solved. In practice, however, it suffices to calculate only the leading moments of the chain-length distribution, where the kth moment of the live-radical (λk) and dead-radical (μk) chaink length distributions are defined as λk = ∑∞ n=1n Rn, and μk = ∞ k ∑n=1n Pn, k = 0, 1, 2, respectively. The weight-average molecular weight (Mw) is calculated from the leading moments through the following equation: Mw = Mm w [λ2 + μ2/(λ1 + μ1)], where Mm is the monomer molecular weight. w In deriving the mathematical model, we have assumed that the diffusional mass effect was negligible when compared against diffusional thermal effects. The assumption is reasonable for this system in which competing transport processes occur. This may be demonstrated by examining the magnitude of the ratio of the thermal to mass diffusivities (as measured by the Lewis number). If the dimensionless Lewis number (Le) is defined as Le = thermal diffusivity/mass diffusivity = αavg/Davg, and the physical properties, like density, heat capacity, thermal conductivity and diffusion coefficient are taken as average values for MMA and PMMA at temperature range from 298 to 400 K, then the Lewis number value ranges from about 6 to 3200. For the present case study the temperature range was taken from 298 to 350 K. Therefore, the Lewis number range varies from 6 to 2000. Hence, it becomes clear that thermal diffusivity has a stronger influence on the system behavior, and this is the main reason why mass diffusivity effects were neglected. The above-described mathematical model features wide variation in the numerical values of some system states spanning several orders of magnitude. To improve the numerical robustness of the model, and to reduce potential difficulties related to the numerical integration and optimization, the model was properly scaled. The dimensionless variables are defined as follows: τ = αt/H2, ξ = x/H, ζ = z/ L, θ = T/T0, θa = Ta/T0, I ̅ = I/I0, λ̅k = ln(λk), μ̅k = ln(μk), k = 0, 1, 2. It should be stressed that due to significant volume contraction of PMMA with respect to MMA (around 20%), strictly speaking the polymerization process features true moving boundary behavior. However, introducing moving boundaries might lead to a problem complex to solve. One of our aims was to keep both the problem formulation and solution as simple as possible. The justification for using fixed boundaries was the realization that even using fixed boundaries our numerical results seem to match the available experimental industrial information. In fact, in a previous publication15 we have shown a comparison of experimental pilot plant results versus simulation ones for two different thicknesses (3 and 6 mm) of the plastic sheet. Comparison against experimental data shows that the prediction capabilities of the proposed model seem to be satisfactory for industrial application, especially for small thickness plastic sheets. From Figure 12 of the abovementioned publication, it can be observed that the shrinkage is homogeneous for 3 and 6 mm thickness plastic sheets. However, higher deviation are observed for 12 and 18 mm thickness plastic sheets. Even in this last case, deviations are not too large (around 4%). These observations support our decision of addressing the polymer contraction behavior as a fixed boundary problem. Moreover, volume contraction was taken into account by

Table 1. PMMA Kinetic Rate Constants12 apparent propagation coefficient

⎛ 2.3ϕm ⎞ 1 1 = + θp(T )λ0 exp⎜− ⎟ kp k p0(T ) ⎝ A + B(T ) ⎠

apparent termination coefficient

⎛ 2.3ϕm ⎞ 1 1 = + θt(T ,I0)λ0 exp⎜− ⎟ kt k t0(T ) ⎝ A + B(T ) ⎠

monomer diffusion time for propagation stage

⎛ 14000 ⎞⎟ θp = exp⎜− 35.167 + ⎝ T ⎠

monomer diffusion time for termination stage

θt =

parameter A in the Fujita− Doolittle equation

A = − 8 × 10−6(T − Tg)2 + 0.1678

parameter B

B = 0.03

true propagation coefficient

⎛ 4353 ⎞ ⎟ k p0 = 2.95 × 107 exp⎜− ⎝ 1.987T ⎠

true termination coefficient

⎛ 701 ⎞ ⎟ k t0 = 5.88 × 109 exp⎜− ⎝ 1.987T ⎠

initiator coefficient

⎛ 14855 ⎞ ⎟ kd = 4.25 × 1016 exp⎜− ⎝ T ⎠

volume fraction of the monomer

ϕm =

initiator concentration [mol/dm3] kinetic coefficient for termination by disproportionation kinetic coefficient for termination by addition

I0 = 0.0258

k tc = k t − k td

kinetic coefficient for chain transfer to monomer

⎛ − 13880 ⎞ ⎟ k f = 9.48 × 103k p exp⎜ ⎝ 1.9872T ⎠

⎛ 3.55481 × 103 ⎞ 0.48139 × 10−22 exp⎜ ⎟ I0 1.987T ⎠ ⎝

1−X 1 + εX

⎡ ⎛ 4090 ⎞⎤ ⎟⎥ k td = k t /⎢1 + 3.956 × 10−4 exp⎜ ⎝ 1.987T ⎠⎦ ⎣

we explain with detail the meaning of each boundary condition. In BC1, we have assumed that the system temperature reaches its maximum value because of symmetry conditions. Hence, ∂T/∂x = 0, ∀z ∈ [0, L]. In BC2, at the system surface, the temperature is governed by interface transport and it is given by the Fourier heat conduction law and the well-known Newton’s law of cooling. Therefore, −k∂T/∂x = ha(T − Ta), ∀z ∈ [0, L]. Regarding BC3, this equation simply states the continuity of the flux of heat. Before reaching this point, the heat flux is governed by convective transport; when heat hits the glass surface, the heat flux equation is governed by conductive transport. At the system interphase (z = 0) a steady-state energy balance states the conservation of the system energy: heat transported by convection ought to be equal to the energy transported by conduction. Hence, −k∂T/∂z = ha(T − Ta), ∀x ∈ [0, H]. Finally, in BC4, the glass surface is insulated from the environment, so no heat is transported from this point into the surroundings. Hence, ∂T/∂z = 0, ∀x ∈ [0, H].

V = Vo(1 + εX )

(6)

where ε is the volume contraction factor and computed from 14356


Industrial & Engineering Chemistry Research ε = (ρm − ρp )/ρp

Article

Using the simultaneous approach for solving the dynamic optimization problem, given by the objective function (eq 8) and the model and operating constraints (eqs 3−5 in their dimensionless form), amounts to transform the original dynamic optimization problem into a NLP. The transformation involves two steps: (a) The method of lines34 was used to take care of the spatial coordinates transforming the original PDAE system into a system of ordinary differential algebraic equations (ODAE). (b) Orthogonal collocation on finite elements was used to transform the ODAE system into a set of nonlinear algebraic equations.35 The integral term in the objective function was approximated using Radau cuadrature. Hence, the discretized NLP problem reads as follows:

(7)

where ρm and ρp are the monomer and polymer densities, respectively.

■

MULTIOBJECTIVE DYNAMIC OPTIMIZATION PROBLEM FORMULATION Physical properties of polymers strongly depend upon the molecular weight distribution. Typical goals in operating a polymerization reactor are to minimize polydispersity or the number molecular weight distribution to produce a polymer featuring target mechanical and rheological properties.32 Nunes et al.32 found that thermal properties, stress−strain properties, impact resistance, strength and hardness of films of poly(methyl methacrylate), and polystyrene were all improved by narrowing the MWD. Hence, the development of a methodology for adjusting the MWD in the manufacturing of MMA plastic sheets should lead to a polymer with improved target properties. Depending upon the intended use of the polymer product, polymers with various types of MWD are required.33 Experimentally, it has been observed that homogeneous polymer properties in the MMA cell cast process can be obtained by proper selection of warming air temperature profiles. Because the MMA cell cast process is essentially a batch process, most of the related polymer properties strongly depend upon the time variation of operating parameters (i.e., warming air temperature). Therefore, the dynamic optimal operation of the reactor should lead to a better and improved MMA homogeneous polymerization system. We propose a multiobjective dynamic optimization formulation when dealing with conflicting performance objectives in polymerization systems. Our goal is to determine dynamic optimal heating patterns leading to homogeneous and improved plastic sheet properties such that trade-off solutions are obtained trough the computation of the Pareto set. To achieve this aim, we use the following objective functions: min fM = θa , Y

w

min fX = θa , Y

∫0

∫0

τ

τ

2 ⎡ 2 ⎛ θa mk ⎞ ⎛ θmki ⎞ ⎢ ⎜ ⎟ min ∑ Δτm ∑ ∑ ⎜1 − d ⎟ + ⎜1 − d ⎟ ⎢⎝ θamk , Ymki θ ⎠ θa ⎠ ⎝ m=1 k=1 i=1 ⎣ Ne

Nc

Nz

2 ⎛ M wmki ⎞ ⎤ ⎥ ⎜ ⎟ + ⎜1 − ⎟ M wd ⎠ ⎥⎦ ⎝

(9a)

2 ⎡ 2 ⎛ θamk ⎞ ⎛ θmki ⎞ ⎢ min ∑ Δτm ∑ ∑ ⎜1 − d ⎟ + ⎜⎜1 − d ⎟⎟ ⎢⎝ θamk , Ymki θ ⎠ θa ⎠ ⎝ m=1 k=1 i=1 ⎣ Ne

Nc

Nz

2⎤ ⎛ X mki ⎞ ⎥ + ⎜1 − d ⎟ ⎝ X ⎠ ⎥⎦

(9b)

where θmki denotes the value of dimensionless temperature at the mth finite element, k is the collocation point within that element, and i is the position point for the sheet thickness (z axis). Ne is the number of finite elements, Nc is the number of collocation points, and Nz is the number of points for discretization along the z coordinate, whereas Δτm is the length of mth finite element. The dynamic optimization formulation is subject to the following set of discretized constraints. The dimensionless states approximation are given by the following Lagrange polynomial interpolation:

⎡ 2 ⎞2 ⎛ ⎛ ⎞2 ⎤ ⎢⎛⎜1 − θ ⎞⎟ + ⎜1 − θa ⎟ + ⎜1 − M w ⎟ ⎥ dτ ⎜ ⎟ ⎟ ⎜ ⎢⎝ M wd ⎠ ⎥⎦ θd ⎠ θad ⎠ ⎝ ⎝ ⎣ (8a)

Nc

0 ̇ Ymki = Y mi + Δτmτtr ∑ A nk Ymni

⎡ 2 2⎤ ⎞2 ⎛ ⎢⎛⎜1 − θ ⎞⎟ + ⎜1 − θa ⎟ + ⎛⎜1 − X ⎞⎟ ⎥ dτ ⎟ ⎜ ⎢⎝ ⎝ X d ⎠ ⎥⎦ θd ⎠ θad ⎠ ⎝ ⎣ (8b)

n=1

(10)

where are the states, Y = [θ, X, I,̅ λ̅0, λ̅1, λ̅2, μ̅0, μ̅1, μ̅2] , at the starting point of each element, τtr is the operating time. The states continuity between finite elements is enforced through the following constraint: T

Y0mi

In the above equation θd, θda , Mdw, and Xd stand for the desired values of dimensionless monomer temperature, dimensionless air temperature, weight-average molecular weight, and monomer conversion, respectively, whereas Y is the vector of system states (θ, X, λ̅k, μ̅k, k = 0, 1, 2). The ratio between the rate of propagation and the rates of propagation and termination of polymerization was used to obtain Mw. In these objective functions we have incorporated three terms, including θ, θd, θa, θda , because our aim is to compute the air temperature as time function (θa) so to drive the plastic sheet temperature θ, and monomer conversion and/or Mw to its desired target values as soon as possible. We can use either f Mw or f X as the main objective function and move the remainder one as a constraint. In any case the resulting MO problem is subject to the partial differential and algebraic equations (PDAE) and the initial and boundary conditions (eqs 3−5), of course, these must be expressed in their dimensionless form.

Nc

0 Y mi = Y m0 − 1i + Δτm − 1τtr ∑ A nNc Yṁ − 1ni n=1

(11)

and, finally, the Ẏmni vector, representing the right-hand side of the dynamic model, is evaluated at the operating conditions, according to the eqs 3 and 4, in dimensionless form.

■

RESULTS AND DISCUSSION In this section we show and discuss the dynamic response of the plastic sheet polymerization reactor when three different operating policies are applied: (a) simple operating policy, (b) optimal control policy resulting from the multiobjective compromise solution and (c) optimal control policy resulting from solving a single objective optimization problem. We also discuss the way the Pareto front was obtained. In all cases, the 14357



Article

Table 2. Design Data for Sheet Reactor monomer initial concentration [mol/dm3] initiator initial concentration [mol/dm3] heat of reaction [kJ/mol] polymer density [kg/m3] polymer thermal conductivity [W/(m K)] polymer heat capacity [J/(kg K)] initial air temperature [K] initial monomer temperature [K] sheet length [m] sheet thickness [m] initiator efficiency [dimensionless]

M0 = 9.98 I0 = 0.0258 ΔH = −58.19 ρ = 1200 k = 0.09 Cp = 1674 Ta0 = 338.15 T0 = 298.15 L = 1.8 2H = 0.018 f = 0.58

Figure 4. Sketches corresponding to the location of the four places where computed results are reported.

dynamic two-dimensional mathematical model for a plastic sheet thickness of 18 mm was used. In Table 2, design and optimization data for the PMMA plastic sheet reactor are shown. We use the method of lines34 for approaching the spatial discretization of the PDAE model of the polymerization reactor. The discretized equations thus obtained form a set of differential algebraic equations that were numerically integrated as an initial value problem by using the dassl solver.36 After some trials, we found that 10 and 6 discretization points suffice to represent the system dynamic response along the longitude and thickness of the plastic sheet, respectively. Monomer temperature, monomer conversion, and weight-average molecular weight dynamic responses were computed in around 7 min CPU time. The CPU times of the calculations carried out in this section are summarized in Table 3. The results of the calculations are only presented for certain points along the sheet. The reason for doing so is because experimentally it is easier to measure the plastic sheet temperature in those places. Accordingly, the results are presented for four points along the plastic sheet (Figure 4): at the center of the figure and the first extreme along the plastic sheet [x = 0, z = 0], as well as the last one [x = 0, z = L], and the surfaces at the first and the last points along the plastic sheet, [x = H, z = 0] and [x = H, z = L], respectively. Moreover, it is worth stressing that the experimental verification of the nonlinear reactor model is crucial in a practical implementation, if the simulation results are to be trusted. In a previous work,15 we have shown an experimental validation of the two-dimensional mathematical model by comparing such a response against industrial pilotplant experiments. Hence, the results presented in this section should provide a clear indication of the polymerization reactor response when the optimal control policies are implemented in the real experimental facility. Computation of the Utopia Point and Pareto Set. To compute the utopia point coordinates (denoted as f LMw and f LX points in Figure 5) we solved the following optimization problem. The optimal solution of the NLP problem given by

Figure 5. Pareto set (○), utopia point (□), compromise solution (▽), and two additional solutions arbitrarily chosen (◁, ▷) using the L and f XL are the normal boundary intersection technique. f M w coordinates of the utopia point, whereas f UMw and f UX stand for the maximum values of the f Mw and f X objective functions, respectively.

eqs 9a, 10, and 11 provides the f LMw and f UX coordinates, whereas by solving the eqs 9b, 10, and 11 the f UMw and f LX points are obtained. In all cases, three internal collocation points, 20 finite elements and 6 and 10 discrete points along the x and z directions, respectively, were used. The optimization formulations were solved using the IPOPT NLP solver37 available in the AMPL Mathematical Programming environment.38 In all the addressed multiobjective optimization problems the number of equations is 33 662 and number of decision variables is 240. The Pareto curve considering both f X and f Mw objective functions was computed using the NBI method previously described. It must be stressed that all the reported values of the f X and f Mw objective functions discussed in this section refer to unscaled values. According to Das and Dennis,9 the problem was solved for a single objective function optimization with β ∈ [0, 1]. The Pareto front was obtained considering the solution

Table 3. Results Obtained on an Intel Core i5, 2 GB RAM, and 2.67 GHz Laptop Computer

a

no.

case

method or algorithm

fX

f Mw

points

iter

1 2 3 4 5 6

simulation 18 mm utopia utopia compromise Pareto set single optimizationa,17

IPOPT IPOPT IPOPT NBI IPOPT

0.0329 0.0446 0.0341

3.2871 0.1508 0.3335

1 1 1 47 1

867 510 314

18.524

CPU time [min] 6.62 10.01 2.03 1.86 261.0 23.85

This case was computed on a Pentium IV computer with 2.8 GHz clock speed and 512 KB RAM. 14358



Article

Figure 6. Control results obtained from applying a “simple operating policy”: The temperature of the plastic sheet was kept constant at 50 °C. The results shown in the first and second columns refer to z = 0 and z = L, respectively. The solid line stands for results at the center of the plastic sheet (x = 0), whereas the dashed line represents results at the surface of the same figure.

computing the full Pareto front consists in obtaining only the compromise solution (i.e., a single point on the Pareto front) as recently demonstrated in ref 18 for the online model predictive control of systems based on multiobjective optimization criteria. However, regardless of the method being used, computing the full Pareto optimal set clearly can be an important task, especially when looking for the effect of different process constraints on the shape and position of the Pareto curve, as well as to propose strategies to select suitable operating policies in an open-loop environment. In Figure 5, two additional solutions has been arbitrarily chosen (◁, ▷), these are located on the Pareto front and the first one is on the left side of the compromise solution and the second one is on the right side of it. Their dynamic optimal polices will be explained later in the “multiobjective optimal control operating policies” section.

of 50 multiobjective optimization problems, 47 of which reached optimal solutions. This procedure allowed us to obtain nearly homogeneous distribution of points, which can be observed in Figure 5 (black dots). Running this algorithm took almost four and half hours for the construction of the whole Pareto curve. The compromise solution (the one on the Pareto front closest to the utopia region) is given by the f X = 0.0341, f Mw = 0.3335 coordinates. These results are in agreement with Arora and Marler,10 in the sense that this technique was able to produce an even spread of points on the Pareto front with a systematic change in their parameters. In a previous work,39 it has been shown that the NBI method is superior against εconstraint and weighted sum techniques because the NBI subproblem does not require exact Hessian values, meaning that it can use any NLP technique. An alternative to the idea of 14359



Article

Figure 7. Dynamic optimal results for the compromise solution. The results shown in the first and second columns refer to z = 0 and z = L, respectively. The solid line stands for results at the center of the plastic sheet (x = 0), whereas the dashed line represents results at the surface of the same figure. In the two top figures, the gray line denotes the computed dynamic optimal temperature profile.

Optimal Control against Simple Operating Policies. a. Simple Operating Policy. To obtain target polymer properties, industrial pilot plant tests have shown that the plastic sheet temperature should be kept constant around 50°C during most of the polymerization process.24 It should be stressed that there is not any kind of support for using this heuristic rule beyond experimental measurements. In fact, target polymer properties are complex functions of other operating variables besides the operating temperature. As a matter of fact, improved average homogeneous polymer properties can be obtained if dynamic nonconstant temperature profiles are used instead.17 In this section we use the term simple operating policy to denote a flowing air temperature profile resulting in the target polymerization temperature (50 °C). The open-loop results obtained from applying this operating policy are shown in Figure 6. As depicted in Figure 6A, because of the warm air flow

direction, at z = 0 the temperature rises quicker compared to the temperature increase at z = L, as displayed in Figure 6B. Indeed, the internal temperature profile (x = 0) is larger than the external one (x = H) because of the cooling effect of air flowing over the plastic sheet. Moreover, as the polymerization reaction goes on, polymer viscosity increases, leading to heat transfer limitations. Hence, heat dissipation becomes harder to achieve. After 20 h, the sheet temperatures acquired the steadystate value. Here, steady state means that, after the reaction has been essentially completed, polymer properties (i.e., conversion, molecular weight averages, etc.) remain constant. However, at z = L, the temperature dynamic response is slower, but finally it also attains the steady-state value when the polymerization reaction has been completed. In both figures, a temperature peak can be clearly seen. The peak is related to the onset of the gel effect, and it is larger at z = L because, at this 14360



Article

average molecular weight profiles. It is well-known that the material properties (stress−strain, toughness, fracture, and fatigue behavior) depend on the molecular weight distribution of a given polymer.32 The molecular weight distributions for the plastic sheet reactor are plotted in Figure 6E,F. In these figures the runaway effect on the molecular weight behavior, associated with the gel effect, is clearly seen. As a consequence, this effect suffices to obtain nonuniform mechanical properties. It can be seen that Mw decreases along the plate, between initial and final edges in agreement with similar reported results.40 The net effect of temperature rising is to decrease the molecular weight distribution. In summary, the results indicate that the use of simple operating policies can result in polymer products with nonhomogeneous properties. Therefore, in this work we claim that the use of optimal nonconstant operating policies could result in polymer products with nearly homogeneous properties. This statement will be validated in the following items. b. Multiobjective Optimal Control Operating Policies. In this section we present the dynamic behavior of the plastic sheet reactor when optimal flowing air temperature profiles, obtained using the Pareto compromise solution, are used. It should be pointed out that the warm air temperature has been arbitrarily initialized at 65°C. From Figure 7A,B we note that dynamic optimal operating policies given by compromise solution lead to a reduction of the time required to manufacture PMMA plastic sheets, in comparison when the above stated simple operating policy is used. This saving in processing time translates into a reduction of heating utilities resulting in increased process profit. Moreover, optimal temperature policies, based on the compromise solution, contribute to obtain polymer homogeneous properties as discussed in this section. Using the optimal warming temperature profile, the polymerization process has almost finished in around 350 min at z = 0, whereas using the simple operating policy temperature profile the gel effect has not yet emerged. This behavior is also clearly seen at z = L. Moreover, the optimal warming air temperature profile gives rise to a smaller gel effect temperature peak. As depicted in Figure 7, due to the warm air flow direction, the monomer temperature increases gradually along the plastic sheet; at z = 0 the dynamic temperature behavior is almost the same as the air temperature, but at z = L, large temperature peaks occur, because the heat removal capacity is decreased. We notice from Figure 7C that optimal temperature profiles (based on the compromise solution) clearly contribute to speedup the polymerization process. For instance, 40% monomer conversion was obtained in around 65 min, whereas 240 min was required using the simple operating policy. The molecular weight gradually increases with the polymerization temperature, but suddenly, a temperature peak appears due to contribution from the heat of reaction. From this moment on, the air temperature gradually fell down to remove the heat generated by the polymerization reaction. From a practical point of view, the main advantage of using the multiobjective optimal control policies consists in reducing the onset of undesired effects, such as bubble formation in the plastic sheet. Because of the higher initial air temperature, larger temperature peaks emerge. Moreover, the slope of the classical s-shaped gel effect curve changes rather quickly explaining the reduction in polymerization times. It is also interesting to note that this time Mw variations are less notorious. In fact, the increase in the polymerization temperature partially suppresses the large increase in Mw due to the gel effect. Although the warm air

Figure 8. (A) Optimal control profiles for two additional points on the Pareto set and percentage relative error of mean values for (B) monomer conversion and (C) molecular weight.

point, thermal effects are stronger due to the heat produced by the reaction contribution. The plastic sheet thermal behavior affects the monomer conversion. Indeed, we can see in Figure 6C that, for the first 350 min, the conversion is almost 35%, but suddenly, it becomes too high, 75% and 85% at the inner and surface plate positions, respectively. When the temperature approach (defined as T(x=0) − T(x=H)) is zero, the monomer conversion is 84%. Parts C and D of Figure 6 display the monomer conversion curves before and after the onset of the gel effect. These curves are concave upward preceding the sharp rise. If the gel effect were suddenly switched on at a certain critical conversion level, the portion of the conversion curve under examination would be concave downward due to the gradual depletion of monomer. Beyond this critical conversion, the gel effect algorithm is suddenly switched on in the model, causing a discontinuous slope in the conversion curve. Onset of the gel effect is only a convenient phrase denoting the region of rapid rise in the conversion. There exists in fact no sharp demarcation in terms of molecular processes before and after the occurrence of this autoacceleration region. The rapid rise in conversion is merely a natural consequence of the increasing importance of mass transfer limitations.12 The influence of the process variables on the characteristics of the product can be understood from the computed weight14361



Article

Figure 9. Dynamic optimal results for single objective. The results shown in the first and second columns refer to z = 0 and z = L, respectively. The solid line stands for results at the center of the plastic sheet (x = 0), whereas the dashed line represents results at the surface of the same figure. In the two top figures, the gray line denotes the computed dynamic optimal temperature profile.

optimal temperature profile features a step-like shape, the switching times may be difficult to estimate by heuristic or trial and error approaches. In summary, it can be observed that using optimal operating policies contributes to obtain homogeneous polymer product in comparison to the use of a simple operating policy. In this section we discuss about multiple solutions on the Pareto set, and how the optimal control profile changes depending on its location on the Pareto set. From Figure 5, we have chosen arbitrarily two points given by the (f X = 0.0320, f Mw = 2.708) (◁) and (f X = 0.0342, f Mw = 0.1515) (▷) coordinates. In Figure 8A three different control actions are shown. The first one is the dynamic optimal policy from the compromise solution (gray continuous line), the second (black continuous line) and third (dotted black line) ones are the optimal control policies on the Pareto set located on the left and the right side of the compromise solution, respectively; these trajectories are called “L-Compromise” and “R-Com-

promise”. In Figure 8B,C the percentage of relative errors for monomer conversion and Molecular weight against the mean behavior of the compromise solution are shown for the LCompromise and R-Compromise solutions. The mean behavior is defined as the sum of the state (monomer conversion or Mw) for four points [(z = 0, x = 0), (z = 0, x = H), (z = L, x = 0), (z = L, x = H)] divided by 4. From Figure 8A we observe that, for some operating periods, the compromise solution is not located between the L- and Rcompromise solutions. In fact, during the first 800 min, the compromise solution closely resembles the R-compromise solution. Also from this figure, we could conclude that the Rcompromise solution could provide similar quality product results in comparison when the compromise solution is used. However, we should stress that there is no way of obtaining this information with computing the whole Pareto front. Clearly, the L-compromise solution could demand larger amounts of auxiliary services because it results in larger temperature profiles 14362



Article

Notes

than the other two solutions. From this short discussion, it seems to be clear than the compromise optimal solution can lead to better and efficient optimal profiles in comparison when only single objective optimal control problems are considered. c. Single-Objective Optimal Control Operating Policies. Finally, in this section we compare the compromise solution against a case study reported in a previous work.17 In ref 17 the problem was solved using the weighted sum method for a single optimization function as follows: min∫ [w1(1 − θ/θd)2 + w2((1 − θa/θda )2) + w3((1 − Mw/Mdw)2) dτ] subject to the constraints given by eqs 3−5. The results are summarized in Table 3 and displayed in Figure 9. As can be seen from Figure 9E,F, improved Mw uniformity was achieved by using the warm air optimal profile computed from the single objective formulation. However, the monomer temperature picks are higher than the ones computed using the multiobjective formulation, as displayed in Figure 9B. It is well-known that careful control of the temperature polymerization is necessary to obtain a bubble-free product of good optical clarity.23 Therefore, to avoid bubble formation, we should not operate above the monomer boiling point (100 °C). In summary, although singleobjective optimal control policies contribute to improve polymer homogeneous properties, the results indicate that polymer properties can be maximized using the multiobjective compromise solution.

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors thank the referees, whose insightful comments helped us improve the quality of the paper.

■

CONCLUSIONS A multiobjective dynamic optimization framework has been developed to simultaneously optimize molecular weight-average dispersion and monomer conversion, which are conflicting targets. The optimization is performed by means of the normal boundary intersection technique. The Pareto set results clearly indicate that there is a trade-off between monomer conversion and Mw. We have demonstrated that dynamic optimal warming temperature profiles help to improve the control on the Mw and monomer conversion simultaneously and lead to a polymer sheet with better homogeneity characteristics, especially when compared against simple isothermal operating policies. Moreover, results based on the multiobjective optimization approach turned out to produce a plastic sheet product with better homogeneous properties, when compared to similar calculations done using a single-objective optimization framework. Even when the set of optimal solutions are off-line computed, presently we are extending this work to compute and implement online optimal dynamic solutions by using nonlinear model predictive control strategies. The optimal policies computed in this paper are relatively easy to implement and maintain. Unfortunately, the experimental validation of the optimized control was not performed because of resource limitation on experimental pilot plant. However, we recognize the importance of the experimental validation for our dynamic optimal policies computed that drives the process to ensures the satisfaction of product quality specifications. Even when no practical implementation of the optimal operating policies was done, we have used a previously validated dynamic mathematical model to compute the results shown in this work. Hence, the results presented should closely follow those observed in a real facility.

■

AUTHOR INFORMATION

Corresponding Author

*M. Rivera-Toledo. E-mail: [email protected]. Phone/fax: +52(55)59504074. URL: http://200.13.98.241/∼martin. 14363

NOMENCLATURE a = thickness-to-length ratio for plastic sheet [dimensionless] Akm = matrix element for orthogonal collocation method Bi = Biot number [dimensionless] Cp = heat capacity [J/(kg K)] e = column vector of ones in the NBI technique F = vector of objective functions f = initiator efficiency g = inequality constraint ha = heat transfer coefficient [W/(m2 K)] h = equality constraint H = sheet thickness [m] I = molar concentration of initiator [mol/dm3] I ̅ = dimensionless molar concentration of initiator k = thermal conductivity [W/(m K)] kd = kinetic coefficient for initiator [1/min] k̅d = dimensionless kinetic coefficient for initiator ki = kinetic coefficient for initiation reaction [dm3/(min mol)] kp = kinetic coefficient for propagation [dm3/(min mol)] k̅p = dimensionless kinetic coefficient for propagation kf = kinetic coefficient for chain transfer to monomer [dm3/ (min mol)] k̅f = dimensionless kinetic coefficient for chain transfer to monomer ktc = kinetic coefficient for termination by addition [dm3/ (min mol)] k̅tc = dimensionless kinetic coefficient for termination by addition ktd = kinetic coefficient for termination disproportionation [dm3/(min mol)] k̅td = kinetic coefficient for termination disproportionation L = sheet length [m] M = molar concentration of monomer [mol/dm3] Mn = number-average molecular weight [kg/kmol] Mw = weight-average molecular weight [kg/kmol] Mm w = monomer molecular weight [kg/kmol] Nc = number of collocation points [dimensionless] Ne = number of finite elements [dimensionless] Nx = number of points for discretization on x-direction [dimensionless] Nz = number of points for discretization on y-direction [dimensionless] Q = heat released by the polymerization reactions [W/m3] Q̅ = dimensionless heat released by the polymerization reactions R = universal gas constant [kJ/(mol K)] T = polymer temperature [K] Ta = air temperature [K] Tg = glass transition temperature of the pure polymer [K] T0 = initial monomer temperature [K] t = reaction time [min] v = quasi-normal vector in the NBI technique W = sheet width X = monomer conversion Y̅ = dimensionless Y state dx.doi.org/10.1021/ie5014162 | Ind. Eng. Chem. Res. 2014, 53, 14351−14365


Article

Ẏ = first-order derivative of the Y state x = x-axis for the sheet thickness [m] w = characterizing parameter of the NBI subproblem z = z-axis for the sheet length [m]

(5) Sun, G.; Li, G.; Gong, Z.; Cui, X.; Yang, X.; Li, G. Multiobjective robust optimization method for drawbead design in sheet metal forming. Mater. Des. 2010, 31, 1917−1929. (6) Gupta, S. K.; Garg, S. Multiobjective optimization of a free radical bulk polymerization reactor using genetic algorithm. Macromol. Theory Simul. 1999, 8, 46−53. (7) Zhou, F.; Gupta, S. K.; Ray, A. K. Modeling of the Sheet-Molding Process for Poly(methyl metacrilate). J. Polym. Sci. A: Polym. Chem. 2001, 81, 1951−1971. (8) Rangaiah, G. P.Multi-objective optimization: techniques and applications in chemical engineering; World Scientific: Singapore, 2009. (9) Das, I.; Dennis, J. Normal-boundary intersection: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems. SIAM J. Control 1998, 8 (3), 631−657. (10) Arora, J. S.; Marler, R. T. Survey of multi-objective optimization methods for engineering. Struct. Multidisc. Optim. 2004, 26, 369−395. (11) Bonilla-Petriciolet, A.; De-Alba-Pérez-de-Gracia, G. G.; Vallecillo-Gómez, S. V.; Tapia-Picazo, J. C. Free Radicals Copolymerization Optimization, System: Acrylonitrile-Vinyl Acetate in CSTR. 21 European Symposium on Computer Aided Process Eng. 2011, 20, 849−854. (12) Chiu, W. Y.; Carrat, G. M.; Soong, S. A Computer Method for the Gel Effect in Free Radical Polymerization. Macromolecules 1983, 16, 348−357. (13) Fonteix, C.; Pla, F.; Benyahia, B.; Latifi, M. A. Multicriteria dynamic optimization of an emulsion copolymerization reactor. Comput. Chem. Eng. 2011, 35, 2886−2895. (14) Altarawneh, I. S.; Srour, M. H.; Gomes, V. Optimal operating strategies for emulsion terpolymerisation. Chem. Eng. Sci. 2008, 63, 4257−4268. (15) Flores-Tlacuahuac, A.; Rivera-Toledo, M.; García-Crispín, L. E.; Vílchis-Ramírez, L. Dynamic modeling and experimental validation of the MMA cells cast process for plastic sheet production. Ind. Eng. Chem. Res. 2006, 45 (25), 8539−8553. (16) Ray, A. K.; Zhou, F.; Gupta, S. K. Multiobjective Optimization of the Continuous Casting Process for Poly(methyl methacrylate) Using Adapted Genetic Algorithm. J. Appl. Polym. Sci. 2000, 78, 1439− 1458. (17) Rivera-Toledo, M.; Flores-Tlacuahuac, A.; Vílchis-Ramírez, L. Dynamic optimization of the mma cell-cast process for plastic sheet production. AIChE J. 2009, 55, 1464. (18) Zavala, V. M.; Flores-Tlacuahuac, A. Stability of multiobjective predictive control: A utopia-tracking approach. Automatica 2012, 48 (10), 2627−2632. (19) Yang, J. B.; Liu, G. P.; Whidborne, J. F.Multiobjective optimization and control. International Journal of Adaptive Control and Signal Processing; John Wiley and Sons, Ltd.: New York, 2003. (20) Yann, C. Multiobjective optimization: principles and case studies. Decision engineering; Springer: Berlin, New York, 2003. (21) Das, I.; Dennis, J. F. A closer look at drawbacks of minimizing weighted sums of objectives for pareto set generation in multicriteria optimization problems. Struct. Optim. 1997, 14 (1), 63−69. (22) Wismer, L. S.; Haimes, D. A.; Lasdon, Y. Y. On a bicriterion formulation of the problems of integrated system identification and system optimization. EEE Trans. Syst. Man Cybern., SMC 1971, 1, 296−297. (23) Rosetti, C.Apparatus for Casting Plastic Sheet. U.S. Patent 3,689,022, 1973. (24) Wanwanichai, P.; Junkasem, J.; Saimaneewong, B.; Mahasan, R.; Tantivess, P.; Magaraphan, R.; Vanichvarakij, Y.; Petiraksakul, P.; Supaphol, P. Sheet-cast Poly(methyl methacrylate): One-step (Water) versus Two-step (Water-Air) Isothermal Processes. Iranian Polym. J. 2005, 14 (1), 61−69. (25) Odian, G.Principles of Polymerization; Wiley Interscience: New York, 1991. (26) Dubé, M. A.; Soares, J. B. P.; Penlidis, A.; Hamielec, A. E. Mathematical Modeling of Multicomponent Chain-Growth olymerizations in Batch, Semibatch, and Continuous Reactors: A Review. Ind. Eng. Chem. Res. 1997, 36, 966−1015.

Greek Letters

α = thermal diffusivity [m2/min] γ = maximization variable of NBI formulation ΔHr = heat of polymer reaction [kJ/mol] θ = dimensionless temperature θa = dimensionless air temperature θg = dimensionless glass transition temperature Θp = monomer diffusion time for propagation stage Θt = monomer diffusion time for termination stage ε = volume expansion factor [dimensionless] ζ = dimensionless variable for z coordinate λ0 = zeroth moment of the growing radicals [mol/dm3] λ̅0 = dimensionless zeroth moment of the growing radicals λ1 = first moment of the growing radicals [mol/dm3] λ̅1 = dimensionless first moment of the growing radicals λ2 = second moment of the growing radicals [mol/dm3] λ̅2 = dimensionless second moment of the growing radicals μ0 = zeroth moment for the dead polymer [mol/dm3] μ̅0 = dimensionless zeroth moment for the dead polymer μ1 = first moments for the dead polymer [mol/dm3] μ̅1 = dimensionless first moment for the dead polymer μ2 = second moments for the dead polymer [mol/dm3] μ̅2 = dimensionless second moment for the dead polymer ρ = density [kg/m3] τ = dimensionless operating time Φ = pay-off matrix in NBI technique Subscripts

a = air i = position in x-direction j = position in z-direction k, n = number of collocation point m = number of finite element tr = transition time rxn = reaction Superscripts

cs = compromise solution d = desired value 0 = starting value L = lower bound so = single optimization sop = simple operating policies U = upper bound u = utopia point

■

REFERENCES

(1) Tsoukas, A.; Tirrell, M.; Stephanopoulos, G. Multiobjective Dynamic Optimization of Semibatch Copolymerization Reactors. Chem. Eng. Sci. 1982, 7, 1785−1795. (2) Karjust, K.; Kuttner, R.; Pohlak, M.; Majak, J. Multi-criteria optimization of large composite parts. Compos. Struct. 2010, 92, 2146− 2152. (3) Alonso, A. A.; Banga, J. R.; Sendin, J. O.; Otero-Muras, I. Improved Optimization Methods for the Multiobjective Design of Bioprocesses. Ind. Eng. Chem. Res. 2006, 45, 8594−8603. (4) Weuster-Botz, D.; Torres-Darias, N.; Franco-Lara, E.; Link, H.; Vera, J. Multi-objective steady state optimization of biochemical reaction networks using a constrained genetic algorithm. Comput. Chem. Eng. 2008, 32, 1707−1713. 14364



Article

(27) Vivaldo-Lima, A.; Hamielec, A. E.; Wood, P. E. AutoAcceleration Effect in Free Radical Polymerization. A Comparison of the CCS and MH Models. Polym. React. Eng. 1994, 2, 17. (28) Mourikas, G.; Kiparissides, C.; Seferlis, P.; Morris, A. J. Online optimizing control of molecular weight properties in batch free-radical polymerization reactors. Ind. Eng. Chem. Res. 2002, 41, 6120. (29) Marten, F. L.; Hamielec, A. E.High-conversion diffusion-controlled polymerization; Henderson, H. N., Bouton, T. C., Eds.; ACS Symposium Series, No. 104; American Chemical Society: Washington, DC, 1979. (30) Achilias, D. S.; Kiparissides, C. Development of a General Mathematical Framework for Modeling Diffusion-Controlled Free Radical Polymerization Reactions. Macromolecules 1992, 25, 3739− 3750. (31) Bird, R. B.; Stewart, W. E.; Lighfoot, E. N.Transport Phenomena, 2nd ed.; John Wiley: New York, 2002. (32) Nunes, R. W.; Martin, J. R.; Johnson, J. F. Influence of Molecular Weight and Molecular Weight Distribution on Mechanical Properties of Polymers. J. Polym. Sci., Part A: Polym. Chem. 1982, 22, 205. (33) Takamatsu, T.; Shioya, S.; Okada, Y. Molecular Weight Distribution Control in a Batch polymerization. Ind. Eng. Chem. Res. 1988, 27, 93−99. (34) Schiesser, W. E.The Numerical Method of Lines. Integration of Partial Differential Equations; Academic Press: New York, 1991. (35) Kameswaran, S.; Biegler, L. T. Simultaneous Dynamic Optimization Strategies: Recent Advances and Challenges. Comput. Chem. Eng. 2006, 30, 1560−1575. (36) Ascher, L. R.; abd Petzold, U. M.Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations; SIAM: Philadelphia, PA, 1998. (37) Wachter, A.; Biegler, L. T.Algorithm for Large-Scale Nonlinear Programming; Technical report, CAPD Technical Report B-00-06; Carnegie-Mellon University: Pittsburgh, PA, 2000. (38) Gay, D. M.; Fourer, R.; Kernighan, B. W.AMPL A Modeling Language for Mathematical Programming, 2nd ed.; Thomson, Brooks/ Cole: Independence, KY, 2003. (39) Rivera-Toledo, M.; Flores-Tlacuahuac, A. A Multiobjective Dynamic Optimization Approach for a Methyl-Methacrylate Plastic Sheet Reactor. Macromol. React. Eng. 2014, 8 (4), 358−373. (40) Baillagou, P. E.; Soong, D. S. MolecularWeight Distribution of Products of Free Radical Nonisothermal Polymerization with Gel Effect Simulation of Poly(Methyl methacrylate). Chem. Eng. Sci. 1985, 40 (1), 87−104.

14365


Multiobjective Dynamic Optimization of the Cell ... - ACS Publications

Recommend Documents