Reverse Engineering of Radical Polymerizations by Multi-objective Optimization

1. Introduction

Radical polymerizations are known to be very robust and to provide access to a wide range of polymers with largely differing properties, which are defined by the process conditions. The strong correlation between production process and material properties is due to the complex reaction mechanism consisting of a large number of elemental reactions even for homopolymerizations with a single monomer [1]. The kinetics of the elemental reactions are strongly dependent on the process conditions. Therefore, the prediction of a suitable radical polymerization process to obtain a polymer with targeted properties is challenging. To allow for on demand polymer synthesis, at first sight it appears highly attractive to apply simulations of polymerization processes, e.g., employing differential equations [2,3,4] or kinetic Monte-Carlo (kMC) methods [5,6,7,8]. Simulations are particularly valuable, because detailed information on polymer microstructure at each time moment is available, which is not accessible from polymerization processes. However, this type of simulations cannot be run backwards and on demand suggestion of polymerization conditions to obtain a pre-defined polymer is not feasible. To overcome this issue, reverse engineering has the potential to provide several solutions, as opposed to a single one to one relation between polymerization variables and microstructural properties [9]. From the input variables, a polymerization process model predicts the concentration vs. time profiles and the polymer properties. Inverse modeling, on the other hand, is more difficult and calls for optimization strategies. In order to determine the optimal input values for systems with complex reaction mechanisms that provide pre-defined reaction outputs (such as pre-set conversion, yield, and/or other product properties), it was suggested to intelligently explore the reaction condition search space [10]. Further, it was proposed to solve reverse engineering problems using machine learning (ML) based prediction [11], in which ML regression models based on the random forest algorithm and a multivariate and multi-target regression problem [12] were applied. The model took a targeted MMD and predicted the initial polymerization conditions to produce polymer with this targeted property, minimizing the errors in predicted recipes. However, the prediction of the recipe was based on the MMD only, and not optimized with respect to multiple objectives like reaction time or conversion.

Pareto or multi-objective optimization (MOO) [13] is the process of maximizing or minimizing many objective functions while taking a set of constraints into consideration. Numerous scientific domains, such as engineering [14], economics [15], and logistics [16], need the use of MOO when making optimum decisions facing trade-offs between two or more competing objectives. An increasing use of MOO has been seen in chemical engineering [17,18]. In 2009, Fiandaca et al. [19] used genetic algorithm-based MOO to optimize a pressure swing adsorption process based on the maximization of two objectives: nitrogen recovery and nitrogen purity. In 2013, Ganesan et al. [20] carried out MOO of the combined carbon dioxide reforming and partial-oxidation of methane with respect to three objectives. MOO utilized a gravitational search algorithm and particle swarm optimization to tackle the problem. With respect to technical application it appears highly important to solve the reverse engineering problem as an optimization task with multiple contradicting objectives [13]. Various ML-based optimization strategies were addressed for the purpose of reverse engineering polymerization processes [21,22,23]. A genetic algorithm-based optimizer is proposed by Mohammadi et al. [9] to generate a variety of polymerization recipes at random and to send them to the kMC simulator for error evaluation.

This study provides a MOO-based reverse engineering approach, not only ensuring that the targeted MMD is obtained by means of minimizing the mean squared error (MSE), but also providing minimal reaction time and maximal monomer conversion. Frequently, there are several recipes for obtaining similar MMDs, which are referred to as candidate recipes in the following. The solution of the MOO approach is referred to as polymerization recipe, which includes temperature, reaction time as well as initial monomer and initiator concentrations. In the proposed MOO approach, the data for the selection of the optimal recipes from the search space is based on kMC simulations, as previously reported for training ML models [11]. The simulation provides the dependence of monomer conversion and the corresponding MMD on the polymerization conditions and reaction time. The proposed reverse engineering algorithm consists of several steps. First, the kMC simulations are run for the selected recipe search space to derive the MMDs and monomer concentrations as an input for the MOO step. Then, MOO is applied for a given target MMD and the Pareto optimal space is found on the base of the search space. In the last step, a weighted sum optimization function is used to calculate the weighted score of each candidate recipe, which is used for evaluating the solutions. The best candidate has the smallest score. To accelerate the MOO procedure an additional search space clustering on the basis of MMDs is considered. The approach pursued is illustrated in Figure 1. The method is tested for the following model system: chemically initiated vinyl acetate (VAc) radical polymerization model system at 60 °C using tert-butyl peroxypivalate as initiator.

2. Reverse Engineering Modeling Approach

2.1. Model Development

