Least Squares Support Vector Machine-Based Multivariate Generalized Predictive Control for Parabolic Distributed Parameter Systems with Control Constraints

This manuscript addresses a new multivariate generalized predictive control strategy using the least squares support vector machine for parabolic distributed parameter systems. First, a set of proper orthogonal decomposition-based spatial basis functions constructed from a carefully selected set of data is used in a Galerkin projection for the building of an approximate low-dimensional lumped parameter systems. Then, the temporal autoregressive exogenous model obtained by the least squares support vector machine is applied in the design of a multivariate generalized predictive control strategy. Finally, the effectiveness of the proposed multivariate generalized predictive control strategy is verified through a numerical simulation study on a typical diffusion-reaction process in radical symmetry.


Introduction
Modeling and control issues in the process industry involving distributed parameter systems (DPSs), which are infinite-dimensional systems, are challenging. The identification of unknown parameters of DPSs is much more complex than that of the lumped parameter systems (LPSs). In general, an infinite-dimensional problem is required to be approximated by finite-dimensional LPSs. Then, the model reduction methods are used to reduce the dimension of the approximate LPSs. For further details, see [1].
A system governed by parabolic partial differential equations (PDEs) is a typical DPS. It can be approximated by a set of low-dimensional ordinary differential equations through the use of Proper Orthogonal Decomposition (POD) technique [2], which is also called Karhunen-Loève expansion (K-L expansion) [3] or Principal Component Analysis (PCA) [4]. The POD technique is a popular technique for carrying out the task of model order reduction of parabolic systems. The spatial basis functions constructed from those slow modes of the system can capture the dominant dynamics while those fast modes can be ignored [5]. For the identification of the unknown parameters of the approximate LPSs constructed from the low-dimensional time-series, conventional approaches, such as projection transformation, can be used. The Singular Value Decomposition and K-L expansion is applied to construct a low-dimensional LPS with a simple structure. However, the accuracy of the approximation is low for a nonlinear DPS [6,7]. The K-L expansion combined with the neural network method can be applied to construct a sufficiently accurate approximation for a nonlinear DPS [8][9][10][11][12][13][14]. The support vector machine (SVM) is a general identification approach, and the least-squares support vector machine (LS-SVM) is its modified version. The LS-SVM, which simplifies the calculation by replacing inequality constraints with equality constraints, inherits the characteristics of the global optimization property and high generalization capabilities of SVM [15]. The model reduction approach based on the K-L expansion combined with LS-SVM is applied to regression problems that do not have a large amount of sample data. It is potentially applicable to a wide range of DPSs [16,17]. Several improved PCA methods are available in the literature, such as the incremental POD method [18], the adaptive POD method [19], the adaptive K-L/Fuzzy method [20], and the nonlinear-PCA method [21], for dealing with more complex nonlinear DPSs. For the nonlinear parabolic DPSs with time-dependent spatial domains, a temporallylocal model order-reduction technique [22] and sparse POD-Galerkin methodology [23] are proposed and have been successfully applied to handle the moving boundary problem for a hydraulic fracturing process.
Generalized predictive control (GPC) strategy is introduced in [24,25]. Due to its simplicity for implementation, and excellent performance and robustness, it has been successfully implemented in a diversity of industrial applications from single-variable to multi-variable, and linear to nonlinear [26].
In this paper, we present a new multivariate GPC strategy used in conjunction with the POD technique to reduce the model order of weakly nonlinear DPSs of parabolic type. Then, the LS-SVM method is used to identify the unknown parameters of the lowdimensional time-series. Detailed methodology used to construct the generalized predictive temperature control in the tubular chemical reactor [27] will be discussed, and the control constraints and stability analysis will also be considered in the paper. The GPC based on SVM method has been successfully applied to LPSs with nonlinear characteristics [28][29][30]. The predictor obtained using the SVM method combined with the POD technique is found to have higher accuracy. It is applicable to a wider range of nonlinear DPSs. On the other hand, the GPC strategy used in conjunction with POD and recursive least square method has been applied to DPSs [31][32][33]. Thus, it is advantageous to combine GPC with SVM for the modeling and control of DPSs.
The remainder of this paper is structured as follows. First, the system equations are described in Section 2. Then, the POD and LS-SVM for regression are briefly introduced, and the POD and LS-SVM-based predictive model is derived in Section 3. The POD and LS-SVM-based multivariate GPC strategy is presented in Section 4. In Section 5, the numerical simulation of a diffusion-reaction process in radical symmetry is carried out to examine the effectiveness of the proposed algorithm in Section 3. The summary and future work of our paper are included in Section 6.

