Next Article in Journal
Application of Raman Spectroscopy to Evaluate the Structure Changes of Lubricating Grease Modified with Montmorillonite after Tribological Tests
Previous Article in Journal
Hollow Fiber Membrane Modification by Interfacial Polymerization for Organic Solvent Nanofiltration
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Improved On-Line Recursive Subspace Identification Method Based on Principal Component Analysis and Sliding Window for Polymerization

1
School of Biological and Chemical Engineering, Zhejiang University of Science & Technology, Hangzhou 310023, China
2
Hangzhou Sinan Intellitech Co., Ltd., Hangzhou 310016, China
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Processes 2024, 12(3), 562; https://doi.org/10.3390/pr12030562
Submission received: 28 January 2024 / Revised: 7 March 2024 / Accepted: 11 March 2024 / Published: 13 March 2024
(This article belongs to the Section Chemical Processes and Systems)

Abstract

:
Polymerization products are indispensable for our daily life, and the relevant modeling process plays a vital role in improving product quality. However, the model identification of the related process is a difficult point in industry due multivariate, nonlinear and time-varying characteristics. As for the conventional offline subspace identification methods, the identification accuracy may be not satisfying. To handle such a problem, an enhanced on-line recursive subspace identification method is presented on the basis of principal component analysis and sliding window (RSIMPCA-SW) in this paper to obtain the state space model for polymerization. In the proposed on-line subspace identification approach, the initial L-factor is acquired by the LQ decomposition of the sampled historical data, firstly, and then it is updated recursively through the bona fide method after the new data have been handled by the sliding window rule. Subsequently, principal component analysis (PCA) is introduced to calculate the extended observation matrix, and finally the on-line model parameters are extracted. Compared with the traditional subspace schemes, smaller computation complexity and higher identification precision are anticipated in the proposed method. A case study on the modeling of the ethylene polymerization verifies the effectiveness of the developed approach, in which the related statistical indexes of the obtained identification model are better.

1. Introduction