This section provides a formal description of solving a reverse engineering problem with an MOO approach. State of the art MOO methods [24] like genetic algorithms require to generate and evaluate new recipes in each step. Due to the high number of required steps, it is very demanding to perform on-line kMC simulations in loop [9].

Instead, here a recipe search space R is selected, where each recipe r ∈ R consists of reaction time t, the initial monomer concentration c_m,0, and the initial initiator concentration c_ini,0. Then, the corresponding monomer concentration c_m(r) and molar mass distribution MMD(r) are obtained for each r ∈ R via kMC simulation. In the MOO approach, the input is a target MMD, MMD^target, and the output is a set of optimal candidate recipes R^*: ${\mathbf{r}}_1^{*}$,${\mathbf{r}}_2^{*}$, …,${\mathbf{r}}_n^{*}$. R^* is a subset of the recipe search space R, R^*⊂ R with MMD(${\mathbf{r}}_i^{*}$) being close to MMD^target as evaluated on the basis of the MSE as well as the maximal conversion and minimal time.

The optimization variables are presented in Table 1. The lower and upper limits for the variables c_mon,0, c_ini,0, and t are defined by the simulated data. In Table 2 for a specific recipe r, the simulated values c_m(r) and MMD(r) are used to calculate the values of three optimization objective functions: objective for the reaction time, f_t(r), and objective for the mean squared error (MSE) between MMD^target and the predicted MMD(r), f_MSE. f_cm(r) is used to turn the maximization problem of monomer conversion f_conv(r) = (c_m,0 − c_m(r)) / c_m,0 into a minimization problem, in which 1 − f_conv(r) is minimized (Table 2).

The final decision takes user preferences into account by assigning specific weights to the objectives applying the weighted sum method [24,25]. Thus, the multi-objective function can be represented in a single-objective way. Then, the values of this function are calculated for each candidate recipe and the recipes with minimal values of the objective function are selected as a set of optimal solutions. A weight w_i is assigned to each normalized objective function f_i as follows:

$$\underset{\mathit{r}}{\min }f\left(\mathit{r}\right)=\underset{\mathit{r}}{\min }\sum _i{w}_i{\bullet f}_i\left(\mathit{r}\right)$$

(1)

where $\sum _i{w}_i=1$, i ∈ {MSE, cm, t}, r ∈ R, and R is a polymerization recipe space. For clarity of presentation, the weight w_cm of the objective f_cm is also referred to as the weight of the conversion objective. If required the number of the objectives can be increased, e.g. the conversion of initiator can be added.

The steps of the proposed algorithm are presented in Figure 2 as direct approach. First, a search space R^S ⊂ R is selected (for details see Section 2.2) and for each r ∈ R^S MMD(r) and c_m(r) are obtained via kMC simulations. Then, MOO is performed over R^S. First, the objective function values are calculated as follows: the values for the objective f_t(r) are already included in r as t, the values of the objective f_cm(r) are calculated in advance for all possible r according to Table 2, and the values of the objective f_MSE are specified by MMD^target. Further, based on the calculated objective function values, the Pareto front points R^par ⊂ R^S leading to MMD^target are identified. For this, the points from R^S are represented in the Pareto optimal space, with coordinates specified by the three objectives. The Pareto front points are found in this Pareto optimal space, such that one value of the objective function cannot be improved without downgrading the value of another objective function. Finally, for the Pareto front points R^par, the weights of each objective function are defined and a set of the best recipe candidates R^*⊂ R^par according to Eq. 1 is selected.

The above-described algorithm is improved with respect to the optimization time by means of clustering the search space, which is illustrated in Figure 2 as the clustering-supported approach. Clustering divides the search space R^S into a number of clusters as illustrated in Figure 3. First, the search space R^S is clustered on the base of MMD(r), which allows for selecting a cluster R^target ⊂ R^S containing the MMDs, which are the closest to MMD^target. In general, a larger number of clusters leads to a smaller number of MMDs per cluster gaining higher similarity of the distributions. However, there is less space for optimization regarding other objectives, e.g., such as polymerization time and monomer conversion. For this reason, an appropriate trade-off between the number of clusters and its size has to be identified. Upon appropriate clustering a cluster for the target MMD (Figure 2, red arrows) R^target is found. Then, by MOO the search space is reduced to the number of Pareto front points R^par ⊂ R^target. Finally, after defining objective weights, the best recipe candidates R^* ⊂ R^par are found according to Eq. 1. Since MOO is applied to a single cluster R^target ⊂ R^S, which is considerably smaller than R^S, the optimization time is significantly reduced.

