Identifiability Analysis and Parameter Estimation of Microgel Synthesis

Feb 6, 2019 - We present a methodology for a parameter identifiability analysis, which approximates the feasible parameter set as a box by solving a s...
0 downloads 0 Views 783KB Size
Article pubs.acs.org/IECR

Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Identifiability Analysis and Parameter Estimation of Microgel Synthesis: A Set-Membership Approach Falco Jung,† Franca A. L. Janssen,† Agnieszka Ksiazkiewicz,‡ Adrian Caspari,† Adel Mhamdi,† Andrij Pich,‡,¶ and Alexander Mitsos*,† †

Aachener Verfahrenstechnik-Process Systems Engineering, RWTH Aachen University, 52074 Aachen, Germany DWI Leibniz Institute for Interactive Materials e.V., 52074 Aachen, Germany ¶ Institute of Technical and Macromolecular Chemistry, RWTH Aachen University, 52074 Aachen, Germany Downloaded via WEBSTER UNIV on February 14, 2019 at 18:05:37 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



S Supporting Information *

ABSTRACT: Functional microgels with tailored structure and specific properties are required for medical and technical applications, thus motivating model-based optimization of their fabrication processes. An important step in the creation of accurate models is parameter estimation. We present a methodology for a parameter identifiability analysis, which approximates the feasible parameter set as a box by solving a series of constrained dynamic optimization problems. The method is applied to the synthesis of microgels based on two monomers, N-vinylcaprolactam and N-isopropylacrylamide, and the cross-linker N,N-methylenebis(acrylamide). The results show that kinetic parameters corresponding to the reaction of the monomers are identifiable as are a subset of the kinetic parameters involving the cross-linker. The reaction kinetics of the cross-linking are faster in comparison to the main polymerization reaction for N-vinylcaprolactam; this allows for an improved understanding of the occurring reaction phenomena. The reaction kinetics of the cross-linking are not identifiable for N-isopropylacrylamide for the given experimental setup; model-based experimental design for parameter precision might enable their identification. The results also indicate potential for model simplification and allow us to make suggestions toward the enhancement of Raman spectroscopy measurements.



on models of the synthesis of microgels. Hoare and McLean19 presented a kinetic model to analyze the local composition of thermoresponsive microgels based on the monomer Nisopropylacrylamide (NIPAM). They determined separate kinetic model parameters based on Price-Alfrey coefficients. Janssen et al.20,21 made a step toward a model-based approach by reviewing current work on the modeling of microgel synthesis and presenting a detailed mechanistic model. They investigated the radical precipitation polymerization of monomers Nvinylcaprolactam (VCL) and NIPAM with the cross-linker N,N-methylenebis- (acrylamide) (BIS) and initiator 2,2-azobis(2-methylpropionamidine) dihydrochloride (AMPA) in aqueous solution and estimated the respective kinetic parameters numerically based on a combination of transition-state theory calculations22 and experimental data. They compared the predictions of the developed model to calorimetric and Raman spectroscopy measurements.23 We simplify the model

INTRODUCTION Functional microgels respond to environmental stimuli, for example changes in temperature or pH level,1 by reversible alteration of their size and shape.2 This property of microgels enables diverse applications, including (i) foliar fertilizer delivery systems that enable the controlled release of nutrients to plants,3 (ii) membrane coatings that adjust their permeability to the temperature of the environment,4 and (iii) drug delivery systems that enable the controlled release of drugs in specific organs.5 The range of applications gives rise to the need for microgels with varying functionality. The properties of microgels strongly depend on the operating conditions of the fabrication process.2 Examples of influencing conditions are the employed reactants,6 the process temperature and pressure,7 and the reactor hydrodynamics.8,9 The understanding of the quantitative effect of operating conditions on microgel functionality is necessary to produce custom-designed microgels. In a model-based approach, this requires the development of appropriate models to improve process understanding and to enable optimization of the process design and its integrated design and operation.10,11 While polymerization modeling has been in the focus of research for several decades,12−18 only a few publications focus © XXXX American Chemical Society

Special Issue: Dominique Bonvin Festschrift Received: October 24, 2018 Revised: January 10, 2019 Accepted: January 17, 2019

A

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research presented by Janssen et al.20 based on their findings. Additionally, we conduct a rigorous parameter identifiability analysis to reduce the number of parameters to be estimated and to improve the confidence in the determined parameter values. Parameter identifiability analysis determines which unknown model parameters can be uniquely estimated on the basis of available experimental data.24 If all model parameters are identifiable, the confidence in the predictive capabilities of the model may be improved by choosing appropriate measurement configurations. In case parameters are not identifiable, model reduction methods have the potential to simplify the model until all remaining parameters are identifiable. Various definitions of parameter identifiability are available in the literature. Parameter identifiability can be analyzed by a range of methods, which are reviewed in detail by Miao et al.25 Commonly, structural and practical identifiability are distinguished. Structural (a priori) identifiability establishes whether the structure of the model equations enables the unique determination of model parameters given a set of measurable variables.24 It does not account for the precision or frequency of the available observations.26 Therefore, it only poses a necessary condition for parameter identifiability. In contrast, practical identifiability accounts for flaws in the measurement sensors and aims to predict confidence intervals for the unknown parameters.27,28 It can be evaluated locally, that is in the neighborhood of a specific parameter set, or globally, which is over the entire defined parameter space.27 A well-established method for the local sensitivity analysis is the estimation of parameter confidence intervals based on the Fisher information matrix.29 This method requires only a sensitivity integration and is therefore numerically cheap in comparison to global methods. Thus, it has been frequently applied in the domain of process systems engineering.28,30 However, results are exact only for linear models,31 as estimates are determined by linearizing the model. Additionally, required estimates of the measurement precision over the course of the experiment are typically hard to obtain. First presented in the scope of set-membership parameter estimation,32 set-inversion methods enable a global identifiability analysis by an approximation of the feasible parameter set, which is the set of parameters for which respective model outputs remain within defined bounds of the available experimental data.33 However, the application of the required set-inversion algorithms is numerically expensive for all but very small process models.34 To lower the computation cost, the feasible set can be approximated via a box by solving a constrained optimization problem.35 The respective constraints of the optimization problem can be formulated as the absolute or relative error of single measurement points, or based on a statistical analysis of the model deviation from all available measurements.36 A box-approximation of the feasible parameter set enables an identifiability analysis of the model parameters at reasonable numerical cost. However, results are guaranteed to be the global solution only if respective optimization problems are solved with a method guaranteeing global optimality. Set-inversion methods have already been applied in the context of parameter identifiability by multiple authors. Paulen et al.34 investigated set-inversion methods in the context of guaranteed parameter estimation of nonlinear dynamic systems. They investigated the effect of various domain and CPU-time reduction strategies to enable the fast solution of set-inversion problems. Gottu Mukkula and Paulen33 applied the method in the context of optimal experimental design. They studied

