Dynamic optimization of batch free radical polymerization with conditional modeling formulation through the adaptive smoothing strategy

Accepted Manuscript

Dynamic Optimization of Batch Free Radical Polymerization with Conditional Modeling Formulation through the Adaptive Smoothing Strategy Jing Kong , Xi Chen PII: DOI: Reference:

S0098-1354(18)30516-7 https://doi.org/10.1016/j.compchemeng.2018.09.023 CACE 6237

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

24 May 2018 12 August 2018 22 September 2018

Please cite this article as: Jing Kong , Xi Chen , Dynamic Optimization of Batch Free Radical Polymerization with Conditional Modeling Formulation through the Adaptive Smoothing Strategy, Computers and Chemical Engineering (2018), doi: https://doi.org/10.1016/j.compchemeng.2018.09.023

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights: Formulated the conditional model for diffusion-controlled FRP



Developed smoothing methods for gel, glass and cage effects



Proposed adaptive smoothing strategy for dynamic optimization of FRP



Presented numerical case studies for efficiency demonstration

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT

Dynamic Optimization of Batch Free Radical Polymerization with Conditional Modeling Formulation through the Adaptive Smoothing Strategy Jing Konga, Xi Chena,b 

a: State Key Laboratory of Industrial Control Technology, College of Control Science

CR IP T

and Engineering, Zhejiang University, Hangzhou, Zhejiang, 310027, China b: National Center for International Research on Quality-targeted Process

AN US

Optimization and Control

Abstract

M

With the increase of conversion and viscosity, a batch free radical polymerization (FRP) may encounter complicated phenomena, including the gel, glass, and cage effects. These effects

ED

normally lead to a conditional modeling formulation described by a set of differential-algebraic

PT

equations (DAEs). Discontinuous gradient exists at model switches triggered by occurrence of the effects, which cause difficulty for dynamic optimization on the batch process. To deal with the

CE

discontinuity, interior-point smooth approximations are proposed to describe the model switches at

AC

the occurrence of the effects. An adaptive smoothing strategy, which can efficiently balance the efficiency of optimization and the accuracy of the approximated model, is proposed for the dynamic optimization. To show the effectiveness of this smoothing strategy, case studies for optimal operating profiles in a polymerization process of methyl methacrylate are conducted. Good computational performance is achieved with the proposed strategy.



Correspondence concerning this article should be addressed to Xi Chen at [email protected].

2

ACCEPTED MANUSCRIPT

Key words: free radical polymerization; conditional modeling; adaptive smoothing; dynamic optimization. Introduction The reaction mechanism of a batch free radical polymerization (FRP) consists of initiation,

CR IP T

propagation, termination, and chain transfer (Xie et al., 1994). Besides these chemical reactions, three diffusion-controlled processes are also involved at high conversion, including the gel, glass, and cage effects, which are related to the termination, propagation, and initiation reactions,

AN US

respectively. Numerous attempts have been made to describe these effects. Marten and Hamielec (1979, 1982) proposed a kinetic conditional model (MH model) based on the free volume theory and demonstrated accuracy of the model on the polymerization of methyl methacrylate (PMMA)

M

and styrene. Stickler et al. (1984) proposed the Stickler-Panke-Hamielec (SPH) model which made a modification on the termination rate expression of the MH model. Tefera et al. (1997)

ED

presented a selection strategy of different FRP models and showed that the parameters in the MH

PT

model and SPH model are robust against changes in the initiator concentration. Vivaldo-Lima et al. (1994) added a description of cage effect to the MH model to get rid of the poor predictions of

CE

weight average chain length at high conversions. The above-mentioned MH models described the

AC

diffusion-controlled batch FRP in a conditional modeling framework. With a first principle model, dynamic optimization can be done for a batch polymerization process to seek for the optimal profiles of operation variables. In these dynamic optimization problems, the objective could be economic, such as maximization of monomer conversion. The constraints are usually sets of differential algebraic equations (DAEs) of the process model and requirements on product qualities. A number of applications have been conducted for polymerizations. Zavala et

3

ACCEPTED MANUSCRIPT

al. (2005) optimized polyurethane copolymerization in semi-batch reactor to achieve high-molecular-weight values for a fixed reaction time using cooling water, steam flowrate, addition of 1,4-butanediol and diamine. Nie et al. (2013) optimized the active catalyst distribution in a catalytic packed bed reactor of ortho-xylene oxidation to achieve higher product rate. Salhi et

CR IP T

al. (2004) studied the temperature profile in a batch copolymerization reactor of styrene and ɑ-methylstyrene to reduce the reaction time without considering gelation. Rivera-Toledo et al. (2009) researched on the PMMA reaction process for plastic sheet production and used warming

AN US

air temperature to optimize homogeneous and plastic sheet properties. Lin et al. (2010) proposed a mathematical model of a seeded suspension polymerization and studied an optimum feeding policy to minimize reaction time with product quality requirement. Gentric et al. (1999) discussed

M

the optimal temperature profile for a batch emulsion polymer reactor to minimize the reaction time and used a non-linear geometric controller to track the temperature profile. Among the numerous

ED

reports on dynamic optimization for polymerization, there are few papers on FRP with the

PT

non-smooth conditional models. At the switch points of these conditional models, state variables are continuous, but the gradients are discontinuous. This phenomenon significantly increases the

CE

difficulty in dynamic optimization, making it difficult to converge.

AC

The batch FRP can be treated as a kind of hybrid system. Optimization of a hybrid system with a conditional model can be conducted by several methods, which can be categorized into six types. The first type of method to deal with such a problem is derivative-free optimization (Audet and Dennis, 2002). These can avoid the difficulty in deriving accurate gradient information, but most of them can be inefficient compared to gradient-based methods. The second type of method is the extension of classical sensitivity theory. Galán et al. (1999) and

4

ACCEPTED MANUSCRIPT

Hiskens and Pai (2000) defined the parametric sensitivity trajectories of hybrid systems and developed sufficient conditions for the existence and uniqueness. Maly and Petzold (1996) and Feehery et al. (1997) exploited the special structure of extended parametric sensitivity. Based on these studies, Galán and Barton (1998) classified hybrid systems into three classes and pointed out

CR IP T

that the first two with fixed mode sequence could be solved by gradient-based methods with extended parametric sensitivity but the third class with changed mode sequence cannot be solved by the same way. Lee (2006) solved two special cases of the third class using branch-and-cut (BC)

AN US

algorithm and control parameterization enhancing transform (CPET), respectively. This extended sensitivity theory can be more accurate and efficient than the derivative-free methods; but not all problems can be solved with it and need special solution tools.

M

The third type of method is based on the generalized derivative of non-smooth functions. Clarke (1990) extended the concept of the derivative to non-smooth functions and developed a

ED

non-smooth analysis system. Based on this, Sussmann (1999) has developed maximum principles

PT

that are analogous to Pontryagin’s Maximum Principle for hybrid systems whose mode sequence does not change with manipulated variables. Piccoli (1999) studied the necessary conditions of

CE

optimality for special cases of hybrid systems with changed mode sequence. Using the necessary

AC

conditions, Caines and Shaikh (2006) developed a two-stage approach for dynamic optimization problems of hybrid system. This method is systematic and has great theoretical value, but special simulation and optimization solution tools based on this theory need to be developed. The fourth type of methods is mathematical programs with equilibrium constraints (MPEC). Baumrucker and Biegler (2009) introduced MPEC to replace conditional modeling, in which real-valued variables satisfying complementarity conditions were used. To use MPEC, dynamic

5

ACCEPTED MANUSCRIPT

optimization problems need to be solved by the simultaneous method, in which the time horizon is divided into several intervals called finite elements. The mode of system is determined by the complementarity conditions at the ends of these intervals. For the simultaneous method, the number of intervals is difficult to determine and may need to be rather large to obtain accurate

CR IP T

solutions, and the convergence of the simultaneous method cannot be guaranteed. Thus, introducing complementarity constraints to the simultaneous method can be a huge challenge.

The fifth type of methods is the mixed-integer program (MIP) (Till et al., 2004). In MIP method,

AN US

by introducing binary (and continuous) auxiliary variables and discretizing the time horizon, discrete part of hybrid system is transformed into equations set. The large numbers of the equations and variables always lead to complex optimization problem.

M

The sixth type of methods is the smoothing method. Zang (1980) introduced a smooth approximation which smoothed out the non-smooth problems only in their neighborhood. After

ED

that, Smale (1986), Chen and Harker (1993), and Kanzow (1996) proposed other smoothing

PT

functions. Chen and Mangasarian (1996) unified the smoothing functions by a family of smoothing functions. Gopal and Biegler (1999) introduced this kind of smoothing functions to

CE

solve complementarity problems in chemical process. After replacing the non-smooth model by a

AC

