Molecular Weight Distribution Control for Polymerization Processes Based on the Moment-Generating Function

The molecular weight distribution is an important factor that affects the properties of polymers. A control algorithm based on the moment-generating function was proposed to regulate the molecular weight distribution for polymerization processes in this work. The B-spline model was used to approximate the molecular weight distribution, and the weight state space equation of the system was identified by the subspace state space system identification method based on the paired date of B-spline weights and control inputs. Then, a new performance criterion mainly consisting of the moment-generating function was constructed to obtain the optimal control input. The effectiveness of the proposed control method was tested in a styrene polymerization process. The molecular weight distribution of the styrene polymers can be approximated by the B-spline model effectively, and it can also be regulated towards the desired one under the proposed control method.


Introduction
Polymers are high molecular compounds consisting of macromolecules with different chain lengths [1]. The properties of polymers are affected by their molecular weight distribution (MWD) [2][3][4][5][6], which can be measured directly or modelled mathematically [7][8][9]. The production of polymers is widespread and long-standing in industrial processes; different applications need different specifications for the polymers [10]. To meet these demands, many performance indices consisting of polydispersity and average molecular weight have been proposed. Under these indices, the control objective is commonly dealt with as an optimization problem, which has been investigated in many studies [11][12][13][14][15][16][17][18][19][20][21]. Though satisfactory control performance could be obtained in these studies, these indices cannot fully characterize the entire molecular weight distribution of polymers.
However, the shape of the molecular weight distribution of polymers is essential in many polymerization processes, such as paints and paper coatings [22,23], and the average molecular weight and polydispersity index is unable to reflect the characteristics of the MWD when the distribution is non-Gaussian. To solve this problem, the stochastic distribution control (SDC) method can be proposed for molecular weight distribution shaping in polymerization processes. The SDC method has been applied to many industrial processes, such as distribution control of the flame temperature in furnace systems [24,25], particulate process control in powder industries [26,27], distribution control of bubble size in flotation processes [28], distribution control of crystal size in crystallization processes [29][30][31] and power probability density function control in nuclear research reactors [32].
Unlike the methods based on performance indices consisting of polydispersity and average molecular weight, the entire molecular weight distribution (MWD) of the polymerization process is able to be regulated towards the target MWD under the SDC method.
Since the MWD cannot be easily used in the MWD shape control of polymerization processes directly, the dynamic modelling of the output MWD as the control input is necessary. B-spline models can be utilized to approximate the output MWD [33]. The main advantage of the B-spline MWD model is the decoupling of the time domain, the MWD definition domain in the formulation and the obtained weights model of the MWD, which can be made mathematically equivalent to any existing physical models for MWD systems subjected to a pre-specified, small modelling error. Although several B-spline-based MWD models have been developed [34][35][36][37], the linear B-spline MWD model is applied in this paper due to its computational simplicity, which can reduce the computational time [38].
Next, a dynamic weights model can be established to make the subsequent control problem easier to manage. Inspired by Greś and colleagues [39], the subspace identification method is used to obtain the state space model between the weight vector of the B-spline model and the control input. Based on the input and weight vector data of the B-spline model, the state space model parameters of the system can be obtained through the row subspace and column subspace mapped by the Hankel matrix. Canonical variate analysis, multivariable output error state space and numerical subspace state space system identification are three influential subspace identification methods. Compared to the traditional identification method, numerical subspace state space system identification (N4SID) requires less calculation and satisfies numeric stability [40]; therefore, the N4SID was adopted in this paper.
Based on the weights model, control algorithms can be derived by minimizing performance criteria, which indicate the difference between the output MWD and the target MWD. As the characteristic of the MWD can be reflected by its moment-generating function (MGF), a new performance criterion using MGF based on the state space model was proposed in this paper. Compared to the traditional performance criteria [33][34][35][36][37][38], the new criterion needs no integral operation on the quadratic error and has less dependence on the regulation of the criterion weights. The effectiveness of the proposed control method was introduced in a styrene polymerization process.
In this paper, a novel control algorithm based on the moment-generating function for polymerization processes was presented. The B-spline model was used to approximate the MWD, and N4SID was introduced for system modelling. Then, a performance criterion mainly consisting of the moment-generating function of the MWD was constructed, and the control law could be obtained by minimizing the performance criterion. This control algorithm was applied to a styrene polymerization process and simulation experiment to verify its effectiveness.