Polymerization refers to the process of converting low molecular weight monomers into high molecular weight polymers, and it is quite common in many industries [1]. There are two kinds of polymerization: addition polymerization and condensation polymerization. For addition polymerization, the unsaturated bonds in monomers (typically carbon–carbon double bonds) are opened and reacted with reactive groups in other monomer molecules, which forms polymer chains. Condensation polymerization is a process whereby small molecular units within monomers (typically dimers or oligomers) are linked together to form larger polymer molecules. These two types of polymerization reactions play a crucial role in the preparation of polymer materials, and the choice between them depends on the structure, properties and applications of the desired materials. Addition polymerization is usually employed for the synthesis of linear polymers, while condensation polymerization is more suitable for producing branched or crosslinked polymers.
During the past few decades, polymer manufacturing technology has made great progress, with the development of dynamic principles, process mathematical models and so on [2]. It is worth mentioning that relevant modeling is one of the significant points for polymerization, due to the fact that it is of great significance for improving the related industrial production efficiency and product quality [3]. As for polymerization, finding the corresponding modeling is not an easy task because of complex reaction mechanisms, uncertain kinetics processes, etc. [4]. Many scholars contributed to the development of the modeling for polymerization. Early mathematical models were built on the basis of first principles, with researchers focusing on reaction kinetics, thermodynamics, mass transfer and particle science. In [5], polymerization was modeled from the perspectives of polymerization kinetics (micro-scale), system thermodynamics (meso-scale) and reactor performance (macro-scale). The Sanchez–Lacombe equation of state model was utilized to simulate the solubility of industrial alkane penetrants in polyethylene, and the analysis of the results showed that the model could provide an accurate prediction of density in olefin polymerization reactors and further guide the production of linear low-density polyethylene in [6]. Zhang et al. [7] derived a polymer thermal reaction kinetic model in the form of an integral, and two undetermined functions were used to compensate for the changes in the reaction rate in the model, and two correction methods were given simultaneously, which made the kinetic model more applicable and effective. In [8], a multi-scale model of the radical polymerization process was effectively built. The employment of the finite element method solved the difficulty of calculating the cost and accuracy of the Fokker–Planck Equation, which can predict the change in particle size distribution on the mesoscale. The above models, established from first principles, can be utilized as substitutes for some polymerization processes after being validated by experimental data. However, the application of these models in practice may be difficult because of inconvenient and expensive maintenance costs, and these mechanism models are often set up under various assumptions, which lead to inconsistency with the actual process [9].
With the development of computer technology and intelligent algorithms, many vital system identification methods have been put forward for modeling polymerization. In order to solve the constraint of the traditional mechanism model for the polymer extrusion process, a system of fuzzy rules was introduced into the extrusion model to fit the sensitive parameters in [10]. Meanwhile, the system was combined with a genetic algorithm to optimize the model structure, and the relevant simulation verified the good adaptability of the built model. In [11], fuzzy logic was applied to the modeling of a discontinuous polymerization reactor to solve the problem of the mechanism model being unable to track the reactor system with transient behavior effectively, and the obtained nonlinear model could effectively predict the change in the average molecular weight of lactic acid and nylon-6 in the pilot scale polymerization experiments. A method based on an expectation maximization algorithm was utilized to intelligently identify batch polymerization reactor process models by Liu et al. [12]. The developed method was mainly concerned with the unknown inputs of the process system, and it was divided into two steps to calculate the model state and the analytic solution of the inputs.
As an important branch of system identification, the subspace method combines linear algebra, statistics and system theory, and it is aimed at establishing the appropriate state space model for processes [13,14]. Comparing with first-principles modeling and other system identification approaches, the subspace method, which adopts algebraic tools such as QR decomposition and singular value decomposition (SVD), possesses the merits of faster computation, better robustness and so on [15]. On the basis of the obtained state space model, the subsequent advanced control strategies can improve the corresponding key indexes in the process further [16]. And many researchers have focused on the study of off-line subspace identification methods and put forward significant results. In order to achieve the effective identification of linear time-invariant systems, a class of output error methods were developed in [17]. Van et al. [18] derived the subspace state system identification algorithm for mixed deterministic stochastic systems. A subspace method in which the bias can be eliminated was developed, and high identification precision was acquired for industries with colored noises in [19]. In order to realize multivariate predictive control in a continuous polymerization reactor of methyl methacrylate, Song et al. [20] developed a subspace algorithm to recognize the Wiener model of the process and then designed a data-based predictive controller. In [21], the state space model based on the subspace identification method was applied to historical process data and nonlinear ultrasonic data to predict the quality, and the results were verified in the experimental data of polyethylene rotary forming process.
However, the identification performance of the off-line subspace method will be reduced for the time-varying process, which promotes the development of the on-line subspace approaches. And there are also many vital related viewpoints. In order to reduce the computation complexity of the on-line identification method, many researchers focused on the selection of approximate methods to solve the updating problem of SVD. In [22], a propagator method which was commonly utilized in array signal processing was introduced to track the subspace of an observable matrix. Meanwhile, many other researchers designed subspace identification methods through updating QR decomposition or SVD directly. The first-order perturbation theory was employed to implement the recursive form of SVD by Shang et al. [23] and the application of the designed recursive canonical variable analysis method in a continuous stirred tank reactor system for polymerization validated the adaptability of this approach to time-varying processes. In order to realize the high-precision identification of time-varying systems, the Heisenberg QR decomposition and eigenvalue decomposition were united to extract the subspace matrix effectively by Wu et al. [24] Meanwhile, the forgetting factor was also introduced in the developed scheme.
It is a fact that there are various complex characteristics in polymerization and inevitable uncertainties in the related production, which presents a big challenge for the relevant accurate modeling process. As for these on-line subspace methods mentioned above, there is still room for the improvement of the recursive computation speed and the identification performance. In order to acquire the desired model for polymerization, it is essential to build an on-line model identification approach which provides faster recursive speed and improved identification precision. To achieve this study objective, an enhanced on-line recursive subspace identification method on the basis of PCA and sliding window is investigated in this paper. By employing PCA, the estimation accuracy of the extended observable subspace is improved, then the enhanced model identification precision is expected. Meanwhile, the utilization of LQ decomposition and the bona fide method simplify the related calculation steps, and thus a faster recursive speed is also anticipated. To validate the effectiveness of the proposed recursive subspace identification method, a case study on the identification experiment of the dynamic model for polyethylene process is demonstrated.

2. Traditional Offline Subspace Identification Method

2.1. Matrix Definition

For simplicity, the following state space model is introduced for polymerization.
x ( k + 1 ) = A x ( k ) + B u ( k ) + w ( k ) y ( k ) = C x ( k ) + D u ( k ) + v ( k )
where u ( k ) m and y ( k ) l are the input and output vectors of the process at time instant k , respectively. x ( k ) n is the state vector, and w ( k ) n and v ( k ) l denote the corresponding zero-mean noise vectors. A n × n , B n × m , C l × n and D l × m are the relevant parameter matrices of the state space model.
We assume that the length of past horizon is p , the length of the future horizon is f , and the length of the dataset used for model identification is N . We can define the past output stack vector y p ( k ) and the future output stack vector y f ( k ) as
y p ( k ) = y ( k p ) y ( k p + 1 ) y ( k 1 ) l p
y f ( k ) = y ( k ) y ( k + 1 ) y ( k + f 1 ) l f
On this basis, then the above stack vectors can be stacked into the following past and future output Hankel matrices, and their columns are defined as J , J = N p f + 1 .
Y p = y p ( k )   y p ( k + 1 ) y p ( k + J 1 ) l p × J
Y f = y f ( k )   y f ( k + 1 ) y f ( k + J 1 ) l f × J
Note that the definitions of the input Hankel matrices U p m p × J , U f m f × J and the noise Hankel matrices W p n p × J , W f n f × J , V p l p × J , V f l f × J are similar to the above formulas and omitted here.
The following input–output Hankel equation is then obtained by the iterative computation of the state space model in Equation (1).
Y f = Γ f X f + H f U f + G f W f + V f
where Γ f , X f are the extended observation matrix and the state sequence with the following definitions.
Γ f = C T   ( C A ) T ( C A f 1 ) T T l f × n
X f = x ( k )   x ( k + 1 ) x ( k + J 1 ) n × J
The remaining two Toeplitz matrices H f , G f in Equation (6), are then defined as
H f = D 0 0 C B D 0 C A f 2 B C A f 3 B D l f × m f
G f = 0 0 0 C 0 0 C A f 2 C A f 3 0 l f × n f

