Abstract

This paper proposes a recursive delta-operator-based subspace identification method with fixed data size. The majority of existing subspace identification methods are constrained to discrete-time systems because of the disparity in Hankel matrices. Additionally, due to the storage cost, LQ-factorization and singular value decomposition in identification methods are best suited for batch processing rather than online identification. The continuous-time systems are transformed into state space models based on the delta-operator to address these issues. These models approach the original systems when the sampling interval approaches zero. The size of the data matrices is fixed to reduce the computing load. By fading the impact of past data on future data, the amount of data storage can be decreased. The effectiveness of the proposed method is illustrated by the continuous stirred tank reactor system.

1. Introduction

In recent decades, subspace identification method has become an important part in the area of identification and control [13]. The discrete-time systems have been the main focus of the identification approaches. Indeed, the mathematical models of continuous-time systems are differential equations, which constitute the majority of existing systems. The numerical techniques like LQ-factorization and singular value decomposition form the foundation for subspace identification. These tools are suitable for time-invariant systems, but due to the high computing weight, they are not suitable for the identification of time-varying systems. Therefore, a new recursive identification approach for continuous-time systems is very desirable.

The direct method and the indirect method are the two basic methods for identifying continuous-time systems. The direct methods determine the system model from the data, and the indirect methods obtain the system model after obtaining the discrete-time systems. Actually, the key to the identification methods is how to deal with the time-derivative issues of the continuous-time systems. The time-derivative problems have been addressed in a variety of ways [4]. The continuous-time model was obtained by utilizing the bilinear connection rather than determining the input-output operator [5]. In [6], the subspace identification approach and Poisson moment functionals (PMFs) were used to handle the time-derivative problems of the signals. The input-output matrix equation in the time-domain was found by using random distribution theory to characterize the time-derivative [7]. The identification approaches of a toroidal continuously variable transmission were investigated in [8], and the study introduced linear integral filter methods and PMF to obtain the submodel parameters. The Laguerre filter-based subspace identification method for continuous-time systems was presented in [9]. The identification problems of fractional commensurate order systems from the nonuniformly data were solved [10]. The continuous-time state space models identification method was obtained by using generalized orthogonal basis functions [11]. The continuous-time system identification approach with missing outputs was investigated by using PMF and nuclear norm minimization [12]. A method of nuclear norm subspace identification based on Kalman filter for the stochastic continuous-time system is proposed in [13]. In [14], the recursive subspace identification method of continuous-time systems via generalized Poisson moment functional is presented.

In the subspace identification process, LQ-factorization and singular value decomposition (SVD) are reliable tools [4, 15]. However, they are computationally expensive in online estimation. To handle the issues with online identification, recursive subspace identification is widely investigated. The traditional recursive identification algorithm was presented by the propagator method after analyzing the connection between array signal processing and subspace identification [16]. By using the least squares support vector machines, the recursive Hammerstein system identification was proposed [17]. The constrained recursive least squares identification method with the forgetting factor was provided [18]. The recursive identification approach was obtained by integrating canonical correlation analysis [19]. In [20], the load disturbance response and the observer Markov parameter matrices were obtained by using the recursive least squares methods, and the identification methods were derived. The local polynomial modeling was used to investigate the two-dimensional recursive least squares identification approach [21]. The recursive tracking approach that identified the parameters of the spacecraft was described [22]. The recursive combined-deterministic stochastic subspace identification algorithm was provided, which used a system observability matrix recursively that avoids the need for SVD [23]. In [24], a new recursive least squares identification algorithm with variable-direction forgetting was proposed for multioutput systems. Based on sparse adaptive hybrid integration, Zhao and Lam [25] demonstrated the general stiffness and redundancy issue in the fast subspace. Rump and Lange [26] proposed the fast computation algorithm of error bounds for all eigenpairs of a Hermitian and all singular pairs of a rectangular matrix. Bhowmik et al. [27] presented the identification approach for structural modal parameters via recursive canonical correlation analysis.

It is clear from the aforementioned references that discrete-time systems are the research focus of the majority of the recursive identification methods. However, identification of continuous-time systems has real value in a variety of circumstances [28, 29]. In addition, the references [16, 18, 22] that mention subspace tracking include the assumption that the system order is known beforehand. Instead, this paper proposes a recursive subspace identification algorithm with the priori order unknown. Based on delta-operator-based method, it fixes the size of the data matrices and proposes a recursive identification algorithm obtained by performing the LQ-factorization recursively. The proposed method has the main contributions as follows: (1) for the online continuous-time identification problems, the state space model that is convergent to the original systems is obtained by using the delta-operator-based method and (2) to reduce the computational burden and storage cost, the size of the data matrices is fixed a priori to fade the influence of old data to the updated data.

