<figure><img src="image_1.png" /></figure>

HAL
open science

# Subspace System Identification with Unknown Disturbance Rejection

Qinghua Zhang

► To cite this version:

Qinghua Zhang. Subspace System Identification with Unknown Disturbance Rejection. IOMAC 2025 - 11th International Operational Modal Analysis Conference, Institut National de Recherche en Informatique et en Automatique [Inria], May 2025, Rennes, France. pp.1-8. hal-05232754

HAL Id: hal-05232754

https://hal.science/hal-05232754v1

Submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Subspace System Identification with Unknown Disturbance Rejection

Qinghua Zhang¹

¹ Inria, Université Gustave Eiffel, 35042 Rennes Cedex, France. qinghua.zhang@inria.fr

## ABSTRACT

System identification usually assumes that the considered system is driven by known inputs and/or stationary random noises. This paper considers the case involving unknown inputs, which are arbitrary disturbances. A typical example is a mechanical structure naturally excited by wind, which changes direction from time to time. To address such disturbances, subspace methods for system identification will incorporate techniques of unknown input observers, which are state estimators with the ability to reject arbitrary unknown disturbances, provided a mathematical model of the system under consideration is available. Simulation results are reported to illustrate the proposed method for system identification while rejecting unknown disturbances.

Keywords: system identification, subspace methods; unknown inputs; disturbance rejection.

## 1. INTRODUCTION

Subspace system identification is a well developed approach to building linear state-space models from sensor data [1, 2]. In general, the sensor data correspond to the input and the output of the underlying system. In some applications, only the output data are available, not the input. This is the case with operational modal analysis, for which the input corresponds to natural excitations. If the unavailable input is stationary white noises, then some subspace system identification methods, known as Stochastic Subspace Identification (SSI), is applicable [3, 4]. The purpose of this paper is to propose a system identification method to address the case involving unknown inputs, which are arbitrary disturbances. The only assumption on the unknown inputs is that they are restricted to some unknown subspace of the state-space.

In order to avoid confusion between SSI, which processes only output data, and more general subspace system identification handling both input and output data, in this paper the latter will be referred to as
N4SID, which is a short name of specific numerical algorithms for subspace state-space system identification [5], though other similar algorithms could also be used for the purpose of this paper.

To address unknown inputs, the design of the new system identification method is inspired by some state estimation methods, which have the ability to reject arbitrary unknown inputs [6, 7], also referred to as disturbances in this paper. These methods, known as Unknown Input Observers (UIO), require that a mathematical model of the system under consideration is available.

In this paper, the main ideas of subspace system identification and UIO are associated in some way, resulting in an identification procedure rejecting unknown inputs. The ability of UIOs to reject unknown inputs is based on output injection. Usually, state estimation and disturbance rejection are jointly considered in UIO design methods, leading to relatively sophisticated design procedures [6, 8]. Some recent advances simplify UIO design by separating it into two steps: disturbances are first rejected by means of output injection, then a standard state estimation method, typically the Kalman filter is applied, as if there was no unknown disturbance [9]. Based on this simplification, in this paper, disturbance rejection is incorporated into system identification in order to estimate state-space models despite unknown inputs.

The association of system identification and UIO has to overcome some difficulties. First, when disturbance rejection is made with output injection, the injected outputs act as inputs. When a model is built from data by applying a system identification method, the fact that the outputs are injected as inputs may lead to a trivial model. Second, when disturbances are rejected, the information in the state equations corresponding to the subspace of the disturbances is canceled out. As a consequence, part of the state equations cannot be found from data by system identification. These difficulties will be explained in detail in the next sections, together with the proposed solutions.

## 2. STATE-SPACE FORMULATION AND SUBSPACE SYSTEM IDENTIFICATION

This paper considers state-space systems in the form of

$$x_{k+1} = Ax_k + w_k \tag{1a}$$

$$y_k = Cx_k + v_k \tag{1b}$$