B-Spline Model and Subspace Identification Method
For a polymerization process, the monomer is transformed into a polymer through a series of chemical reactions. The free radical polymerization process includes chain initiation, chain growth, chain transfer and chain termination. The reaction forming monomer free radicals is called the chain initiation reaction, and it is the key to controlling the rate throughout the polymerization reaction. An initiator added to the monomer is heated to generate initiating radicals, and once the free radical monomer is generated, it will immediately be added to the second monomer to generate a free radical chain containing two monomer units. The activity of the free radical chain does not decay; it will immediately be added to the third monomer and subsequent monomers in the chain. Then, the degree of polymerization of chain radicals increases rapidly; these reactions that increase the degree of polymerization are called chain growth reactions. The binding reaction of free radicals is called a termination reaction, which has two forms: coupling and disproportionation. In addition, the chain free radicals undergo a transfer reaction with monomers, initiators or formed macromolecules, whereby a polymer chain is terminated, and a new radical capable of propagating is formed. The schematization of proposed reaction mechanisms is shown as Figure 1. Entropy 2022, 24, x FOR PEER REVIEW 3 radical capable of propagating is formed. The schematization of proposed reaction m anisms is shown as Figure 1. The dynamic MWD of the polymerization process can be obtained by the M model described in [33] or measured directly through online techniques, and it cann easily used in MWD shape control directly. Therefore, a B-spline approximation mo established through the obtained MWD data. Through the weight data of the B-s model and the control input, the output MWD model can be identified by the N method, which is used in MWD control described later. The schematic diagram o proposed control algorithm is shown in Figure 2. By using the moment-generating tion (MGF), the pseudo-state vector of the target MWD and approximated MWD c obtained as ref z and z , respectively. The control input, u , is regulated by the MWD troller for the shaping of output MWD.    Since the integration of ( , ) k yu  is equal to 1, the number of independent expa weights is n − 1, and the remaining expansion weight can be expressed by other expa weights.
Defining The dynamic MWD of the polymerization process can be obtained by the MWD model described in [33] or measured directly through online techniques, and it cannot be easily used in MWD shape control directly. Therefore, a B-spline approximation model is established through the obtained MWD data. Through the weight data of the B-spline model and the control input, the output MWD model can be identified by the N4SID method, which is used in MWD control described later. The schematic diagram of the proposed control algorithm is shown in Figure 2. By using the moment-generating function (MGF), the pseudo-state vector of the target MWD and approximated MWD can be obtained as z re f and z, respectively. The control input, u, is regulated by the MWD controller for the shaping of output MWD. The dynamic MWD of the polymerization process can be obtained by the MWD model described in [33] or measured directly through online techniques, and it cannot be easily used in MWD shape control directly. Therefore, a B-spline approximation model is established through the obtained MWD data. Through the weight data of the B-spline model and the control input, the output MWD model can be identified by the N4SID method, which is used in MWD control described later. The schematic diagram of the proposed control algorithm is shown in Figure 2. By using the moment-generating function (MGF), the pseudo-state vector of the target MWD and approximated MWD can be obtained as ref z and z , respectively. The control input, u , is regulated by the MWD controller for the shaping of output MWD.

B-Spline Approximation for Output MWD
Consider a polymerization process where

B-Spline Approximation for Output MWD
Consider a polymerization process where u k ∈ d×1 is the manipulated variable to control the shape of MWD, defined as γ(y, u(k)). This can be approximated by the B-spline neural network where u k represents the control input at time instant k, B i (y)(i = 1, · · · , n) stands for the basis functions which have been pre-designed on the domain of [a, b], n represents the number of basis functions and ω i (u k ) denotes the expansion weight. e 0 is the approximation error which can be ignored for expression simplification.
Since the integration of γ(y, u k ) is equal to 1, the number of independent expansion weights is n − 1, and the remaining expansion weight can be expressed by other expansion weights.
Defining Then, the MWD γ(y, u k ) can be reformed as Multiplying C(y) T to both sides of Equation (2) and carrying out an integral operation over the definition domain, the expansion weights vector v k can then be calculated as Since the dynamics of the output MWD can be represented by Equation (1), we assume that the input u k and expansion weights vector v k satisfy the following state space form: where x k ∈ (n−1)×1 denotes the state vector and A, B, C and D stand for the coefficient matrices of the system, which can be identified by the subspace state space system identification method in the next part.

Subspace State Space System Identification
The relationship between control input u k and expansion weight vector v k can be represented by Equation (4). This paper employs the N4SID method to identify the coefficient matrices A, B, C and D [40]. Compared with traditional methods, N4SID requires less calculation and has a high modelling accuracy.
To save computing space and speed up computer processing, some Hankel matrices are constructed as follows: where i is the row number that should be larger than the order of the identified system, j is the column number, and 2i + j − 2 should be smaller than the length of the data.
To simplify the derivation of the N4SID method, some matrices of input and expansion weights are defined as Similarly, the state matrices can be defined as The generalized observability matrix (Γ i ) and generalized controllability matrix (∆ i ) of the identified system are shown as Entropy 2022, 24, 499