2.2. Double-Orthogonal-Projection-Based Subspace Identification Method (2ORT-SIM) [19]

In order to simplify the derivation process, here the noise term in Equation (6) is defined as E f = G f W f + V f , and then Equation (6) can be converted into the following form.
Y f = Γ f X f + H f U f + E f
To obtain the parameter matrices A and C , the extended observation matrix Γ f needs to be computed, firstly, and then the two required parameter matrices can be extracted from it. Taking the orthogonal projection of Equation (11) onto the orthogonal complement of U f , the relevant effects on the dataset can be eliminated.
Y f / U f = ( Γ f X f + H f U f + E f ) / U f
where / denotes the projection operator, and ( ) is the orthogonal complementary space of the corresponding matrix. By adopting U f / U f = 0 , Equation (12) can be rewritten as
Y f / U f = ( Γ f X f + E f ) / U f
In order to eliminate the effect of E f in Equation (13), another orthogonal projection needs to be performed. And the following equation is obtained by projecting Equation (13) onto the past input U p .
Y f / U f / U p = ( Γ f X f + E f ) / U f / U p
If the model input is a persistent excitation, then the following condition is satisfied.
lim N E f / U f / U p = 0
Here, the proof of Equation (15) can be found in reference [19] and the relevant contents are omitted here for brevity. Then, Equation (14) can be transformed into
Y f / U f / U p = Γ f X f / U f / U p
By taking an SVD for the left-hand side of Equation (16), the following formula is acquired.
Y f / U f / U p = Ψ ^ 1 Ψ ^ 1 Σ ^ 1 0 0 Σ ^ 2 ϒ ^ 1 T ( ϒ ^ 1 ) T
The extended observation matrix Γ f is estimated as
Γ f = Ψ ^ 1
and then the parameter matrices A and C are extracted as
A = ( K 1 Γ f ) K 2 Γ f C = K 3 Γ f
where K 1 = I ( f 1 ) l   0 ( f 1 ) l × l , K 2 = 0 ( f 1 ) × l   I ( f 1 ) l , K 3 = I l   0 l × ( f 1 ) l . ( ) denotes the pseudo-inverse of the matrix.
As for the identification of parameter matrices B and D , they can be obtained by introducing Γ f and U f to both sides of Equation (11). And the following equation is obtained by utilizing Γ f Γ f = 0 and U f U f = I .
Γ f Y f U f = Γ f H f + Γ f E f U f
Similarly, lim N Γ f E f U f = 0 will be tenable if model inputs are persistent excitations. That is, the noise in Equation (20) will be eliminated accordingly. For the convenience of subsequent calculations, the remaining items in Equation (20) are defined as follows.
M = Γ f Y f U f = M 1 M f F = Γ f = F 1 F f
where the corresponding matrix dimensions are M ( l f n ) × m f and F ( l f n ) × l f .
Then, the parameter matrices B and D can be obtained by solving the linear combination in Equation (22) through least squares.
D B = M 1 M f [ F 1 F f 1 F f F 2 F f 0 F f 0 0 I 0 0 K 1 Γ f ]

3. Proposed On-Line Subspace Identification Method