System Description
We consider a class of weakly nonlinear PDEs of parabolic type as given below: where t is the time variable, x ∈ Ω is the spatial variable, and Ω is the spatial domain, y(x, t) is the controlled output variable, u(t) is the control input variable, A and B are the linear operators in Hilbert space H , and f (·) is an unknown weakly nonlinear function. System (1) is applicable to model a variety of physical and biochemical processes, such as catalytic reaction rods, steel casting, and the snap oven [5].
For system (1), it contains an unknown weakly nonlinear function f (·). To obtain accurate information of such a system, a sufficient number of sensors need to be placed along the spatial location. Due to actual physical conditions, only a small number of actuators is allowed to be mounted to observe the state. The input-output data are obtained from the actual production process under random signal excitation. The control algorithm is developed in three stages: In Section 3, the output {y(x m , t n )} M,N m=1,n=1 excited by the input {u(k, t n )} K,N k=1,n=1 are used by the POD technique to reduce the dimensionality of the approximate model, and M, N, and K are, respectively, the number of the sensors, the sample duration, and the actuators. A low-dimensional time-series obtained after model order reduction is identified by LS-SVM. In Section 4, the POD used in conjunction with the LS-SVM-based multivariate GPC strategy is utilized to design the required controller through solving a receding horizon optimization problem at every sample period.

Model Reduction by POD and LS-SVM
The POD and LS-SVM-based model will be established in this section. Since a high dimensional output is not ideal for the design of the controller, the POD technique will be applied first to reduce the dimension of the model. Let the input u(·, t) ∈ K×N , and let y(x m , t n ) ∈ M×N . The POD basis of rank can be calculated by Algorithm 1 proposed in [2].

Algorithm 1 Proper Orthogonal Decomposition (POD) basis of rank
where is the rank of the POD basis. Step Step 2. Determine R = YY T ∈ M×M .
Step 4. Set Step Choose the minimum such that E ( ) ≥ 99%. Then, where Ψ (i) (x) and Y (i) (t), i = 1, · · · , , represent, respectively, the i-th basis function and the corresponding low-dimensional time-series, Through the use of the LS-SVM approach proposed in [15], the i-th estimated value of the time-series is: where is a function which maps the vector X (i) from function space into feature space. The optimization problem can be expressed as: where γ (i) is the regularization factor, and ξ (i) is the relaxation factor. According to Lagrange multiplier method, the Lagrange function is: where α (i) is the Lagrange operator. Let ∂L /∂w (i) = ∂L /∂b (i) = ∂L /∂e (i) = ∂L /∂α (i) = 0, the coefficients α (i) and b (i) are obtained through solving the following linear algebraic equations:  The time-series Y (i) can be obtained as: where The matrix form of the time-series autoregressive exogenous (ARX) model obtained by the LS-SVM approach is: where

POD and LS-SVM-Based Multivariate GPC
GPC strategy proposed in [23,24] is a class of predictive control algorithm developed for adaptive control due to its excellent performance and simplicity for implementation. It has been applied in a diversity of actual industrial plants during the past three decades. A multivariate GPC strategy will be proposed for parabolic DPSs with unknown coefficients in this section.

Multivariate GPC Strategy
Consider the following ARX model containing a zero-mean white noise: where the output Y(t) ∈ × 1, the input U(t) ∈ K × 1, and the white noise with zero mean e(t) ∈ × 1. ∆ = 1 − z −1 , C(z −1 ) is the identity matrix. Minimizing the quadratic cost function on a finite horizon to obtain the optimal control increment as follows: whereŶ(t + j|t) is the output predictive j-step ahead value based on previous state up to time t. λ is the weighted coefficient of ∆u. N p and N u are, respectively, the prediction horizons and control horizons. The reference trajectory w(t + j) is defined as: where y sp is the set point value. Consider the Diophantine equation: where A = A∆. The order of the polynomial matrix E j (z −1 ) is j − 1, and that of the By solving E j (z −1 ) recursively, the N p -step prediction valueŶ can be written as: . . .
. . . where The matrix form of the prediction valueŶ is: Rewrite the quadratic criterion (11) as: If there are no constraints, the optimum ∆u without constraint can be obtained as: where H = (G T G + λI) −1 . It is difficult to avoid control constraints in practical systems. The constrained incremental control law directly truncated at the boundary is not optimal. Considering the control constraints Θ∆u ≤ θ, the quadratic criterion using a Lagrange multiplier is given as follows [31]: where ε is the weighting marix of control constraint which is positive definite. The optimum ∆ũ with constraint can be obtained as: The quadratic criterion is solved repeatedly at the subsequent sampling times, only the first term of ∆ũ is applied to the actual system, so the control action is: where k, k = 1, · · · , K, is defined in Section 2. The constrained control action is handled by the truncated algorithm: where u min ≤ u(k, t j ) ≤ u max . The proposed POD and LS-SVM-based multivariate GPC algorithm can be summarized in Algorithm 2.

Algorithm 2 The POD and least-squares support vector machine (LS-SVM)-based multivariate generalized predictive control (GPC).
Require: A set of output y(x m , t n ) is derived by appropriate excitation signals.
Step 2. The predictorŶ and the ARX model are built by LS-SVM with a linear kernel.
Step 4. Use the control law u(k, t j ) to the actual system at instant t j .