Different methods can be applied for clustering the search space on the basis of MMD. Clustering of distributions and their representation in histograms is an important topic, which attracted a lot of attention, because of specific metrics, which should be used to compare the distributions. One of the most popular and fast clustering methods is the kMeans method [26]. A modified kMeans clustering algorithm was applied to the clustering of histograms [27]. Further, a novel non-parametric clustering algorithm of empirical probability distributions was proposed [28]. Here, the classical kMeans clustering method was used. This algorithm starts with a random separation of the MMDs into clusters. At each step, it recalculates the centroids of each cluster and relocates the data points to the new centroids. The clustering process finishes when the clusters are stable or the given number of iterations is reached. In this study, the simplest Euclidean distances are used for calculation of the distances between multi-dimensional data points, while specific metrics for clustering of distributions are also available [27,29].

There are different strategies for the data generation for the MOO procedure: use of exclusively in advance generated data from kMC simulation (Figure 2, blue arrows), on demand kMC-generated data, ML-generated data, kMC-based and ML-generated hybrid data sets, etc. Currently, as a first step, the focus is exclusively on the use of kMC-simulated data.

2.2. Data acquisition and Processing

The in-house developed kMC simulator mcPolymer was used to carry out the simulations required for the generation of polymerization data according to the search space [30]. The simulator allows for exporting the concentration profiles of all reactants and products as well as microstructural data like MMD, chain composition and branching of all polymeric species involved in the process. The simulator output was adapted to be well-machine readable. The data were filtered, further abstracted, logically connected, and stored in the well-structured no-SQL database MongoDB. The kMC simulated MMDs and monomer concentrations were obtained for the selected search space R^S, allowing for MOO and subsequent finding of the weighted optimal solution.

The kMC simulations were performed for radical polymerizations with VAc as monomer, tert. butyl peroxypivalate as initiator, and methanol as solvent. The simulations are based on a full kinetic model for the VAc radical polymerization containing all elemental reactions [31]. The following polymerization conditions were used: constant temperature of 60 °C, c_ini,0 in the range of 1.0 to 20.0 mmol·L⁻¹, and c_m,0 in the range of 2.0 to 5.0 mol·L⁻¹ with uniformly distributed grid size of c_ini,0 (geometrically scaled grid points) and c_m,0 (arithmetic scaled grid points), resulting in 225 simulations of the process. The geometric scale was selected for c_ini,0 to put more attention on the small values of this parameter. The polymerization process was simulated for a constant reaction time of 6 hours and the properties of interest were recorded every 20 minutes, thus, obtaining in total 18 data points at different time moments for each investigated property. Thus, the data set contains 4050 different MMDs.

For training the ML prediction models for single-objective reverse engineering the obtained data set was divided into training and test set in proportion of 80:20. The same data was used for MOO, again taking 80 % of the data as training set for the search space R^S and 20 % of the data as test set R^test. The test set R^test contains 810 recipes, which correspond to a set MMD^test consisting of 810 kMC-simulated MMDs. The evaluation of all optimization approaches was performed with MMD^test with each element serving as MMD^target. In order to test the MMO approach a single MMD is selected from MMD^test and used as MMD^target for the MMO approach. The performance of the MMO is evaluated by testing it with every MMD from MMD^test.

3. Results and Discussion

In order to compare results obtained via the previously described ML modeling-based reverse engineering strategy [11] with data from the MOO approach introduced in this work, the technique reported for butyl acrylate was adopted for the polymerization of vinyl acetate. Previously, it was described how an ML regression model provides a recipe r = (c_m,0, c_ini,0, T) for a fixed time for a given MMD^target. c_m,0 and c_ini,0 are initial concentrations of the monomer and the initiator, respectively, and T is the polymerization temperature. The prediction was performed with the random forest method for a target MMD to minimize the errors in the predicted recipe r(MMD). In the case of VAc polymerizations, the temperature is kept constant at T = 60 °C and the reaction time t was predicted by the model: r = (c_m,0, c_ini,0, t). The reverse engineering result for a sample MMD^target presented in Figure 4 corresponds to the recipe (c_m,0; c_ini,0; t)/ (2.0 mol/L; 1.9 mmol/L; 140 min). Only very small differences between the target and predicted MMD as well as the predicted recipe are seen. The overall performance of the developed model was evaluated using the R² determination coefficient, which is a measure for the quality of the prediction by the ML model. Keeping in mind that R²=1 indicates perfect prediction, it is remarkable to note that a rather low value of 0.78 for R² is obtained while visual inspection of the MMDs for the representative example given in Figure 4 indicates only minor differences. However, this ML-based approach does not consider multiple contradicting optimization objectives and does not provide multiple alternative solutions as proposed in the Pareto optimization.