Here, the state space model of the polymerization in Equation (1) is also considered, and the relevant equations before Equation (14) can be obtained easily by referring to Section 2.
In order to improve the corresponding identification precision for the polymerization in the presence of various uncertainties in practice and to obtain a smaller computation complexity for the on-line subspace algorithms, here, the recursive technique, which is the integration of LQ decomposition with the bona fide method, is employed to update the model parameters in real time.
The initial L-factor at time instant k needs to be obtained, firstly, and the second orthogonal projection performed for Equation (14) is replaced by the following LQ decomposition, in which Z f = Y f / U f .
U p ( k | k J + 1 ) Z f ( k | k J + 1 ) = L 11 ( k ) 0 L 21 ( k ) L 22 ( k ) Q 1 T ( k ) Q 2 T ( k )
where the variables k J + 1 and k in parentheses denote the starting modeling moment and current modeling moment, respectively.
Then, the left-hand side of Equation (14) is equivalent to
Z f / U p = L 21 ( k ) Q 1 T ( k )
Define Z f 1 = L 21 ( k ) Q 1 T ( k ) , and the following formula for the initial L-factor can be constructed by employing the property of Q 1 T ( k ) Q 1 ( k ) = I .
Z f 1 ( Z f 1 ) T = L 21 ( k ) L 21 T ( k )
Then, the bona fide method is utilized to obtain a recursive update expression for L 21 ( k ) L 21 T ( k ) .
When the input–output data pairs { u ( k + 1 ) , y ( k + 1 ) } at time instant k + 1 are collected, the following equation can be obtained by adding the new dataset to the left-hand side of Equation (23).
U p ( k | k J + 1 )     U p ( k + 1 ) Z f ( k | k J + 1 )     Z f ( k + 1 ) = U p ( k J + 1 )     U p ( k + 1 | k J + 2 ) Z f ( k J + 1 )     Z f ( k + 1 | k J + 2 )
The formulas as follows are satisfied according to the LQ decomposition shown in Equation (23).
U p ( k | k J + 1 ) = L 11 ( k ) Q 1 T ( k ) Z f ( k | k J + 1 ) = L 21 ( k ) Q 1 T ( k ) + L 22 ( k ) Q 2 T ( k )
On the basis of Equation (27), the following equations can be acquired.
U p ( k + 1 | k J + 2 ) = L 11 ( k + 1 ) Q 1 T ( k + 1 ) Z f ( k + 1 | k J + 2 ) = L 21 ( k + 1 ) Q 1 T ( k + 1 ) + L 22 ( k + 1 ) Q 2 T ( k + 1 )
Combining Equations (27) and (28) with Equation (26), and redefining the relevant block matrices as
Θ 1 = L 11 ( k ) Q 1 T ( k ) L 21 ( k ) Q 1 T ( k ) + L 22 ( k ) Q 2 T ( k ) Θ 2 = U p ( k + 1 ) Z f ( k + 1 ) Θ 3 = U p ( k J + 1 ) Z f ( k J + 1 ) Θ 4 = L 11 ( k + 1 ) Q 1 T ( k + 1 ) L 21 ( k + 1 ) Q 1 T ( k + 1 ) + L 22 ( k + 1 ) Q 2 T ( k + 1 )
the relationship of Θ 1 Θ 2 = Θ 3 Θ 4 reveals that Θ 1 Θ 1 T + Θ 2 Θ 2 T = Θ 3 Θ 3 T + Θ 4 Θ 4 T holds, and the recursive updating equations for the L-factor from time instant k to time instant k + 1 can be gained according to Q i T Q j = 0 , i j .
L 11 ( k + 1 ) L 11 T ( k + 1 ) = L 11 ( k ) L 11 T ( k ) + U p ( k + 1 ) U p T ( k + 1 ) U p ( k J + 1 ) U p T ( k J + 1 )
L 11 ( k + 1 ) L 21 T ( k + 1 ) = L 11 ( k ) L 21 T ( k ) + U p ( k + 1 ) Z f T ( k + 1 ) U p ( k J + 1 ) Z f T ( k J + 1 )
L 21 ( k + 1 ) L 11 T ( k + 1 ) = L 21 ( k ) L 11 T ( k ) + Z f ( k + 1 ) U p T ( k + 1 ) Z f ( k J + 1 ) U p T ( k J + 1 )
In order to construct the form of L 21 ( k + 1 ) L 21 T ( k + 1 ) , the following equation is obtained by uniting Equations (30)–(32).
L 21 ( k + 1 ) L 21 T ( k + 1 ) = L 21 ( k + 1 ) L 11 T ( k + 1 ) ( L 11 ( k + 1 ) L 11 T ( k + 1 ) ) L 11 ( k + 1 ) L 21 T ( k + 1 )
where
L 11 T ( k + 1 ) ( L 11 ( k + 1 ) L 11 T ( k + 1 ) ) L 11 ( k + 1 ) = I
In order to improve the estimation precision of the extended observation matrix, PCA, which has a good performance in data feature extraction, is introduced after the L-factor at time instant k + 1 has been acquired. And the left singular vector Z f 1 can be obtained by carrying out the following PCA on the left-hand side of Equation (33).
L 21 ( k + 1 ) L 21 T ( k + 1 ) = P s T s T + P r T r T
where P , T are the load matrix and score matrix, respectively. Subscripts s and r denote the main subspace and residual subspace.
When the residual load matrix P r is acquired by estimation, then the extended observation matrix at time instant k + 1 can be obtained by the following formula.
Γ f ( k + 1 ) = ( ( P r ( k + 1 ) ) T )
Finally, the model parameters A ( k + 1 ) , B ( k + 1 ) , C ( k + 1 ) , D ( k + 1 ) at time instant k + 1 can be gained by referring to Equations (19)–(22).
A ( k + 1 ) = ( K 1 Γ f ( k + 1 ) ) K 2 Γ f ( k + 1 ) C ( k + 1 ) = K 3 Γ f ( k + 1 )
D ( k + 1 ) B ( k + 1 ) = M 1 ( k + 1 ) M f ( k + 1 ) [ F 1 ( k + 1 ) F f 1 ( k + 1 ) F f ( k + 1 ) F 2 ( k + 1 ) F f ( k + 1 ) 0 F f ( k + 1 ) 0 0 I 0 0 K 1 Γ f ( k + 1 ) ]
Remark 1.
Note that the size of the data matrix is fixed in the bona fide method, and the related operation is performed by sliding the data window, and thus the computation complexity of the proposed subspace method is reduced greatly.
Remark 2.
The choice of model order plays a crucial role for the final modeling results. There are two main methods for model order estimation in most subspace methods. One is determined by checking the intervals in the singular value profile of the Hankel matrix, which is called the singular value method [25], and another is determined by the spectrum of the AIC criterion [26]. However, the optimal order of the model may not be obtained because of various unavoidable noises. An effective method for estimating the model order is given for this polyethylene reaction process in this paper. Here, an optional range of the system order is assumed, firstly, and then the errors between the model and the actual value for different orders are calculated. Finally, the order with smaller errors and acceptable computation complexity will be selected as the optimal order.
In order to enhance the readability of the paper, the related detailed steps are summarized as follows for the proposed subspace strategy (Algorithm 1).
Algorithm 1. Steps for the Proposed RSIMPCA-SW
(1) Initialize inputs of the algorithm at time instant k , including past horizon p , future horizon f , number of samples N , input dataset u , output dataset y .
(2) Iterative computation of Equation (1) to construct the input–output Hankel formula in Equation (11).
(3) For the first orthogonal projection, construct Equation (13) by projecting Y f onto the orthogonal complement of U f .
(4) Perform the LQ decomposition of Equation (23) instead of the second orthogonal projection, and construct L 21 ( k ) L 21 T ( k ) for the initial L-factor at time instant k .
(5) Add the new data at time instant k + 1 to the last column of the original data matrix to construct Equation (26).
(6) Recursive update of L-factor by bona fide method according to Equations (30)–(34).
(7) Perform PCA on the L-factor at time instant k + 1 according to Equation (35).
(8) Calculate the extended observation matrix Γ f by referring to Equation (36).
(9) Extract the model parameters A ( k + 1 ) , B ( k + 1 ) , C ( k + 1 ) , D ( k + 1 )  according to Equations (37) and (38).
(10) Slide the data window to remove the first column of the data matrix in Equation (26) and return to step 5.