smoothing approximation, this problem can be solved by efficient classical solvers for smooth simulation and optimization. Though a number of efforts have been proposed to deal with non-smooth models in process optimization, as we know none of them have been designed for the batch FRP. In this project, the smoothing method with an adaptive strategy is proposed for the dynamic optimization of FRP with gel, glass, and cage effects. In the following sections, we first describe the reaction

6

ACCEPTED MANUSCRIPT

mechanism of FRP and address these phenomena that occur during the reaction process. Then, we present a dynamic optimization to seek the optimal operating profiles. To deal with the non-smoothness caused by these three effects, the adaptive smoothing strategy is proposed. The performance of the proposed method is validated through case studies.

CR IP T

2 Modeling and simulation of batch FRP with gel, glass, and cage effects

2.1 Reaction mechanism

AN US

As shown in Table 1, the FRP process can be described by the following reactions: initiator decomposition, chain initiation, propagation, termination by combination and disproportion, and chain transfer to monomer and polymer. According to the reaction mechanism, the mass balance equations of a well-mixed batch reactor are as follows:

M

For initiator 𝐼2 : 𝐼2

PT

For initial radical 𝐼 • :

(1)

ED

𝐼2 𝑡

𝐼• 𝑡

𝐼• 𝑀

(2)

CE

𝐼2

For monomer 𝑀:

AC

𝑀 𝑡

𝐼• 𝑀

𝑝

𝑀

(3)

Where

(4)

𝑝

The jth moment of free radicals can be expressed as: (

) 𝑡

𝐼• 𝑀

𝐼2

𝑀(

)

7

(

)

(5)

ACCEPTED MANUSCRIPT

𝑝𝑀 ∑

where (6)

𝑜

(7)

The jth moment of dead polymer species is expressed as:

(

∑∞ 𝑛 𝑅𝑛• and 𝑛

where

)

∑

(9)

AN US

𝑀

𝑡

CR IP T

(8)

𝑝

∑∞ 𝑛 𝑃𝑛 are the jth moments of the chain length distribution 𝑛

of free radicals and dead polymers, respectively.

𝑀

M

For a batch process, the monomer conversion can be expressed as: 𝑀

ED

𝑀

(10)

And the dynamic variation of reaction mixture volume is:

PT

(11)

𝑀

AC

𝑀𝑛

CE

The number and weight average molecular weights are described as:

𝑀

𝑀

2

(12)

2

(13)

Equations (1) to (13) form the frame model described in differential and algebraic equations (DAEs) with reaction kinetic coefficients. In an FRP, the system viscosity keeps growing with the reaction proceeding. Therefore, diffusion effects must be taken into consideration, which make the reaction kinetic coefficients time-varying in the model. The details are elaborated as follows.

8

ACCEPTED MANUSCRIPT

2.2 Gel effect Under normal circumstances, at high conversion the concentrations of monomer and initiator decrease, then the conversion rate slows down. For an FRP, however, the conversion rate starts to increase when the monomer conversion goes to 15%~20% (Marten and Hamielec, 1979). This

CR IP T

phenomenon is the so-called gel effect, which can be explained by diffusion control. With increasing viscosity, long-chain radical becomes more difficult to move and diffuse. Thus, the termination reaction is restrained and slows down. The termination rate

can drop 2 to 3 orders

AN US

of magnitude when the conversion goes to 40%~50%, resulting in the increase of the conversion rate.

Let K=𝑀𝑚 𝑒𝑥𝑝(𝐴/ 𝑓 ), and the critical value of the gel effect (Stickler et al., 1984) is denoted as

M

𝐾3 . At the beginning of the reaction, the termination rate is defined as the processing of the reaction, 𝐾 increases. When 𝐾

becomes diffusion controlled (Marten and Hamielec, 1979; Scorah et al.,

PT

2006): /

𝑀

where 𝑀

𝑛

/𝑀

and