3.1. Direct Pareto Optimization

Figure 5A shows all Pareto front points with three objectives f_i, i ∈ {MSE, cm, t}. After definition of objective weights, the optimal solution from the Pareto front points is selected, which satisfies the weighted score calculated with Eq. 1 best. Several combinations of objective weights are considered: time focused, conversion focused, MSE focused. For example, in the time focused case, the weight of time is higher than the other weights. To identify the optimal recipe leading to the target MMD, the MSE should have a higher weight than the other objectives. To reduce the number of Pareto front points an MSE limit is chosen. As an example, the data in Figure 3 was obtained with an MSE limit of 2∙10⁻³.

The optimization procedure is illustrated in Figure 5 for an example target MMD, which corresponds to the recipe (c_m,0; c_ini,0; t) / (2.0 mol/L; 1.9 mmol/L; 140 min). Rather than presenting the objective f_cm, which represents the minimized monomer concentration, the conversion is shown in Figure 5A,B. Figure 5B demonstrates the reduction of the number of Pareto front points from 340 points (A) to 46 points (B) by filtering with the MSE limit of 2∙10⁻³. The MMDs of all Pareto front points and of the filtered Pareto front points are given in Figure 5C,D, respectively. The filtering approach allows for considering only MMDs of similar shape. Figure 6 shows the impact of different MSE limits on the number of Pareto front points. Even with an MSE limit of 10⁻³ the average number of points is around 30, which is a relevant number for consideration of other objectives. The minimal number of filtered Pareto points for the MSE limit of 10⁻³ is about 10.

Figure 7 allows for comparison of MSEs obtained via ML prediction-based (left) and MOO-based approaches (right) for reverse engineering calculated for all elements of MMD^test. Minimal MSEs are preferable. For the ML-based prediction (left), the MSE values of almost all MMD^target from MMD^test are less than 10⁻³, while most MSEs obtained by the MOO-based approach with an MSE limit of 2∙10⁻³ are by two orders of magnitude lower. Most values are less than 10⁻⁵.

The next step is to evaluate the candidate recipes based on the weighted sum of the objective values according to Eq. 1. Figure 8 and Table 3 present the candidate recipes for the different combinations of objective weights indicated. In Figure 8, a color code is assigned to the individual candidate recipes according to the score of the weighted sum of the values of the three objective functions f_MSE, f_cm, f_t calculated according to Eq. 1. For the weights w_MSE = 0.9, w_cm = 0.02, and w_t = 0.08 in Figure 8A the focus is on the MSE. As a minimization problem is solved, a minimal score value (visualized by big pink points) determines the best candidates, given by the recipes with IDs 0, 15, 16 in Table 3. Other weight combinations with focus on the monomer conversion in Figure 8B and with equal weights in Figure 8C demonstrate, which recipe candidates are selected depending on specific requirements.

3.2. Clustering Supported Pareto Optimization

The clustering-supported optimization is performed for the same MMD^target as in the previous section. It is assumed that all points contained in the target cluster R^target have satisfactory MSEs. For this reason, only the two objectives f_cm and f_t have been chosen to be optimized. Figure 9 shows R^target in the coordinate space, which corresponds to the objective functions. The Pareto front points R^par ⊂ R^target are indicated by colored markers. Table 4 presents all Pareto front points with the corresponding recipes and objective weights.

To find the weighted solution for a time focused result with w_cm = 0.2 and w_t = 0.8 the candidate with ID 1 reaching a conversion of 47 % at a reaction time of 40 min is best. If conversion is considered to be more important than reaction time as in the case of w_cm = 0.8 and w_t = 0.2, the candidate with ID 4 is the best solution leading to a conversion of 77 % with a reaction time of 100 min.

MMD clustering focuses only on a single objective function f_MSE when identifying the most similar cluster for the target MMD. The accuracy is determined on the basis of the MSE, which represents the deviation of the selected solution from MMD^target and depends on the cluster size, and consequently, on the number of clusters. For a small number of clusters, every cluster contains a broader spectrum of MMDs as compared to a large number of clusters and leads to an increase in the MSEs of MMD^target and the MMDs contained in this cluster. However, the cluster covers a larger value range of the objective function, allowing for a broader range of predicted recipes. Figure 10 depicts how the MMD^target can be reached for clustering the search space into 20, 40, and 60 clusters, respectively. Figure 10A–C shows the points of the target clusters in the objective function space. The shape of the clusters is similar. As the number of clusters grows, the number of candidates per cluster decreases, because there are fewer points fitting the required accuracy with respect to the MMD. Therefore, the size of the whole cluster shrinks. The red points in Figure 10 mark the Pareto front points, from which the optimal candidates are selected.