where $x_k \in \mathbb{R}^n$ is the state vector, $y_k \in \mathbb{R}^m$ is the vector of sampled sensor measurements, $w_k \in \mathbb{R}^n$ is due to external excitations and modeling uncertainties, $v_k \in \mathbb{R}^m$ is the measurement noises, and $A \in \mathbb{R}^{n \times n}$, $C \in \mathbb{R}^{m \times n}$ are matrices characterizing the state-space model [10]. Any linear transformation of the state vector $x_k$ leads to another equivalent state-space formulation of the same structure. It means that, viewed from sensor data, the state-space model of a given structure is not unique, and subspace system identification methods build models corresponding to arbitrary linear transformations. Nevertheless, all these equivalent state-space models contain the same modal parameters, notably, the eigenvalues of the matrix A, which are invariant regarding the linear transformations of $x_k$.

When SSI is applied to system (1), it is assumed that $w_k$ is a stationary white noise vector [3, 4]. There are situations where this assumption is by far not true, when $w_k$ is mainly due to external excitations, which may be strongly non-stationary and unknown. Such situations are considered in the next sections.

## 3. DISTURBANCE REJECTION

To address unknown inputs or disturbances, let us add an extra term $Ed_k$ into the state equation (1a), so that $w_k$ represents only a stationary white noise. Accordingly, the state-space model (1) becomes

$$x_{k+1} = Ax_k + Ed_k + w_k \tag{2a}$$

$$y_k = Cx_k + v_k, \tag{2b}$$

in which $d_k \in \mathbb{R}^p$ is most often referred to as disturbance, to avoid possible confusion with the known input or the control variable, and $E \in \mathbb{R}^{n \times p}$ is a matrix specifying the subspace to which the disturbance
is restricted. It is assumed that the dimension p of $d_k$ is smaller than the dimension n of $x_k$. In practice, it is related to the fact that, typically, disturbances affect only part of the degrees of freedom of the underlying system.

The disturbance rejection by output injection [9] was originally studied for the purpose of state estimation, i.e., compute the trajectory of the state vector $x_k$ from sensor measurements $y_k$, assuming that the system matrices A, C, E are known (as considered in [9], but this assumption will be discarded in the present paper). To recall the main results of [9], let

$$G = E[(CE)^T(CE)]^{-1}(CE)^T. \quad (3)$$

At the discrete time instant $k+1$, equation (2b) leads to

$$0 = G(y_{k+1} - Cx_{k+1} - v_{k+1}) \quad (4)$$

$$= Gy_{k+1} - GC(Ax_k + Ed_k + w_k) - Gv_{k+1}. \quad (5)$$

Inserting it into (2a) yields

$$x_{k+1} = (I - GC)Ax_k + (I - GC)Ed_k + Gy_{k+1} + (I - GC)w_k - Gv_{k+1}. \quad (6)$$

With G defined in (3), it is straightforward to check that

$$(I - GC)E = 0. \quad (7)$$

Then the unknown disturbance $d_k$ is eliminated from the state equation (6), which, together with (2b), constitutes the transformed state-space model

$$x_{k+1} = (I - GC)Ax_k + Gy_{k+1} + (I - GC)w_k - Gv_{k+1}. \quad (8a)$$

$$y_k = Cx_k + v_k. \quad (8b)$$

The state estimation problem considered in [9] treats $Gy_{k+1}$ as an input term in the state equation. Then the Kalman filter is applied to estimate the trajectory of $x_k$, assuming that the matrices A, C, E are known. In this paper, these matrices are unknown. The purpose is to estimate an equivalent model, while rejecting the unknown disturbance $d_k$.

## 4. SUBSPACE IDENTIFICATION WITH DISTURBANCE REJECTION

Subspace system identification methods, such as N4SID [5], are widely applied for building state-space models of the form

$$x_{k+1} = \bar{A}x_k + Bu_k + \bar{w}_k \quad (9a)$$

$$y_k = Cx_k + v_k \quad (9b)$$

by estimating the matrices $\bar{A}, B, C$ form data sequences of $u_k$ and $y_k$. Usually $u_k$ is the input of the underlying system, $y_k$ is its output. If such a method was directly applied to (8), then it would be related to (9) through