The rest of this paper is structured as follows. The problem formulation is described in Section 2. In Section 3, the proposed recursive subspace identification algorithm is derived. Section 4 provides the application to the continuous stirred tank reactor system. Section 5 presents the conclusions.

2. Problem Statement

Consider a continuous-time system as follows:where is the state vector, is the output vector, is the input vector, is system matrix, is input matrix, is output matrix, is direct transmission matrix, and and are white Gaussian processes. The corresponding covariance matrices are as follows:where is the Dirac delta function and is the expectation operator. Identifying the matrices , , , and via the input-output data is our purpose.

It can be clearly seen that the input and output derivatives are at least equal to -. However, the output vector and the state vector have no time-derivative, and the noises and are not differentiable. Under these circumstances, it is also impossible to adapt traditional subspace identification techniques, which are initially studied in the context of offline discrete-time models, to the identification of online continuous-time systems. Therefore, to deal with the issue of online identification of continuous-time stochastic systems, we provide a novel approach RSI-DOSM (recursive subspace identification based on delta-operator state space model).

3. Recursive Subspace Identification via Delta-Operator State Space Model

3.1. Deriving the State Space Model by Delta-Operator

We initially establish the state space model based on delta-operators to handle the identification issues of systems (1). The sampled input-output behavior of system (1) can be derived by using a straightforward antialiasing filter and zero-order hold.wherewhere is the shift operator with , is the sampling period, denotes the both sides are nearly equal, and and are stochastic disturbances; they satisfy the following equation:where cov is the abbreviation for covariance.

Theorem 1. According to equation (3) described above, the state space model of system (1) via delta-operator can be obtained.

Proof. As the sampling period goes to zero, it is found from equation (4) that . Hence, the shift operator causes a model degeneracy. The delta-operator in [30, 31] is introduced as follows:For a sufficiently small , we havewhere is the differentiation operator.
According to [32], the model based on delta-operator is given.whereThe covariance matrices are as follows:It is found that are both approximately equal to .
The spectral densities are defined to achieve sampling independence as follows:, and when . Hence, the state space model (6) via delta-operator converges to the original continuous-time system (1).
In terms of [33], the optimal filter for the model (8) is as follows:where is the Kalman gain, is the one-step predicted estimate of , and is the delta form of Riccati equation in [34, 35].
In addition, the new innovation process is described as follows:The innovation form of system (8) is given by the following equation:Hence, the state space model based on delta-operator is derived.

3.2. Obtain Input-Output Matrix Equation

Under the assumption that is larger than , the extended observability matrix for the model (14) is defined as follows:

The input data matrix is given as follows:where the matrices are the same form as .

The state vector sequence is described as follows:

According to equation (14), the extended model is given by the following equation:where

To handle the high frequency noises of input-output signals, the stable prefilter with order is first defined as follows:where are the prefilter parameters.

Based on filter (20), system (14) is transformed as follows:where , , , and are similar to .

The input-output matrix equation (18) becomes as follows:whereand are defined in a similar way.

Definewhere is the -th component of , which is similar for .

In view of prefilter (20), we get the following equation:

Equation (25) can be converted to the following equation:where

In terms of equation (27), the input-output data matrices can be obtained.

Remark 2. It should be pointed out that, compared with the identification methods in [3638], the computation of matrix logarithm is avoided by using the delta-operator state space model (8).

Remark 3. It is important to note that the delta-operator method has fewer design parameters than the more general continuous-time subspace identification methods in [7, 11, 12].

3.3. Recursive Estimation of the Extended Observability Matrix

Let and the instrumental variable . The matrix is similarly defined as .

Given a set of new output data, construct data vectors , , and as updated data. The relation is described as follows:

Theorem 4. Based on equation (28), the projection matrix constructed by input-output data can be obtained by the LQ-factorization and matrix transformation.

Proof. Perform the LQ-factorization on the left submatrix in the left-hand side of equation (28), it gives the following equation:Then, perform the LQ-factorization on the right submatrix in the right-hand side of equation (28), it gives the following equation:Substituting equations (29) and (30) into the two sides of equation (28) and according to the relation , it is found thatThen, we have the relations as follows:In terms of equation (33), we get the following equation:According to equation (31), we have the following equation:According to equations (39)–(41), we have the following equation:In view of equation (31), we obtain the following equation:Perform the eigenvalue decomposition on , we get the following equation:According to equations (43) and (44), the matrices , , , and can be estimated recursively. We can determine the system order by inspecting the eigenvalue value in equation (46). The order of the system is equal to the number of eigenvalue values different from zero.
In terms of the literature [39], we have the following equation:where denotes the matrix of without the last rows.
According to equation (47), the state sequence is obtained. The matrices , , , and can be solved by the following equation:where is a block Hankel matrix with only one row of outputs, is the state sequence, and and are the Kalman filter residuals. Intuitively, since the Kalman filter residuals and are uncorrelated with the state, it seems natural to solve this set of equations in a least squares sense.
For clarity, the proposed method RSI-DOSM is shown in Table 1.