4. Simulation

4.1. Polyethylene Process Flow Description

The process described in this section runs on a full-density polyethylene plant, which adopts typical Unipol low pressure gas phase fluidized bed polymerization technology. Here, ethylene is used as the main monomer and butene-1 or hexene-1 is utilized as the co-monomer, with a content range of 1.2–3 wt%. In addition, hydrogen and isopentane are employed as the molecular weight regulator and the condensation inducer, respectively. Note that the operation is carried out in the condensation mode, and it should be guaranteed that the dew point of the circulating gas is above 45 degrees Celsius. During polymerization, the condensate content in the circulating gas is kept at about 10 wt% to improve the reaction efficiency. The products are linear low-density polyethylene and partial medium- and high-density polyethylene particle resin, and the relevant density and melt index ranges are 0.915–0.965 g/cm3 and 0.3–6 g/10 min. This process is mainly divided into the raw material refining system and polymerization system. In terms of the comprehensive product quality and the production cost, this technological process is still the current mainstream technology for polyethylene production. Figure 1 shows the detailed flow chart of the polyethylene process.
The polymerization system mainly consists of the reactor (R-4000), cooler (E-4000) and compressor (C-4000). The raw materials and catalysts are sent to the fluidized bed reactor continuously from the refining system of polymerization to produce polyethylene resin, and the storage tank (D-4090) stores the intermediate products temporarily. The gas phase materials are circulated under the action of a centrifugal circulating gas compressor and a water-cooled circulating gas cooler incessantly. Circulating gases improve the fluidized bed’s back-mixing ability and carry the reaction material to the reaction zone. At the same time, the circulating gases also take away the heat generated by the exothermic reaction. As for the reactor, it has two discharge systems which operate alternately in general, but can also be operated separately. Each discharge system consists of a product disposal tank and a product ejection tank. The discharge system starts automatically when the reactor bed reaches a predetermined level, and the powder products are dropped continuously from the disposal tank to the ejection tank, and then the power products are transported to the product degassing bin of the degassing system directly. Here, a termination system is set up for the polymerization system, and the polymerization will be stopped partially or completely when the reversible catalyst poison is injected into the reactor. Note that removing the poison from the bed will reactivate the catalyst, and the reaction can be restarted quickly.

4.2. Data Acquisition and Preprocessing

It is not an easy task to choose the appropriate variables for the effective modeling of complex processes. Meanwhile, many variables may not be measured because of the limitations of the relevant measuring sensors. Moreover, variables in the processes are commonly coupled due to the fact that there are multiple sub-reaction processes. All these factors affect the corresponding precision of the process modeling greatly.
In this section, four typical variables are selected for model identification. Among these, ethylene flow rate (C2f (kg/h)) and catalyst feed amount (Catf (kg/h)) are chosen as inputs, and ethylene partial pressure (ppc2 (ppm)) and polyethylene production (prod (kg/h)) are selected as outputs. Their relationships are shown in Figure 2.
Here, 2000 sets of real sample data collected from a full-density polyethylene plant of a chemical company in 33.3 h (sampling frequency is 60 s) were selected for the simulation. In order to illustrate the effectiveness and generalization ability of the proposed on-line algorithm preferably, a little more data are required for validation, and the ratio of 3:2 was chosen for splitting the 2000 sample data. On the basis of this ratio, the first 1200 samples of the dataset were classified as the identification set, and the remaining 800 samples were chosen as the validation set in the model identification. Since the difference of scale between data is too large, the normalization preprocessing method is utilized, which takes the following form.
x ˜ i , j = ( b a ) x i , j min ( j ) max ( j ) min ( j ) + a
where x i , j and x ˜ i , j denote the data before and after the normalization. min ( j ) , max ( j ) are the minimum and maximum values of the corresponding attributes in data. a is the relevant lower limit of the magnitude after normalization, and b is the related upper limit. The preprocessed data of the process for modeling are shown in Figure 3 and Figure 4.