The Stability Analysis
We assumed that the following assumptions are satisfied throughout: Assumption 1. The reference trajectory w(t + j) is bounded, and the prediction error of model (8) is bounded. (1) is open-loop stable, so that the influence of the control action exerted at the past time tends to zero with the increase of time.

Assumption 2. System
Assumption 3. The fast modes of the dynamic of system (1) are ignored, so that it suffices to consider the stability conditions of the nominal system (9). Assumption 4. The future disturbances are equal to that at time t. The mismatch and any uncertainty of model (15) can be estimated as: where Y(t + 1) = Ψ T (x)y(x, t + 1), and 0 ≤ ≤ 1.
Lemma 2. (Constrained case) [6] If the eigenvalue λ i of the matrix A satisfy | λ i | < 1, then the closed-looped system is asymptotically stable, where Proof. The proof is similar to that given for Lemma 1.

Case Study
The proposed multivariate GPC strategy is applied in the simulation study of an exothermic catalytic reaction in a long thin rod in radical symmetry, which has the form A → B [34] of transport reaction as shown in Figure 1. The goal is to produce as much product B as possible while consuming the least amount of energy by adjusting the axial cooling medium temperature. Assume that the density and the heat capacity of the rod are constant, and that the reaction rate of the rod is uniform. Due to radial symmetry, the dynamical distribution of the axial temperature on the rod is governed by the following PDE of parabolic type: subject to the first kind of boundary conditions, and initial condition: where t, x, and L are, respectively, the dimensionless time, axial position, and the length of reactor, and the controlled variables T(x, t) and the control variables T c (x, t) are, respectively, the axial temperature and the cooling medium temperature of the reactor. In the simulation, let the length of the rod L = π, the initial temperature T 0 (x) = 0.5, the heat of reaction β T = 50, the heat transfer coefficient β u = 2, and the activation energy γ = 4. The 21 temperature sensors and 4 cooling medium temperatures T c (x, t) are equally located along the long rod. Each sample period ∆t = 0.01, and the total period is t s = 2.
All the actual simulation values are calculated using Matlab PDE toolbox. The random excitation input signal that varies randomly between −4 to 1 are shown in Figure 2 and the corresponding output is shown in Figure 3.
Different structures of the low-dimensional time-series model are considered for our numerical experiments, and the orders of ny = 4 and nu = 3 are finally selected for which the desired accuracy of the corresponding reduced dimensional model is achieved. The parameters N = 160, and K = 4. The 10-fold cross-validation is used to obtained the regularization factor γ. The LS-SVM-based model can approximate the low-dimensional time-series satisfactorily which is shown in Figure 5. The error of the predictive model is between [−0.02, 0.015] which is shown in Figure 6, and the mean of the root mean square error is 0.0001. The average training times of the model are [16.4, 5.5, 20.5, 2.8] ms, respectively.   For expressions (11) and (12), choose λ = 1, α w = 0.3, N p = 3, N u = 2 as the parameters of the controller. The unconstrained control increment and control action are calculated by expressions (19) and (22), and the first column of the control variable T c is used in the actual system. The process is repeated at every sample period. The control variables T c and the output response curves for the unconstrained case are shown in Figures 7 and 8, respectively.
In Figures 7 and 8, the controller can quickly drive the system dynamics back to the steady-state in 30 sample periods, and the steady-state error is below 1 × 10 −5 . From the simulation results, we can see that the proposed controller has excellent tracking performance and can be effectively applied to parabolic DPSs involving unknown coefficients. The fluctuation range of control variables between −8 to 1.5 is large.
We now consider the case of control constraints. More specifically, −0.5 ≤ ∆u ≤ 0.5, and −4 ≤ u ≤ 1. The constrained control increment and control action are updated by expressions (21) and (23). The control variables T c and the output response curves under the constrained case are shown in Figures 9 and 10, respectively.   In Figures 9 and 10, we can see that the controller is working well in the case of control constraints. Due to the presence of the control constraints, the settling time is now extended to 50 sample periods. The time taken for control action calculation in each sampling period is given in Figure 11.
In Figure 11, we can see that the time taken for unconstrained case is in milliseconds, and it is a bit longer for the constrained case.

Conclusions
A new control strategy, named POD and LS-SVM-based multivariate GPC, was developed for weakly nonlinear parabolic DPSs with unknown coefficients. The POD was used for model order reduction, resulting in low-dimensional lumped parameter systems. The desired control was designed based on this low-dimensional model. The LS-SVM with linear kernel function is found effective in the identification issues for linear and weakly nonlinear plants. An interesting question to ask is: Are there kernel functions for general nonlinear systems? The POD used in conjunction with the LS-SVM-based multivariate GPC strategy is easy to implement and is applicable to weakly nonlinear parabolic systems with unknown coefficients. As an illustration, a diffusion-reaction process was considered and solved using the proposed methods . The controller obtained has achieved the desired good performance. For future research, the online LS-SVM and radial basis kernel function will be applied to improve the modeling accuracy and to a more complex problem with nonlinearity. Furthermore, the method will also be extended to handle moving boundary problems, where the single POD basis function is replaced by a temporally-local basis function.