of 13
A Toeplitz matrix is defined as Then, the state space Equation (4) can be reformulated by the Hankel matrices as Substituting Equations (12) and (13) to Equation (14), X f can be represented as and Γ † i is the Moore-Penrose pseudo-inverse of the generalized observability matrix.
Substituting Equation (15) to (13), V f can be denoted as Making an oblique projection of V f to W p along the direction of U f where O i is the oblique projection. Based on the assumption that System (4) is observable and controllable, rank(O i ) = n, and then the singular value decomposition of O i can be shown as where S 1 denotes a diagonal matrix consisting of singular values far greater than 0 and S 2 denotes a diagonal matrix composed of singular values close to 0. From Equations (17) and (18), Γ i and X f can be estimated as Then, the coefficient matrices of System (4) can be obtained by least squares (LS) Therefore, the estimated expansion weights vector (v k ) can be obtained by wherev k is the estimated expansion weights vector, andx k is the estimated state variable vector of the dynamic weights model.

MWD Shaping Control Algorithm Using MGF for Polymerization Processes
We can note that the MWD is uniquely determined by its moment-generating function. We can carry out a sampling operation on the moment-generating function of the MWD, and the complete information of the MWD can be reflected if the sampling points are sufficiently large. The shape of the MWD approaches the target MWD when the sampling points of the output MWD approach those of the target MWD. Then, a performance criterion using the moment-generating function can be constructed based on the identified dynamic weights model, and the control algorithm can be derived by minimizing the criterion.

Moment-Generating Function
For the traditional performance criterion [33][34][35][36][37][38], if the integral value of the quadratic error between the output MWD and target MWD is very small, then the regulation of criterion weights is critical. In order to reduce the impact of criterion weights on control performance, a moment-generating function is proposed in this section to construct the performance criterion of the proposed control method. The moment-generating function can be represented as where ξ ∈ R 1 and γ(y, u k ) is the output MWD at time k. By expanding (23), the following equation can be obtained as follows: where m s denotes the sth moment. For the convenience of the following calculation, the cumulant-generating function φ k (ξ) can be defined as φ k (ξ) = logM k (ξ) (25) The estimated value of φ k (ξ), denoted asφ k (ξ), can be calculated utilizing the obtained output MWD. By selecting ξ 1 , ξ 2 , . . . , ξ m , a pseudo-state vector (z k ) can be denoted as follows: Note that enough characteristics of the output MWD can be reflected if m is sufficiently large.
Similarly, the pseudo-state vector of the target MWD can be denoted as Then, the following performance criterion can be constructed for the controller design: where Q and R denote the weights, and the second term denotes energy constraints on the control input. The first term on the right side of Equation (28) can be reformed as which indicates the difference between the output MWD and target MWD.

Control Algorithm
As shown in Performance Criterion (28), it is obvious that the output MWD approaches the target MWD when the performance criterion decreases. In practical polymerization processes, the constraint on the control input is usually taken into consideration in the case of actuator saturation.
Then, the optimal control input can be solved as follows where U min and U max stand for the lower and upper bounds of the control input, respectively.
Since the analytical solution of the nonlinear programming optimization problem, (30), cannot be easily solved, the control input can be obtained by quadratic programming in this paper.

Simulation Study
In this section, the performance of the proposed controller was tested in a styrene polymerization process. Styrene was the monomer for polymerization, and azobisisobutyronitrile was used as the initiator. The monomer and the initiator were blended and pumped into a tank reactor where the polymerization occurred, after which the styrene polymers were produced. Since the ratio of the monomer to the initiator is essential and critical to the polymerization process in the reaction tank, the ratio of the flow rate of the monomer to the sum flow rate, defined as c, is regarded as the manipulated variable in the MWD control system.