$$\bar{A} = (I - GC)A \quad (10)$$

$$B = G \quad (11)$$

$$u_k = y_{k+1} \quad (12)$$

$$\bar{w}_k = (I - GC)w_k - Gv_{k+1}. \quad (13)$$

Such a direct application cannot work correctly, as it would lead to a trivial model, due to the presence of the same sequence $y_k$ both as input and as output. The trivial model would be

$$x_{k+1} = u_k \quad (14a)$$

$$y_k = x_k \quad (14b)$$
corresponding to

$$ \bar{A} = 0, B = I, C = I, \bar{w}_k = 0, v_k = 0. \qquad (15) $$

In order to avoid such trivial results, the idea is to *inject part of $y_k$ corresponding to a subset of the available sensors*. Let a partition of $y_k$ and the corresponding partitions of $C$ and $v_k$ be

$$ C = \begin{bmatrix} C^a \\ C^{\bar{a}} \end{bmatrix}, \quad y_k = \begin{bmatrix} y_k^a \\ y_k^{\bar{a}} \end{bmatrix}, \quad v_k = \begin{bmatrix} v_k^a \\ v_k^{\bar{a}} \end{bmatrix}, \qquad (16) $$

where the exponent “a” indicates a subset of the rows of $C$, $y_k$ and $v_k$, whereas “$\bar{a}$” indicates the complementary rows. The sub-vector $y_k^a$ is then injected instead of the whole vector $y_k$. Accordingly, the matrix $G$ becomes

$$ G^a = E^{(a)}[(C^a E^{(a)})^T (C^a E^{(a)})]^{-1} (C^a E^{(a)})^T. \qquad (17) $$

where $C^a$ is part of $C$ as in (16). The matrix $E$ is *not* partitioned. Since $E$ may change due to non-stationary external force excitations, the notation $E^{(a)}$ is used to indicate the matrix $E$ when the data are collected and then processed by injecting $y_k^a$. Then (9) becomes

$$ x_{k+1} = \bar{A}^{(a)} x_k + G^a y_{k+1}^a + \bar{w}_k \qquad (18a) $$

$$ y_k = C^{(a)} x_k + v_k. \qquad (18b) $$

The matrices $\bar{A}^{(a)}$, $G^a$, $C^{(a)}$ are to be estimated by applying N4SID (or a similar identification method) from sensor data, with $u_k = y_{k+1}^a$ as input and $y_k$ as output. In this result, $C^{(a)}$ is the full matrix $C$ estimated by injecting $y_k^a$.

Subspace system identification with disturbance rejection assumes that some unknown disturbance $d_k$ is actually present. In practice, if the presence of disturbances (or non-stationary external excitations) is uncertain, the standard SSI should be first tried, without trying to reject any disturbance. If the result is satisfactory by checking the model residul, then it is useless to try the new method with output rejection. Otherwise different partitions of $y_k$ should be tried. Because the matrices $C^a$, $E^{(a)}$ are unknown, and (17) involves the inverse of $[(C^a E^{(a)})^T (C^a E^{(a)})]$, it is not possible to know in advance if a chosen partition of $y_k$ corresponds to a invertible $[(C^a E^{(a)})^T (C^a E^{(a)})]$.

With an appropriate partition of $y_k$, N4SID is then applied to (18) with $u_k = y_{k+1}^a$ as input. It yields an estimate of the matrix $\bar{A}^{(a)}$, which is related to the $A$ matrix of (2) in a way similar to (10). More specifically,

$$ \bar{A}^{(a)} = (I - G^a C^a) A. \qquad (19) $$

It is this matrix $\bar{A}^{(a)}$ which is mainly estimated, together with $G^a$, $C^{(a)}$ ($C^a$ is part of $C^{(a)}$), by the application of N4SID to (18) with $u_k = y_{k+1}^a$ injected.