4.3. Model Order Estimation

As is mentioned in Remark 2, the first 200 sets of the total sampled dataset are utilized for the model order estimation here, and the relevant error spectrum is shown in Figure 5.
Here, the theoretical range of the optimal model order is chosen between [1,9] in this section. The circle scatters indicate the distribution of relative squared error (RSE) for ppc2 at each order, and the asterisk scatters indicate the distribution of RSE for prod. The right triangular scatters indicate the average of the two outputs, i.e., the mean relative squared error (MRSE). It can be easily seen that the model possesses the minimum value of MRSE at n = 3 . Here, the MRSE values for some higher orders may be also acceptable; however, modeling the polymerization with a lower order can reduce the relevant complexity for the subsequent controller design, and thus the optimal model order is selected as 3 here.

4.4. Identification Results and Validation Analysis

The initial input parameters of the proposed algorithm are set as follows: the length of the past horizon is p = 10 , the length of the future horizon is f = 10 , and the number of total samples is N = 2000 . In order to verify the effectiveness of the proposed algorithm further, several existing off-line subspace identification algorithms (MOESP [17], N4SID [18], 2ORT-SIM [19]) and the on-line subspace identification algorithm (OSIMPCA-E [24]) are introduced as the comparison. Here, the 1200 sets of sample data in the identification set are used for the relevant training.
Then, the 800 sets of sample data in the validation set are utilized to test the prediction accuracy of these built models in this section, and the RSE, mean absolute error (MAE) and variance accounted for (VAF) are calculated to evaluate the corresponding model accuracy.
The calculation methods of RSE and MAE are shown in Equations (40) and (41), and their values range from 0 to 1. When their values are close to 0, a high matching degree between the prediction model and the actual process will be obtained. The VAF is calculated as shown in Equation (42); note that an ideal identification precision is expected when its value is near 1.
R S E = j = 1 N v d ( y ( j ) y ^ ( j ) ) 2 j = 1 N v d ( y ( j ) ) 2
M A E = 1 N v d j = 1 N v d ( y ( j ) y ^ ( j ) )
V A F = 1 v a r i a n c e ( y y ^ ) v a r i a n c e ( y )
where N v d denotes the number of samples in the validation set. j is the sampling serial number, and y ( j ) and y ^ ( j ) represent the actual value and the model predicted value of the j th data, respectively.
Note that fixed identified models are used for these off-line methods, and recursive update models are adopted for on-line approaches. And the related statistical results of these algorithms are given in Table 1.
In order to test the performance of the RSIMPCA-SW model, the relevant predicted values and the real values are displayed in Figure 6 and Figure 7, further, where the OSIMPCA-E approach is introduced as the comparison.

5. Discussions

There are some general guidelines for the selection of the parameters for the proposed on-line subspace identification approach. In terms of the past and future horizons p and f , they will affect the dimensions of the relevant Hankel matrices, and their values must be bigger than the maximum order of the pre-set process model. For the selection of the optimal model order n , it influences the relevant model identification precision greatly. As is mentioned in Section 4.3, the value of the corresponding MRSE can be chosen as the criterion. Note that although the MRSE values of some other orders may also be acceptable, the selection of the optimal order as 3 will offer a better compromise between model precision and the subsequent computation complexity.
By observing the statistical results for both approaches in Table 1, it can be easily seen that better identification performance is achieved for two on-line methods in which smaller RSE, MAE and bigger VAF are obtained. It is also implied that the recursion technology in on-line identification algorithms is beneficial for coping with various uncertainties in practice and complex characteristics in polymerization. Among the two on-line subspace identification schemes, the proposed RSIMPCA-SW algorithm shows superiority with improved statistical results in which RSE and MAE are smaller and VAF is larger. The corresponding responses of two on-line approaches are shown in Figure 6 and Figure 7. It can be found that the predicted values of the proposed strategy are closer to the actual values, comparing with those of the OSIMPCA-E approach, by and large, which also indicates that enhanced model identification precision is achieved for the proposed RSIMPCA-SW method. In a word, the proposed on-line scheme provides better ensemble identification performance.
Note that it is common to fill in the missing data using interpolation methods (e.g., linear interpolation and polynomial interpolation) during the data preprocessing stage. When a variable which needs to be modeled is not available, the other highly correlated variable will be chosen as the alternate. Then, the result of the alternate variable can be regarded as the related reference. For a polymerization processes with more variables, faster recursive computation speed and improved identification accuracy are also expected for the proposed algorithm because of the adoption of PCA and LQ decomposition with the bona fide method.

6. Conclusions