Remark 5. It is worth mentioning that the main advantages of the proposed method RSI-DOSM are twofold: (i) Optimal filter (12) converges to the original filter for system (1) when the sampling period goes to zero, which circumvents the identification difficulty from the time-derivative. (ii) The size of the data matrices is fixed a priori to perform the LQ-factorization recursively, which is a key to reduce computational burden of the identification algorithm.

4. Application to Continuous Stirred Tank Reactor

A simulated continuous stirred tank reactor (CSTR) system is used to validate the efficacy and accuracy of RSI-DOSM [40, 41]. The process flow diagram is depicted in Figure 1.

The reactor’s main purpose is to deliver the concentration of the outlet effluent at a prescribed value, by manipulating the coolant flow rate circulating in the reactor’s jacket. A first-order irreversible reaction in the CSTR system is assumed. One component , which is an outflow stream, is produced by the reactant and solvent’s incoming flow. The heat produced by the reaction is transferred via the jacket’s cooling flow. The temperature is maintained by adjusting the amount of heat transferred through the reactor jacket. To control the reactor’s liquid level and temperature, the cascade control approach is used. It is a challenge for the process of system identification because of the complex nonlinear dynamics and the temperature-dependent rate constant as well as the coupled dynamics between temperature and concentration.

The dynamic model of the CSTR system can be written by using the fundamental idea of the balances of components, mass, and energy as a foundation.where is the reactor’s cross-sectional area, is the activation energy, is pre-exponential factor, is the content’s density, is the gas constant, is the heat capacity of coolant, is the heat capacity of the contents, is the reaction heat, is the total heat transfer area, is the heat transfer coefficient, and is the coolant’s density.

Table 2 lists the system’s remaining specifications. It lists the associated standard deviations of the Gaussian noise utilized in the simulation process. Data sets are produced using the remaining variables in Table 2 as inputs and as the measured outputs in the CSTR simulation system’s usual operating mode. Table 2 also presents the corresponding Gaussian noise standard deviations. The first 200 samples of the data set are used to build the model. The scaled data have a unit variance and zero mean.

Two forms of natural changes between and are used to describe the time-varying properties: (1) a slow variation in with a slope of and (2) a slow variation of catalyst deactivation with 3 . We set each number of future and past block rows .

The prediction outputs of RSI-DOSM are shown in Figure 2. The measured data from CSTR and the predicted data from RSI-DOSM are displayed, respectively, by the red and black lines. It can be found that the predicted trajectory obtained by RSI-DOSM and the measured trajectory are quite consistent.

In contrast, the estimated results are also obtained using the recursive subspace identification based on Laguerre filters (RSILF) approach [42]. The predicted outputs of RSILF are shown in Figure 3. The red and black lines represent the CSTR measured data and RSILF predicted data, respectively. Figure 4 displays the prediction error between the measured and predicted values for each of the output figures from RSI-DOSM and RSILF. The RSI-DOSM and RSILF prediction errors are represented by the black and blue lines, respectively. It can be demonstrated that RSI-DOSM has more consistency.

Two distinct performance metrics, the variance accounted for (VAF) criterion and the mean square error (MSE), are taken into consideration in order to compare the RSI-DOSM and RSILF models [39, 43]. The definitions of the MSE and VAF are as follows:where is the prediction output of the system and is the data length.

The results are summarized in Table 3.

Clearly, for the two separate criteria taken into consideration, the RSI-DOSM model’s predictive ability exceeds that of the RSILF model. In particular, the improvements of 26%, 40%, 35%, and 18% in the MSE of output variables are obtained by resorting to the RSI-DOSM model as an alternative to RSILF model, respectively. It demonstrates that RSI-DOSM’s prediction performance is superior to RSILF’s prediction performance.

5. Conclusions

A recursive subspace identification method based on delta-operator state space model (RSI-DOSM) is proposed. The state space model that is converged to the original systems is derived utilizing the delta-operator. The size of the data matrices is chosen a priori to minimize the impact of historical data and reduce the calculation cost of recursive algorithms. The results of the simulation suggest that the proposed identification algorithm has the improvements of 26%, 40%, 35%, and 18% in the MSE of output variables. It is important to extend the proposed method to eigen perturbation techniques, which will be the future work.

Data Availability

The simulation data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (NSFC) (no. 62003082), the Natural Science Foundation of Hebei Province (no. F2021501018), and the Science and Technology Project of Hebei Education Department (no. ZD2022148).