It turns out that, with $G^a$ defined in (17), $(I - G^a C^a)$ is a singular (non invertible) square matrix. It is then impossible to deduce $A$ solely from $\bar{A}^a$ through (19), nor its eigenvalues. This is the second difficulty when applying subspace system identification methods with disturbance rejection. Because the arbitrary unknown disturbance $d_k$ cancels out the information in the subspace of the state-space related to $E^{(a)}d_k$, the collected data do not contain sufficient information about $A$. This is an intrinsic difficulty caused by arbitrary unknown disturbances, which are typically non-stationary natural excitations. To overcome this difficulty. different datasets will be combined together so that they compensate each other.

In the example of a mechanical structure naturally excited by wind, the disturbance $d_k$ is in fact due to the non-stationary part of the unknown external wind force. Its non-stationary nature implies that it may affect different degrees of freedom of the underlying mechanical structure, corresponding to different subspaces of the state-space associated to the state-space model (2). In other words, over the time, the
matrix E in (2) may change. After the first dataset processed by injecting $y_k^a$, some time later, a second dataset is collected and then processed by injecting $y_k^b$, which corresponds to another partition of $y_k$. Assume that the matrix E is $E^{(b)} \neq E^{(a)}$ when this dataset is collected. Then N4SID is applied, with $u_k = y_{k+1}^b$, to

$$x_{k+1} = \bar{A}^{(b)}x_k + G^b y_{k+1}^b + \bar{w}_k \qquad (20a)$$

$$y_k = C^{(b)}x_k + v_k, \qquad (20b)$$

yielding estimates of $\bar{A}^{(b)}$, $G^b$ and $C^{(b)}$, similar to the estimates of $\bar{A}^{(a)}$, $G^a$, $C^{(a)}$ corresponding to (18).

This result is then related to A through an equation similar to (19). Combining the new equation with (19) will help to solve for A. However, there is one more difficulty.

The extra difficulty is related to the fact that the state representation of a linear system is not unique. By applying a linear transformation $\bar{x}_k = Tx_k$ to (9) with any square invertible matrix $T$, an equivalent state-space model of the same system is then obtained with the new state vector $\bar{x}_k$ replacing $x_k$. Therefore, the same input-output data of (sequences of $u_k$ and $y_k$) can be generated either with $x_k$ or with $\bar{x}_k$ as state vector. Therefore, the result of system identification is not unique. Nevertheless, after the state transformation $\bar{x}_k = Tx_k$, the matrix $\bar{A}$ becomes $T\bar{A}T^{-1}$, which has the same eigenvalues as $\bar{A}$. Therefore, though the result corresponds to $T\bar{A}T^{-1}$ with some unknown matrix $T$, the eigenvalues of $\bar{A}$ can still be estimated from this result.

Due to this non-uniqueness of the result of system identification, after the second identification corresponding to (20) yielding the estimates of $\bar{A}^{(b)}$, $G^b$ and $C^{(b)}$, to link the result with the matrix A, the equation similar to (19) is

$$\bar{A}^{(b)} = (I - G^b C^b) T A T^{-1}, \qquad (21)$$

with some unknown matrix $T$. Moreover, $C^{(a)}$ and $C^{(b)}$ are related by

$$C^{(a)} = C^{(b)}T. \qquad (22)$$

The proposed solution is then to extract A from (19), (21) and (22). In these equations, $\bar{A}^{(a)}$, $G^a$, $C^a$, $C^{(a)}$ and $\bar{A}^{(b)}$, $G^b$, $C^b$, $C^{(b)}$ are already estimated by the two applications of N4SID. The unknowns A and T are then solved for in the least squares sense.

## Summary of the whole identification procedure

1. Collect sensor data consisting of a sequence of $y_k$.

2. Apply standard stochastic subspace identification to the data by assuming a model in the form of (1), in which $w_k$ is a stationary white noise. If the result is satisfactory by checking the model residual, then return the estimated model and stop the procedure, otherwise go to the next step.