In this paper, a new on-line recursive subspace identification method based on PCA and sliding window is developed to obtain a more precise model with smaller computation complexity for polymerization. Compared with the conventional subspace methods, PCA is employed to improve the estimation precision of the extended observation matrix, and LQ decomposition with the bona fide method are also utilized for reducing the computation complexity of the recursion in the on-line identification. Based on the introduction of the above two technologies, improved identification precision and lower computation complexity are anticipated. At the same time, an optimal estimation method for the process model order under a recommended criterion is given. The case study of the model identification of the polyethylene process shows that smaller ensemble prediction errors were obtained for the proposed RSIMPCA-SW approach. Meanwhile, the VAF results and the related responses also imply that better fitting performance is gained for the proposed on-line algorithm. And all these above results verify the superiority of the proposed RSIMPCA-SW approach compared with the introduced off-line and on-line subspace identification algorithms.

Author Contributions

Conceptualization, J.Q., T.L. and G.H.; methodology, J.Q.; software, J.Q. and J.Z.; validation, J.Q., S.L. and C.S.; formal analysis, C.S. and S.L.; investigation, J.Q.; resources, J.Z.; data curation, J.Q. and S.L.; writing—original draft preparation, J.Q., T.L. and C.S.; writing—review and editing, J.Q., T.L. and G.H.; visualization, J.Q.; supervision, G.H. and B.W.; project administration, C.S. and G.H.; funding acquisition, J.Z., G.H. and B.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Zhejiang Provincial Science and Technology Plan Project, China (Grant No. 2023C01117).

Data Availability Statement

Data are available upon request due to privacy restrictions. The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest

Authors Jubin Zhang and Bin Wen were employed by the company Hangzhou Sinan Intellitech Co., Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

  1. Ray, W.H.; Villa, C.M. Nonlinear dynamics found in polymerization processes—A review. Chem. Eng. Sci. 2000, 55, 275–290. [Google Scholar] [CrossRef]
  2. Kreitser, T.; Tumina, S.; L’Vovskii, V. Effect of the parameters of the mathematical model of polymerization kinetics, on estimation of the molecular weight of polyethylene. Polym. Sci. U.S.S.R. 1978, 20, 2989–2998. [Google Scholar] [CrossRef]
  3. Vivaldo-Lima, E.; Mohammadi, Y.; Penlidis, A. Modeling and simulation of polymerization processes. Processes 2021, 9, 821. [Google Scholar] [CrossRef]
  4. Fiosina, J.; Sievers, P.; Drache, M.; Beuermann, S. Polymer reaction engineering meets explainable machine learning. Comput. Chem. Eng. 2023, 177, 108356. [Google Scholar] [CrossRef]
  5. Touloupidis, V. Catalytic Olefin Polymerization Process Modeling: Multi-Scale Approach and Modeling Guidelines for Micro-Scale/Kinetic Modeling. Macromol. React. Eng. 2014, 8, 508–527. [Google Scholar] [CrossRef]
  6. Bashir, M.A.; Kanellopoulos, V.; Ali, M.A.; McKenna, T.F. Applied Thermodynamics for Process Modeling in Catalytic Gas Phase Olefin Polymerization Reactors. Macromol. React. Eng. 2020, 14, 1900029. [Google Scholar] [CrossRef]
  7. Zhang, Z.; Guo, F.; Song, W.; Jia, X.; Wang, Y. Empirical correction of kinetic model for polymer thermal reaction process based on first order reaction kinetics. Chin. J. Chem. Eng. 2021, 38, 132–144. [Google Scholar] [CrossRef]
  8. Urrea-Quintero, J.-H.; Marino, M.; Hernandez, H.; Ochoa, S. Multiscale modeling of a free-radical emulsion polymerization process: Numerical approximation by the Finite Element Method. Comput. Chem. Eng. 2020, 140, 106974. [Google Scholar] [CrossRef]
  9. Mueller, P.A.; Richards, J.R.; Congalidis, J.P. Polymerization reactor Modeling in Industry. Macromol. React. Eng. 2011, 5, 261–277. [Google Scholar] [CrossRef]
  10. Tan, L.P.; Lotfi, A.; Lai, E.; Hull, J. Soft computing applications in dynamic model identification of polymer extrusion process. Appl. Soft Comput. 2004, 4, 345–355. [Google Scholar] [CrossRef]
  11. Lima, N.M.; Linan, L.Z.; Melo, N.C.; Manenti, F.; Filho, R.M.; Embirucu, M.; Maciel, R.W. Nonlinear fuzzy identification of batch polymerization processes. Comput. Aided Chem. Eng. 2015, 37, 599–604. [Google Scholar]
  12. Liu, Z.; Luan, X. Intelligent Modeling for Batch Polymerization Reactors with Unknown Inputs. Sensors 2023, 23, 6021. [Google Scholar] [CrossRef]
  13. Qin, S.J. An overview of subspace identification. Comput. Chem. Eng. 2006, 30, 1502–1513. [Google Scholar] [CrossRef]
  14. Corbett, B.; Mhaskar, P. Subspace identification for data-driven modeling and quality control of batch processes. AIChE J. 2016, 62, 1581–1601. [Google Scholar] [CrossRef]
  15. Ghosh, D.; Hermonat, E.; Mhaskar, P.; Snowling, S.; Goel, R. Hybrid Modeling Approach Integrating First-Principles Models with Subspace Identification. Ind. Eng. Chem. Res. 2019, 58, 13533–13543. [Google Scholar] [CrossRef]
  16. Lou, H.; Su, H.; Gu, Y.; Hou, W.; Xie, L.; Rong, G. Nonlinear predictive control with modified closed-loop subspace identifi-cation-piecewise linear model for double-loop propylene polymerization process. Control. Theory Appl. 2015, 32, 1040–1051. [Google Scholar]
  17. Verhaegen, M. Identification of the deterministic part of MIMO state space models given in innovations form from in-put-output data. Automatica 1994, 30, 61–74. [Google Scholar]
  18. Van Overschee, P.; De Moor, B. N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica 1994, 30, 75–93. [Google Scholar]
  19. Hou, J.; Liu, T.; Chen, F. Orthogonal projection based subspace identification against colored noise. Control. Theory Technol. 2017, 15, 69–77. [Google Scholar] [CrossRef]
  20. Song, I.H.; Rhee, H.K. Nonlinear control of polymerization reactor by Wiener predictive controller based on one-step sub-space identification. Ind. Eng. Chem. Res. 2004, 43, 7261–7274. [Google Scholar]
  21. Gomes, F.P.C.; Garg, A.; Mhaskar, P.; Thompson, M.R. Data-Driven Advances in Manufacturing for Batch Polymer Processing Using Multivariate Nondestructive Monitoring. Ind. Eng. Chem. Res. 2019, 58, 9940–9951. [Google Scholar] [CrossRef]
  22. Mercère, G.; Bako, L.; Lecœuche, S. Propagator-based methods for recursive subspace model identification. Signal Process. 2008, 88, 468–491. [Google Scholar] [CrossRef]
  23. Shang, L.; Liu, J.; Zhang, Y.; Wang, G. Efficient recursive canonical variate analysis approach for monitoring time-varying processes. J. Chemom. 2017, 31, e2858. [Google Scholar] [CrossRef]
  24. Wu, P.; Chen, L.; Zhou, W.; Guo, L. Online subspace identification based on principal component analysis and noise estima-tion. J. Zhejiang Univ. 2018, 52, 1694–1701. [Google Scholar]
  25. Misra, S.; Nikolaou, M. Adaptive design of experiments for model order estimation in subspace identification. Comput. Chem. Eng. 2017, 100, 119–138. [Google Scholar] [CrossRef]
  26. Sugiyama, M.; Ogawa, H. Subspace Information Criterion for Model Selection. Neural Comput. 2001, 13, 1863–1889. [Google Scholar] [CrossRef]