𝑒𝑥𝑝[ 𝐴( /

where

{

𝐾 is first time satisfied.

5

𝑓

𝑝

𝑝( +𝜀 𝑋 +𝜀𝑋

𝐾 𝑀

/

𝑓

𝑓

)]

(14)

are the weight-average molecular weight and free volume fraction at the

𝑓

CE

AC

time when 𝐾3

𝐾3 , the gel effect is triggered, after which,

ED

the termination rate

, if 𝐾 < 𝐾3 . With

𝑝) 𝑝

𝑓

(

and

𝑝.

is free volume fraction with the general form as: )

(15)

Hence the gel effect model is described as follows:

𝑀𝑚 𝑒𝑥𝑝(𝐴/ 𝑓 ) < 𝐾3 /𝑀

𝑛

𝑒𝑥𝑝[ 𝐴( /

𝑓

/

𝑓

)]

𝐾

𝑀𝑚 𝑒𝑥𝑝(𝐴/ 𝑓 )

𝐾3

(16)

Among these equations, 𝐴, 𝑛, 𝑚 are general model parameters, which can be found the value in

9

ACCEPTED MANUSCRIPT

, 𝐾3 ,

Table 2;

𝑝,

,

𝑝,

are specific model parameters, which have different

values for different systems.

2.3 Glass effect

CR IP T

When conversion continues to increase, the reaction system gradually becomes solid-like. In this state, polymer performs like glass, so this phenomenon is called the glass effect. When this occurs, diffusion of monomer is restrained. As a result, the chain propagation rate dramatically. This effect happens when the free volume 𝑓

2.

Before the occurrence of this effect, propagation rate

the initial value of 𝑝

𝑝

𝑒𝑥𝑝[

𝑝.

( /

Then,

𝑝

/

𝑓

𝑝

equals to

𝑝

, which is

varies as (Marten and Hamielec, 1979; Stickler et al., 1984): 2 )]

𝑓

starts to decrease

is below the critical value of the

AN US

monomer

𝑓

𝑝

(17)

𝑝

𝑝

𝑝

𝑝

𝑓

𝑒𝑥𝑝[

𝑓

( /

2 𝑓

/

𝑓

2 )]

ED

,

M

Therefore, the model of glass effect has the conditional formulation:

𝑓

<

𝑓

(18) 2

PT

In these equations, B is general model parameter referred to Table 2; kp0, Vfcr2 are specific

CE

parameter referred to Table 3.

AC

2.4 Cage effect

After the decomposition of initiator, only some of initial radicals react with monomer, while the others fail to diffuse from the cage made by reactants and are wasted by side reactions, e.g., transfer to initiator. Hence, the cage effect reduces the efficiency of the initiator expected to happen when the critical free volume of initiator,

𝑓

3,

. Cage effect is

is reached0 (Scorah et al.,

2006; Vivaldo-Lima and Hamielec, 1994). Hence the model of cage effect can be described as

10

ACCEPTED MANUSCRIPT

follows: 𝑓

{

𝑓

𝑒𝑥𝑝 *

(

𝑓

3

𝑓

3

)+

𝑓

<

𝑓

(19)

3

In this equation, C, Vfcr3, f0 are specific parameters, which can be found in Table 3.

CR IP T

2.5 Simulation

In this article, we use the bulk homopolymerization of methyl methacrylate (PMMA), which is an important polymer for its good properties (Yan et al., 2013), initiated by 2,2’-azoisobutyronitrile

AN US

(AIBN) as an example. The mathematical model including Equations (1) to (19) can be applied for the PMMA and styrene. The values of general model parameters for all monomer-initiator systems (Marten and Hamielec, 1979) are shown in Table 2. The numerical values of specific model

M

parameters describing the FRP of MMA initiated by AIBN (Achilias and Kiparissides, 1992; Chen

ED

et al., 2009; Marten and Hamielec, 1979; Matyjaszewski and Davis, 2002; Stickler et al., 1984; Scorah et al., 2006) are listed in Table 3. To observe how the gel, glass, and cage effects influence

PT

the reaction, a simulation is conducted at the operation conditions listed in Table 4. Figure 1 shows

CE

the changes of the free volume, conversion, number molecular weight (Mn) and weight molecular weight (Mw) during the reaction process. The activations of the gel effect, glass effect and cage

AC

effect are marked using blue, black and red dots, respectively. Figure 2 illustrates the change of termination rate coefficient. Intersection of the dashed line and dash-dotted line indicates that

starts to decrease at this point, i.e., the start of gel effect, which

in this example is at 15.61 min. At this conversion, viscosity is not high enough to influence the diffusion of the monomer. Hence, the chain propagation rate makes the polymerization rate increase (in proportion to 11

𝑝/

𝑝 ⁄2

does not change much, which

). As the diffusion of free radicals

ACCEPTED MANUSCRIPT

is influenced by the gel effect, the radical life is prolonged, which increases average molecular weight. As shown in Figure 1, the conversion increases quickly after this point. Consequently, the average molecular weight also increases, especially when the conversion increase rate is high. This process lasts until the occurrence of the next effect at 37.07 min. During this process, the

average molecular weight increase from 5

5

, respectively.

Figure 3 shows the change of

𝑝

5

to 7 3

4

and 9 7

to

with the glass effect. Intersection of the dashed line and

AN US

6

94

CR IP T

conversion increase from 0.195 to 0.732, the weight average molecular weight and number

dash-dotted line indicates the start of the glass effect. At that point (37.07 min),

𝑝

reduces

sharply to a small value, and discontinuous gradient occurs at this jump point. As a result, the rate 𝑝/

⁄2

decreases.

M

of polymerization proportional to

As seen in Figure 4, the intersection of the dash line and dash-dotted line indicates the start of the

ED

cage effect (37.24 min), which is shortly after the start of glass effect. After this point,

reduces the rate of

PT

sharply to a small value. As shown in Equations (2) and (3), the decrease of

reduces

chain initiation, and the reaction process slows down. Together with the glass effect, this effect

CE

slows down the increase of the conversion rate and the average molecular weight, as shown in

AC

Figure 1. From the start of the glass effect to the end of this reaction, the conversion change from 0.732 to 0.872, and weight average molecular weight and number average molecular weight change from 7 3

5

to 9 3

5

and

6

5

to

75

5

, which are rather slow

compared with the gel effect. The previous description shows that the batch FRP is event-trigger system with several effects. They can be described with conditional modelling. The occurrence of the effects obviously

12

ACCEPTED MANUSCRIPT

changes the system features. Discontinuity will occur at the switching points, which may cause trouble for optimization.

3 Dynamic optimizations with adaptive smoothing strategy

CR IP T

In this project, we seek the optimal operation profile to maximize the conversion of monomer or minimize the reaction time with specified product quality requirement. The control variables are assumed to be piecewise linear functions of time over a specified number of time intervals. The general mathematical formulation of this optimization problem could be expressed in the

AN US

following form: 𝑚 𝑥/𝑚 𝑛

𝑀 𝑀𝑛

𝐴 𝑚𝑜 𝑒 𝑒𝑞𝑢 𝑡 𝑜𝑛𝑠 𝑞𝑠 𝑢𝑚 𝑛 𝑢 𝑢𝑚 𝑟 𝑀 (𝑡𝑓 ) 𝑀 𝑟 𝑀𝑛 𝑡𝑓 𝑀𝑛 𝑛

9

𝑟 𝑟

(20)

M

𝑠𝑡

ED

where n denotes the number of stages, Mw,target and Mn,target means the target product quality, r

PT

represents the acceptable constraint range. The sequential optimization strategy is adopted for this problem. The DAE solver is DASPK, and the NLP solver is SNOPT.

CE

The optimization problem is very difficult to solve. For most initial values, the solver cannot

AC

obtain any solution. These are caused by the discontinuity of gradients at the switch points in the conditional model. Discontinuous gradients make continuous-time integration algorithms either fail or become inefficient when trying to integrate over these discontinuities; and the convergence property of classical optimization algorithms is no longer guaranteed, which brings numerical difficulties in the computation. To solve this problem, a smoothing strategy is introduced as follows.

13

ACCEPTED MANUSCRIPT

3.1 Smoothing functions Chen and Mangasarian (1996) proposed a class of smoothing functions with some good properties. Conditional problems can be formulated by using the fundamental function 𝑥

be written as, 𝜎 𝑥

∫

∞

+

𝛿 𝑦

∫

∞

𝜎 𝑦

𝑦, where 𝜎 ∙ is the step function. Similarly, it can

𝑦, where 𝛿 𝑥 is the Dirac delta function, which has the same

properties as probability density functions in that: +∞

𝛿 𝑥

∫

𝛿 𝑦

𝑚 𝑥{𝑥 }.

CR IP T

For this function, we have 𝑥

+

(21)

𝑦

∞

AN US

Therefore, using a piecewise-continuous probability density function

𝑥 to smooth the Dirac

delta function can be reasonable. Then parameterize the density function: 𝑡̂ 𝑥

/

𝑥/

is a positive parameter. Then we have:

𝑠̂ 𝑥

∫ 𝑡̂ 𝑥

𝑡

𝜎 𝑥

𝑝̂ 𝑥

∫ 𝑠̂ 𝑥

𝑡

𝑥

(24)

+

PT

∞

(23)

ED

∞

M

where

(22)

As

approaches zero, smoothing functions approach their non-smooth counterparts. This

CE

smooth formulation has the key theoretical properties as follows.

AC

Proposition 1 (Chen and Mangasarian, 1996): Let

𝑥 be a probability density function and 𝑡̂ 𝑥

parameter. Let (A1)

/

𝑥/

, where

is a positive

𝑥 satisfy the following:

𝑥 is piecewise continuous with a finite number of pieces and satisfies Equation (21).

(A2) E[|𝑥|]

∫

+∞ |𝑥| ∞

𝑥

𝑥<

∞

Then

14

ACCEPTED MANUSCRIPT

(P1) 𝑝̂ 𝑥

is continuously differentiable.

(P2)

𝑝̂ 𝑥

2

𝑥

, where

+

+∞

∫ |𝑥|

𝑥

𝑥

∞

2

𝑚 𝑥 ,∫

𝑥

𝑥

𝑥

-

(25)

∞

Now we obtain a class of smooth plus functions, and the smooth function 𝑝̂ 𝑥

goes to 0. In this article, we use a specific smoothing function,

CR IP T

plus function as the parameter

approaches the

Interior-Point smooth approximation, to smooth the three effects, where:

𝑥2

(26)

3/2

AN US

𝑥

The fundamental smoothing function is as following:

𝑥

/

and for this function 𝑝̂ 𝑥

𝑥2

𝑥

𝑝̂ 𝑥

𝑚 𝑥

𝑥

⁄2

2

and

(27)

. Hence the approximation error can be evaluated:

2

M

𝑚 𝑥

/

(28)

ED

In the following sections, we use this smoothing function to approximate the three effects.

PT

3.2 Smoothing of three effects

CE

3.2.1 Smoothing of glass and cage effects The conditional model of the glass and cage effects are described in Equations (18) and (19).

AC

These two effects have similar form. The smoothing process is discussed using glass effect as an example. To smooth this conditional equation, we need re-express them as max-formula. From section 2.3, when when

𝑓

<

𝑓

2

,

𝑝

𝑓

𝑓

𝑝

𝑒𝑥𝑝[

2,

𝑝

( /

𝑝

𝑓

, which equals to /

re-expressed as:

15

𝑓

2 )]

𝑝

𝑝

𝑒𝑥𝑝[

]. And

. Hence the glass effect can be

ACCEPTED MANUSCRIPT

𝑝

𝑝

𝑒𝑥𝑝[

𝑚 𝑥(

/

/

𝑓

2 )]

𝑓

(29)

In the same way, the cage effect can be re-formulated as: 𝑒𝑥𝑝[

𝑚 𝑥(

/

/

𝑓

𝑓

3 )]

(30)

Using the interior point function (Equation (27)) and introducing smoothing parameters 𝜉𝑘𝑝 and

CR IP T

𝜉𝑓 , the max-formulas of the glass and cage effects can be approximated as: Glass effect: 𝑝

𝑝

𝑒𝑥𝑝

(31)

𝑘𝑝

where /

𝑓

𝑓

2

√ /

/

𝑓

𝑓

2

Cage effect: 𝑒𝑥𝑝 where 𝑓

/

𝑓

3

√ /

𝑓

/

𝑓

3

ED

/ 𝑓

𝜉𝑘𝑝

(32)

(33)

M

𝑓

2

AN US

/ 𝑘𝑝

2

𝜉𝑓

(34)

PT

3.2.2 Smoothing of gel effect

CE

The conditional model of gel effect is described in Equation (16). To smooth this effect, we need re-express it as max-formula. Compared with the glass and cage effects, the max-formula model

AC

of gel effect is not that obvious. To obtain the max-formula model, we first introduce a binary variable 𝑢

{

} to indicate the activation of the gel effect, such that when 𝑢

and when 𝑢

𝑀

/𝑀

𝑛

𝑒𝑥𝑝[ 𝐴( /

/

𝑓

𝑓

,

)]. Then the conditional model

can be expressed as: 𝑢

𝑢{

𝑀

/𝑀

𝑛

𝑒𝑥𝑝[ 𝐴( /

Then, introduce two intermediate variables,

and 16

𝑓

/

𝑓

)]}

, to denote the effect condition:

(35)

ACCEPTED MANUSCRIPT

𝑚 𝑥

𝐾3

𝐾

(36)

𝑚 𝑥

𝐾

𝐾3

(37)

Hence when 𝐾 < 𝐾3 (before the gel effect, 𝑢 the gel effect, 𝑢 ,

𝐾

; when 𝐾

𝐾3 (after

𝐾3 . Therefore, we can obtain the relationship between 𝑢

as:

𝑢

CR IP T

and

𝐾

),

𝐾3

),

(38)

Under these definitions, we obtain the max-formula model of the gel effect: 𝑢{

𝑀

/𝑀

𝑛

𝑒𝑥𝑝[ 𝐴( /

𝑢 𝑚 𝑥 𝑚 𝑥

𝑓

/

𝑓

)]}

AN US

𝑢

𝐾3 𝐾 𝐾 𝐾3

(39)

using the interior point function: 𝑢{

𝑢

𝐾3

𝐾

√ 𝐾3

𝐾

2

𝐾3

2

PT

𝐾3

√ 𝐾

/𝑀

𝑛

𝑒𝑥𝑝[ 𝐴( /

𝜉𝑘

𝑓

/

𝑓

)]}

(40)

𝜉𝑘

CE

𝐾

𝑀

ED

𝑢

M

The smooth formulation for the gel effect is derived by introducing smoothing parameter 𝜉𝑘 , and

AC

3.3 Adaptive smoothing strategy for optimization According to section 3.1, the smoothing parameter can control the smoothness and accuracy of the function. Using this property, we design an adaptive smoothing strategy for dynamic optimization. 3.3.1 Unified smoothing of the three effects To smooth the conditional model, three smoothing parameters 𝜉𝑘 𝜉𝑘𝑝 and 𝜉𝑓 have been introduced. Taking the gel effect in Equation (40) as an example, if 𝜉𝑘 is too small in 17

ACCEPTED MANUSCRIPT

comparison to 𝐾

𝐾3 2 , the influence of 𝜉𝑘 will be neglected, hence the smoothing effect in

Equation (40) hardly works. Smoothing parameters for different effects should have different value ranges. But if we adjust these three smoothing parameters one by one, the computation efficiency will be low. Therefore, a unified smoothing parameter ξ is used to balance the

CR IP T

complexity of optimization and the accuracy of the approximate models. At the same time, we need to ensure that when the smoothing parameter 𝜉 changes, the smoothing errors of different effects are kept at similar orders of magnitude. For this reason, three smoothing parameter 𝜉𝑘𝑝 and 𝜉𝑓 are introduced as:

AN US

coefficients 𝜉𝑘 𝜉𝑘 𝜉𝑘 𝜉 𝜉 { 𝑘𝑝 𝜉𝑘𝑝 𝜉 𝜉𝑓 𝜉𝑓 𝜉

The values of three smoothing parameter coefficients 𝜉𝑘

(41)

𝜉𝑘𝑝 𝜉𝑓 in this article are set as 4.0e7,

M

60 and 400, respectively. Table 5 compares the end results of the simulations on the smooth model

ED

under different smoothing parameter settings with the original model. We can clearly see that the

PT

difference is increased with the increased smoothing parameter. Figure 5 further illustrates the comparative results of conversion and average molecular weight. We can conclude that the model

CE

becomes smoother and more biased with the increase of the smoothing parameter. In the dynamic

AC

optimization process, therefore, the smoothing parameter should be carefully tuned to balance the smoothness and the model accuracy. 3.3.2 Constraint slackness In the previous sections, we introduce the basic smoothing function and conclude that if the smoothing parameter is too small, the smooth model would approximate the original model accurately, but the smoothness of the function is decreased. On the contrary, if the smoothing

18

ACCEPTED MANUSCRIPT

parameter is too large, the optimization problem with a smoother model becomes easier to solve, but the model reliability will be affected due to the inaccuracy; particularly, meeting constraints will be difficult due to the model mismatch. Thus, besides the smoothing parameter, the acceptable constraint range r as shown in Equation (20) is also adjusted so that the smooth model can be

CR IP T

solved with a large feasible range. In this project, the original constraints are set at r = 5%. During the optimization process, the parameter is adjusted between [5%, 80%]. With the increase of r, a wider feasible region is applied, which makes the optimization easier to solve. During the

AN US

smoothing process, whenever the original constraints cannot be met, r is increased to make it feasible. When the smoothing parameter 𝜉 is reduced, r is also reduced. Eventually, it should be restored to the original setting, 5%.

M

3.3.3 Adaptive smoothing strategy

An adaptive strategy is proposed to tune the smoothing parameter and the acceptable constraint

ED

range to balance the optimization difficulty and accuracy. As a smoother model usually takes less

PT

computational time to converge, the solution process starts by a model with a large smoothing parameter and then gradually decrease it. Meanwhile, the acceptable constraint range is also

CE

adaptively changed to enlarge the feasible region of the optimization problem at a fixed smoothing

AC

parameter. The process will increase the acceptable constraint range when the smoothing parameter is too large to meet the constraints. With a decrease of the smoothing parameter, the acceptable constraint range is also decreased and eventually reaches the original setting. However, optimization of the original model should be avoided until the smooth model is close to it; otherwise, optimization of the original model tends to fail and the calculation cost for the original model can be rather long. Therefore, we introduce a so-called “test model” with a small smoothing

19

ACCEPTED MANUSCRIPT

parameter (1.0e-4) and the original acceptable constraint range (5%) to bridge the gap. Each time an optimization result is obtained, it will be passed to the test model as an initial guess. If optimization of the test model is successful, optimization on the original problem will be activated with a good initialization obtained from the test model. If either the test model or the original

CR IP T

model fails to converge during the solution process, the adaptive smoothing strategy will be restored and optimization on a new smooth model will be conducted. This process will be repeated until optimization on the original model succeeds, or, a failed end is reported when the process

AN US

fails to find any solution during the smoothing process. As shown in Figure 6, the flowchart of the adaptive strategy can be summarized as follows. Strategy steps:

Step 1. Attempt to optimize the initial problem with 𝜉

𝜉 and 𝑟

step 2.

ED

M

(a) If a solution is found, record these two parameters, 𝜉𝑠

𝑟∙𝑝 , 𝑝

𝑠𝑠

𝜉, 𝑟𝑠

𝑠𝑠

𝑟, and go to

𝜉, 𝑟𝑠

𝑠𝑠

𝑟, and go to

.

PT

(b) Otherwise, loosen the constraints by 𝑟

𝑟.

If 𝑟 > 𝑟𝑚 , go to the failed end.

2.

Otherwise, go back to step 1.

CE

1.

𝜉

𝑠

AC

Step 2. Attempt to optimize the test problem with 𝜉

, 𝑟

(a) If a solution is found, record these two parameters, 𝜉𝑠

𝑟

𝑠

𝑠𝑠

.

step 3.

(b) Otherwise, go to step 4. Step 3. Attempt to optimize the original conditional model. (a) If a solution is found, accept the optimization result, go to the successful end.

20

ACCEPTED MANUSCRIPT

(b) Otherwise, simulate the original conditional model. Test the accuracy of the smooth problem. 1.

If the accuracy requirement is satisfied, accept the result and go to the successful end. Otherwise, go to step 4.

Step 4. Attempt to optimize the problem with 𝜉

𝜉𝑠

𝑠𝑠 /𝑝𝜉 ,

(a) If a solution is found, record these two parameters, 𝜉𝑠 If 𝜉

𝜉

and 𝑟

2.

Otherwise, go to step 5.

𝑠

(b) Otherwise, let 𝑝𝜉

𝑟

𝑠

, go to step 3.

𝑝𝜉 / .

If 𝑝𝜉 < 𝑝𝜉 𝑚 𝑛 , go to the failed end.

2.

Otherwise, go back to step 4.

M

1.

ED

Step 5. Attempt to optimize the problem with 𝜉

𝜉𝑠

𝑠𝑠 ,

𝑟

PT

(a) If a solution is found, record these two parameters, 𝜉𝑠 1.

If 𝜉

2.

Otherwise, go to step 2.

𝑠

CE

𝜉

AC

(b) Otherwise, let 𝑝

and 𝑟

𝑟

𝑟𝑠

𝑠𝑠

𝑠𝑠 .

𝜉, 𝑟𝑠

𝑠𝑠

𝑟.

𝑠𝑠

𝑟.

AN US

1.

𝑟

CR IP T

2.

𝑠

𝑟𝑠

𝑠𝑠 /𝑝

𝑠𝑠

.

𝜉, 𝑟𝑠

, go to step 3.

𝑝/ .

1.

If 𝑝 < 𝑝

2.

Otherwise, go back to step 5.

𝑚 𝑛,

go to step 2.

4 Results and discussion In this section, the bulk homopolymerization of methyl methacrylate (PMMA) initiated by 2,2’-azoisobutyronitrile (AIBN) is taken as an example. Three case studies are presented and 21

ACCEPTED MANUSCRIPT

discussed. In all cases, the PMMA model and smoothing strategy are coded in FORTRAN and SNOPT is used as the optimization solver. All programs are executed on a 3.30 GHz Intel Xeon CPU E3-1230 v3.

CR IP T

4.1 Case 1: Conversion Maximization with Changing Initiator Feed The first case study aims to determine the optimal profile of initiator feed to maximize monomer conversion with specified final weight average molecular weight and number average molecular weight in a given batch operation time (𝑡𝑓 ). Based on the control parameterization method, the

AN US

profile of the initiator feed is parameterized for different time stages. With increasing number of time stages, better optimal results can be obtained. The specified final average molecular weights are set according to the simulation results as listed in Table 6. After conducting the simulation and 5% tolerance on the final number/weight average molecular

ED

weights as the end point constraints.

M

obtaining the results, we set 𝑟

The mathematical formulation of this problem is shown in Equation (42):

𝑠 𝑡

𝐴 𝑚𝑜 𝑒 𝑒𝑞𝑢 𝑡 𝑜𝑛𝑠

AC

𝑀𝑛

𝑚

𝑞𝑠 𝑛

9 (42)

𝑟

𝑀 (𝑡𝑓 )

𝑀

𝑟

𝑟

𝑀𝑛 (𝑡𝑓 )

𝑀𝑛

𝑟

CE

𝑚𝑛

𝑀

𝑡𝑓

PT

𝑚 𝑥

In this case, we take 𝑡𝑓

4 𝑚 𝑛 as a batch reaction time limitation, and use 1, 5 and 10 stages

to solve this optimization problem respectively. As the optimal results rely heavily on initialization, each problem is calculated for 100 times with different randomly generated initializations. The results of directly optimizing the original model and using the smoothing strategy are shown in Table 7. In the table, the success number denotes how many successful optimization results are

22

ACCEPTED MANUSCRIPT

obtained among the 100-time calculation; and the total time counts the overall time cost of the 100 optimization problems. We can see that the optimization cannot converge with a single stage. The key reason is that the specified quality constraints cannot be met with a batch reaction time set as 40 min and only one-stage initiator profile. As shown in Table 6, the quality constraints were

CR IP T

obtained from the simulation for 100 min reaction at a constant initiator feed. However, when the number of profile stages is increased to 5, the problem can be solved, which means that increasing control variables helps to extend feasible region and thus meet the quality constraints. By using

AN US

the direct optimization on the original problem, most runs fail to find a solution owing to the difficulty in solving the conditional model. However, the results of using the smoothing strategy show that success rate can be greatly enhanced to over 95%, and the total calculation time is also

M

greatly decreased. Additionally, increasing the stage number results in the increasing complexity of the optimization problem. The direct optimization becomes difficult to converge. But the

ED

smoothing method keeps a very high success rate.

PT

To further analyze how the optimization results are affected by the manipulated variables, we extend the profile stage number to 20 and 40. The conversion results are listed in Table 8. We can

CE

see that with the increase of the profile stage number and thus the degrees of freedom, the

AC

maximized conversion rate is also slightly increased and eventually reaches to a limit. The dynamic results of key variables are illustrated in Figure 7, in which the solid line shows the results of the 5-stage profile, the dashed line shows the result of the 10-stage profile, the dotted line shows the result of the 20-stage profile, and the chain-dotted line shows the result of the 40-stage profile. Basically, the results follow similar trends, indicating the consistency of the optimal results. The first subplot in Figure 7 shows the change of initiator feed with time. The

23

ACCEPTED MANUSCRIPT

second figure represents the conversion over time. The third subplot shows the change of free volume over time, where the color solid dots on the lines represent the occurrence of the gel effect, glass effect and cage effect, respectively, in order of the time sequence. The fourth subplot is the change in number average molecular weight and weight average molecular weight with time.

CR IP T

From the first subplot in Figure 7, we can see that the most amounts of initiators tend to be added around the middle stage of the batch reaction, close to the occurrence of glass and cage effects as shown in the third subplot. Considering that the cage effect decreases the efficiency of initiator,

AN US

adding initiator can make up the influence of cage effect and benefits the maximization of conversion rate.

4.2 Case 2: Conversion Maximization with Changing Initiator Feed and Isothermal

M

Condition

ED

To increase the conversion, reaction temperature is added as a control variable in dynamic optimization problem. The mathematical formulation of this problem can be shown in Equation

PT

(43):

𝑡𝑓

𝐴 𝑚𝑜 𝑒 𝑒𝑞𝑢 𝑡 𝑜𝑛𝑠

CE

𝑠 𝑡

𝑚𝑛

𝑚

𝑚𝑛

AC

𝑀 𝑀𝑛

𝑞𝑠 𝑛

9 (43)

𝑚

𝑀 (𝑡𝑓 )

𝑀

𝑀𝑛 (𝑡𝑓 )

𝑀𝑛

Because the MH model was established under isothermal condition, the temperature profile is set at isothermal condition varying between 4 ℃ to 9 ℃. The other constraints of this problem are the same as the previous optimization. The reaction time is also set at 40 min. The numbers of stages are first set as 1, 5, and 10,

24

ACCEPTED MANUSCRIPT

respectively. The results are shown in Table 9. For the direct optimization method, the approximating success rates of these problems decreases with the number of stages increases. Especially for the 10-stage case, the success rate is only 7%. However, the smoothing strategy can raise the success rate to above 80%.

CR IP T

Additionally, the optimization is run with the profile stage number set as 20 and 40. The conversion results are show in Table 8, and the profile of the optimized initiator feed and temperature profile, conversion, free volume (indicating the three effects) and average molecular

AN US

weight are presented in Figure 8. As shown in Table 9, adding temperature as a new manipulated variable improves the monomer conversion. Besides, the conversion increases as the number of stages increases and eventually reaches a limit. The representation of Figure 8 is very similar to

M

Figure 7, except that the new optimization variable, temperature, is added into the first subplot. Moreover, this dynamic optimization problem becomes feasible for a 1-stage profile of the

ED

initiator feed with the reaction temperature optimized at 76 87 ℃. For other cases with more

PT

stages, the optimized temperature is at its upper bound, 9 ℃, which means that high temperature helps to raise the conversion. Like case 1, initiator tends to be added around the occurrence of

CE

glass and cage effects.

AC

4.3 Case 3: Reaction Time Minimization with Changing Initiator Feed and Isothermal Condition

Previously we use fixed reaction time (40min) to maximize monomer conversion and show the advantages of the smoothing strategy. In this section, we change the optimization formulation by minimizing the reaction time subject to 90% conversion rate. This formulation is also very practical in some cases. The optimization formulation is as follows: 25

ACCEPTED MANUSCRIPT

𝑡𝑓 𝑠 𝑡

𝐴 𝑚𝑜 𝑒 𝑒𝑞𝑢 𝑡 𝑜𝑛𝑠 𝑚𝑛

𝑚 𝑚𝑛

𝑀

9 𝑀 (𝑡𝑓 )

𝑀𝑛

𝑀𝑛 (𝑡𝑓 )

𝑞𝑠 𝑛

9 (44)

𝑚

𝑀 𝑀𝑛

CR IP T

The numbers of stages are set as 1, 5, 10, 20 and 40 respectively. The optimal results are presented in Table 10 and Figure 9. Because temperature reaches the upper bound of 90℃ in all of the optimal solutions, the temperature profiles are not drawn in Figure 9. Compared with the result in

AN US

section 4.2, the reaction time of the 1-stage problem is increased and the temperature reaches 90℃ to meet the constraints on the conversion rate and average molecular weight at the same time. With the increasing number of stages, the minimum reaction time becomes shorter. When the time intervals increased from 1 to 5, the reaction time decreases almost 40%, but when the number of

M

time intervals increased from 5 to 10 and 10 to 20, time saving becomes small. When the stages

ED

change from 20 to 40, the reaction time almost keeps the same. With increasing number of stages,

PT

the minimum reaction time will tend to a fixed value. Comparing the 40-stage results in Tables 8

CE

and 10 shows that we can save 9.95 min (24.88%) at a cost of only 0.0041 decrease of conversion.

4.4 Detailed illustration of the adaptive smoothing strategy

AC

The above case studies have shown the advantage of the proposed smoothing strategy. Hereafter, we will illustrate in detail how the smoothing strategy can work for dynamic optimization. Basically, a smooth model is easier to solve than the original model. Table 11 presents the comparison of how different values of the smoothing parameter and the acceptable constraint range affect the dynamic optimization formulated in Equation (43). To keep consistency of the comparison, 100 different initial guesses of the optimized variables are generated randomly, and 26

ACCEPTED MANUSCRIPT

the tests under each parameter settings are all calculated for 100 times with those 100 initializations. The following conclusions can be made from the data in Table 11. (1) A smoother model with a larger smoothing parameter takes less calculating time; (2) The smaller a smoothing parameter is, the better it approximates the original model; (3) The success rate of an optimization

CR IP T

problem depends on both the smoothness of the model and tightness of the constraints; (4) The average SNOPT iterations do not have a clear relation with the smoothing parameter; but there is a clear trend that the smoother the model is, the less time is required for one run of computation.

AN US

Based on these properties, we propose the smoothing strategy as mentioned in Section 3.3.3, which balances the smoothing parameter and acceptable constraint range at the same time. Next, we present the calculating process with the adaptive smoothing strategy by using the case study in

M

Section 4.2 with a 10-stage profile. The process of changing the smoothing parameter and the acceptable constraint range is shown in Figure 10. First, Step 1 in Section 3.3.3 is taken. The

𝑟

𝜉

𝑒

and acceptable constraint

% is optimized. It fails to find any feasible solution. We mark it as run 1 with

PT

range 𝑟

ED

initial problem with a large smoothing parameter 𝜉

an empty circle which is used to denote a failure solution as shown in Figure 10. The switch (b) is

CE

taken and the acceptable constraint range 𝑟 is automatically extended to 20%. As the acceptable 8 %, it goes back to Step 1. This procedure is repeated until

AC

constraint range is less than 𝑟𝑚

the 3rd run with the acceptable constraint range extended to 40%. This problem is successfully solved and represented by the solid circle as shown in the figure. Then the switch (a) in Step 1 is taken, and it goes to Step 2. In this step, a “test problem” with 𝜉

𝑒

4 and r = 5% is

calculated using the successful result of the 3rd run as the initial value. This optimization fails to find a feasible solution as denoted by the empty circle marked “4” in the figure. Then the process

27

ACCEPTED MANUSCRIPT

goes to Step 4. In this step, the smoothing parameter is decreased to

𝑒

, and the acceptable

constraint range stays at the value of the 3 rd run. This problem is solved successfully and represented by the solid circle marked “5”. Then the switch (a)-2 in Step 4 is taken, and it goes to the Step 5. A new optimization with the same smoothing parameter 𝜉

𝑒

and a tighter

CR IP T

acceptable constraint range r =10% is calculated. This problem fails and is represented by the empty circle marked “6” in the figure. Hence, the switch (b)-2 in Step 5 is taken. A new problem with the acceptable constraint range extended to 20% is solved successfully, and shown by the

AN US

solid circle marked “7”. After obtaining a successful optimization result, it goes to Step 2 and optimize the “test problem”. This time, the “test model” is solved successfully as represented by the solid circle marked “8” in the figure. Then the process moves to Step 3. In this step, the

M

original model is optimized successfully with the initial value provided by the 8 th run, thus the final result is obtained. As Figure 10 is drawn in logarithmic coordinates, it cannot represent the

ED

calculation of the original conditional model whose smoothing parameter is 0. Therefore, the last

PT

successful point is not represented in Figure 10. From this example, we can see that this adaptive smoothing strategy starts from the smooth model

CE

with a large smoothing parameter and decreases the smoothing parameter gradually, which can

AC

save a lot of calculating time. Changing the acceptable constraint range together with the smoothing parameter can guarantee successful results. With the decrease of the smoothing parameter and the acceptable constraint range, the problem approximates to the original optimization problem. In this way, better initial values are provided step by step which makes the original model more likely to converge. At the same time, introducing the “test problem” can avoid unnecessary calculation on the time-consuming original problem and thus save the

28

ACCEPTED MANUSCRIPT

calculating time. For these reasons, comparing with direct optimization on the original model, using the adaptive smoothing strategy can improve the success rate and decrease solution time at the same time.

CR IP T

5 Conclusions By taking the gel, glass, and cage effects in FRP into consideration, a conditional model is first established. Owing to the non-smooth nature of the three effects, the original unsmoothed dynamic optimization problem is difficult to converge. To overcome this difficulty, an adaptive

AN US

smoothing strategy is proposed by balancing the complexity of optimization and the accuracy of the approximated model. To show the effectiveness of this smoothing strategy, dynamic optimization for a PMMA process is studied. The numerical results show that this strategy solves

M

the dynamic optimization problem effectively with high rate of success. The results of dynamic

PT

Acknowledgements

ED

optimization also give advice on operating the PMMA reaction.

The authors gratefully acknowledge the financial support of National Key Research and

CE

Development Program (No. 2017YFE0106700), NSFC-Zhejiang Joint Fund for the Integration of

AC

Industrialization and Informatization (No. U1509209), and the Fundamental Research Funds for the Central Universities (No. 2018QNA5013).

Nomenclature 𝑡

time, min volume of reaction mixture, 𝐿

29

ACCEPTED MANUSCRIPT

initial volume of reaction mixture, 𝐿 volumetric contraction coefficient ∙ 𝑚𝑜

molecular weight of monomer,

𝐼2

concentration of initiator, 𝑚𝑜 ∙ 𝐿

𝐼•

concentration of primary radical, 𝑚𝑜 ∙ 𝐿

𝑀

concentration of monomer, 𝑚𝑜 ∙ 𝐿

𝑀

initial concentration of monomer, 𝑚𝑜 ∙ 𝐿

𝑅𝑛•

concentration of free radical with n repeat units, 𝑚𝑜 ∙ 𝐿

𝑃𝑛

concentration of dead polymer with n repeat units, 𝑚𝑜 ∙ 𝐿

AN US

CR IP T

𝑀

jth moments of the chain length distribution of free radicals, 𝑚𝑜 ∙ 𝐿

monomer conversion

M

jth moments of the chain length distribution of dead polymers, 𝑚𝑜 ∙ 𝐿

cumulative number average molecular weight,

𝑀

cumulative weight average molecular weight,

∙ 𝑚𝑜 ∙ 𝑚𝑜

PT

ED

𝑀𝑛

decomposition rate coefficient for initiator, refer to Table 3, 𝐿 ∙ 𝑚𝑜 rate coefficient for the initiation reaction, equals to

𝑝

CE

𝐿 ∙ 𝑚𝑜

,

∙𝑚 𝑛

rate coefficient of propagation reaction before the glass effect, in Table 3, 𝐿 ∙ 𝑚𝑜

AC

𝑝

𝑝

∙𝑚 𝑛

rate coefficient of propagation reaction, 𝐿 ∙ 𝑚𝑜

∙𝑚 𝑛

rate coefficient of chain transfer to monomer, 𝐿 ∙ 𝑚𝑜

∙𝑚 𝑛

rate coefficient for chain transfer to initiator, in this model, equals to 0, 𝐿 ∙ 𝑚𝑜 rate coefficient of chain transfer to polymer, 𝐿 ∙ 𝑚𝑜 rate coefficient of termination by combination, 𝐿 ∙ 𝑚𝑜

30

∙𝑚 𝑛

∙𝑚 𝑛 ∙𝑚 𝑛

∙𝑚 𝑛

ACCEPTED MANUSCRIPT

rate coefficient of termination by disproportion, 𝐿 ∙ 𝑚𝑜

∙𝑚 𝑛

rate coefficient of termination before the gel effect, 𝐿 ∙ 𝑚𝑜

∙𝑚 𝑛

rate coefficient of termination by combination and disproportion, 𝐿 ∙ 𝑚𝑜 𝑜

equals to

𝑘𝑡𝑑 𝑘𝑡

∙𝑚 𝑛

, is a specific parameter, refer to Table 3

initiator efficiency specific model parameter, refer to Table 3

𝑀

AN US

specific model parameter, refer to Table 3

CR IP T

initiator efficiency before the cage effect, refer to Table 3

weight average molecular weight where the gel effect starts, free volume where the gel effect starts

𝑓

free volume

M

𝑓

∙ 𝑚𝑜

the thermal expansion coefficient of polymer, refer to Table 3

𝑝

ED

the thermal expansion coefficient of monomer, refer to Table 3 glass transition temperature of polymer, refer to Table 3, 𝐾

PT

𝑝

glass transition temperature of monomer, refer to Table 3, 𝐾 volume fraction of polymer, refer to Table 3

CE

𝑝

AC

volume fraction of monomer, refer to Table 3

𝐾3

model parameter determining the onset of the gel effect, refer to Table 3

𝐾

variable determining the onset of the gel effect

𝑓

2

free volume where the glass effect starts, refer to Table 3

𝑓

3

free volume where the cage effect starts, refer to Table 3

31

ACCEPTED MANUSCRIPT

References: Achilias, D. S., Kiparissides, C., 1992. Development of a general mathematical framework for modeling diffusion-controlled free-radical polymerization reactions. Macromolecules 25(14), pp. 3739-3750.

CR IP T

Audet, C., Dennis, JR. J. E., 2002. Analysis of generalized pattern searches. SIAM Journal on Optimization 13(3), pp. 889-903.

Baumrucker, B. T., Biegler, L. T., 2009. MPEC strategies for optimization of a class of hybrid

AN US

dynamic systems. Journal of Process Control 19(8), pp. 1248-1256.

Clarke, F. H., 1990. Optimization and nonsmooth analysis, second ed. Society for Industrial and Applied Mathematics, Philadelphia.

M

Chen, B., Harker, P. T., 1993. A non-interior-point continuation method for linear complementarity problems. SIAM Journal on Matrix Analysis and Applications 14(4), pp. 1168-1190.

ED

Chen, C., Mangasarian, O. L., 1996. A class of smoothing functions for nonlinear and mixed

PT

complementarity problems. Computational Optimization and Applications 5(2), pp. 97-138. Caines, P. E., Shaikh, M. S., 2006. Optimality zone algorithms for hybrid systems: Efficient

CE

algorithms for optimal location and control computation, in: Hespanha J., Tiwari A. (Eds.),

AC

Hybrid Systems: Computation and Control. Springer, Berlin, Heidelberg, pp. 123-137. Chen, X., Chen, L., Feng, J., Yao, Z., Qian, J. X., 2009. Toward polymer product design. I. Dynamic optimization of average molecular weights and polydispersity index in batch free radical polymerization. Industrial & Engineering Chemistry Research 48(14), pp. 6739-6748. Feehery, W. F., Tolsma, J. E., Barton, P. I., 1997. Efficient sensitivity analysis of large-scale differential-algebraic systems. Applied Numerical Mathematics 25(1), pp. 41-54.

32

ACCEPTED MANUSCRIPT

Galán, S., Barton, P. I., 1998. Dynamic optimization of hybrid systems. Computers & Chemical Engineering 22, pp. 183-190. Galán, S., Feehery, W. F., Barton, P. I., 1999. Parametric sensitivity functions for hybrid discrete/continuous systems. Applied Numerical Mathematics 31(1), pp. 17-47.

CR IP T

Gentric, C., Pla, F., Latifi, M. A., Corriou, J. P., 1999. Optimization and non-linear control of a batch emulsion polymerization reactor. Chemical Engineering Journal 75(1), pp. 31-46.

Gopal, V., Biegler, L. T., 1999. Smoothing methods for complementarity problems in process

AN US

engineering. AIChE Journal 45(7), pp. 1535-1547.

Hiskens, I. A., Pai, M. A., 2000. Trajectory sensitivity analysis of hybrid systems. IEEE Transactions on Circuits and Systems Part I: Regular Papers 47(2), pp. 204-220.

M

Kanzow, C., 1996. Some noninterior continuation methods for linear complementarity problems. SIAM Journal on Matrix Analysis and Applications 17(4), pp. 851-868.

ED

Lee, C. K., 2006. Global optimization of hybrid systems. The department of Chemical

PT

Engineering, Massachusetts Institute of Technology. Lin, W., Biegler, L. T., Jacobson A M., 2010. Modeling and optimization of a seeded suspension

CE

polymerization process. Chemical Engineering Science 65(15), pp. 4350-4362.

AC

Li, Y., Chen, Z. H., Zeng, C. C., 2013. Poly (Methyl Methacrylate) (PMMA) nanocomposite foams, in: Mittal, V. (Ed.), Polymer Nanocomposite Foams. CRC Press, Boca Raton, pp. 1-29.

Marten, F. L., Hamielec, A. E., 1979. High conversion diffusion-controlled polymerization, in: Henderson and Bouton, Polymerization Reactors and Processes ACS Symposium Series. American Chemical Society, Washington, D.C., pp. 43-70. Marten, F. L., Hamielec, A. E., 1982. High‐ conversion diffusion‐ controlled polymerization of

33

ACCEPTED MANUSCRIPT

styrene. I. Journal of Applied Polymer Science 27(2), pp. 489-505. Maly, T., Petzold, L. R., 1996. Numerical methods and software for sensitivity analysis of differential-algebraic systems. Applied Numerical Mathematics 20(1), pp. 57-79. Matyjaszewski, K., Davis, T. P., 2002. Handbook of radical polymerization. Wiley-Interscience,

CR IP T

New York. Rivera-Toledo, M., Flores-Tlacuahuac, A., Vílchis-Ramírez, L., 2009. Dynamic optimization of the methylmethacrylate cell-cast process for plastic sheet production. AIChE journal 55(6), pp.

AN US

1464-1486.

Nie, Y., Witt, P. M., Agarwal, A., Biegler, L. T., 2013. Optimal Active Catalyst and Inert Distribution in Catalytic Packed Bed Reactors: ortho-Xylene Oxidation. Industrial &

M

Engineering Chemistry Research 52(44), pp. 15311-15320.

Piccoli B., 1999. Necessary conditions for hybrid optimization. The 38th IEEE Conference on

ED

Decision and Control 1, pp. 410-415.

PT

Stickler, M., Panke, D., Hamielec, A. E., 1984. Polymerization of methyl methacrylate up to high degrees of conversion: Experimental investigation of the diffusion‐controlled polymerization.

CE

Journal of Polymer Science: Polymer Chemistry Edition 22(9), pp. 2243-2253.

AC

Smale, S., 1986. Algorithms for solving equations. The International Congress of Mathematicians, pp. 172-195.

Sussmann, H. J., 1999. A maximum principle for hybrid optimal control problems. The 38th IEEE Conference on Decision and Control 1, pp. 425-430. Salhi, D., Daroux, M., Gentric, C., Corriou, J. P., Pla, F., Latifi, M. A., 2004. Optimal temperature-time programming in a batch copolymerization reactor. Industrial & engineering

34

ACCEPTED MANUSCRIPT

chemistry research 43(23), pp. 7392-7400. Scorah, M. J., Dhib, R., Penlidis, A., 2006. Modelling of free radical polymerization of styrene and methyl methacrylate by a tetrafunctional initiator. Chemical Engineering Science 61(15), pp. 4827-4859.

CR IP T

Tefera, N., Weickert, G., Westerterp, K. R., 1997. Modeling of free radical polymerization up to high conversion. I. A method for the selection of models by simultaneous parameter estimation. Journal of applied polymer science 63(12), pp. 1649-1661.

AN US

Till, J., Engell, S., Panek, S., Stursberg, O., 2004. Applied hybrid system optimization: An empirical investigation of complexity. Control Engineering Practice 12(10), pp. 1291-1303.

Vivaldo-Lima, E., Hamielec, A. E., Wood, P. E., 1994. Auto-acceleration effect in free radical

M

polymerization. A comparison of the CCS and MH models. Polymer reaction engineering 2(1-2), pp. 17-85.

ED

Xie, T. Y., McAuley, K. B., Hsu, J. C. C., Bacon, D. W., 1994. Gas phase ethylene polymerization:

PT

production processes, polymer properties, and reactor modeling. Industrial & engineering chemistry research 33, pp. 449-479.

CE

Zang, I., 1980. A smoothing-out technique for min-max optimization. Mathematical Programming

AC

19(1), pp. 61-77.

Zavala, V. M., Flores-Tlacuahuac, A., Vivaldo-Lima, E., 2005. Dynamic optimization of a semi-batch reactor for polyurethane production. Chemical Engineering Science 60(11), pp. 3061-3079.

35

ACCEPTED MANUSCRIPT

List of figures: Figure 1 Simulation Results of Free Volume, Conversion, Mn and Mw. Figure 2 Simulation Result of Termination Rate Coefficient Kt in Gel Effect. Figure 3 Simulation Result of Chain Propagation Rate Coefficient Kp in Glass Effect. Figure 4 Simulation Result of the Efficiency of Initiator f in Cage Effect.

Parameters. Figure 6 Flowchart of the Adaptive Smoothing Strategy.

CR IP T

Figure 5 Simulation Result of Original Model and Smooth Models under Different Smooth

Figure 7 Calculation Results of Conversion Maximization with Changing Initiator Feed.

Figure 8 Calculation Results of Conversion Maximization with Changing Initiator Feed and

AN US

Temperature.

Figure 9 Calculation Results of Reaction Time Minimization with Changing Initiator Feed and Temperature.

AC

CE

PT

ED

M

Figure 10 Calculating Process of the Adaptive Smoothing Strategy.

36

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 1 Simulation Results of Free Volume, Conversion, Mn and Mw.

37

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 2 Simulation Result of Termination Rate Coefficient Kt in Gel Effect.

38

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 3 Simulation Result of Chain Propagation Rate Coefficient Kp in Glass Effect.

39

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 4 Simulation Result of the Efficiency of Initiator f in Cage Effect.

40

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 5 Simulation Result of Original Model and Smooth Models under Different Smooth Parameters.

41

ACCEPTED MANUSCRIPT

Start

Y

Optimize on problem with ξ=ξ0 r=r0

N r>rmax

N r=r
p0

Optimization success? Y

ξsuccess=ξ rsuccess=r Successful End

Y

Optimize on the original model

Y

Accurate enough?

Optimize on the test problem

Optimization success?

Y

ξsuccess=ξ rsuccess=r

CR IP T

N

Optimization success?

N

N

Y Failed End

AN US

Simulate on the original model

N

pξ < pξ,min

N

pξ= pξ / 2

Optimize on problem with ξ=ξsuccess /pξ r=rsuccess

Optimization success?

Y

ED

M

Y

ξsuccess=ξ rsuccess=r

N

PT

From r adjusting procedure?

Y

N Optimize on problem with ξ=ξsuccess r=rsuccess / pr

CE AC

ξsuccess<ξtest rsuccess=rtest

Y

Optimization success?

Y

pr < pr,min N

N

Figure 6 Flowchart of the Adaptive Smoothing Strategy.

42

pr = pr /2

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 7 Calculation Results of Conversion Maximization with Changing Initiator Feed.

43

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 8 Calculation Results of Conversion Maximization with Changing Initiator Feed and Temperature.

44

ACCEPTED MANUSCRIPT

x 10

-4

1 1 stage 5 stages 10 stages 20 stages 40 stages

7 -

I2/mol*min1

6 5

0.8

Conversion

8

4 3 2

0.6 1 stage 5 stages 10 stages 20 stages 40 stages

0.4 0.2

1 0

0

10

20

30

40

0

50

0

10

1 stage 5 stages 10 stages 20 stages 40 stages

0.2

Free Volume

Avergae Molecular Weight

2.5

0.15 0.1

10

20

30

40

Time/min

x 10

5

2

40

50

1 stages 5 stages

1.5

10 stages 20 stages 40 stages

1 0.5

AN US

0.05

0

30

Time/min

0.25

0

20

CR IP T

Time/min

0

50

0

10

20

30

40

50

Time/min

AC

CE

PT

ED

M

Figure 9 Calculation Results of Reaction Time Minimization with Changing Initiator Feed and Temperature.

45

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 10 Calculating Process of the Adaptive Smoothing Strategy.

46

ACCEPTED MANUSCRIPT

List of Tables: Table 1 Mechanism of Free Radical Polymerization Table 2 General Model Parameters for Monomer-Polymer-Initiator Systems Table 3 Specific Model Parameters for MMA-PMMA-AIBN System Table 4 Operation Conditions of Simulation Table 5 Comparison of the End Results of Simulation between Smooth Model and Original Model

CR IP T

Table 6 Operation Conditions and Target Average Molecular Weight

Table 7 Calculation Results of Conversion Maximization with Changing Initiator Feed Table 8 The Max Conversion Results of Dynamic Optimization

Table 9 Calculation Results of Conversion Maximization with Changing Initiator Feed and

AN US

Constant Temperature Table 10 Optimization Results of Minimize Reaction Time

Table 11 Optimization Results of Smooth Models with Different Smoothing Parameter and

AC

CE

PT

ED

M

Original Model

47

ACCEPTED MANUSCRIPT

Table 1 Mechanism of Free Radical Polymerization reaction types

descriptions 𝐼2 →

𝑘𝑑

𝐼• 𝑘

chain initiation

𝐼•

chain propagation

R•𝑛

𝑀→

bimolecular termination

R•𝑛

• 𝑅𝑚 →

disproportion termination

R•𝑛

• 𝑅𝑚 →

transfer to monomer

R•𝑛

𝑀→

transfer to the polymer

R•𝑛

𝑃𝑚 →

𝑅•

𝑘𝑝

𝑘𝑡𝑐

𝑘𝑡𝑑

𝑘𝑡𝑟 𝑀

𝑘𝑡𝑟 𝑃

• 𝑅𝑛+

𝑃𝑛+𝑚 𝑃𝑛

𝑃𝑛 𝑃𝑛

𝑃𝑚

𝑅•

• 𝑅𝑚

AC

CE

PT

ED

M

AN US

𝑀→

CR IP T

initiator decomposition

48

ACCEPTED MANUSCRIPT

Table 2 General Model Parameters for Monomer-Polymer-Initiator Systems

m

0.5

n

1.75

A

1.11

B

1

R

1.987

Units

cal g-1 mol-1

CR IP T

Values

AC

CE

PT

ED

M

AN US

Parameters

49

ACCEPTED MANUSCRIPT

Table 3 Specific Model Parameters for MMA-PMMA-AIBN System Parameters

Values 7

85

𝑝

e p

Units L mol-1 min-1

4353 /𝑅

9

5 88

e p

L mol-1 min-1

7 /𝑅 44 /𝑅

CR IP T

6 9 6e p /

167 383

𝑝

2

K-1

AN US

3

𝐾3

K

4

48

𝑝

𝑓

K

e p

4

446 /

e p

3

8

K-1

/

63

M

0.58

e p

3 66

ED

6

min-1

/𝑅

PT

C 𝑓

AC

0.069

3

CE

ε

0.42

96647 95 4

9 48

3

3

64 33 e p

4

388 /𝑅

5

50

ACCEPTED MANUSCRIPT

Table 4 Operation Conditions of Simulation on PMMA Operation Variables Operation Values

Units

Initiator

0.01766

mol

Monomer

9.62767

mol mol min-1

Initiator Feed 80

℃

Reaction time

100

min

Initial Volume

1.2

CR IP T

Temperature

AC

CE

PT

ED

M

AN US

L

51

ACCEPTED MANUSCRIPT

Table 5 Comparison of the End Results of Simulation between Smooth Model and Original Model 0 Smooth Parameter 1.0e-1

1.0e-2

1.0e-3

1.0e-4

1.0e-5

1.0e-6

1.0e-7 (Original model)

Mn (105)

1.7098 1.7321

1.7481

1.7530

1.7539

1.7540

1.7540

1.7540

Mw (105)

8.3199 9.1366

9.2818

9.3099

9.3144

9.3151

9.3151

9.3152

0.85587 0.87078 0.87194 0.87209 0.87210 0.87211 0.87211

0.87211

AC

CE

PT

ED

M

AN US

CR IP T

Conversion

52

ACCEPTED MANUSCRIPT

Table 6 Operation Conditions and Target Average Molecular Weight Variables

Values

Units

Initiator

0.01766

mol

Monomer

9.62767

mol

𝑒

mol min-1

5

Temperature

80

Reaction time

100

Initial Volume

1.2

Result of Mw

89

min L

4 5

g mol-1 g mol-1

AC

CE

PT

ED

M

AN US

Result of Mn

4 5679

℃

CR IP T

Initiator Feed

53

ACCEPTED MANUSCRIPT

Table 7 Calculation Results of Conversion Maximization with Changing Initiator Feed Direct Optimizing

Smooth Strategy

Number Success

Total Time

Success

Total Time

Number

(s)

Number

(s)

1

0

7364.1

0

2431.2

5

44

8028.0

99

1942.7

10

33

30144.0

95

6223.3

AC

CE

PT

ED

M

AN US

CR IP T

of Stages

54

ACCEPTED MANUSCRIPT

Table 8 The Max Conversion Results of Dynamic Optimization Number of Case 1

Case 2

1

NA

0.8692

5

0.8830

0.9025

10

0.8844

0.9037

20

0.8850

0.9040

40

0.8853

0.9041

AC

CE

PT

ED

M

AN US

CR IP T

Stages

55

ACCEPTED MANUSCRIPT

Table 9 Calculation Results of Conversion Maximization with Changing Initiator Feed and Constant Temperature Direct Optimizing

Smooth Strategy

Number Success

Total Time

Success

Total Time

Number

(s)

Number

(s)

1

22

1612.3

98

1514.7

5

15

7370.3

83

7812.5

10

7

16089.0

82

CR IP T

of Stages

AC

CE

PT

ED

M

AN US

18118.1

56

ACCEPTED MANUSCRIPT

Table 10 Optimization Results of Minimize Reaction Time 1

5

10

20

40

Reaction Time (min)

51.62

33.07

31.16

30.16

30.05

Conversion

0.9016

0.9000

0.9000

0.9000

0.9000

AC

CE

PT

ED

M

AN US

CR IP T

Number of Stages

57

ACCEPTED MANUSCRIPT

Table 11 Optimization Results of Smooth Models with Different Smoothing Parameter and Original Model Smooth Models with Different Smooth Parameter

Original

Acceptable Constraints Range 1.0e-3

1.0e-4

1.0e-5

1.0e-6

model

Success Number

7

10

17

10

12

7

7

Calculating Time (s)

1825.3

4089.8

5524.7

4834.3

8358.4

10351.5

16089.0

Aver. SNOPT Iterations

122

161

180

146

141

116

139

Success Number

41

54

46

35

29

29

28

Calculating Time (s)

5025.2

6358.7

5277.7

6082.3

9717.5

12580.0

15999.2

Aver. SNOPT Iterations

211

246

207

168

157

187

184

Success Number

61

60

54

30

33

20

26

Calculating Time (s)

2366.3

2264.6

3086.3

3572.0

15993.4

27513.3

166

192

CR IP T

0.4

1.0e-2

AN US

0.05

1.0e-1

10070.

0.8

195

212

181

AC

CE

PT

ED

M

Aver. SNOPT Iterations

58

152

6

170