3. Choose a subset of sensors corresponding a sub-vector of $y_k$, say $y_k^a$. Apply a subspace identification method with the input data $u_k = y_{k+1}^a$ and the output data corresponding to the full vector $y_k$. If the result is satisfactory by checking the model residual, return the estimates of $\bar{A}^{(a)}$, $G^a$, $C^a$, $C^{(a)}$ and go to the next step. Otherwise repeat this step by choosing another sub-vector of $y_k$.

4. Some time after the collection of the previous data set, collect a new data set when the non-stationary external forces may have significantly changed. For example, if the external forces are due to wind, wait so that the wind direction has changed since the last dataset collection.

5. Choose a subset of sensors corresponding a sub-vector of $y_k$, say $y_k^b$. Apply a subspace identification method to the new dataset with the input data $u_k = y_{k+1}^b$ and the output data corresponding to the full vector $y_k$. If the result is satisfactory by checking the model residual, return the estimates of $\bar{A}^{(b)}$, $G^b$, $C^b$, $C^{(b)}$ and go to the next step. Otherwise repeat this step by choosing another sub-vector of $y_k$. If all the sub-vectors of $y_k$ have already been tried, go back to Step 4.
6. Solve for A and T in the least squares sense the equations

$$ \bar{A}^{(a)} = (I - G^a C^a)A \qquad (23) $$

$$ \bar{A}^{(b)} = (I - G^b C^b)TAT^{-1} \qquad (24) $$

$$ C^{(a)} = C^{(b)}T. \qquad (25) $$

Return the eigenvalues of A and stop the procedure.

## 5. SIMULATION EXAMPLE

In this preliminary study, consider a simple simulation example to illustrate the proposed identification procedure. Data are generated with a state-space model in the form of (2). Because of the continuous time nature of most physical systems, matrices are first chosen for a continuous time model:

$$ A_c = \begin{bmatrix} -0.9 & 0.9 & 0 & 0 \\ -0.9 & 0 & 0.90 & 0 \\ 0 & -0.9 & -0.90 & 0 \\ 0 & 0 & 0 & -1 \end{bmatrix}, \quad C_c = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}, \quad E_c = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \qquad (26) $$

The corresponding discrete time model matrices, as in (2), are then computed as

$$ A = e^{A_c \tau}, \quad C = C_c, \quad E = A_c^{-1}(A - I)E_c \qquad (27) $$

with the sampling period $\tau$ equal to 1 ms. These parameters are used with (2) for generating the first dataset. To simulate the non-stationary nature of external excitation, a second dataset is simulated with

$$ E_c = [0 \ 0 \ 1 \ 0]^T, \qquad (28) $$

while the other parameters remain unchanged.

The disturbance $d_k$ is simulated as

$$ d_k = 5 \sin(0.1k) + 5 \sin(0.1\sqrt{2}k) + g_k \qquad (29) $$

where $g_k$ is a white Gaussian noise of unitary variance, and $\sqrt{2}$ is added so that the deterministic part is (in theory) not periodic. Each component of $w_k$ is independently generated as a white Gaussian noise of variance 0.001, and similarly each component of $v_k$ has the variance 0.0005.

The matrix A has two complex eigenvalues $0.2366 \pm 0.5921i$ and two real eigenvalues $0.4066$ and $0.3679$. As the mode corresponding to the last eigenvalue is not observable, the data generated by the 4th order system could also be generated by an equivalent 3rd order system. For this reason system identification leads to a 3rd order model.