Modelling of the Molecular Weight Distribution
Since the molecular weight distribution in the styrene polymerization process is not easily measured directly, dynamic modelling of the molecular weight distribution in the styrene polymerization process is presented in this section.
Under the assumption that the sum flow rate into the tank reactor is invariant, the dynamic model of molecular weight distribution can be established by following the work in [33]. To uphold conciseness of the description, only the major differential equations are presented.
The concentrations in the tank reactor can be described by the following mass balance equations: where C I and C M stand for the concentrations of azobisisobutyronitrile and styrene, respectively. C I0 and C M0 denote the initial concentrations of azobisisobutyronitrile and styrene, respectively. κ is a constant which represents the average residential time in the tank reactor. K d , K i , K p and K ct stand for different rate constants in the polymerization process. φ 0 denotes the radical concentration. Through the use of the generation function technique, the following differential equations can be deduced: where φ 1 and φ 2 are leading moments for radicals, and K t stands for the termination rate constant. Similarly, the differential equations concerning the leading moments for dead polymers can also be derived as where G 0 , G 1 and G 2 stand for the leading moments for dead polymers. Then, the MWD of the polymers can be described by a Schultz-Zimme distribution function, shown as where y denotes the chain length and Therefore, the molecular weight distribution of the produced styrene polymers can be calculated from the leading moments of polymers by Equation (39). The number-average molecular weight M n and weight-average molecular weight M w can be easily obtained from the molecular weight distribution, and the polydispersity of the polymers can be calculated as H I = M w /M n .
Once the molecular weight distribution of the produced styrene polymers is obtained, it can be approximated by a B-spline model. In this simulation, ten third-order basis functions are applied to construct the B-spline model, and the basis functions are formulated as where i stands for the ith basis function of the B-spline model, j denotes the chain length of the dead polymers and w(·) is the width of the basis function, and 9 199, i = 10 The error of the B-spline approximation of the produced styrene polymer MWD is presented in Figure 3. It is obvious that the approximation error keeps within a proper range, indicating that the B-spline model is able to approximate the produced styrene polymer MWD effectively. The initial MWD, as shown in Figure 4, reaches the maximum value around chain length 200 and then gradually decreases; thus, the error in the model decreases with an increasing chain length.
Entropy 2022, 24, 499 9 of 13 presented in Figure 3. It is obvious that the approximation error keeps within a proper range, indicating that the B-spline model is able to approximate the produced styrene polymer MWD effectively. The initial MWD, as shown in Figure 4, reaches the maximum value around chain length 200 and then gradually decreases; thus, the error in the model decreases with an increasing chain length.

MWD Shape Control of the Styrene Polymerization Process
The control objective in this simulation is regulating the molecular weight distribution of the produced styrene polymers towards the target distribution, which is shown in Figure 5. The ratio, c , is the control input of the MWD shape control system for the styrene polymerization process.
The sampling period is set to 1s, the initial input is 0.47 and for the weights of the performance criterion (27), Q is a 10-order identity matrix, and The simulation results are shown in Figures 4-8.

MWD Shape Control of the Styrene Polymerization Process
The control objective in this simulation is regulating the molecular weight distribution of the produced styrene polymers towards the target distribution, which is shown in Figure 5. The ratio, c, is the control input of the MWD shape control system for the styrene polymerization process.

MWD Shape Control of the Styrene Polymerization Process
The control objective in this simulation is regulating the molecular weight distribution of the produced styrene polymers towards the target distribution, which is shown in Figure 5. The ratio, c , is the control input of the MWD shape control system for the styrene polymerization process.
The sampling period is set to 1s, the initial input is 0.47 and for the weights of the performance criterion (27), Q is a 10-order identity matrix, and The simulation results are shown in Figures 4-8.    The sampling period is set to 1s, the initial input is 0.47 and for the weights of the performance criterion (27), Q is a 10-order identity matrix, and R = 0.1.
The simulation results are shown in Figures 4-8.           Figure 4 demonstrates that the performance criterion can be regulated to a small value in a relatively short time by the proposed control method. The MWDs of the polymerization process at certain times are demonstrated in Figure 5. The variations of the ratio (manipulated variable of the control system) are illustrated in Figure 6, and the variations are within a proper range; the ratio keeps relatively stable when the MWD of the polymerization process is regulated close to the target one. The evolution of the MWD is illustrated by a three-dimensional graph in Figure 7. The change in the polydispersity over time is shown in Figure 8, which can be kept within a proper range. It can be seen that the MWD of the polymerization process gradually approaches the target MWD when the performance criterion gradually becomes smaller under the proposed algorithm.
The simulation results shown above illustrate that the moment-generating functionbased control algorithm for polymerization processes' MWD control is effective, and both the regulating time and the evolution of the MWD are satisfactory.

Conclusions
In this work, a control algorithm using the moment-generating function was presented for MWD shaping of polymerization processes. The output MWD was approximated by a discrete-time B-spline MWD model. Based on the input and expansion weight data of the B-spline model (which was used to approximate the output MWD), a system model could be identified through the N4SID method. Then, a performance criterion based on the moment-generating function was proposed, and the control algorithm could be derived by minimizing the performance criterion. The proposed control algorithm was tested in a styrene polymerization process, and the simulation results confirmed its effectiveness. The manipulated variable can be kept within a proper range, and the regulation time is satisfactory. Output MWD can also be adjusted to the target MWD gradually. The contributions can be summarized as:

1.
The moment-generating function was used to construct the performance criterion, which simplifies the regulation of criterion weights so that the application of the proposed control method to the practical polymerization process becomes more convenient; 2.
The subspace identification method was applied to establish the dynamic weights model of the controlled system based on the pair data of the control input and expansion weights, which avoids the numerical ill-conditioning and parameter overlap problems in traditional identification methods; 3.
The proposed moment-generating, function-based control algorithm is effective for the MWD shaping of the styrene polymerization process. Both the evolution of the output MWD and the regulation time are satisfactory.