Abstract
Simulated moving bed (SMB) chromatographic separation is a new type of separation technology based on traditional fixed bed adsorption operation and true moving bed (TMB) chromatographic separation technology, which includes inlet-outlet liquid, liquid circulation, and feed liquid separation. The input-output data matrices were constructed based on SMB chromatographic separation process data. The SMB chromatographic separation process was modeled by utilizing two subspace system identification algorithms: multivariable output-error state-space (MOESP) identification algorithm and numerical algorithms for subspace state-space system identification (N4SID), so as to obtain the 3rd-order and 4th-order state-space yield models of the SMB chromatographic separation process, respectively. The model predictive control method based on the established state-space models is used in the SMB chromatographic separation process. The influence of different control indicators on the predictive control system response performance is discussed. The output response curves of the yield models were obtained by changing the related parameters so that the yield model parameters are optimized set meanwhile. Finally, the simulation results showed that the yield models are successfully controlled based on the each control period and given yield range.
1. Introduction
Due to the complexity of simulated moving bed (SMB) chromatographic separation mechanism and many parameters affecting the separation effect [1, 2], especially the sensitivity to operating parameters and disturbances, the SMB chromatographic separation process is difficult to run at the optimal operating point for a long time. The separation zones of SMB chromatographic separation process were analyzed [3]. A method to optimize the separation was proposed so as to improve yield purity and economic efficiency [4]. It is an effective method to set up the effective and accurate model on the separation process and then implement advanced process control methods in order to solve the existing problems of SMB chromatographic separation process. There are many modeling methods for the SMB model, such as extracting the device features and simplifing mathematical calculation. Many researchers have studied the modeling methods of SMB. The general rate model (GRM) is the most accurate, but the actual model is difficult to solve. The analytical solutions of GRM were obtained by using the Laplace transform by using different boundaries of linear adsorption isotherms [2]. The model optimization sequence and the impact on the separation efficiency under different influencing factors were analyzed [5]. An alternative model method was proposed to obtain a reduced-order model of SMB so as to reduce the computational complexity of the model and achieve the desired optimization results [6]. The front velocity method was proposed for SMB process modeling in order to remove the method that adsorption isotherm was obtained by means of equilibrium experiments in the traditional SMB modeling process. Then, the feasibility of the model was verified by experiments [7]. Aiming at the existing problems of the discretization model by utilizing the finite-dimensional partial differential equations, two Krylov-type reduced-order models were proposed to optimize the computational accuracy and computational efficiency of the model [8]. The main control objective of the SMB chromatographic separation process is to increase the target yield of the SMB separation process. Many papers have made relevant research on process control methods. The genetic algorithm was adopted to optimize the parameters in the SMB separation process [9]. The standing wave design method was proposed to improve the purity and yield of the separated components while reducing the consumption of the analytical solution [10]. The flow rate control strategy was adopted, where two operational variables were added in each switching cycle and the SMB separation performance was improved by the SimCon operation [11]. Because the traditional numerical optimization method has many uncertainties, the system parameters of size exclusion simulated moving bed (SEC-SMB) were studied by batch chromatography [12]. The impact of the fault column on the SMB separation process was analyzed through experiments and theoretical analysis [13]. An adaptive nonlinear model predictive control method was proposed for an SMB device to realize the praziquantel enantiomer separation. Under the set purity levels, the controller can respond quickly [14]. Some scholars have carried out research and application on the control methods of MPC. Markus Bambach and Michael Herty had studied to control the strain rate during the process using model predictive control [15]. Saeed Shamaghdari and Mohammad Haeri had discussed the design of a new robust model predictive control algorithm for nonlinear systems [16]. Wahid A et al. study investigates a multivariable model predictive control (MPC) based on two-point temperature control strategy [17]. Some articles on the study of process control methods are also discussed. Yunjian Hu et al. proposed an optimal tension and thickness control method based on the receding horizon control (RHC) strategy to improve control performance and enhance system robustness [18]. Ren et al. studied a reduced-order extended state observer- (ESO-) based sliding mode control scheme for friction compensation of a three-wheeled omnidirectional mobile robot [19].
In this paper, a predictive control method based on the state-space yield models of SMB chromatographic separation process obtained by utilizing the subspace system identification method is proposed. The subspace system identification is an identification algorithm based on the input-output data matrix by means of mathematical tools, such as geometry, matrix analysis, and probability statistics so as to obtain the state-space matrix in state-space equations. Many papers have proposed many model identification methods for complex process models. van Overschee and de Moor elaborated on the multivariable output-error state-space (MOESP) and the numerical algorithms for subspace state-space system identification (N4SID) [20]. Hachicha et al. compared the morbid conditions of MOESP and N4SID and confirmed that they are robust in the parameter estimation process of the controlled object [21]. The discrete state-space equation was vectorized and scaled to solve the bilinear low-rank minimization problem [22]. The subspace identification method was applied to the fault diagnosis of the wind turbine shaft [22]. Model predictive control is an advanced optimization control method composed of the predictive model, rolling optimization, and feedback correction. Clarke developed a controlled autoregressive integral moving average model (CARIMA) [23]. Garic and Morari proposed the internal model control (IMC) algorithm [24]. Ding made a related study on the robustness of multiple time-delay systems [25]. Kayacan et al. successfully solved the trajectory tracking problem of the self-made tractor trailer system by using the fast distributed nonlinear model predictive control algorithm [26]. Li et al. designed a distributed predictive controller to satisfy the stability conditions by using a distributed implementation [27]. Farina and Scattolini proposed a new distributed predictive control algorithm for linear discrete-time systems [28]. Some newly published papers also propose new methods for model identification. Hou et al. designed a subspace identification method which is proposed for the Hammerstein-type nonlinear systems subject to periodic disturbances with unknown waveforms [29]. Gradient-based and least-squares-based iterative identification algorithms are developed for output-error (OE) and output-error moving average (OEMA) systems, and the simulation results confirm theoretical findings [30]. Ding et al. proposed an algorithm about a least-squares-based and a gradient-based iterative algorithm for the Hammerstein nonlinear systems using the hierarchical identification principle [31]. Ding used an iterative least-squares algorithm to estimate the parameters of output-error systems and enhance computational efficiencies [32]. Xu and Ding discussed the parameter estimation problems for dynamical systems by means of the impulse response measurement data [33]. Xu and Ding used the hierarchical identification principle and the iterative search to study the modeling of multifrequency signals based on measured data [34]. The signal modeling methods based on statistical identification are proposed by means of the dynamical window discrete measured data [35]. Ding proposed a method for identifying the system model parameters and the noise model parameters for stochastic systems described by the CARARMA models [36]. Xu et al. used the hierarchical principle and the parameter decomposition strategy to develop a parameter estimation algorithm for linear continuous-time systems [37].
The structure of the paper is organized as follows: Section 2 introduces the technique of SMB chromatography separation process. Section 3 introduces the two subspace system identification algorithms (MOESP and N4SID). Section 4 introduces the structure of the model predictive control algorithm and the predictive control method based on the state-space model. Section 5 shows the experimental simulation and results analysis. Finally, the conclusion is summarized in Section 6.
2. Technique of Simulated Moving Bed Chromatographic Separation Process
The devices in the SMB chromatographic separation process are composed of a circulation pump, filter, constant temperature system, pressure sensor, flow rate meter, rotary valve, pipeline, control part, and a plurality of separation columns [38, 39]. According to the specific requirements of the separation target compound, the number of separation columns is changed from 4 to 24, and the number of pumps is 3 to 5. The SMB system is mainly divided into the mobile phase circulation system and the stationary phase circulation system. The mobile phase circulation system makes the compounds of mixed components circulate in the columns under the action of the pressure pump through the valve of the entrance port. The corresponding target separation substance is obtained at the valve outlet of extraction solution. The stationary phase circulation system simulates the flow of the stationary phase by the cooperation of two valves. The stationary phase flows in the opposite direction compared to the direction of the mobile phase.
2.1. Import and Export Liquid of Chromatographic Separation Process
The chromatographic separation process is divided into different zones depending on the function, namely, the solid phase regeneration zone, the two separation zones, and the eluent regeneration zone. The zones of chromatographic separation process are shown in Figure 1 [22].
The first zone is located between the mobile phase inlet and the extract liquid outlet, the second zone is located between the extract liquid outlet and the entrance liquid inlet, the third zone is located between the entrance liquid inlet and the raffinate outlet, and the fourth zone is located between the raffinate outlet and the mobile phase inlet. In Section I, the mobile phase detaches the adsorbent and the separation material regeneration, that is, the solid phase regeneration zone. In Section II, the strongly adsorbed component moves relatively slowly in the column due to adsorption and is concentrated in the second region. In Section III, the weakly adsorbed component moves from the second region to the third region at a faster velocity and is aggregated in the third region. In Section IV, the mixture is separated in the second zone and the third zone, and the mobile phase is regenerated. Before the mixture is subjected to a new cycle, the mixed components must be completely desorbed to allow the stationary phase to be purified. The weak components in the separation zone are desorbed, as the liquid phase of the mobile phase moves, and the strong components are adsorbed and move with the stationary phase.
2.2. Cycle of Material Liquid in Chromatographic Column
In the SMB chromatographic separation process, only the mobile phase in the column is circulating, and the stationary phase does not move. By controlling the valves, the position of the liquid in each zone is switched to simulate the flow of the stationary phase. The state switching of chromatographic separation is shown in Figure 2 [22]. At the initial moment, the positions of the eluent, extract liquid, entrance liquid, and raffinate are shown in Switching State 1 in Figure 2. After t time, through the valve switching, the position of each inlet and outlet is shown in Switching State 2 in Figure 2. After the position switching, the result is that the stationary phase forms a relative movement with the mobile phase at a speed of L/t so as to simulate the countercurrent process of the moving bed.
(a)
(b)
2.3. Separation of Mixture
By carrying out switching the location of inlet-outlet continuously, after several cycles, the mixed compound is under the flushing of the eluent so that the sample components are gradually separated. After the separation for a period of time, the mixture in the column can be separated to some extent by adjusting and optimizing the interval of valve switching, the moving speed of the inlet-outlet, and in-out of double inlet-outlet solutions per unit time. Meanwhile, the concentration separation curve in the column is always maintained in a stable distribution state. Under cyclic operation of the SMB chromatographic system, the mixed compound is separated continuously.
The core structure of the SMB device used in this paper is shown in Figure 3. The universal SMB system is built by seven self-designed chromatographic dedicated two-way solenoid valves in each chromatographic column. There are four channels in the mobile phase inlet of each chromatographic column, which are, respectively, introduced into the circulating liquid of the former root column. The pumps 1, 2, and 3 are used to deliver the raw material liquid F, eluent D, and eluent P. There are four pipes in the mobile phase outlet, which, respectively, output the current chromatographic column circulating liquid, the weak adsorption raffinate R, the strong adsorption extract liquid E, and the ordinary adsorption substance or the feedback liquid M. The seven pipes are controlled by solenoid valves. Each pipe of outlet is equipped with a detector. The microprocessor controls the working state of the solenoid valves (“on” or “off”). Under the control of the solenoid valves, the SMB system can work on multiple working modes as required.
2.4. Selection of SMB Process Variables
The variable information of the SMB chromatographic separation process is shown in Table 1. F pump flow rate, D pump flow rate, and switching time are chosen as input variables, and target substance yield and impurity yield are chosen as output variables. Using the subspace identification algorithms, the corresponding state-space model can be established. The actual operational data of the SMB chromatographic separation process are collected, whose input data are pump flow rate, flushing pump flow rate, and switching time data and whose output data are target substance purity, impurity purity, target substance yield, and impurity yield. The number of data is more than 1400 groups. The whole input-output data set is shown in Figure 4.
3. Yield Model Identification of SMB Chromatographic Separation Process Based on Subspace System Identification Methods
The structure of the typical discrete state-space model is shown in Figure 4. The state-space model can be described aswhere , , , , , and .
The constructed input Hankel matrices are as follows:where i denotes the current moment, denotes from first row to ith row, p denotes past, and denote future.
Similarly, the output Hankel matrices can be defined as follows:
The status sequence is defined as
We obtain the following after iterating Formula (1):where V is the noise term, is the extended observability matrix, and is the lower triangular Toeplitz matrix, which are expressed asand the oblique projections is expressed as
In [1], the oblique projections has been explained. Finally, we get
Singular value decomposition for is expressed as
Get the extended observability matrix and the Kalman estimator of as . At the same time, in the process of singular value decomposition, the order of the model is obtained from the singular value of .
In [2], system matrices A, B, C, and D are determined by estimating the sequence of states. Therefore, the state-space model of the system is obtained.
The mathematical expression of the MOESP algorithm is given below.
LQ decomposition of system matrices is constructed by using the Hankel matrices:and singular value decomposition on is as follows:
The extended observability matrix is equal to
Multivariable output-error state-space (MOESP) is based on the LQ decomposition of the Hankel matrix formed from input-output data, where L is the lower triangular and Q is the orthogonal. A singular value decomposition (SVD) can be performed on a block from the L (L22) matrix to obtain the system order and the extended observability matrix. From this matrix, it is possible to estimate the matrices C and A of the state-space model. The final step is to solve overdetermined linear equations using the least-squares method to calculate matrices B and D. The MOESP identification algorithm has been described in detail in [40, 41].
The mathematical expression of the N4SID algorithm is
Given two matrices ,
The extended observability matrix is equal to
The method of numerical algorithms for subspace state-space system identification (N4SID) starts with the oblique projection of the future outputs to past inputs and outputs into the future inputs’ direction. The second step is to apply the LQ decomposition, and then the state vector can be computed by the SVD. Finally, it is possible to compute the matrices A, B, C, and D of the state-space model by using the least-squares method. N4SID has been described in detail in [42].
4. Model Predictive Control Algorithm
4.1. Fundamental Principle of Predictive Control
4.1.1. Predictive Model
The predictive model is used in the predictive control algorithm. The predictive model just emphasizes on the model function rather than depending on the model structure. The traditional models can be used as predictive models in predictive control, such as transfer function model and state-space model. The pulse or step response model, more easily available in actual industrial processes, can also be used as the predictive models in predictive control. The unsteady system online identification models can also be used as the predictive models in predictive control, such as the CARMA model and the CARIMA model. In addition, the nonlinear system model and the distributed parameter model can also be used as the prediction models.
4.1.2. Rolling Optimization
The predictive control is an optimization control algorithm that the future output will be optimized and controlled by setting optimal performance indicators. The control optimization process uses a rolling finite time-domain optimization strategy to repeatedly optimize online the control objective. The general adopted performance indicator can be described aswhere and are the actual output and expected output of the system at the time of , is the minimum output length, is the predicted length, is the control length, and is the control weighting coefficient.
This optimization performance index acts on the future limited time range at every sampling moment. At the next sampling time, the optimization control time domain also moves backward. The predictive process control has corresponding optimization indicators at each moment. The control system is able to update future control sequences based on the latest real-time data, so this process is known as the rolling optimization process.
4.1.3. Feedback Correction
In order to solve the control failure problem caused by the uncertain factors, such as time-varying, nonlinear, model mismatch, and interference in actual process control, the predictive control system introduces the feedback correction control link. The rolling optimization can highlight the advantages of predictive process control algorithms by means of feedback correction. Through the above optimization indicators, a set of future control sequences are calculated at each control cycle of the prediction process control. In order to avoid the control deviation caused by external noises and model mismatch, only the first item of the given control sequence is applied to the output, and the others are ignored without actual effect. At the next sampling instant, the actual output of the controlled object is obtained, and the output is used to correct the prediction process.
The predictive process control uses the actual output to continuously optimize the predictive control process so as to achieve a more accurate prediction for the future system behavior. The optimization of predictive process control is based on the model and the feedback information in order to build a closed-loop optimized control strategy. The structure of the adopted predictive process control is shown in Figure 5.
4.2. Predictive Control Based on State-Space Model
Consider the following linear discrete system with state-space model:where is the state variable, is the control input variable, is the controlled output variable, and is the noise variable.
It is assumed that all state information of the system can be measured. In order to predict the system future output, the integration is adopted to reduce or eliminate the static error. The incremental model of the above linear discrete system with state-space model can be described as follows:where , , and .
From the basic principle of predictive control, it can be seen that the initial condition is the latest measured value. is the setting model prediction time domain. is the control time domain. The measurement value at current time is . The state increment at each moment can be predicted by the incremental model, which can be described aswhere is the prediction of time at time . The controlled outputs from time to time can be described as
Define the step output vector as
Define the step input vector as
The equation of step predicted output can be defined aswhere
The given control optimization indicator is described as follows:
Equation (21) can be rewritten in the matrix vector form:where , , and .
The reference input sequence is defined as
Then, the optimal control sequence at time can be obtained bywhere
From the fundamental principle of predictive control, only the first element of the open loop control sequence acts on the control system; that is to say,
Let
Then, the control increment can be rewritten as
4.3. Control Design Method for SMB Chromatographic Separation
The control design steps for SMB chromatographic separation are as follows: Step 1. Obtain optimized control parameters for the MPC controller. The optimal control parameters of the MPC controller are obtained by using the impulse output response of the yield model under different parameters. Step 2. Import the SMB state-space model and give the control target and MPC parameters, calculating the control action that satisfied the control needs. Finally, the control amount is applied to the SMB model to get the corresponding output.
The MPC controller structure of SMB chromatographic separation process is shown in Figure 6.
5. Simulation Experiment and Result Analysis
5.1. State-Space Model Based on Subspace Algorithm
5.1.1. Select Input-Output Variables
The yield model can be identified by the subspace identification algorithms by utilizing the input-output data. The model input and output vector are described as and , where are, respectively, flow rate of F pump, flow rate of D pump, and switching time, is the target substance yield, and is the impurity yield.
5.1.2. State-Space Model of SMB Chromatographic Separation System
The input and output matrices are identified by using the MOESP algorithm and N4SID algorithm, respectively, and the yield model of the SMB chromatographic separation process is obtained, which are the 3rd-order MOESP yield model, the 4th-order MOESP yield model, the 3rd-order N4SID yield model, and the 4th-order N4SID yield model. The parameters of four models are listed as follows.
The parameters for the obtained 3rd-order MOESP yield model are described as follows:
The parameters for the obtained 4th-order MOESP yield model are described as follows:
The parameters for the obtained 3rd-order N4SID yield model are described as follows:
The parameters for the obtained 4th-order N4SID yield model are described as follows:
5.1.3. SMB Yield Model Analysis and Verification
Select the first 320 data in Figure 4 as the test sets, and bring in the 3rd-order and 4th-order yield models, respectively. Calculate the respective outputs, and compare the differences between the model output values and the actual values. The result is shown in Figure 7.
(a)
(b)
(c)
(d)
It can be seen from Figure 7 that the yield state-space models obtained by the subspace identification method have higher identification accuracy and can accurately fit the true value of the actual SMB yield. The results verify that the yield models are basically consistent with the actual output of the process, and the identification method can effectively identify the yield model of the SMB chromatographic separation process. The SMB yield model identified by the subspace system provides a reliable model for further using model predictive control methods to control the yield of different components.
5.2. Prediction Responses of Chromatographic Separation Yield Model
Under the MOESP and N4SID chromatographic separation yield models, the system prediction performance indicators (control time domain , predicted time domain , control weight matrix , and output-error weighting matrix ) were changed, respectively. The impulse output response curves of the different yield models under different performance indicators are compared.
5.2.1. Change the Control Time Domain
The control time domain parameter is set as 1, 3, 5, and 10. In order to reduce fluctuations in the control time-domain parameter, the remaining system parameters are set to relatively stable values, and the model system response time is 30 s. The system output response curves of the MOESP and N4SID yield models under different control time domains are shown in Figures 8(a) and 8(b). It can be seen from the system output response curves in Figure 8 that the different control time domain parameters have a different effect on the 3rd-order model and 4th-order model. Even though the control time domain of the 3rd-order and 4th-order models is changed, the obtained output curves y1 in the MOESP and N4SID output response curves are coincident. When m = 1, 3, 5, although the output curve y2 has quite difference in the rise period between the 3rd-order and 4th-order response curves, the final curves can still coincide and reach the same steady state. When m = 10, the output responses of MOESP and N4SID in the 3rd-order and 4th-order yield models are completely coincident; that is to say, the output response is also consistent.
(a)
(b)
5.2.2. Change the Prediction Time Domain
The prediction time-domain parameter is set as 5, 10, 20, and 30. In order to reduce fluctuations in the prediction time-domain parameter, the remaining system parameters are set to relatively stable values, and the model system response time is 150 s. The system output response curves of the MOESP and N4SID yield models under different predicted time domains are shown in Figures 9(a) and 9(b).
(a)
(b)
It can be seen from the system output response curves in Figure 9 that the output curves of the 3rd-order and 4th-order MOESP yield models are all coincident under different prediction time domain parameters; that is to say, the response characteristics are the same. When , the output response curves of the 3rd-order and 4th-order models of MOESP and N4SID are same and progressively stable. When , the 4th-order y2 output response curve of the N4SID model has a smaller amplitude and better rapidity than other 3rd-order output curves.
5.2.3. Change the Control Weight Matrix
The control weight matrix parameter is set as 0.5, 1, 5, and 10. In order to reduce fluctuations in the control weight matrix parameter, the remaining system parameters are set to relatively stable values, and the model system response time is 20 s. The system output response curves of the MOESP and N4SID yield models under different control weight matrices are shown in Figures 10(a) and 10(b). It can be seen from the system output response curves in Figure 10 that under the different control weight matrices, although y2 response curves of 4th-order MOESP and N4SID have some fluctuation, y2 response curves of 4th-order MOESP and N4SID fastly reaches the steady state than other 3rd-order systems under the same algorithm.
(a)
(b)
5.2.4. Change the Output-Error Weighting Matrix
The output-error weighting matrix parameter is set as [1 0], [20 0], [40 0], and [60 0]. In order to reduce fluctuations in the output-error weighting matrix parameter, the remaining system parameters are set to relatively stable values, and the model system response time is 400 s. The system output response curves of the MOESP and N4SID yield models under different output-error weighting matrices are shown in Figures 11(a) and 11(b). It can be seen from the system output response curves in Figure 11 that when , the 4th-order systems of MOESP and N4SID has faster response and stronger stability than the 3rd-order systems. When , the 3rd-order and 4th-order output response curves of MOESP and N4SID yield models all coincide.
(a)
(b)
5.3. Yield Model Predictive Control
In the simulation experiments on the yield model predictive control, the prediction range is 600, the model control time domain is 10, and the prediction time domain is 20. The target yield was set to 5 yield regions with reference to actual process data: [0, 0.4], [0.4, 0.5], [0.5, 0.7], [0.7, 0.8], [0.8, 0.9], and [0.9, 0.99]. The 4th-order SMB yield models of MOESP and N4SID are, respectively, predicted within a given area, and the output curves of each order yield models under different input are obtained by using the model prediction control algorithm, which are shown in Figures 12(a) and 12(b). It can be seen from the simulation results that the prediction control under 4th-order N4SID yield models is much better than that under 4th-order MOESP models. The 4th-order N4SID can bring the output value closest to the actual control target value. The optimal control is not only realized but also the actual demand of the model predictive control is achieved. The control performance of the 4th-order MOESP yield models has a certain degree of deviation without the effect of noises. The main reason is that the MOESP model has some identification error, which will lead to a large control output deviation. The 4th-order N4SID yield model system has a small identification error, which can fully fit the output expected value on the actual output and obtain an ideal control performance.
(a)
(b)
6. Conclusions
In this paper, the state-space yield models of the SMB chromatographic separation process are established based on the subspace system identification algorithms. The optimization parameters are obtained by comparing their impact on system output under the model prediction control algorithm. The yield models are subjected to corresponding predictive control using those optimized parameters. The 3rd-order and 4th-order yield models can be determined by using the MOESP and N4SID identification algorithms. The simulation experiments are carried out by changing the control parameters, whose results show the impact of output response under the change of different control parameters. Then, the optimized parameters are applied to the predictive control method to realize the ideal predictive control effect. The simulation results show that the subspace identification algorithm has significance for the model identification of complex process systems. The simulation experiments prove that the established yield models can achieve the precise control by using the predictive control algorithm.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
Zhen Yan performed the main algorithm simulation and the full draft writing. Jie-Sheng Wang developed the concept, design, and critical revision of this paper. Shao-Yan Wang participated in the interpretation and commented on the manuscript. Shou-Jiang Li carried out the SMB chromatographic separation process technique. Dan Wang contributed the data collection and analysis. Wei-Zhen Sun gave some suggestions for the simulation methods.
Acknowledgments
This work was partially supported by the project by National Natural Science Foundation of China (Grant No. 21576127), the Basic Scientific Research Project of Institution of Higher Learning of Liaoning Province (Grant No. 2017FWDF10), and the Project by Liaoning Provincial Natural Science Foundation of China (Grant No. 20180550700).