When system identification with disturbance rejection is applied to the first data set, the agreement between the simulated state sequence $x_k$ and the state sequence $\hat{x}_k$ computed with the estimated model is inspected to evaluate the result. This is possible in a simulation study. Because system identification leads to a model corresponding to an arbitrary basis of the state-space, the linear transformation between $x_k$ and $\hat{x}_k$ is first found by solving a least squares problem. Then the estimated state sequence $\hat{x}_k$ is transformed accordingly to be compared with the simulated state sequence $x_k$. In Figure 1, $\hat{x}_k$ is the estimated state sequence after the linear state transformation. On the left side, its 3 components are superposed on the 3 corresponding components of $x_k$. The simulated system has 4 state components, of which the 4th one is not observable. The 3 trajectories shown in Figure 1 are those of the observable components. Because the difference between $x_k$ (blue curves) and $\hat{x}_k$ (red curves) is barely visible, their difference $x_k - \hat{x}_k$ is plotted on the right side of Figure 1.
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Left: Components of the estimated state sequence $\hat{x}_k$ superposed on the observable components of the simulated state sequence $x_k$. Right: the difference $x_k - \hat{x}_k$. The estimated state sequence has been linearly transformed so that it corresponds to the same state-space basis as the simulated state sequence.</figcaption></figure>

Then system identification is applied to the second data set, and the matrix $A$ is extracted from the two results of system identification, as summarized in the procedure at the end of the previous section. To compare the estimated complex mode with the simulated complex mode. Let $a \pm bi$ denote the pair of complex eigenvalues of the simulated system, and $\hat{a} \pm \hat{b}i$ their estimates. The relative errors for the real part and the imaginary part are respectively defined as

$$ \text{Real part relative error} = \frac{a - \hat{a}}{|a|}, \quad \text{Imaginary part relative error} = \frac{b - \hat{b}}{|b|}. \qquad (30) $$

Their mean values and standard deviations evaluated over 100 random realizations of the simulation and estimation are presented in Table 1.

<table><caption>**Table 1:** Means and standard deviations of the real part relative error and the imaginary part relative error, as defined in (30), evaluated over 100 random realizations.</caption><thead><tr><th></th><th>Relative error mean</th><th>Relative error standard deviation</th></tr></thead><tbody><tr><td>Real part</td><td>-0.0103</td><td>0.0467</td></tr><tr><td>Imaginary part</td><td>-0.0169</td><td>0.0366</td></tr></tbody></table>

## 6. CONCLUSION

In this paper, techniques of subspace system identification and unknown input observers are associated in some way, resulting in a subspace identification method rejecting unknown inputs. A simple simulation example has been presented to illustrate the effectiveness of the proposed method for modal properties characterization. Future work will address numerical efficiency in order to apply the method to larger systems.

## REFERENCES

[1] P. Van Overschee and B. De Moor. *Subspace Identification for Linear Systems: Theory, Implementation Applications*. Springer-Verlag, 1996.

[2] Lennart Ljung. *System Identification – Theory for the User*. Prentice-Hall, 2nd edition edition, 1999.

[3] Michael Döhler and Laurent Mevel. Efficient multi-order uncertainty computation for stochastic subspace identification. *Mechanical Systems and Signal Processing*, 38(2):346–366, 2013.
[4] Fahad Bin Zahid, Zhi Chao Ong, and Shin Yee Khoo. A review of operational modal analysis techniques for in-service modal identification. *Journal of the Brazilian Society of Mechanical Sciences and Engineering*, 42(8):398, 2020.

[5] Peter Van Overschee and Bart De Moor. N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. *Automatica*, 30(1):75-93, 1994.

[6] He Kong and Salah Sukkarieh. An internal model approach to estimation of systems with arbitrary unknown inputs. *Automatica*, 108:1-11, 2019. doi: https://doi.org/10.1016/j.automatica.2019.06.034.

[7] Bernard Delyon and Qinghua Zhang. On the optimality of the Kitanidis filter for state estimation rejecting unknown inputs. *Automatica*, 2021. URL https://hal.inria.fr/hal-03041232. Published on-line https://doi.org/10.1016/j.automatica.2021.109793.

[8] P. K. Kitanidis. Unbiased minimum variance linear state estimation. *Automatica*, 23(6):775-778, 1987.

[9] Qinghua Zhang and Liangquan Zhang. State estimation for stochastic time varying systems with disturbance rejection. In *18th IFAC Symposium on System Identification (SYSID 2018)*, pages 55-59, Stockholm, Sweden, 2018.

[10] Jer-Nan Juang. *Applied system identification*. Prentice-Hall, Inc., 1994.