dynamic and static models with a maximum of two system states. The problem was solved using a multiple-shooting algorithm via CasADi with IPOPT. Walz et al.35 approximated the feasible parameter set by a box also in the context of optimal experimental design. They studied dynamic systems with a maximum of three system states and an analytic solution. They solved the optimization problems to global optimality with the solver BARON37 in GAMS.38 Goerke and Engell39 determined the feasible parameter-set for the kinetic parameters in a radical polymerization case study. They employed an evolutionary algorithm to determine the feasible set. All of the above publications used fictitious data generated based on the available model. Moreover, the bounding errors defining the precision of the generated fictitious data was assumed to be known. Recently, we analyzed the feasible parameter set of a microgel homopolymerization process and compared various formulations for the constraints of the feasible parameter set.36 We computed separate feasible parameter sets by sampling a twodimensional parameter space. This method is feasible only for small examples because of the exponential increase in computational time with the number of parameters to be inspected. Hence, while set-inversion has been applied in several case studies involving artificial data generated based on numerical models, its application to real experimental data has not yet been reported. In this work, we conduct a rigorous parameter identifiability analysis based on the approximation of the feasible parameter set for the synthesis of microgels based on VCL or NIPAM. In contrast to the existing studies, we base our investigation on real measurement data; thus, we have to consider the mismatch between model and experiment. First, we present the model of the microgel synthesis reaction and some details about the available experimental data. We then introduce the methods for parameter estimation and the identifiability analysis in detail. Finally, we present and discuss the obtained results and draw final conclusions.



FUNCTIONAL MICROGEL SYNTHESIS We now present the microgel synthesis under investigation and review available measurements. Additionally, we develop a process model and normalize it for parameter estimation and identifiability analysis. We focus on the synthesis of functional microgels based on the monomer VCL or NIPAM with the cross-linker BIS. Detailed experimental procedures of the synthesis carried out in batch experiments have been published recently.1,20,40 Here, we focus on isothermal synthesis conditions at 70 °C. The experimental procedure involves first the dissolution of monomer and cross-linker in solvent. The radical polymerization reaction is initiated by the thermally decomposing initiator AMPA (I), which forms primary radicals (R0). The radicals propagate with monomer (Mj) to form growing radical chains (Rnj ) with chain length n and the terminal end of type j, with j = 1 denoting a monomer and j = 2 a crosslinker terminal end. Upon reaching a critical chain length, the polymer chains precipitate from the liquid phase and form a disperse phase. The cross-linker has two double bonds. Hence, when one cross-linker is attached to the polymer chain, the second double bond is a pendant double bond (PDB). When the PDB reacts with another radical, a cross-link (X) is formed. The cross-linked polymer chains form the microgel. The synthesis can be carried out with various initial concentrations of the reactants to influence the cross-linker density or particle size. B

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research



MECHANISTIC MODEL OF MICROGEL SYNTHESIS A mechanistic model of the microgel synthesis has already been published by Janssen et al.20,21 The results presented therein allow for additional assumptions and simplifications of the model, which are introduced and discussed below. They enable the reduction of the number of unknown model parameters. Moreover, the reduced model size enables a faster computation for future model-based control applications. We aim to use the same model equations (also called model structure) for the models describing the synthesis of microgels based on the two monomers VCL or NIPAM. The differences of the two systems are then captured only in the values of the parameters describing the kinetics of the various reactions. In the model used here, we account for the decomposition of the initiator, the propagation reactions between an active radical chain and the monomers, and a termination reaction between two active radical chains. Additionally, we consider the crosslinking reaction between a PDB and an active radical chain. The main assumptions we make are the following: • The disperse system can be modeled in a pseudohomogeneous phase because the reaction locus is mainly in the gel phase.20 The shift of the reaction locus is caused by the short critical chain length at which microgels precipitate and the phase equilibrium of monomers that favors the gel phase. • The initiator efficiency, which describes the fraction of initiator that initiates a radical chain, is assumed to be equal to one because of the thorough nitrogen purging conducted here. The efficiency can be lower for other experimental setups. • The kinetics of the cross-linking reaction, which is the reaction of a pendant double bond with an active radical chain, are assumed to be as fast as the kinetics of a propagation reaction with a BIS monomer. • The kinetics of termination reactions do not depend on the terminal end of the active radical chains, as they are mostly diffusion-limited. • Chain-transfer reactions are excluded from the model, as quantum-mechanical calculations determined that they are slow in comparison to the propagation reactions.20 • The volume and total mass of the reactor content is equal to the volume and mass of the solvent, as the system is strongly diluted. • The kinetics of the propagation reactions depend only on the type of radical at the end of the chain.15 • The enthalpy accumulation caused by reactants than other the solvent is negligible, as the system is strongly diluted. The reaction mechanism with associated reaction constants is given in Table 1. The rate of decomposition of the initiator, kd, is provided by the manufacturer.41 The model resulting from the listed assumptions is c İ (t ) = −kdc I(t )

(1a)

c M1 ̇ (t ) = −k p11c R1(t )cM1(t ) − k p21c R2(t )c M1(t )

(1b)

c M2 ̇ (t ) = −k p12c R1(t )cM2(t ) − k p22c R2(t )c M2(t )

(1c)

Table 1. Reaction Mechanism with Associated Reaction Rate Constant (km)a type of reaction

a

decomposition

I → 2R0

propagation

R in + Mj ⎯→ ⎯ Rnj + 1 + (j − 1)PDB

termination by disproportionation

R in + Rmj → Pn + Pm

cross-linking

PDB + R in ⎯⎯⎯⎯→ R n2+ 1 + X

k pij

kt

k p, i 2

km, where m is the type of reaction and the involved reactants.

c R2 ̇ (t ) = −k t(c R2(t )2 + c R1(t )c R2(t )) + k p12c R1(t )c M2(t ) − k p21c R2(t )c M1(t )

(1e)

c PDB ̇ (t ) = ( +k p12c R1(t ) + k p22c R2(t ))(c M2(t ) − c PDB(t )) (1f)

