Continuous-time (CT) system identification plays an important role in many application problems because most dynamical systems encountered in the physical world are CT in nature^{[1, 2, 3, 4]}. In general,the approaches for CT system identification can be classified into two categories: the indirect approach and the direct approach. The indirect approach uses discretization and discrete-time (DT) identification methods to obtain DT model first and then converts it to CT form,while the direct approach identifies the CT model directly from the input and output measurement data. Although the indirect approach has dominated the CT system identification for some time with benefits of well-performing DT system identification methods developed over last decades,there are still difficulties in this approach due to primarily the estimation of equivalent DT model and the conversion of equivalent DT model to its CT counterpart. Difficulties,such as nonminimum phase and undesirable sensitivity problems,are encountered when the sampling time interval is either too small or too large^[2],which makes the indirect approach difficult to implement in some practical applications.

The direct approach can overcome the difficulties associated with the indirect approach,but it needs to compute the time-derivatives of input/output (I/O) variables,which is a nontrivial and delicate task in the presence of measurement noise^{[2, 3, 5]}. Many methods have been devised to circumvent the difficulty in reconstructing the time-derivatives,and a concise survey can be found in ^[2]. The generalized Poisson moment functional (GPMF)^{[1, 6, 7]} is one of the efficient methods to circumvent the difficulty of time-derivative and improve the accuracy of CT identification. GPMF belongs to the category of linear filtering method which reconstructs the time derivative of I/O data with Poisson moment function. Its digital implementation and parameter setup are detailed in ^[1]. GPMF has been extensively applied to system identification,industrial control and micro-electronics areas^{[1, 8, 9]}.

As a convex heuristic of the rank minimization problem,nuclear norm minimization (NNM) has recently been extensively studied^{[10, 11, 12, 13, 14]}. NNM is guaranteed to find the minimum rank matrix under suitable assumptions and can be solved by the iterative reweighted algorithms^[15] or the interior point method^[10]. Due to the convex relaxation of rank minimization^{[13, 16]},NNM has been used in SIM to minimize the rank of Hankel matrix and determine the order of the model^{[10, 12, 17, 18, 19, 20]}. As we all know,subspace identification method (SIM)^[21],is an effective and easy-to-implement identification method for linear time-invariant (LTI) system. In theory,the system order is equal to the rank of Hankel matrix which,in turn,equals the number of nonzero singular values of the matrix. However,the I/O variables in practice often contain noise that perturbs the singular values of Hankel matrix and makes the order determination of SIM unreliable.

The NNM can overcome the difficulty of model order determination in SIM and obtain a correct order identification when the priori order is unknown. However,by now the results of NNM based subspace system identification are confined to the DT models and do not carry over easily to CT case,due to the difference between the DT and CT Hankel matrices. Furthermore,the time-derivative in continuous-time SIM (CSIM) can bring in extra errors due to the measurement noise in the I/O data. How to correctly identify the CT system where the model order is not known a priori is still a big challenge in system identification field.

In light of the above observations,this paper proposes a novel method that incorporates NNM into the GPMF-based subspace method for the CT system identification. The proposed method not only circumvents the difficulty in the computation of time-derivatives of the noise corrupted I/O variables,but also alleviates the difficulty of model order estimation in the CT subspace identification. The algorithm for solving the NNM in CT SIM is deduced by ADMM and NNM is extended to CT cases for the first time.

The paper is organized as follows. Section Ⅱ reviews the basic DT SIM and its NNM. Section Ⅲ introduces the theory and algorithm of the GPMF method. Section Ⅳ analyzes the NNM and presents the framework of the novel hybrid algorithm. In Section Ⅴ,two typical numerical examples are presented to evaluate the GPMF based subspace identification with NNM,and concluding remarks are given in Section Ⅵ.

The notations used in the paper are standard. The estimate of a variable $x$ is denoted by $\hat{x}$. For a matrix $X$,$\|X\|_F$ denotes its Frobenius norm,$\|X\|_*$ its nuclear norm and $X^{\bot}$ its orthogonal projection. $s$ represents the Laplace variable. $R$ denotes the set of all real numbers,${\bf R}^+$ the set of all positive real numbers,and $R^{+*}$ the projectively extended positive real numbers. $I_q \in {\bf R}^{q\times q}$ is the identity matrix and $0_q \in {\bf R}^{q\times q}$ is the null matrix where all elements are zeros.

Consider a linear time-invariant DT state space system^[22]:

where $u(k) \in R^m$ and $y(k) \in R^p$ are respectively the input and output,$x(k) \in R^n$ is the state vector,and $A \in R^{n\times n}$,$B \in {\bf R}^{n\times m}$,$C \in R^{p\times n}$,$D \in {\bf R}^{p\times m}$ are the state,input,output and feed-forward matrices,respectively.

It is well known that for the conventional discrete-time SIM (DSIM)^[22],the basic state space model shown in (1) and (2) can be rewritten in matrix form at the periodic instants $k = 0, 1,...,N-1$,where $N$ is the number of samples.

In (3)

$U_{0|r-1}$is the same as $Y_{0|r-1}$ except its elements are the input variable $u(k)$,$k = 0,1,...,N-1$,and $X_0=[x(0) x(1) \cdots x(N-r)] \in R^{n \times (N-r+1)}$. Both $Y_{0|r-1}$ and $U_{0|r-1}$ are the standard block Hankel matrices with the identical structure for the DSIM. $\mathcal{O}_r \in {\bf R}^{rp\times n}$ is the observability matrix and $\varPsi_r \in {\bf R}^{rp \times rm}$ is the system matrix of (1) and (2) with the structures shown in (13) and (14).

However,SIM may suffer from incorrect estimation of the model order. Theoretically,when $r$ and $N$ approach to infinity,the model order is given by^[22]

where $U_{0|r-1}^{\bot}$ is the projection of $U_{0|r-1}$,whose columns span the nullspace of $U_{0|r-1}$ . The model order is normally determined by the number of singular values of $Y_{0|r-1} U_{0|r-1}^{\bot}$ that are above the median value of the singular values in a logarithmic scale^[17]. However,this approach can be difficult in practice as the singular values may become inseparable when the I/O measurements are noise corrupted. Some examples of bad cases are presented by Fazel et al.^[12] to show this problem.

The NNM of a matrix $X$ is the heuristic convex relaxation of the rank minimization of X^[16]. When NNM is applied to the DSIM, it is usually used to constraint the rank of matrix product $Y_{0|r-1} U_{0|r-1}^{\bot}$ which determines the model order. Denote $\|X\|_*=\sum_{i=1}^{{\rm min}\{m,n\}}\sigma_i(X)$ as the nuclear norm of $X \in R^{m\times n}$ where $\sigma_i(X)$ is the $i\text{th}$ singular value of the matrix $X$.

Based on results of ^[12] ,the DT subspace identification can be cast into the following optimization problem

In (6)

where $\tilde{y}(\cdot)$ is the measured system output containing noise,$y_d \in R^{p\times N}$ is the optimization variable defined similarly as $\tilde{y}_d$ but with $\tilde{y}(\cdot)$ replaced by the computed system output obtained from the optimization process,$y_d-\tilde{y}_d$ is the output error,and $Y_{0|r-1}(y_d)$ is similar to $Y_{0|r-1}$ given in (4) with the elements from $y_d$ .

Extensive analysis and simulation results in ^[12] have shown that the minimization of $\|Y_{0|r-1}(y_d) U_{0|r-1}^{\bot}\|_*$ dramatically improves the performance of SIM,especially in the estimation of model order. However,it is not straightforward to extend the above NNM based DSIM to the CT case,because of the time-derivative of I/O data and the different structures of Hankel matrix in the CT subspace identification. The different structures of Hankel matrix make the solution of NNM in DT case cannot be directly used to solve NNM in CT case.

Consider a linear time-invariant CT state space system

where $u(t)\in R^m$ and $y(t)\in R^p$ are respectively the input and output,$x(t)\in R^n$ is the state vector,and $A\in R^{n\times n}$,$B\in R^{n\times m}$,$C\in {\bf R}^{p\times n}$,$D\in R^{p\times m}$ are the state,input, output and feed-forward matrices,respectively. Taking the time-derivative on (9) and substituting $\dot{x}(t)$ into (8),then repeating the time-derivative of $y(t)$ up to the $(r-1)$-{th} order,the state space representation (8) and (9) can be converted to the matrix form^[2]

Equation (10) is the basis for the CSIM where

where $(\cdot)^{(r-1)}:=\frac{{\rm d}^{r-1}(\cdot)}{{\rm d}t^{r-1}}$ and $U_{0|r-1}^C$ is identical to $Y_{0|r-1}^C$ in structure but with $y(t)$ replaced by $u(t)$,

and

The time sequence $t=\{t_1,t_2,...,t_N\}$ is used to represent different sampling instants,which are not necessarily equidistant. $\mathcal{O}_r^C (r>n)$ is the observability matrix and $\varPsi_r^C$ is the (block Toeplitz) system matrix of (8) and (9). Although the matrices $Y_{0|r-1}^C$ and $U_{0|r-1}^C$ do not have the structure of standard Hankel matrix,they will still be called Hankel matrix hereafter for simplification of terminology and convenience of comparison with the DSIM.

Note that (10) is in the same form as its DT counterpart (3), hence the computation procedures of the DSIM can be readily used here to obtain the CT state-space model for the system (8) and (9). In particular,when $r$ and $N$ are sufficiently large,the model order $ n={\rm rank}(Y_{0|r-1}^C U_{0|r-1}^{C\bot})$ can be determined by calculating the nonzero singular values of $(Y_{0|r-1}^C U_{0|r-1}^{C\bot})$^[20],where $U_{0|r-1}^{C\bot}$ is the projection of $U_{0|r-1}^{C}$. These computations are based on $(Y_{0|r-1}^C U_{0|r-1}^{C\bot})$ which consists of the time-derivatives of the system output and input variables up to the order $r-1$. Since the measured system output is always noise corrupted,its direct differentiations,especially the high order ones,will greatly amplify the noise and hence distort the identification results. To overcome this problem,the generalized Poisson moment functional filtering is often used in the CT subspace identification to obtain the derivatives of system input and output variables.

The $j$-th order GPMF transformation of a scalar signal over $[0,t]$ is given by the following convolution product^[23]

where $g_j(t)=\beta^{j+1} \frac{t^j}{j!} {\rm e}^{-\lambda t} $ is the $j$-th Poisson pulse function,$\lambda \in R^+$ and $\beta \in R^{+*}$ are the Poisson filter constant and gain, respectively (most of the time,$\beta=\lambda$ is set).

The Laplace transform of the Poisson pulse function with the Laplace variable $s$ is

Therefore,the $j$-{th} order GPMF of $h(t)$ can be measured as the output at time-instant $t$ of a cascaded low-pass filter chain of $(j+1)$ identical stages,each with a transfer function given by $\beta/(s+\lambda)$ .

Neglecting the initial condition of the GPMF (which is valid for GPMF filtering the input-output signals over interval $[0,T_0]$^[1]),the most important property of GPMF can be summarized as follows.

where

and $\Omega(u,v)$ represents the $(u,v)$-{th} element of $\Omega \in R^{i\times (j+1)}$ with $j+1>i$.

According to (17),the GPMF filter can be expressed as a single input state-space model. Taking $\pmb{{h}}_{j}(t)$ as the state vector,the state-space model of the $j$-th order GPMF-based filtering for the input signal $h(t)$ can be written as^[1]

where

$B_F=[0~\cdots~0~\beta]^{\rm T} \in R^{j+1}$,$C_F=\Omega \in R^{i\times (j+1)}$. In the case of multi-variable systems,the inputs and outputs are vectors $\left( h(t) \in R^q \right)$ and the state-space description of the filter is extended to

$B_{Fq}=[0_q~\cdots~0_q~\beta I_q]^{\rm T} \in R^{q(j+1)\times q}$,$C_{Fq}=\Omega_q \in R^{qi\times q(j+1)}$,where $I_q\in R^{q\times q}$ is the identity matrix,$0_q\in {\bf R}^{q\times q}$ is the null matrix and $\Omega_q(u,v)$ is the $(u,v)$-th column block matrix of $\Omega_q$.

For digital implementation of the GPMF-based filtering,the CT state space model (22) is discretized with the zero order hold or the Tustin technique,which are both available in Matlab. The influences of the discretization on the GPMF-based filtering have been analyzed and compared extensively in ^[1],and the Tustin technique is recommended in the paper.

Now reconsider the CSIM back in (10). Applying GPMF-based filtering to (10) gives

Define

then (26) can be simply written as

As (28) is in the form of (10),it also admits all the technical tools discussed in Section II-A for DSIM. Thus,the GPMF-based CSIM^[1] simply applies the computation procedure of discrete-time SIM in (28) to determine the model order and estimate the CT state-space model for the system (8) and (9).

Similar to the DSIM,the GPMF-based CSIM^[1] described in Section III-B uses the rank of ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) defined in (27) to determine the model order,under the assumption that the order of time-derivative and the number of samples contained in ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) are large enough. As a result,it shares the same problems with its DT counterpart, particularly,the problem of model order overestimation when the measured system output variables contain noise. Because of the estimated derivatives of $y(t)$ and $u(t)$ contained in ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$),the computed singular values of ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) can be even more problematic than in DT case,resulting in even unreliable estimation of the model order. Therefore,it is necessary to optimize the output variables and their derivatives in the Hankel matrix $Y_{0|r-1}^j$ to improve the robustness and accuracy of subspace identification. To this end, the following new formulation is proposed for the GPMF-based CT subspace identification.

In (29),

where $\tilde{y}(\cdot)$ is the measured CT system output filtered by the GPMF described in Section III-B, $y_c \in {\bf R}^{p\times rN}$ is the optimization variable having the same structure as $\tilde{y}_c$ but with the elements replaced by the computed system output filtered by GPMF,$Y_{0|r-1}^j$ has the same structure as that in (11) but with its elements consisting of $y_c$,$U_{0|r-1}^j$ is identical to $Y_{0|r-1}^j$ with the elements being the GPMF filtered inputs,and $\lambda$ is a non-negative weighting factor.

In the conventional GPMF-based CSIM,the estimation of the model order is mainly based on the truncated SVD of ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$)^[17] which is sensitive to noise corruption. The new formulation (29) seeks a new GPMF filtered $y$ that minimizes the nuclear norm of ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) and is close to the GPMF filtered output measurement $\tilde{y}$. This will improve the robustness in model order estimation and parameter matrix identification.

The generalized NNM problem in (29) can be solved by various first-order methods including alternating direction methods of multipliers (ADMM),proximal point algorithms and gradient projection method. The advantages and the drawbacks of these algorithms have been discussed and analyzed in ^[12]. In the sequel,the ADMM algorithm is adopted to solve the optimization problem (29).

Define

and

where $U_{0|r-1}^{j\bot}$ is known. Then the optimization problem in (29) can be rewritten as

Definition 1. Let $E$ and $F$ be two inner product spaces and $\phi$ be a linear mapping $\phi: E\rightarrow F$. Then the adjoint mapping $\phi^*: F\rightarrow E$ of the linear mapping $\phi$ is defined with the relation^[24]

Assume that $X\in R^{pr\times q}$ is a matrix with the same dimension as that of $\mathcal{H}(y_c)$. According to Definition 1,the adjoint mapping of $\mathcal{H}(\cdot)$ is given by

Define

then the adjoint mapping of Hankel matrix transformation for CT case can be shown as

Denoting $\pmb{Y}=-\mathcal{H}(y_c)$ and using $\|\mathcal{H}(y_c)\|_* = \|-\mathcal{H}(y_c)\|_*$,(34) can be converted to the standard form of minimization with equality constraint that can be solved by ADMM^[25]

Following the same line as that of ^[12],define the following augmented Lagrangian function for ADMM.

where $\beta$ is a positive number and $\Lambda \in {\bf R}^{pr\times q}$ is the Lagrangian multiplier. In ADMM,the optimization variables $\pmb{Y}$,$y_c$,$\Lambda$ are successively updated by minimizing Lagrangian function $L$ with respect to $\pmb{Y}$ then $y_c$ then $\Lambda$ individually,and the values of objective function and optimization variables are used as termination criteria^[12]. The procedures of ADMM for solving (39) are summarized below.

Step 1. Preset $(y_c^0,\Lambda^0)$,$\beta>0$, $\sigma=\frac{\beta}{\beta\eta+1}$$(\eta=\min\{r,N\})$,and $\tau\in (0,\frac{\sqrt{5}+1}{2})$.

Step 2. Take the SVD

where the columns of $U$ and $V$ are left- and right-singular vectors respectively and $\Sigma$ is a diagonal matrix consisting of the singular values,then update the parameters ${\pmb Y}$,$y_c$, and $\Lambda$

Step 3. Check the termination criteria,if all of them fail, go to Step 2.

The above procedure generally terminates when the relative duality gap is less than $10^{-4}$ or the change in Frobenius norm of successive iterations is smaller than $10^{-8}$ or the maximum number of iteration hits 2000^[12].

In this section,the novel hybrid algorithm proposed in Section IV for GPMF-based CT subspace identification with NNM is evaluated. The NNM is used as a pre-processing step to calculate an optimized filtered output sequence,which is then used to construct the Hankel matrix for the basic subspace identification algorithm.

Two typical multi-variable systems given in (47) and (48) are simulated as the true systems in evaluation,one is a 2-input and 1-output system and the other is a 3-input and 1-output system. Both systems are not easy to identify with the existing subspace identification methods when the prior information on the model order is unavailable.

The excitation signals to the true systems are chosen as pseudorandom signals with normal distribution to satisfy the persistent excitation condition^[2],which were generated by the Matlab function "randn". The output signals from the true systems are corrupted by additive noise with the signal-noise-ratio (SNR) of 20 dB,and the input and output signals are sampled with a constant sampling time 0.01 s. In order to demonstrate the advantage and effectiveness of the NNM,the Matlab function "n4sid" is used as the baseline identification algorithm for comparison. The "n4sid" is a build-in function in system identification toolbox of Matlab for the subspace identification of both DT and CT systems. For CT system identification,the "n4sid" uses indirect method,namely, identifies a DT model first and then converts it to CT model. The model order is determined by setting order as "best",which uses the same order determination criterion discussed in Section II-A. The CT identification for "n4sid" is implemented by setting the sampling time "Ts" equal to 0.

Three algorithms are applied to the simulated systems respectively to identify the models which are then compared. The three algorithms are n4sid,GPMF and NNM GPMF algorithms. The n4sid is as described above,the GPMF is the GPMF-based SIM without NNM, and NNM GPMF is the algorithm presented in Section IV with combination of GPMF-based SIM and NNM. In the GPMF and NNM GPMF algorithms,the model orders are determined by truncating the singular values of Hankel matrix with the rule given in ^[17].

In the evaluation results presented in the sequel,$n$ represents the estimated model order,$e_I$ the output error and $e_V$ the validation error,where $e_I$ and $e_V$ are defined as^[10]

In the definitions (44) and (45),$N_I$ is the number of data samples used in the identification experiment,$N_V$ is the number of data samples used for the model validation that also include the $N_I$ data samples used in identification,hence $N_V>N_I$, $\hat{y}(t)$ is the output of the identified model,$\tilde{y}(t)$ is the measured system output,and

Example 1.

Table Ⅰ compares the model order estimations,the output errors and the validation errors of the three algorithms for the example. As seen from the table,the NNM GPMF performs best in all the three performance measures.

TABLE Ⅰ

TABLE Ⅰ THE ORDER ESTIMATIONS, OUTPUT ERRORS AND VALIDATION ERRORS OF THE THREE ALGORITHMS n4sid GPMF NNM GPMF n eI eV n eI eV n eI eV Example 1 1 1.56 1.78 2 4.19 3.09 3 0.24 0.24 Example 2 2 1.48 1.23 2 17.54 16.69 3 0.16 0.17

TABLE Ⅰ THE ORDER ESTIMATIONS, OUTPUT ERRORS AND VALIDATION ERRORS OF THE THREE ALGORITHMS

Fig. 1 gives the first four singular values of the matrix ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) consisting of the GPMF filtered outputs of the estimated models,where the singular values are normalized by the maximum singular value. As seen clearly from the figure,for Example 1,the singular values of ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) from the GPMF method are very close to each other,which makes the determination of model order difficult. In contrast,the singular values from NNM GPMF algorithm show a significant difference at the fourth singular value,which makes the determination of order easier.

Fig. 1

	Download: JPG larger image
Fig. 1. Singular values of the Hankel matrices from GPMF and NNM GPMF algorithms for Example 1.

Fig. 2 compares the magnitudes of frequency response of the true system and the estimated models with algorithms n4sid and NNM GPMF,the general GPMF SIM is ignored due to the non-convergent result. As shown in Fig. 2,NNM GPMF algorithm matches best with the true system no matter the frequency is low or high. Fig. 3 compares the noise free output of the true system (47) and the outputs of the models identified by n4sid and NNM GPMF. As seen from the plots,NNM GPMF algorithm outperforms the n4sid and GPMF algorithms significantly in the precision of system identification.

Fig. 2

	Download: JPG larger image
Fig. 2. The magnitude of frequency responses of Example 1.

Fig. 3

	Download: JPG larger image
Fig. 3. The system outputs of Example 1.

Example 2.

This example is of the same order as Example 1 but with 3-input and 1-output for further assessment of the performance of NNM GPMF. The results of the model order estimations,output errors and validation errors of the three algorithms for this example are shown in Table Ⅰ. As seen from the table,the baseline algorithms n4sid and GPMF give an underestimated model order,2,and larger identification and validation errors than those from the proposed NNM GPMF algorithm. The first four singular values of the matrix ($Y_{0|r-1}^j U_{0|r-1}^{j\bot}$) are given in Fig. 4. The magnitudes of frequency response and outputs of the true system and the models of n4sid and NNM GPMF algorithms are respectively given in Figs. 5 and 6,and the result of GPMF SIM is not convergent as well. Again,the proposed NNM GPMF algorithm performs best.

Fig. 4

	Download: JPG larger image
Fig. 4. Singular values of the Hankel matrices from GPMF and NNM GPMF algorithms for Example 2.

Fig. 5

	Download: JPG larger image
Fig. 5. The magnitudes of frequency responses of Example 2.

Fig. 6

	Download: JPG larger image
Fig. 6. The system outputs of Example 2.

A novel method has been presented for the subspace identification of CT systems. This method incorporates the recently developed NNM technique into the conventional GPMF-based CSIM. It takes the advantage of GPMF algorithm to avoid the time-derivatives of input and output measurements and exploits the power of NNM to minimize the rank of Hankel matrix consisting of the GPMF filtered input and output measurements. An ADMM based computation algorithm has been derived to implement this novel method,and two numerical examples have been presented to verify the validity and evaluate the performance of the proposed method. The examples have demonstrated that the proposed method is superior to the existing methods for CT subspace system identification.