Figure 10D–F shows the MMDs for all Pareto front points. With increasing number of clusters, the MSE of MMD^target and Pareto front MMDs decreases. The objective function value range also decreases as illustrated in Table 5. For the example, given f_MSE of MMD^target improves by a factor of 100 when going from a total number of 20 clusters to a total number of 60 clusters. This improvement is achieved by a loss in conversion, which is reduced from 58 % to 27 %. In this example, reaction time is also lowered with increasing number of clusters. The last two rows in Table 5 provide the scores calculated by Eq. 1 for the best recipes obtained with clustering of the search space into 20, 40, and 60 clusters for two cases with a different number of objectives. Note, that only the scores obtained with the same combination of weights (given in a single row of Table 10) are comparable. For equal weights of two objectives (w_cm = 0.5, w_t = 0.5) the minimal score of 0.39 was obtained with clustering into 40 clusters. Although the optimization has been carried out within one cluster using only two objectives, the results can be related to the direct Pareto optimization, including the MSE objective. It is shown that for equal weight of three objectives (w_MSE = 1/3, w_cm = 1/3, w_t = 1/3) the minimal score of 0.30 was also obtained with the same clustering into 40 clusters.

The comparison of the direct and the clustering-supported Pareto approach is based on the MSE limit. For the direct approach, the MSE limit can be chosen by the operator. For the clustering-supported approach, the individual MSE limit for a given MMD^target, which is evidently depending on the number of clusters, is defined as the maximal MSE value for all MMDs contained in R^target.

Both approaches were compared on the basis of MMD^test. Then, for the clustering supported approach, the MSE limit is defined as the maximal value of all individual MSE limits over MMD^test. Figure 11 shows the dependency of the MSE limit on the total number of clusters. For more than 30 clusters, the MSE limit is almost constant at a level of 0.003. The MSE limit increases rapidly for less than 30 clusters.

Figure 12A illustrates the average number of Pareto front points depending on the MSE limit. Increasing the MSE limit results in a reduction in the number of clusters (see Figure 11) associated with an enlargement in cluster size, and thereby leading to a higher count of Pareto front points. The clustering-supported approach is more selective considering the average number of Pareto front points. The standard deviation is used to describe the value ranges of the objective function and of the MSE function. Wider ranges are advantageous, because they provide more space for optimization. Averaging the value ranges of each objective function over MMD^test yields the average ranges illustrated in Figure 12B–D. The value ranges of the reaction time and conversion are wider (Figure 12C,D) for the clustering supported approach, while the MSE value range is more restricted (Figure 12B). The findings indicate that the clustering approach allows for a better optimization of monomer conversion and reaction time while the direct approach is better suited for optimization of the MSEs.

Table 6 presents the execution time for both approaches for the test set MMD^test consisting of 810 target MMDs. For the clustering-supported approach, the optimization time depends slightly on the number of clusters. However, it is by a factor of ten smaller than for the direct approach.

The rapidness of the MOO algorithm is especially important, for the next step, when in future genetic algorithms together with ML models will be applied to extend the search space for the polymerization parameters, because Pareto optimization should be executed many times for each evolution stage of genetic algorithms. In this situation, the clustering supported Pareto optimization is preferable.

4. Conclusions and Outlook

The study bridges the gap between multi-objective optimization (MOO) and its application for reverse engineering of polymerization processes. The proposed MOO models allow for simulation-supported determination of an optimal recipe for a targeted molar mass distribution, taking multiple objectives, e.g., minimal reaction time, maximal conversion of monomer, and similarity of the found MMD compared to the target MMD into consideration. The proposed approach can be accelerated by an additional clustering of the simulated MMDs. Moreover, it is possible to obtain a number of suitable candidates considering different weights of the selected objectives. A set of alternative optimal solutions is obtained for each specific combination of weights. First insights are provided on the formulation of a polymerization reverse engineering task as a MOO problem. In future, the understanding gained in combination with previously proposed ML models [11] for polymerizations will be applied to facilitate and to speed-up finding of optimal solutions with a limited number of kMC simulations. ML models can be used to set-up the search space and to evaluate the candidate polymerization procedures rather than using kMC simulations. The search for candidate solutions will be performed with the help of genetic algorithms. Moreover, the consideration of more complex microstructural details, e.g., such as branching in high-temperature polymerization of acrylate [32], may require to increase the number of objectives in MOO.