where cI, cM1, and cM2 are the molar concentrations of initiator AMPA (I), monomer VCL or NIPAM (M1), and cross-linker BIS (M2), respectively. cR1 and cR2 are the molar concentrations of active radical chains with a monomer or cross-linker terminal end, respectively; they are equivalent to the zeroth-order moments of the radical chains. The enthalpy-transfer rate ΣR, which is the enthalpy flow rate due to the exothermic polymerization reaction, is modeled based on the assumption that the major contributors to the enthalpy release are the propagation reactions. Hence, the enthalpytransfer rate can be determined as 2

Σ R (t ) =

2

∑ ∑ −ΔHR,ij(k pijcRi(t )cMj(t ) i=1 j=1

+ (j − 1)k p,i2c Ri(t )c PDB(t ))V

(1g)

where ci is the molar concentration of component i ∈ [VCL/ NIPAM, BIS], ΔHR,ij is the enthalpy of reaction of the respective reactions, determined by quantum mechanical calculations,22 and V is the volume of the solvent in the reactor. The weight fraction of the monomer can be determined based on the concentration by wM1(t ) =

c M1(t )VMM1 m H 2O

(1h)

where MM1 is the molecular weight of monomer and mH2O is the mass of solvent in the reactor. The weight fraction of the crosslinker wM2 cannot be measured with Raman spectroscopy because of the low amounts of cross-linker present in the reactor. For the parameter estimation and identifiability analysis, we normalized the model (1a)−(1h). We define the kinetic reaction constant kp11 as a reference value, as it is the main reaction occurring during the synthesis. For the scaling of the concentration of active radical chains, the parameters c0R1 and c0R2 are introduced with values of 10−2 and 10−4 mol/m3, respectively. Normalization yields five unknown ratios between kinetic reaction constants and the reference reaction to be estimated and analyzed in the parameter identifiability analysis. The respective ratios are presented in Table 2. The explanation of the normalized model is provided in the Supporting Information.

c R1 ̇ (t ) = +2kdc I(t ) − k t(c R1(t )2 + c R1(t )c R2(t )) − k p12c R1(t )c M2(t ) + k p21c R2(t )c M1(t )

mechanism kd

(1d) C

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Table 2. Unknown Parameters in the Normalized Model σi̅ =

unknown parameter

rI = rt =

kd

description

0 k p11c R1

ratio of the reaction constants of the initiator decomposition and the reference reaction

kt k p11

ratio of the reaction constants of the termination reaction and the reference reaction

rp12 = rp21 = rp22 =

k p12 k p11 k p21 k p11 k p22 k p11

(3)

j=1 k=1

(4)

where mH2O is the mass of water in the reactor, Ĉ p,ins the heat capacity of the instrumentation, Ĉ pH2O the specific heat capacity of water, TR the reactor temperature, and Q̇ loss the heat loss to the environment via the lid. The derivative of the reactor temperature dTR is provided by the measurement equipment; we dt checked the provided values using finite-differences of the reactor temperature measurements. The enthalpy flow across the reactor wall, Q̇ rtc, is also measured by the measurement equipment based on two sensor bands that are mounted to the reactor. The heat loss to the environment, Q̇ loss, varies over the course of the reaction because of changes of the temperature of the reactor lid. The reactor lid heats up because of the heat released during the reaction. The heat loss to the environment is measured immediately before and after the reaction, as Q̇ loss,0 and Q̇ loss,∞, respectively. During the reaction, it is interpolated as a linear function of the thermal conversion of the reactants ξth

ratio of the reaction constants of the propagation reaction between an active radical chain with BIS terminal end and BIS and the reference reaction

The final model has a total of six differential and two algebraic equations with six model parameters, of which five are unknown. Because of the good scaling achieved by the normalization, the model is less stiff and derivatives are moderate. Hence, it can be integrated and optimized quickly and robustly.



EXPERIMENTAL DATA Experimental data have been published, including calorimetric20 and Raman spectroscopy measurements.23 In addition to the already published data, new data of the syntheses with an improved nitrogen purging in the reaction mixture prior to polymerization are available. Oxygen hinders the initiation of the polymerization by slowing it down. To ascertain minimum levels of oxygen, the solution was purged with nitrogen directly for 30 min, rather than purging the air−solution interface, as in the previously reported procedure.1 Initial concentrations of monomer and initiator were the same for all respective reactions; exact values can be found in the Supporting Information. Aside from the initial amounts of reactants, the only degree of freedom that can be chosen freely is the reactor temperature. All syntheses were carried out at 70 °C, as the synthesis duration at this temperature is long enough to enable Raman spectroscopy measurements, while the heat release is sufficient for detection by calorimetry. Data for initial crosslinker molar fractions of 0, 2.5, and 5.0 mol % are available for the syntheses based on VCL, and data for initial cross-linker molar fractions of 2.5 mol % are available for the syntheses based on NIPAM. Links to the data are included in the Supporting Information. For the parameter estimation and identifiability analysis conducted here, we obtain the precision of each measurement in terms of the standard deviation σ̅ i based on repetitions of measurements ỹi,k in separate experiments, where k = 1, ···, nexp and nexp is the number of repetitions of an individual experiment. First, we determine the mean measurement value, ỹi, as the mean of repetitions of a measurement as

Q̇ loss(t j) = Q̇ loss,0 + (Q̇ loss, ∞ − Q̇ loss,0)ξth(t j)

(5)

The determined enthalpy-transfer rate, Σ̃R, is one measurement that is available for parameter estimation and the identifiability analysis. Additionally, Raman spectroscopy23 measurement of the mass fractions of monomer and polymer were reported20 and are available for parameter estimation and identifiability analysis. The weight fraction measurements of the monomer and polymer content are redundant in a batch experiment, as the sum of monomer and polymer weight fraction measurements should always equal the initial amount of monomer added to the reactor. It was shown that the predictions of the polymer content still suffer from significant measurement noise, while the monomer measurements show much smaller measurement noise.20 Therefore, we consider only the measurements of the monomer mass fraction in this work. However, if the accuracy of polymer measurements can be improved, it is potentially fruitful to include the measurements again in the analysis.



SET-MEMBERSHIP PARAMETER ESTIMATION To analyze parameter identifiability with the set-membership approach presented, here we solve several optimization problems: 1. Parameter estimation problem (7) to estimate the parameter vector and determine the corresponding optimal model outputs. 2. Single-parameter identifiability problem (10) for every parameter to determine the rectangular approximation of the feasible parameter set.

nexp k=1