Figure 1. Gas phase fluidized bed polyethylene process flow.
Figure 1. Gas phase fluidized bed polyethylene process flow.
Processes 12 00562 g001
Figure 2. Variables for subspace identification.
Figure 2. Variables for subspace identification.
Processes 12 00562 g002
Figure 3. Preprocessed process input data.
Figure 3. Preprocessed process input data.
Processes 12 00562 g003
Figure 4. Preprocessed process output data.
Figure 4. Preprocessed process output data.
Processes 12 00562 g004
Figure 5. Error spectrum for different orders.
Figure 5. Error spectrum for different orders.
Processes 12 00562 g005
Figure 6. The response comparison of ppc2.
Figure 6. The response comparison of ppc2.
Processes 12 00562 g006
Figure 7. The response comparison of prod.
Figure 7. The response comparison of prod.
Processes 12 00562 g007
Table 1. Statistical results of different algorithms.
Table 1. Statistical results of different algorithms.
AlgorithmRSEMAEVAF
ppc2prodppc2prodppc2prod
MOESP0.27260.22960.12420.10330.80590.9055
N4SID0.30910.23250.14360.10300.76990.8981
2ORT-SIM0.51720.29540.25670.14150.52050.8853
OSIMPCA-E0.13270.11950.04890.04800.88540.9450
RSIMPCA-SW0.12600.04830.04710.02170.91650.9885
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Qian, J.; Zhang, J.; Lei, T.; Li, S.; Sun, C.; He, G.; Wen, B. An Improved On-Line Recursive Subspace Identification Method Based on Principal Component Analysis and Sliding Window for Polymerization. Processes 2024, 12, 562. https://doi.org/10.3390/pr12030562

AMA Style

Qian J, Zhang J, Lei T, Li S, Sun C, He G, Wen B. An Improved On-Line Recursive Subspace Identification Method Based on Principal Component Analysis and Sliding Window for Polymerization. Processes. 2024; 12(3):562. https://doi.org/10.3390/pr12030562

Chicago/Turabian Style

Qian, Jiayu, Jubin Zhang, Ting Lei, Silin Li, Chen Sun, Guanghua He, and Bin Wen. 2024. "An Improved On-Line Recursive Subspace Identification Method Based on Principal Component Analysis and Sliding Window for Polymerization" Processes 12, no. 3: 562. https://doi.org/10.3390/pr12030562

APA Style

Qian, J., Zhang, J., Lei, T., Li, S., Sun, C., He, G., & Wen, B. (2024). An Improved On-Line Recursive Subspace Identification Method Based on Principal Component Analysis and Sliding Window for Polymerization. Processes, 12(3), 562. https://doi.org/10.3390/pr12030562

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop