2.6: Introduction of Gray- and Black-Box Modeling Using System Identification and Machine Learning

Learning Objectives

After completing this section, you will be able to:

1. Define gray- and black-box modeling approaches.
2. Describe the machine learning modeling method for developing a gray-box model.
3. Use system identification to develop a block-box model.

Sections 2.2 to 2.4 introduced physical based (or physics-based) modeling techniques. There are two additional data-based modeling approaches, called gray- and black-box methods, that construct the system model based on only system input and output data. Because most of the data is collected digitally with a fixed sample rate, data-based modeling is often in the discrete-time domain described by difference equations. For instance, the following general form \ref{2.163} with \(n\ge m>0\) can be used to model a single-input and single-output (SISO) system that is also called a regression model.

\[y\left(k\right)+a_1y\left(k-1\right)+\dots +a_ny\left(k-n\right)=b_0u\left(k\right)+b_1u\left(k-1\right)+\dots +b_mu(k-m) \label{2.163}\]

Note that a linear continuous-time system described by a continuous-time ODE can always be discretized using sample-and-hold devices such as A/D (analog-to-digital) and D/A (digital-to analog) converters (see Chapter 4). For a simple system, a time-derivative in an ODE can be replaced by a difference function as follows when the sample period \(\Delta T\) is small enough.

\[\dfrac{df(t)}{dt}\approx \dfrac{f\left(k\Delta T\right)-f(\left(k-1\right)\Delta T)}{\Delta T} \label{2.164}\]

In this case, for a small enough sample period \(ΔT\), the time derivatives can be replaced by the above difference formula to convert a continuous-time ODE into a discrete-time difference equation. Assuming that the model structure is known (that is, \(n\) and \(m\) are known), there are multiple methods in this case that can be used to find the difference equation’s coefficients (\(a_1,a_2,\dots ,\ a_n\) and \(b_0,b_1,\dots ,b_m\)). One commonly used method is called the Least-Squares (LS) optimization method. We will provide a simple example to demonstrate it below in Example 2.18.

Example 2.18

Let \(n=2\) and \(m=1\). Equation \ref{2.164} can then be reorganized into the following.

\[y\left(k\right)=b_0u\left(k\right)+b_1u\left(k-1\right)-a_1y\left(k-1\right)-a_2y\left(k-2\right) \label{2.165}\]

Using information at steps \(k,k+1,\ \dots ,\ k+l\) leads to \(l+1\) equations. We can organize these equations into the following matrix equation

\[\underbrace{\left[ \begin{array}{c} \begin{array}{c} y(k) \\ y(k+1) \end{array} \\ \begin{array}{c} \vdots \\ y(k+l) \end{array} \end{array} \right]}_{Y(k)}=\underbrace{\left[ \begin{array}{cc} \begin{array}{cc} u(k) & u(k-1) \\ u(k+1) & u(k) \end{array} & \begin{array}{cc} -y\left(k-1\right) & \ -y(k-2) \\ -y\left(k\right) \ & \ -y(k-1) \end{array} \\ \begin{array}{cc} \vdots & \vdots \\ u(k+l) & u(k+l-1) \end{array} & \begin{array}{cc} \vdots & \vdots \\ -y(k+l-1) & -y(k+l-2) \end{array} \end{array} \right]}_{D(k)}\underbrace{\left[ \begin{array}{c} \begin{array}{c} b_0 \\ b_1 \end{array} \\ \begin{array}{c} -a_1 \\ -a_2 \end{array} \end{array} \right]}_{B} \label{2.166}\]

where vector \(Y(k)\) and matrix \(D(k)\) can be formed using the collected experimental control input \(u(k+j)\) (\(j=0,\ 1,\ \dots ,\ l\)) and correspondent system output \(y(k+j)\) (\(j=0,\ 1,\ \dots ,\ l\)). Normally, parameter \(l\) is selected to be much greater than \(n+m\) or, in other words, \(l\gg m+n\). With this assumption, matrix \(D(k)\) is a tall matrix with the number of rows much higher than the number of columns. Due to measurement noises in experimental data and system nonlinearity, linear equation \ref{2.166} will generally not be satisfied. As a result, there is no feasible solution for \(B\). Only an optimal approximated solution can be obtained that minimizes the error defined below

\[{\left\|\left[ \begin{array}{c} \begin{array}{c} y(k) \\ y(k+1) \end{array} \\ \begin{array}{c} \vdots \\ y(k+l) \end{array} \end{array} \right]-\left[ \begin{array}{cc} \begin{array}{cc} u(k) & u(k-1) \\ u(k+1) & u(k) \end{array} & \begin{array}{cc} -y\left(k-1\right) & \ -y(k-2) \\ -y\left(k\right) \ & \ -y(k-1) \end{array} \\ \begin{array}{cc} \vdots & \vdots \\ u(k+l) & u(k+l-1) \end{array} & \begin{array}{cc} \vdots & \vdots \\ -y(k+l-1) & -y(k+l-2) \end{array} \end{array} \right]\left[ \begin{array}{c} \begin{array}{c} b_0 \\ b_1 \end{array} \\ \begin{array}{c} -a_1 \\ -a_2 \end{array} \end{array} \right]\right\|}_2 \label{2.167}\]

where \({\left\|\cdot \right\|}_2\) is a 2-norm in the Euclidian space. That is, if \(x={\left[x_1,\ x_2,\dots ,x_n\right]}^T\), then

\[{\left\|x\right\|}_2=\sqrt{x^2_1+x^2_2+\dots +x^2_n} \label{2.168}\]

With the assumption that \(l\gg m+n\) and matrix \(D(k)\) is full column rank, the optimal solution minimizing the 2-norm in equation \ref{2.167}, called the Least-Squares (LS) solution, is shown below

\[B={\left(D^T\left(k\right)D\left(k\right)\right)}^{-1}D^T\left(k\right)Y(k) \label{2.169}\]

where the obtained parameter vector \(B\) represents the identified model defined in equation \ref{2.165}.

Remark 2.1

In order to have a high-quality identified model, the control input \(u(k)\) should contain rich information so that the dynamic system is fully excited. One indication of a good excitation is that the condition number of matrix \(D^T\left(k\right)D\left(k\right)\) is small; in this case, the condition number is defined as the maximum eigenvalue of \(D^T\left(k\right)D\left(k\right)\) divided by the minimum one. Note that matrix \(D^T\left(k\right)D\left(k\right)\) is positive definite, assuming that it is invertible. Here, the condition number of matrix \(D^T\left(k\right)D\left(k\right)\) can be used to justify the quality of the identified gray-box system model.

Remark 2.2

The formulation in equation \ref{2.166} can also be applied to nonlinear difference equations for gray-box system identification. That is, the method described in Example 2.18 is not limited to linear system identification.

Example 2.19 uses the LS method to identify a nonlinear system model defined by a nonlinear difference equation.

Example 2.19

Consider a nonlinear system described by the following nonlinear difference equation with its \(n=2\) and \(m=1\).

\[y\left(k\right)+a_1y^3\left(k-1\right)+a_2y\left(k-2\right)=b_0u\left(k\right)+b_1u^2\left(k-1\right) \label{2.170}\]

Letting \(l\gg m+n\) and using a similar method to that in Example 2.18, we can reorganize equation \ref{2.170} into a matrix equation as follows with nonlinear terms in \(D(k)\).

\[\underbrace{\left[ \begin{array}{c} \begin{array}{c} y(k) \\ y(k+1) \end{array} \\ \begin{array}{c} \vdots \\ y(k+l) \end{array} \end{array} \right]}_{Y(k)}=\underbrace{\left[ \begin{array}{cc} \begin{array}{cc} u(k) & u^2(k-1) \\ u(k+1) & u^2(k) \end{array} & \begin{array}{cc} -y^3\left(k-1\right) & \ -y(k-2) \\ -y^3\left(k\right) \ & \ -y(k-1) \end{array} \\ \begin{array}{cc} \vdots & \vdots \\ u(k+l) & u^2(k+l-1) \end{array} & \begin{array}{cc} \vdots & \vdots \\ -y^3(k+l-1) & -y(k+l-2) \end{array} \end{array} \right]}_{D(k)}\underbrace{\left[ \begin{array}{c} \begin{array}{c} b_0 \\ b_1 \end{array} \\ \begin{array}{c} -a_1 \\ -a_2 \end{array} \end{array} \right]}_{B} \label{2.171}\]

With the assumption that \(l\gg m+n\) and matrix \(D(k)\) is full column rank, we can obtain the LS optimal solution minimizing the 2-norm of \({\left\|Y\left(k\right)-D\left(k\right)B\right\|}_2\). In this case, the entries in columns 2 and 3 of matrix \(D(k)\) are nonlinear.

There is a real-time system identification scheme using the iterative LS method with a forgetting factor so that the identified model will adapt to current physical system status [1] and [9]. This approach can also be used for predicting system future response (see [19] and [20] for its application to vehicle control systems).

Machine learning is now widely used for optimization based on the model obtained from the learning process. To obtain a high-quality model, the excitation needs to be designed to fully excite the physical system. The machine learning approach can fill this gap by automatically selecting the new test points to improve the current model accuracy [21] and [22]. In this section, we will discuss one type of surrogate model—the Kriging model—that is widely used in machine learning applications.

A surrogate model (static gray-box model) is cheap to evaluate because it does not include the physical dynamics that maps from inputs to outputs. Reference [23] reviews how to train a gray-box surrogate model, based on a collection of input-output data, from either simulating a detailed physical model or conducting experiments on a physical system. Among these methods, the Kriging (Gaussian process) model is now widely used. It predicts the mean value; it also predicts variance that includes both the intrinsic uncertainty inherent in a stochastic process and the extrinsic uncertainty of the unknown response surface [24].

The standard Kriging model is deterministic as it includes only model uncertainty, also called extrinsic uncertainty; it does not consider system intrinsic noises. Much literature has discussed Kriging modeling in detail [25] with a few modified for stochastic processes [24], [25], [26]. In a stochastic Kriging model, by allowing the Kriging model to regress the data instead of interpolating it, the original deterministic approach was extended to a broad one. This was accomplished by adding a regression constant, also called a regularization constant, to the leading diagonal of the Kriging correlation matrix, which will be defined in equation \ref{2.172}. This section briefly introduces both deterministic and stochastic Kriging models. They can be modeled based on the extension of the Design and Analysis of Computer Experiments (DACE) MATLAB toolbox by Nielsen [27] and the extended SK code package developed by Nelson’s group [24]. The reasons for using both deterministic and stochastic surrogate models will be explained in the next paragraph.

As seen in equation \ref{2.172}, a scalar stochastic Kriging model essentially includes three terms. The first term, \({\boldsymbol{f}}^T\left(\boldsymbol{x}\right)β\), is a function of input vector \(\boldsymbol{x}\) to capture the general trend. The second, \(z\left(\boldsymbol{x}\right)\), is a realization of a zero mean random field so that the response of unknown design sites is modeled within a statistical framework. The third, \(ε\left(\boldsymbol{x}\right)\), is introduced for system and measurement noise consideration to model the nature of stochastic processes. The probabilistic format of a scalar response surface model is

\[y\left(\boldsymbol{x}\right)={\boldsymbol{f}}^T\left(\boldsymbol{x}\right){\boldsymbol{β}}+z\left(\boldsymbol{x}\right)+ε\left(\boldsymbol{x}\right) \label{2.172}\]

where \(\boldsymbol{x}\) is a \(p\) dimensional vector of design variables; \(\boldsymbol{f}\left(\boldsymbol{x}\right)\) is a \(q\) dimensional vector of known (given) functions; and \(β\) is a constant parameter vector to be determined with compatible dimension. In addition, \(z\left(\boldsymbol{x}\right)\) is a statistical model of random process with zero mean and covariance based on spatial correlation functions between two design points \({\boldsymbol{x}}^i\)and \({\boldsymbol{x}}^j\), where \(i,\ j\in \left[1,\ n\right]\) and \(n\) is the total number of evaluation points Furthermore, \(ε\left(\boldsymbol{x}\right)\) is an intrinsic noise that is assumed to be Gaussian with its variance \(λ^2\left(\boldsymbol{x}\right)=Ver(ε\left(\boldsymbol{x}\right))\) called the regularization constant.

Without noise \(ε\left(\boldsymbol{x}\right)\), the expression for elements of \(n\) sample covariance matrix \(R_n\) can be expressed in equation \ref{2.173}

\[\left[R_n\left(i,j\right)\right]=\left[Cov\left(z\left({\boldsymbol{x}}^i\right),z\left({\boldsymbol{x}}^{\boldsymbol{j}}\right)\right)\right]\boldsymbol{=}σ^2R_d\left(θ,{\boldsymbol{x}}^i,{\boldsymbol{x}}^{\boldsymbol{j}}\right) \label{2.173}\]

where \(σ^2\) is the variance of the process to be determined, and \(R_d\) is an \(n\times n\) correlation function matrix based on the spatial distance for all pairs \({\boldsymbol{x}}^i\) and \({\boldsymbol{x}}^{\boldsymbol{j}}\) (\(i,j\in \left[1,n\right]\)). For the Kriging model, \(R_d\) is a Gaussian correlation function and can be calculated using equation \ref{2.174} based on the spatial distance. Note that the correlation between points decreases as their distance increases.

\[\left[R_d\left(i,j\right)\right] =\left[Cov\left(z\left({\boldsymbol{x}}^i\right),z\left({\boldsymbol{x}}^{\boldsymbol{j}}\right)\right)\right] = \left[\mathrm{e}\mathrm{xp}\left(\boldsymbol{-}\sum^p_{l=1}{θ_l{\left|x^i_l-x^j_l\right|}^2}\right)\right] \label{2.174}\]

where \(θ_l\ \)(\(1\le l\le p\)) is the \(l^{th}\) element of parameter vector \(θ\boldsymbol{=}{\left[θ_1,\ θ_2,\cdots ,θ_p\right]}^T\) associated with the width of the Gaussian function and \(x^i_l\) (\(1\le l\le p\)) is the \(l^{th}\) element of design variable vector \({\boldsymbol{x}}^i\boldsymbol{=}{\left[x^i_1,\ x^i_2,\cdots ,x^i_p\right]}^T\).

The effect of parameter vector \(θ\) can be thought of as determining how quickly the function changes as \({\boldsymbol{x}}^{\boldsymbol{j}}\) moves away from \({\boldsymbol{x}}^{\boldsymbol{i}}\). A large \(θ_l\) implies a significant activity, or low correlation, in dimension \(l\). From this, with the known input-output data set \(D\) formed by \(n\) input points, an \(n\times n\) covariance matrix \(R_n=σ^2R_d\) can be constructed.

For the stochastic case with independent noises, the corresponding covariance matrix \(R_s\) is obtained by adding regulation constant \(λ^2\) to the diagonal components of \(σ^2R_d\) such that \(R_s=σ^2R_d+λ^2I\). With all these, the model has three unknowns to be solved for the Kriging model. They are scalar \(σ\), and parameter vectors \(β\) and \(θ\), determined by maximizing the likelihood function in equation \ref{2.175} over the observed points. To simplify the likelihood maximization, the natural logarithm will normally be taken to have the following,

\[{\mathrm{l}\mathrm{n} L\left(σ,β,\mathrm{\ }θ\right)\ }=-\dfrac{1}{2}\left[n{\mathrm{l}\mathrm{n} \left(2πσ^2\right)+{\mathrm{l}\mathrm{n} \left({\mathrm{d}\mathrm{et} \left(R_s\right)\ }\right)\ }\ }\right]\boldsymbol{+}\dfrac{1}{σ^2}{\left(Y-Fβ\right)}^TR^{-1}_s\left(Y-Fβ\right) \label{2.175}\]

where \(R_s\) is a covariance matrix for these observed points, \(Y={\left[y\left({\boldsymbol{x}}^1\right),y\left({\boldsymbol{x}}^2\right),\dots ,y({\boldsymbol{x}}^N)\right]}^T\). In addition, \(F={\left[{\boldsymbol{f}}^T\left({\boldsymbol{x}}^1\right),{\boldsymbol{f}}^T\left({\boldsymbol{x}}^2\right),\dots ,{\boldsymbol{f}}^T({\boldsymbol{x}}^N)\right]}^T\) is a vector formed by \(n\) observed function responses and a function matrix evaluated at \(n\) design vectors.

After the likelihood function is maximized, the stochastic Kriging model can be obtained. Accordingly, the predicted mean \(μ_y\) and mean square error \(s^2_y\) for \(y(\boldsymbol{x})\) can be calculated as

\[μ_y\left(\boldsymbol{x}\right)={\boldsymbol{f}}^T\left(\boldsymbol{x}\right)β+{\boldsymbol{r}}^T\left(\boldsymbol{x}\right)R^{-1}_s\left(Y-Fβ\right) \\ s^2_y\left(\boldsymbol{x}\right)=σ^2+{\boldsymbol{u}}^T\left(\boldsymbol{x}\right){\left({\boldsymbol{F}}^TR^{-1}_s\boldsymbol{F}\right)}^{-1}\boldsymbol{u}\left(\boldsymbol{x}\right)-{\boldsymbol{r}}^T\left(\boldsymbol{x}\right)R^{-1}_s\boldsymbol{r}\left(\boldsymbol{x}\right) \label{2.176}\]

where \({\boldsymbol{r}}^T\left(\boldsymbol{x}\right)=σ^2\left[R_d\left(\boldsymbol{x},{\boldsymbol{x}}^1\right),R_d\left(\boldsymbol{x},{\boldsymbol{x}}^2\right),\dots ,R_d\left(\boldsymbol{x},{\boldsymbol{x}}^n\right)\right]\) and \(\boldsymbol{u}\left(\boldsymbol{x}\right)=\boldsymbol{f}\left(\boldsymbol{x}\right)-{\boldsymbol{F}}^TR^{-1}_sr(\boldsymbol{x})\).

Note that the exact stochastic property of \(z\left(\boldsymbol{x}\right)\) is not modeled in this process. Instead, the stochastic property of \(y\left(\boldsymbol{x}\right)\), required for optimization, is modeled. The general structure of deterministic and stochastic Kriging models is summarized in Figure 2.44 (solid line path). For the deterministic case, variance \(λ^2\) is set to zero to obtain the predicted mean and mean square errors. In this case, each point provides an exact correlation with itself, forcing the predictor to interpolate evaluated points. Detailed derivation is presented in reference [25].

Existing research has shown that the regression function used in the Kriging model has less impact on the model performance [28]. However, the zero-order function in the ordinary Kriging model has proven to be erroneous for certain nonlinear applications [29]. The universal Kriging model adopts polynomial functions, mostly linear or quadratic, as the regression part to capture the simple spatial trend [29], [30]. A polynomial function with order higher than two should be avoided because complicated data trends are better described by the stochastic model component [30] (Gaussian process). With all these details, a surrogate Kriging model can be developed.

In order to solve the optimization problem efficiently, a surrogate-model-assisted optimization (SMAO) algorithm can be adopted to improve the model accuracy in the region near the optimum using proposed new test observation data points; see the dotted-line path in Figure 2.44. For detailed application information, see references [21] and [22].

For the black-box model, there is no assumption on the model structure. Usually, the detailed model information such as model order is determined during the system identification process based purely on system input and output data. As discussed earlier; in order to have a high-quality system identified model, the input excitation needs to contain signals over a wide bandwidth so that all system modes can be excited. If a system mode were not excited, the associated mode would not be contained in the identified model. Therefore, white noise, which has a constant frequency spectrum over the entire frequency range, is ideal for system identification (black-box modeling). However, every practical system has its physical limitations such as limited actuator magnitude. In practice, ideal white noise may not be feasible to be used because saturations on the input signals (actuators) have to be applied to avoid violating actuator constraints.

An impulse signal, such as the impulse function \(δ\)(t) introduced earlier in Section 2.1 about the Laplace transform, is also a choice for system excitation. However, its large magnitude makes it infeasible for accurate implementation. One of the input signals widely used in practice for system identification is PRBS (pseudo-random-binary-signal); its absolute magnitude remaining constant makes it feasible to be accurately implemented in practice. In the remainder of this section, we will introduce a system identification method using PRBS as the excitation to obtain a black-box model.

The most commonly used PRBS signals are based on maximum length sequences (called \(m\)-sequences) [31] for which the length of the PRBS signals is \(m=2^n-1\), where \(n\) is an integer called the order of PRBS. Let \(z^{-1}\) represent a delay operator and define \(\hat{p}\left(z^{-1}\right)\) and \(p\left(z^{-1}\right)\) to be the polynomials shown below.

\[p\left(z^{-1}\right)=a_nz^{-n+1}\bigoplus \ \cdots \bigoplus a_2z^{-1}\bigoplus a_1=\hat{p}\left(z^{-1}\right)z^{-1}\bigoplus 1 \label{2.177}\]

where \(a_i\) is either \(1\) or \(0\), and \(\bigoplus\) obeys binary addition, i.e.,

\[1\bigoplus 1=0=0\bigoplus 0 \ ~\ \& \ 1\bigoplus 0=1=0\bigoplus 1 \label{2.178}\]

and the non-zero coefficients \(a_i\) of the polynomial are defined in Table 2.5 and in [31].

Table 2.5: Nonzero coefficients of PRBS polynomial.

Polynomialorder (\(n\))	Period ofsequence (\(m\))	Non-zeroCoefficient
2	63	\(a_5,a_6\)
3	127	\(a_4,a_7\)
4	255	\(a_2,a_3,a_4,a_8\)
5	511	\(a_5,a_9\)
6	1023	\(a_7,a_{10}\)
7	2047	\(a_9,a_{11}\)

The PRBS can then be generated by the following formula

\[\hat{u}(k+1)=\hat{p}\left(z^{-1}\right)\hat{u}(k) \label{2.179}\]

where \(\hat{u}\left(0\right)=1\) and \(\hat{u}\left(-1\right)=\hat{u}\left(-2\right)=\cdots =\hat{u}\left(-n\right)=0\). Define the following sequence.

\[s(k)=\left\{ \begin{array}{cc} a; & \mathrm{i}\mathrm{f}\ k\ \mathrm{i}\mathrm{s\ even} \\ -a; & \mathrm{i}\mathrm{f}\ k\mathrm{\ }\mathrm{is\ odd} \end{array} \right.\label{2.180}\]

The signal shown in the following equation

\[u(k)=s(k)\bigotimes \left[-a+2a\hat{u}(k)\right] \label{2.181}\]

is called the inverse PRBS, where \(\bigotimes\) obeys the following rules.

\[a\bigotimes a=-a=\left(-a\right)\bigotimes \left(-a\right), a\bigotimes (-a)=a=(-a)\bigotimes a \label{2.182}\]

It is clear after some analysis that \(u(k)\) has a period of \(2m\) and that \(u\left(k\right)=-u(k+m)\). The mean of the inverse PRBS is shown below.

\[m_u=E_{2m}\left[u(k)\right]\triangleq \dfrac{1}{2m}\sum^{2m-1}_{k=0}{u\left(k\right)}=0 \label{2.183}\]

The autocorrelation of \(u(k)\), which is \(R_{uu}\left(τ\right)=E_{2m}\left[u(k+τ)u(k)\right],\) is shown below.

\[R_{uu}\left(τ\right)=\dfrac{1}{2m}\sum^{2m-1}_{k=0}{u\left(k+τ\right)u\left(k\right)}=\left\{ \begin{array}{c} \begin{array}{cc} a^2, & τ=0 \ \ -a^2, & τ=m \end{array} \\ \begin{array}{cc} -a^2/m, & τ\ \mathrm{ev}\mathrm{en} \\ a^2/m, & τ\ \mathrm{o}\mathrm{dd} \end{array} \end{array} \right. \label{2.184}\]

Note that the first and second order information of the inverse PRBS signal are very close to those of the white noise for a large enough \(m\). The inverse PRBS is used in the \(q\)-Markov Cover (\(q\)-Markov COVariance Equivalent Realization) identification algorithm, discussed in the next [10]. For convenience, in the rest of this section, the term “PRBS” is used to represent the inverse PRBS.

A discrete-time asymptotically stable system can be represented in a standard discrete-time state-space model format.

\[\begin{array}{ccc} x(k+1) & = & Ax\left(k\right)+Dw(k) \ y(k) & = & Cx\left(k\right)+Hw(k) \end{array} \label{2.185}\]

where \(x\in {\mathbb{R}}^{n_x}\), \(w\in {\mathbb{R}}^{n_w}\), and \(y\in {\mathbb{R}}^{n_y}\) are state, input, and output vectors, respectively. Let \(X\) be the positive definite solution of the following Lyapunov equation

\[X=AXA^T+DWD^T \label{2.186}\]

where \(W\) is the covariance matrix of white noise input \(w(k)\). The system Markov parameter matrices (\(H_i,\ i=0,\ 1,\ \dots ,\ q-1\)) and the associated covariance matrices (\(R_i,\ i=0,\ 1,\ \dots ,\ q-1\)) are defined below.

\[\begin{array}{c} H_0=H, \ H_i=CA^{i-1}D, \ i=1,\ 2,\ \dots ,\ q-1 \\ R_i=CA^iXC^T+H_iWH^T_i, \ i=0,\ 1,\ \dots ,\ q-1 \end{array}\label{2.187}\]

Note that these parameter matrices are dependent on system matrices \(A\), \(D\), \(C\), and \(H.\) They can be calculated based on the input-output data and its autocorrelations. Thus inversely, the results will be used further to generate the system state-space model.

Because PRBS signal \(u(k)\) is a periodic signal with period \(2m\), its system steady-state response is also periodic with the same period as defined in equations \ref{2.188} and \ref{2.189}.

\[x^j_s\left(k+i\right)=A^{k+i}x^j_s\left(0\right)+\sum^{k+i-1}_{n=0}{A^{k+i-1-n}De_j}u(n) \label{2.188}\]

where \(e_j={\left[0,\ \dots ,\ 1,\ 0,\ \dots ,0\right]}^T\)is a unit base vector with its \(j^{th}\) element equal to \(1\) and the rest of elements equal to zero. In addition, \(x^j_s\left(k+j\right)\) is the steady-state periodic response at the beginning of the period. The associated system output responses are shown below.

\[y^j_s\left(k+i\right)=Cx^j_s\left(k+i\right)+HDe_ju(k+i) \label{2.189}\]

for \(i=0.\ 1,\ \dots ,\ q-1\), where \(e_j\) selects the \(j^{th}\) input channel of \(w(k)\) to be excited by PRBS. Now with all equations from \ref{2.185} to \ref{2.189}, we have the following.

\[Q_qy^j_q\left(k\right)=Q_qO_qx^j_s\left(k\right)+Q_qH_qw^j_q(k) \label{2.190}\]

where \(Q_q=block\ diag\left[Q,\ Q,\dots ,\ Q\right]\) is a block diagonal matrix with \(q\) blocks of positive definite output weighting matrix \(Q\); vectors \(y^j_q\left(k\right)\) and \(w^j_q\left(k\right)\) are defined in equation \ref{2.191}; and matrices \(O_q\) and \(H_q\) are defined in equation \ref{2.192}.

\[y^j_q\left(k\right)=\left[ \begin{array}{c} \begin{array}{c} y^j_s\left(k\right) \\ y^j_s\left(k+1\right) \end{array} \\ \begin{array}{c} \vdots \\ y^j_s\left(k+q-1\right) \end{array} \end{array} \right], w^j_q\left(k\right)=\left[ \begin{array}{c} \begin{array}{c} e_ju\left(k\right) \\ e_ju\left(k+1\right) \end{array} \\ \begin{array}{c} \vdots \\ e_ju\left(k+q-1\right) \end{array} \end{array} \right] \label{2.191}\]

\[O_q=\left[ \begin{array}{c} \begin{array}{c} C \\ CA \end{array} \\ \begin{array}{c} \vdots \\ CA^{q-1} \end{array} \end{array} \right], H_q=\left[ \begin{array}{cc} \begin{array}{cc} H_0 & 0 \\ H_1 & H_0 \end{array} & \begin{array}{cc} \cdots & 0 \\ \cdots & 0 \end{array} \\ \begin{array}{cc} \vdots & \vdots \\ H_{q-1} & H_{q-2} \end{array} & \begin{array}{cc} \ddots & \vdots \\ \cdots & H_0 \end{array} \end{array} \right] \label{2.192}\]

Taking the expectation over one PRBS period from both sides of equation \ref{2.190} yields the following equations from \ref{2.193} to \ref{2.198} for any \(q>0\).

\[Q_qR_qQ^T_q=Q_qO_qXO^T_qQ^T_q+Q_qH_qWH^T_qQ^T_q+Q_q{\Delta }_qQ^T_q \label{2.193}\]

where \({\Delta }_q\ \)is the error term caused by finite period number \(2m\) of the PRBS and based upon Theorem 3.1 in reference [10]. In addition, the Frobenius norm of \({\Delta }_q\) goes to zero as \(m\) goes to infinity. The rest of the matrices are defined below.

\[R_q=\left[ \begin{array}{cc} \begin{array}{cc} R_{yy}(0) & R_{yy}\eqref{GrindEQ__1_} \\ R_{yy}\eqref{GrindEQ__1_} & R_{yy}(0) \end{array} & \begin{array}{cc} \cdots & R_{yy}(q-1) \\ \cdots & R_{yy}(q-2) \end{array} \\ \begin{array}{cc} \vdots & \vdots \\ R_{yy}(q-1) & R_{yy}(q-2) \end{array} & \begin{array}{cc} \ddots \ & \vdots \\ \cdots \ & R_{yy}(0) \end{array} \end{array} \right] \label{2.194}\]

\[W_q=block\ diag\left[W,\ W,\dots ,W\right] \label{2.195}\]

where

\[R_{yy}(i)=\sum^{n_w}_{j=1}{\left(\dfrac{1}{2m}\sum^{2m-1}_{k=0}{y^j_s\left(k+i\right){\left[y^j_s\left(k\right)\right]}^T}\right)} \label{2.196}\]

\[X=\sum^{n_w}_{j=1}{\left(\dfrac{1}{2m}\sum^{2m-1}_{k=0}{x^j_s\left(k\right){\left[x^j_s\left(k\right)\right]}^T}\right)} \label{2.197}\]

\[W=\sum^{n_w}_{j=1}{\left[\dfrac{1}{2m}\sum^{2m-1}_{k=0}{e_ju^2(k)e^T_j}\right]} \label{2.198}\]

For a sufficiently large \(m\), we have the following.

\[Q_qR_qQ^T_q=Q_qO_qXO^T_qQ^T_q+Q_qH_qWH^T_qQ^T_q \label{2.199}\]

This equation allows the use of PRBS to develop a \(q\)-Markov Cover for system identification of these system matrices defined in equation \ref{2.185}. As shown in Figure 2.45, we use the five steps of system identification to complete the entire identification process. The detailed \(q\)-Markov Cover system identification algorithm is provided in the next subsection.

---