nexp

∼ dT (t ) y ji ̂ ̂ O R j zzzz − Q̇ (t j) Σ R (t j) = jjjjCp,ins + m H2OCp,H rtc 2 dt z{ k − Q̇ loss(t j)

ratio of the reaction constants of the propagation reaction between an active radical chain with BIS terminal end and VCL/NIPAM and the reference reaction

∑ yi ,̃ k (t j)

N

∑ ∑ [yi ̃ (t j) − yi ,̃ k (t j)]2

Real-time calorimetric measurements of the enthalpy flow across the reactor wall Q̇ rtc and measurements of the reactor temperature TR are available at an interval of two seconds. They enable the determination of the measurement of the enthalpytransfer rate, Σ̃R, as

ratio of the reaction constants of the propagation reaction between an active radical chain with VCL/NIPAM terminal end and BIS and the reference reaction

1 yi ̃ (t j) = nexp

1 nexpN − 1

(2)

where N is the number of measurement time points and tj is the time of measurements with j = 1, ···, N. We then compute the measurement precision as the standard deviation of the residuals of individual measurements from the mean measurement value as D

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

In a second step, we determine the feasible parameter set, which is the set of parameters Pe for which the mismatch between model predictions and optimal model predictions is within specified constraints, that is

The model of the microgel polymerization reaction is a semiexplicit index-1 differential-algebraic equation system (DAE) of the form ẋ(t ) = f(x(t ), z(t ), u(t ), p)

(6a)

0 = g(x(t ), z(t ), u(t ), p)

(6b)

y(t ; p) = m(x(t ), z(t ), u(t ), p)

(6c)

0 = h(x(t0), z(t0), p)

(6d)

Pe: = {p ∈ P: c(y(t j ; p), y(t j ; popt), δ(t j)) ≤ 0}

f: nx × nz × nu × np → nx

g: nx × nz × nu × np → n y m: nx × nz × nu × np → n y

δiLB(t j) ≤ (yi (t j ; popt) − yi (t j ; p)) ≤ δiUB(t j), ∀ i , t j

define the process model with the initial conditions h: nx × nz × np → nx . For improved readability, we use the notation y(t; p) to indicate that the output trajectories y are computed by integration of eq (6a) for a given parameter vector p. To access the identifiability of the unknown model parameters based on available experimental data ỹ, optimal values of the unknown parameters popt are determined in a first step by minimization of the weighted residual between model predictions and available measurements, that is n ỹ

N

p∈P

i=1 j=1

1 |y (t j ; p) − yi ̃ (t j)| σi , j i

(8b)

where Pe is the set of feasible model parameters. The function c constrains the deviation between a measurement i at time tj and the respective model prediction to acceptable errors δ. Their formulation depends on the assumptions made about the statistical distribution of the measurement error.36 Here, we make the assumption of a normally distributed measurement error, which can be bounded by three standard deviations 3σ̅ i. A parameter vector p is considered infeasible if the deviation between any model prediction yi from the optimal model predictions yi at any time is greater than the set lower (δLB i ) or upper measurement error bounds (δUB ) i

where x(t ) ∈ nx are differential states, z(t ) ∈ nz are algebraic states, u(t ) ∈ nu are inputs, y(t ) ∈ n y are measured variables, and p ∈ np are time-invariant parameters. The functions

min ∑ ∑

(8a)

(9a)

The lower and upper measurement error bounds are determined based on the bounding mean measurement precision 3σ̅ i and the mismatch between the model predictions with the optimal parameter set yi and the available experimental data ỹi as δiLB(t j) = min[yi (t j ; popt) − yi ̃ (t j), −3σi̅ ]

(9b)

δiUB(t j) = max[yi (t j ; popt) − yi ̃ (t j), +3σi̅ ]

(9c)

We take the mismatch between the model predictions with the optimal parameter set and the available experimental data into account to avoid a mismatch with a greater absolute value than the bounding error of 3σ̅ i reducing the size of the feasible parameter set. For a two-dimensional case, Figure 1 illustrates a nonconvex feasible parameter set Pe, marked in a dark-gray color. It contains all parameter vectors for which the respective model predictions satisfy the constraints defined by (8b). The exact feasible parameter set is difficult to determine. Therefore, it is approximated by a rectangular box, whose edges are defined

(7)

In problem formulation (7), P is an a priori estimate of the parameter space and ỹi(tj) are available measurements with respective standard deviations σi,j at times tj ∈ [t0,tend] with j = 1,···,N. A solution of problem (7) yields an optimal parameter popt, which leads to the optimal model outputs y (t; popt) by integration of the model (6a). The optimization problem should be solved globally to ensure that a parameter vector for optimal agreement with available measurements is found. Therefore, we solve it with the global enhanced scatter-search method implemented in the open-source software for dynamic optimization MEIGO42 with the dynamic hill climbing method43 as solver. However, global convergence cannot be guaranteed as the solver is nondeterministic. As introduced, the available experiments have been performed with varying initial molar fractions of the cross-linker BIS of 0, 2.5, and 5 mol % for the syntheses based on VCL, and crosslinker molar fractions of 2.5 mol % are available for the syntheses based on NIPAM. For the parameter estimation and identifiability analysis, we use only the experiments with an initial concentration of 5 mol % for the VCL case and leave other experimental results for validation. The lower and upper bounds of the parameters are set to 10−3 and 103, respectively. This choice is based on the assumption that if a reaction is more than 1000 times faster or slower than the main propagation reaction, its effects would not be measurable with the available measurement apparatus, as respective reactions are too fast or too slow. We run the parameter estimation for a total of 1000 s; this choice is based on our prior experience with the model and the used solver.

Figure 1. Feasible parameter set with an approximation by a box. E

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research by the parameter vectors (p1,min, p2,min)T and (p1,max, p2,max)T for the two-dimensional case. The minimum and maximum parameter values defining the box approximation can be obtained by solving a series of constrained dynamic optimization problems. For each parameter pk to be analyzed, a maximum pk,max and a minimum pk,min are determined by solving the respective optimization problems34 max pk,max − pk,min

quarter of the evaluations. This increases the confidence that the result is close to the global optimum. The numerical results of the parameter estimation are shown in Table 3 for the syntheses Table 3. Parameter Estimation Result for VCL- and NIPAMBased Syntheses parameter rI rt rp12 rp21 rp22

[p k,min, p k,max ]∈ P × P

s.t. δiLB(t j) ≤ (yi (t j ; popt) − yi (t j ; pk,min)) ≤ δiUB(t j), ∀ i , ∀ t j δiLB(t j) ≤ (yi (t j ; popt) − yi (t j ; pk,max )) ≤ δiUB(t j), ∀ i , ∀ t j

estimated value VCL −2

1.79 × 10 6.37 × 10−1 3.43 × 101 2.81 2.56 × 101

estimated value NIPAM 1.78 × 10−2 2.48 × 10−1 1.46 × 10−1 1.60 × 10−3 1.00 × 10−3

based on VCL and NIPAM. The results for the parameters rI and rt are comparable in their order of magnitude. Kinetic parameters describing the propagation reactions between monomers and cross-linker deviate by several orders of magnitude; reaction kinetics in the NIPAM system is slower than in the VCL system. The difference in the kinetics of the cross-linking reactions is expected to cause a difference in the distribution of the cross-linker within the microgel. Because of the fast reaction kinetics, the cross-linker is expected to accumulate at the center of microgels based on VCL, an expectation that was already verified experimentally.40 For NIPAM-based microgels, the accumulation of cross-linker in the core of the microgel is expected to be less dominant because of the slower cross-linking reaction kinetics. Figure 2 shows the model predictions and experimental data after the parameter estimation step. The error bars indicate the standard deviation of the respective measurements determined based on the available repetitions of the experimental setups; see eq (3) for details. Overall, the agreement of the model predictions and experimental data is good. For the syntheses based on VCL, especially the two peak-structure of the enthalpytransfer rate is predicted well by the model. Some deviations in the enthalpy-transfer rate occur between 50 and 100 s. For the syntheses based on NIPAM, the measurements of the mass fraction of NIPAM agree well with model predictions. The enthalpy-transfer rate is overpredicted by the model in the last third of the synthesis. A reason for the difference between model predictions and experimental data could be an invalid modeling assumption. However, a relaxation of the discussed modeling assumptions does not lead to an improved agreement. Another explanation is a minor heat stream that is not yet considered in the energy balance, cf. eq (4) or imprecisions in the numerical values of the enthalpies of reaction ΔHR,ij, which were determined based on quantum-mechanical calculations.22 Set-Membership-Based Identifiability Analysis. The results of the approximation of the feasible parameter set and the parameter identifiability are presented in Figure 3. The respective measurement error bounds are presented in Table 4. The values for measurement standard deviation of the enthalpy-transfer rate, σ̅ ΣR, and the monomer weight fraction, σ̅ wM1, are of the same order of magnitude for the syntheses based on VCL and NIPAM. σ̅ ΣR is higher for the synthesis based on VCL; this is a result of the faster dynamics of the VCL synthesis in comparison to the synthesis based on NIPAM that impact the precision of the measurements (cf. Figure 2). σ̅ wM1 is slightly higher for the synthesis based on NIPAM; several effects may cause the worse performance, as will be discussed below. The standard deviations were determined based on the available data

(10)

where the vectors pk,min and pk,max are the vectors of model parameters for which parameter pk is at its respective minimum or maximum value, that is for the case of the minimum pk,min := (p1,k,min, p2,k,min, ···, pk,min, ···, pnp,k,min)T. For the example of Figure 1, the parameter vectors p1,min = (p1,min, p2,1,min) and p1,max = (p1,max, p2,1,max) are the solution of problem (10) for parameter p1. The size of the feasible parameter set and the identifiability of model parameters are closely connected. An identifiability analysis determines whether a model parameter can be estimated based on available experimental data. A parameter is identifiable if the feasible parameter set covers only a limited space in the direction of one parameter pk. A parameter is not identifiable if the feasible parameter set covers a wide space in the direction of one parameter pk, because the parameter cannot be determined to sufficient accuracy; a cutoff value for the extent of the feasible parameter set is required to distinguish between identifiable and non-identifiable parameters. However, in practice, a cutoff value is often difficult to define and a qualitative discussion of the results of the set-membership parameter estimation is necessary. The feasible parameter set allows for the determination of confidence ranges, which are equal to the extent of the feasible parameter set in the direction of the parameter under investigation. The dynamic optimization problems in (10) are challenging to solve because of the high number of constraints and the nonlinearity of the model. We solved them using the framework for dynamic optimization DyOS.44 The system of differential equations is integrated with the integrator sLIMEX45 with absolute and relative tolerances set to 10−6. The model is imported as a functional mock-up unit (FMU)46 from the commercial software Dymola.47 The enhanced scatter search method implemented in the toolbox MEIGO42 is used for global optimization, while the dynamic hill climbing method is run for local searches.43 As the solver is not deterministic, global optimality cannot be guaranteed.



RESULTS AND DISCUSSION This section presents the results of the set-membership parameter estimation. Additionally, we discuss the numerical results to show the benefits of carrying out model normalization and highlight potential improvements to the used measurement methods. Parameter Estimation. The enhanced scatter search method for the parameter estimation evaluated the model over 105 times, while the result is already found after less than a F

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 2. Agreement between model predictions and experimental data for copolymerization at 70 °C of VCL and BIS 5 mol % (top) and NIPAM and BIS 2.5 mol % (bottom). Whiskers indicate one standard deviation.

order of magnitude. Parameters rp12, rp21, and rp22 all cover multiple orders of magnitude. The results for both monomers indicate that parameters corresponding to the reaction of an active chain with a crosslinker end rp21 and rp22 have rather wide confidence ranges. This result agrees with our expectations, as the initial amount of crosslinker is low in comparison to the initial amount of monomer. Hence, reactions involving the cross-linker occur rarely; therefore, modeling predictions are less sensitive to changes in their kinetics. Parameters rI and rt affect the main polymerization reaction of the monomer and are therefore determined to a good degree of accuracy, as changes in the kinetic parameters impact the model predictions considerably. The values for these parameters are similar in the case based on VCL and the case based on NIPAM, which is in line with the assumption that the kinetics of the termination are mostly dependent on diffusive limitations and do not depend on the involved reactants. Parameter rp12 is estimated to good accuracy for the VCL case and to bad accuracy for the NIPAM case. This can be explained by the measurement of the enthalpy-transfer rate for the system based on VCL. It shows a two-peak structure and the height of the first peak depends on the amount of cross-linker added to the system. Hence, the measurement allows us to distinguish between the reactions of the cross-linker and the main monomer; hence, the kinetics of the reaction involving the cross-linker can be identified. The results of the analysis show that the considered model with the estimated model parameters can be assumed to be valid only for molar cross-linker fractions below 5%, which is the highest molar fraction that experimental data is available for. For higher cross-linker fractions, additional experimental analysis would be necessary to ensure that parameters corresponding to the cross-linker reactions rp21 and rp22 are estimated to a better accuracy. Experimental Validation. On the basis of the results of the set-membership parameter estimation, we expect the proposed model with the estimated parameters to be valid as long as the molar fraction of cross-linker remains below 5%, because the

Figure 3. Results of the parameter identifiability analysis. Crosses mark the optimal value determined in the parameter estimation, and solid lines indicate the confidence range of the parameter values. Blue are results for syntheses based on VCL, and red are those based on NIPAM.

Table 4. Measurement Standard Deviation for VCL- and NIPAM-Based Syntheses bounding error σ̅ ΣR [W] σ̅ wM1

value VCL

value NIPAM

−1

1.1 × 10−1

−4

9.0 × 10−4

5.1 × 10 2.4 × 10

from four repetitions of the syntheses with VCL as the main monomer and three repetitions of the syntheses with NIPAM as the main monomer. Each experiment provided approximately 400 data points from calorimetry and 30 data points from Raman spectroscopy measurements. Figure 3 presents also the confidence ranges of the parameter estimates. The results show a variation in the degree of accuracy to which the single parameters can be determined. For microgels based on the monomer VCL, parameters rI, rt, and rp12 are determined to good accuracy as their confidence range covers approximately 1 order of magnitude. However, the accuracy of parameters rp21 and rp22 is worse, as their confidence ranges cover multiple orders of magnitude. For microgels based on the monomer NIPAM, parameters rI and rt are accurate within 1 G

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 4. Agreement between model predictions and experimental data for copolymerization of VCL and BIS at 70 °C for BIS molar fractions of 0 mol % (top) and 2.5 mol % (bottom).

Figure 5. Active constraints for parameter vectors prt,max for microgels based on VCL (left) and NIPAM (right). Top: Enthalpy-transfer rate ΣR with ashed lines indicating point constraints at measurement times every 2 s. Bottom: Mass fraction of monomer wM1 with crosses indicating the time points of measurements.

Requirements on Measurement Techniques. Figure 5 presents the deviation ϵi of model predictions of measurement i with the parameter vector prt,max from the model predictions with the optimal parameter set popt

propagation reaction of an active radical chain with a cross-linker end and a cross-linker could not be estimated with satisfying precision. Therefore, higher fractions of cross-linker may lead to false predictions of the model. The agreement between model predictions and experimental data for varying initial concentrations of cross-linker for microgels based on the monomer VCL are presented in Figure 4. No additional parameters were adjusted to obtain the results; only the initial amount of crosslinker was changed to the respective experimental conditions. The agreement in both cases is good, which is an indicator for the validity of the modeling assumptions and the estimated parameter values.

ϵi(t j) = yi (t j ; popt) − yi (t j ; pr ,max ) t

(11)

Additionally, it shows the respective measurement error bounds δi (see eq (9b)). If the absolute value of the deviations is equal to the measurement error bound, the respective constraint is active. The vector prt,max is shown here as an example; it is the vector that contains the maximum value of the parameters rt. H

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article



CONCLUSION We presented an approach for estimation and identifiability analysis of the model parameters in a microgel synthesis process. We approximated the feasible parameter set as a box by solving a series of constrained dynamic optimization problems. The obtained results show that the parameters describing the main propagation reactions as well as a subset of the parameters describing the reactions of the cross-linker can be identified. However, parameters corresponding the propagation reaction of the pure cross-linker cannot be identified. For microgels based on VCL, the results reveal that the available measurements of the enthalpy flow rate via calorimetry and the monomer weight fraction via Raman spectroscopy both actively constrain the feasible parameter set. Hence, to achieve maximum parameter precision, both measurement techniques should be applied during the synthesis. To further constrain the feasible parameter set and hence to increase parameter precision, longer measurement times of the Raman spectrometer should be accepted to achieve a higher accuracy. For microgels based on NIPAM, only enthalpy-transfer rate measurements constrain the feasible parameter set. It should be further investigated why Raman spectroscopy measurements lack precision in comparison to measurements taken during the synthesis of microgels based on VCL. The model with the determined identifiable parameters was validated by comparison with other experimental setups. The agreement is good without further adjustment of the model or the parameters. Nevertheless, parameter precision could potentially be improved by model-based design of experiment. The model can be applied in the model-based design of microgels by dynamic optimization of synthesis experiments. The experimental validation of the optimized experiments is a challenging next step on the path to custom-designed microgels.

For microgels based on VCL, Figure 5 shows that the constraint on the enthalpy-transfer rate ΣR is active in the beginning of the reaction. Also, the constraint on the monomer weight fraction wM1 is active approximately half way through the synthesis. For microgels based on NIPAM, Figure 5 shows that only the constraint on the enthalpy-transfer rate ΣR is active, while the constraint on the monomer weight fraction wM1 is not active. On the basis of these results, the importance of the experimental measurement procedures for the determination of the kinetic parameters can be analyzed and possible suggestions toward the improvement of measurement frequencies and accuracy can be made. First, Figure 5 shows that both measurement sensors may constrain the feasible parameter set in the case of microgels based on VCL. Therefore, both sensors should be used simultaneously to improve the precision of parameters estimated based on the measurements. For microgels based on NIPAM, only the enthalpy-transfer rate sensor constrains the feasible parameter set. This is caused by the low accuracy of the monomer weight fraction measurements in comparison to the good accuracy of the enthalpy-transfer rate measurements for the case of NIPAM (cf. Table 4). A possible reason for the low accuracy of the Raman spectroscopy measurements is attachment of growing microgels to the Raman probe. Another reason may be rooted in the evaluation method used for the determination of the monomer weight fraction based on the raw Raman spectra. Monomer weight fractions are determined by fitting a summation of the spectra of pure components to the Raman spectra of the mixture.23 If the spectra of the pure components are similar, it is more difficult to differentiate them in the raw Raman spectra and the fitting procedure may lead to a lower accuracy. Second, Figure 5 shows that the measurements’ sensitivity to the parameter values changes over the course of the synthesis. For the system based on VCL, the enthalpy-transfer rate is close to the bounding error during the first third of the experiment and the constraint on the weight fraction of monomer is violated only after approximately half of the reaction. This is caused by the correlation of the enthalpy-transfer rate with the reaction rate of the propagation reactions. Hence, a change in a reaction rate constant has a more significant effect on the enthalpytransfer rate, while concentrations of reactants are high in the beginning of the reaction. In contrast, the weight fraction of monomer correlates to the integral of the reaction rates over time. A change in a reaction rate constant that results in a change of the reaction rate is integrated over time and only after some time the respective constraint may be violated. Figure 5 allows for suggestion for the improvement of Raman spectroscopy measurements based on the results for VCL. For Raman spectroscopy, a trade-off between measurement frequency and precision exists.23 Figure 5 suggests that, to further limit the feasible parameter set, the frequency should be lowered to improve the measurement precision, as several constraining measurements are close to being violated at the same time. Methods for the online determination of the measurement duration for the application in microgel synthesis enable the adjustment of the measurement time depending on the signal-to-noise ratio48 and can also be employed with the aim of further constraining the feasible parameter set. However, it should be noted that the results of the above analysis depend on the individual experiment setting and the respective realization of the measurement error and cannot be generalized for all experimental setups.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.8b05274. Explanation of the normalized model (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +49 (0)241 80 97717. Fax: +49 (0)241 80 92326. ORCID

Andrij Pich: 0000-0003-1825-7798 Alexander Mitsos: 0000-0003-0335-6566 Notes

The authors declare no competing financial interest. The model in the Modelica language and experimental data can be accessed at http://permalink.avt.rwth-aachen.de/?id= 720172.



ACKNOWLEDGMENTS We dedicate this article to our friend and colleague Prof. Dominique Bonvin. Among his many accomplishments, he has convincingly argued for appropriate modeling and he has extensively modeled batch systems, including polymerization. We hope he will approve of our model. This work was funded by the German Research Foundation (DFG) within project B4 Microgel synthesis: Kinetics, particle formation and reactor I

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Microgels by Precipitation Polymerization: Process Modeling and Experimental Validation. Ind. Eng. Chem. Res. 2017, 56, 14545−14556. (21) Janssen, F.; Ksiazkiewicz, A.; Kather, M.; Kröger, L.; Mhamdi, A.; Leonhard, K.; Pich, A.; Mitsos, A. Kinetic Modeling of Precipitation Terpolymerization for Functional Microgels. Comput.-Aided Chem. Eng. 2018, 43, 109−114. (22) Kröger, L.; Kopp, W.; Leonhard, K. Prediction of Chain Propagation Rate Constants of Polymerization Reactions in Aqueous NIPAM/BIS and VCL/BIS Systems. J. Phys. Chem. B 2017, 121, 2887− 2895. (23) Meyer-Kirschner, J.; Kather, M.; Pich, A.; Engel, D.; Marquardt, W.; Viell, J.; Mitsos, A. In-line Monitoring of Monomer and Polymer Content During Microgel Synthesis Using Precipitation Polymerization via Raman Spectroscopy and Indirect Hard Modeling. Appl. Spectrosc. 2016, 70, 416−426. (24) Walter, É .; Pronzato, L. Identification of Parametric Models from Experimental Data; Springer: London, 1997. (25) Miao, H.; Xia, X.; Perelson, A.; Wu, H. On Identifiability of Nonlinear ODE Models and Applications in Viran Dynamics. SIAM Rev. 2011, 53, 3−39. (26) Raue, A.; Kreutz, C.; Maiwald, T.; Klingmuller, U.; Timmer, J. Addressing Parameter Identifiability by Model-based Experimentation. IET Syst. Biol. 2011, 5, 120−130. (27) Gabor, A.; Villaverde, A.; Banga, J. Parameter identifiability analysis and visualization in large-scale kinetic models of biosystems. BMC Syst. Biol. 2017, 11, 54. (28) Quaiser, T.; Monnigmann, M. Systematic Identifiability Testing for Unambiguous Mechanistic Modeling - Application to JAK-STAT, MAP Kinase, and NF-kappaB Signaling Pathway Models. BMC Syst. Biol. 2009, 3, 50. (29) Franceschini, G.; Macchietto, S. Model-based Design of Experiments for Parameter Precision: State of the Art. Chem. Eng. Sci. 2008, 63, 4846−4872. (30) Nogueira, I. B.; Pontes, K. V. Parameter Estimation with Estimability Analysis applied to an Industrial Scale Polymerization Process. Comput. Chem. Eng. 2017, 96, 75−86. (31) Villaverde, A. F.; Barreiro, A.; Papachristodoulou, A. Structural Identifiability of Dynamic Systems Biology Models. PLoS Comput. Biol. 2016, 12, e1005153. (32) Jaulin, L.; Walter, E. Set Inversion via Interval Analysis for Nonlinear Bounded-error Estimation. Automatica 1993, 29, 1053− 1064. (33) Gottu Mukkula, A.; Paulen, R. Model-based Design of Optimal Experiments for Nonlinear Systems in the Context of Guaranteed Parameter Estimation. Comput. Chem. Eng. 2017, 99, 198−213. (34) Paulen, R.; Villanueva, M.; Chachuat, B. Guaranteed Parameter Estimation of Non-linear Dynamic Systems using High-order Bounding Techniques with Domain and CPU-time Reduction Strategies. IMA Journal of Mathematical Control and Information 2016, 33, 563−587. (35) Walz, O.; Djelassi, H.; Caspari, A.; Mitsos, A. Bounded-error Optimal Experimental Design via Global Solution of Constrained minmax Program. Comput. Chem. Eng. 2018, 111, 92−101. (36) Jung, F.; Mhamdi, A.; Mitsos, A. 13th International Symposium on Process Systems Engineering; pp 925−931. (37) Tawarmalani, M.; Sahinidis, N. V. A Polyhedral Branch-and-cut Approach to Global Optimization. Mathematical Programming 2005, 103, 225−249. (38) GAMS. https://www.gams.com (accessed 2018-09-01). (39) Goerke, T.; Engell, S. Application of Evolutionary Algorithms in Guaranteed Parameter Estimation. CEC 2016, 5100−5105. (40) Schneider, F.; Balaceanu, A.; Feoktystov, A.; Pipich, V.; Wu, Y.; Allgaier, J.; Pyckhout-Hintzen, W.; Pich, A.; Schneider, G. Monitoring the internal structure of Poly(N-Vinylcaprolactam) Microgels with Variable Cross-link Concentration. Langmuir 2014, 30, 15317−15326. (41) Industries, Wako Pure Chemical, Ltd. V-50 (AMPA) Properties, 2017; https://www.wako-chemicals.de/en/product/v-50. (42) Egea, J.; Henriques, D.; Cokelaer, T.; Villaverde, A.; MacNamara, A.; Danciu, D.; Banga, J.; Saez-Rodriguez, J. MEIGO: An Open-source

modelling of the Collaborative Research Center SFB 985 Functional Microgels and Microgel Systems. The authors also gratefully acknowledge the financial support of the Kopernikus project SynErgie by the Federal Ministry of Education and Research (BMBF). We thank Preet Joy, Leif Kröger, Johannes Faust, Julian Meyer-Kirschner, Jörn Viell, and Michael Kather for valuable discussions.



REFERENCES

(1) Pich, A.; Richtering, W. Microgels by Precipitation Polymerization: Synthesis, Characterization, and Functionalization; Springer, 2011; pp 1− 37. (2) Plamper, F.; Richtering, W. Functional Microgels and Microgel Systems. Acc. Chem. Res. 2017, 50, 131−140. (3) Meurer, R.; Kemper, S.; Knopp, S.; Eichert, T.; Jakob, F.; Goldbach, H.; Schwaneberg, U.; Pich, A. Biofunctional Microgel-Based Fertilizers for Controlled Foliar Delivery of Nutrients to Plants. Angew. Chem., Int. Ed. 2017, 56, 7380−7386. (4) Lohaus, T.; de Wit, P.; Kather, M.; Menne, D.; Benes, N. E.; Pich, A.; Wessling, M. Tunable Permeability and Selectivity: Heatable Inorganic Porous Hollow Fiber Membrane with a Thermo-responsive Microgel Coating. J. Membr. Sci. 2017, 539, 451−457. (5) Guerzoni, L.; Bohl, J.; Jans, A.; Rose, J.; Koehler, J.; Kuehne, A.; de Laporte, L. Microfluidic Fabrication of Polyethylene Glycol Microgel Capsules with Tailored Properties for the Delivery of Biomolecules. Biomater. Sci. 2017, 5, 1549−1557. (6) Ramos, J.; Imaz, A.; Forcada, J. Temperature-sensitive Nanogels: Poly(N-Vinylcaprolactam) versus Poly(N-Isopropylacrylamide). Polym. Chem. 2012, 3, 852−856. (7) Kather, M.; Ritter, F.; Pich, A. Surfactant-free Synthesis of Extremely Small Stimuli-responsive Colloidal Gels using a Confined Impinging Jet Reactor. Chem. Eng. J. 2018, 344, 375−379. (8) Varga, Z.; Wang, G.; Swan, J. The Hydrodynamics of Colloidal Gelation. Soft Matter 2015, 11, 9009−9019. (9) Varga, Z.; Swan, J. Hydrodynamic Interactions Enhance Gelation in Dispersions of Colloids with Short-ranged Attraction and Longranged Repulsion. Soft Matter 2016, 12, 7670−7681. (10) Bonvin, D. Optimal Operation of Batch Reactors - A Personal View. J. Process Control 1998, 8, 355−368. (11) Bonvin, D.; Georgakis, C.; Pantelides, C. C.; Barolo, M.; Grover, M. A.; Rodrigues, D.; Schneider, R.; Dochain, D. Linking Models and Experiments. Ind. Eng. Chem. Res. 2016, 55, 6891−6903. (12) Odian, G. Principles of Polymerization, 4th ed.; WileyInterscience, 2004. (13) Kiparissides, C. Polymerization Reactor Modeling: A Review of Recent Developments and Future Directions. Chem. Eng. Sci. 1996, 51, 1637−1659. (14) D’hooge, D. R.; van Steenberge, P. H.; Reyniers, M.-F.; Marin, G. B. The Strength of Multi-scale Modeling to Unveil the Complexity of Radical Polymerization. Prog. Polym. Sci. 2016, 58, 59−89. (15) Mayo, F.; Lewis, F. Copolymerization. I. A Basis for Comparing the Behavior of Monomers in Copolymerization. J. Am. Chem. Soc. 1944, 66, 1594−1601. (16) Smith, W. V.; Ewart, R. H. Kinetics of Emulsion Polymerization. J. Chem. Phys. 1948, 16, 592−599. (17) François, G.; Srinivasan, B.; Bonvin, D.; Hernandez Barajas, J.; Hunkeler, D. Run-to-Run Adaptation of a Semiadiabatic Policy for the Optimization of an Industrial Batch Polymerization Process. Ind. Eng. Chem. Res. 2004, 43, 7238−7242. (18) Kadam, J. V.; Marquardt, W.; Srinivasan, B.; Bonvin, D. Optimal Grade Transition in Industrial Polymerization Processes via NCO Tracking. AIChE J. 2007, 53, 627−639. (19) Hoare, T.; McLean, D. Multi-Component Kinetic Modeling for Controlling Local Compositions in Thermosensitive Polymers. Macromol. Theory Simul. 2006, 15, 619−632. (20) Janssen, F.; Kather, M.; Kröger, L.; Mhamdi, A.; Leonhard, K.; Pich, A.; Mitsos, A. Synthesis of Poly(N-Vinylcaprolactam)-Based J

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Software Suite based on Metaheuristics for Global Optimization in Systems Biology and Bioinformatics. BMC Bioinf. 2014, 15, 136. (43) Yuret, D.; de La Maza, M. Dynamic Hill Climbing: Overcoming the Limitations of Optimization Techniques; 1993. (44) DyOS RWTH Aachen University (AVT.SVT). http:// permalink.avt.rwth-aachen.de/?id=295232 (accessed 2018-09-01). (45) Schlegel, M.; Marquardt, W.; Ehrig, R.; Nowak, U. Sensitivity Analysis of Linearly-implicit Differential-algebraic Systems by One-step Extrapolation. Applied Numerical Mathematics 2004, 48, 83−102. (46) FMI standard - Functional mock-up interface for model exchange and co-simulation. http://fmi-standard.org (accessed 2018-05-01). (47) Dassault Systems. https://www.3ds.com/products-services/ catia/products/dymola/ (accessed 2018-06-01). (48) Meyer-Kirschner, J.; Kather, M.; Ksiazkiewicz, A.; Pich, A.; Mitsos, A.; Viell, J. Monitoring Microgel Synthesis by Copolymerization of N-isopropylacrylamide and N-vinylcaprolactam via In-Line Raman Spectroscopy and Indirect Hard Modeling. Macromol. React. Eng. 2018, 12, 1700067.

K

DOI: 10.1021/acs.iecr.8b05274 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX