{'fraction_non_alphanumeric': 0.0430297671676982, 'fraction_numerical': 0.012967875036840553, 'mean_word_length': 4.169408897014016, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving Linear State-Space Models with Additional Iterations

An estimated state-space model can possibly be improved by further iterations with estimation data. This contribution specifically studies if models obtained by subspace estimation can be improved by subsequent re-estimation of the B, C, and D matrices (which involves linear estimation problems). Several tests are performed, which shows that it is generally advisable to do such further re-estimation steps using the maximum likelihood criterion. Stated more succinctly in terms of MATLAB functions, ssest generally outperforms n4sid.


INTRODUCTION
Linear state-space models are perhaps the most common model structure used in system identification. They can be estimated in several different ways. Among the most common techniques are so called subspace methods, such as MOESP, Verhaegen and Dewilde (1992), and N4SID, Overschee and Moor (1994). These were also developed to work with frequency domain data, McKelvey et al. (1996). Whichever way the model was obtained, an interesting the question is if it can be improved in some way by further polishing using observed data.
In this contribution we investigate how further iterations on the B, C and D matrices can possibly improve the model quality and also whether it is worthwhile to reestimate the A matrix.
We argue that such further estimation work is well motivated. The model properties are improved most of the time (but not always), and in some cases, the improvements are significant.
MATLAB R and Simulink R are registered trademarks of The MathWorks, Inc. See mathworks.com/trademarks for a list of additional trademarks. Other product or brand names may be trademarks or registered trademarks of their respective holders.
Data for Section 7, Fig. 5,and Fig. 6 presented in this paper is provided courtesy of Volvo Car Corporation. Presentation of the plot and the study results do not grant authorization to extracting and reusing the data for any purpose. Extraction and reuse of the data require authorization from Volvo Car Corporation.

THE STATE-SPACE MODEL
A linear state-space model in output error form (no noise model) is given by By assuming Gaussian noise distribuion for e, the maximum likelihood estimate (MLE) of the matrices is given by which can be rewritten as We note that for fixed A, C the predictionŷ is linear in B, D and for fixed A, B the prediction is linear in C, D. Remark 1. The matrices A, B, C, D may be parameterized by some paremeter vector θ. The estimation and minimization is carried out with respect to θ. It will be assumed that any such parameterization is linear in θ.

REFINING THE ESTIMATE
With a given state-space model (1) and access to data [y, u] from the system (either data that was used to estimate (1) or a fresh dataset), it is natural to ask if the maximum likelihood (ML) method in Eq.
(2) can be used to refine the estimate.
y(t|A, B, C, D) is linear in B, D for fixed A, C, and therefore reestimating B, D is a simple linear regression problem. Once this is done, one could fix A, B and reestimate C, D which again is a linear regression problem. Repeating that gives some simple B, C, D iterations, which is a way of minimizing (2a) for a fixed A, which is a bilinear problem. See, for example, Eq. (10.68) in Ljung (1999). Simple experiments show that the B, C, D iterations typically lead to the same (possibly local) minimum of Eq. (2a) as applying the Gauss-Newton method to minimize Eq. (2a) for fixed A.
The aforementioned B, C, D iterations is a block coordinate descent approach. It is faster than classical coordinate descent since in two steps all coefficients in B, C, D are updated and at each iteration more accurate values are used from previous iteration: updated C in B, D iteration and updated B in C, D iteration. The same approach is called separable least squares for a slightly different cost function appearing in identification of Hammerstein systems, Bai and Li (2004). The convergence of block coordinate descent is guaranteed when the objective function has a unique minimum in each coordinate block, see pages 261-262 in Luenberger and Ye (2008). This condition corresponds to the least squares problem in Eq. 2 not being rank deficient during each B, D and C, D iteration. This is frequently satisfied in practice per our experiments.
After the refinement of B, C, D is done, one can indeed fix these parameters and reestimate A using Eq. (2a) -or alternatively minimize the ML criterion with respect to all matrices simultaneously.

SISO, MISO, SIMO, AND MIMO
A typical implementation of subspace methods extract an estimate of the A, C matrices first, followed by a linear regression for the B and D matrices with fixed A, C. The result is not necessarily a (local) optimum in terms of A, C for MIMO systems per the cost function Eq. (2a), and refinements are possible.
However, there are two special cases: • SISO and MISO systems: The solution of the linear regression for B, D is also a (local) optimum in terms of C. In other words, unless A is refined, further iterations for B, C, D matrices will not reduce the cost function. • SIMO systems: The solution of the linear regression for B, D is not a local optimum in terms of C. However, solving one additional linear regression problem for C, D matrices with fixed A, B corresponds to a (local) optimum with respect to the B matrix. After this additional linear regression, further iterations on B, C, D matrices are not needed unless A is changed.
To clarify the statements above we have the following results: be an order n state-space realization of a MISO linear system. Let C 0 be any row vector of the same size as C such that (A, C 0 ) is an observable pair. Then there exists matrices B 0 and D 0 such that the realization (A, B 0 , C 0 , D 0 ) is input-output equivalent with (A, B, C, D).
Proof: Trivially we have the choice D 0 = D and we will show that there exists a B 0 such that the rest of the impulse response matrices coincide. Pick any matrix Since the pair (A, C 0 ) is observable the corresponding observability matrix has full rank which imply that the row vectors C 0 A i to the right in (3) are all linearly independent. Hence a solution exists. Note that by construction T A = AT . Finally we have the identity (4) which directly shows that B 0 = T B. Theorem 1. Consider the case when the A matrix is fixed and we define (A, B * , C * , D * ) to be a state-space realization which minimizes (2b), w.r.t. A, B, C with A fixed. Let C 0 be any row vector such that (A, C 0 ) is observable. Define B 0 and D 0 to be the minimizers to (2b) w.r.t. B and D, when A and C 0 are fixed. Then V (y, u, A, B * , C * , D * ) = V (y, u, A, B 0 , C 0 , D 0 ).
Proof: By Lemma 1 we can conclude that for the realization (A, B * , C * , D * ) there exists an input-output equivalent realization (A, B , C 0 , D ) which imply V (y, u, A, B * , C * , D * ) = V (y, u, A, B , C 0 , D ) since the loss function only depend on the input output properties and not on the specific realization. Since the minimization of V w.r.t. B and D only has a subset of the free parameters compared to the minimization w.r.t. B, C, D it follows that B and D are the minimizer to V , i.e. B 0 = B and D 0 = D .
The results above can also be formulated and proved for SIMO systems.
Consider the following result. Theorem 2. Given a matrix A ∈ R n×n and assume there exists a vector v such that (A, v) is observable. Then for a matrix B the following two statements are equivalent i=0 forms a basis for R n . Then it follows that Since B commutes with A by assumption it also commutes with all powers of A. Hence for all integers r ≥ 0 Since the set {A r v} n−1 i=0 span R n we have shown that for any vector x we have Bx = ( The result in Theorem 2 shows a direct limitation for MIMO systems. The set of matrices that commutes with A is only an n-dimensional subspace of all n × n matrices. Hence we can only fix n parameters in C 0 , e.g. one row as shown in the proof of Lemma 1. The remaining rows in C 0 must be free parameters to optimize over. For MIMO systems it is thus beneficial to iterate between minimizing V w.r.t. C, D and B, D respectively until convergence. Theorem 1 and Theorem 2 solely rely on the observability assumption on the (A, C) pair. One extra simplifying assumption can be made to illustrate the underlying structure of the least squares problem Eq. (2a). Assume A has distinct eigenvalues. Then A admits the decomposition A = P ΛP −1 with a diagonal matrix Λ. Rewrite Eq. (2c) with B ∈ R n×nu , C ∈ R 1×n , D ∈ R 1×nu exactly as: . . , n} are the diagonals of Λ. For real-valued A, B, and C matrices the columns ofC and the rows ofB come in conjugatepairs. Since Λ t−k is diagonal,CΛ t−k = L t−k diag(C) where diag(C) is the diagonal matrix with elements of the vectorC on the diagonals. Therefore: This form shows that the eigenvalues of A, contained in L t−k , form a set of basis functions for the linear regression problem for B (equivalently,B). Fixed C, and in turn diag(C), is a scaling of the basis functions. Note thatC ∈ C 1×n having a zero element is equivalent to (A, C) pair not being observable when A has distinct eigenvalues. When an element ofC is zero, the corresponding basis function is multiplied by zero and not utilized in regression.
A less formal view of Theorem 1 and its proof can also be seen from this form. Note that the unknowns in B, C can be combined into a single matrix X = diag(C)B, and the linear regression can be performed for X, D for fixed A, where the rows of X are in conjugate-pairs if A, B, C are real-valued. diag(C) is nonsingular when all elements ofC are nonzero. Any fixed choice of C that correspond to all nonzero elements inC = CV (i.e. (A, C) pair is observable when A has distinct eigenvalues) can be used to extract a B matrix estimate from X, without affecting the cost function value observed in Eq. (2a).
The remarks made through distinct eigenvalues in A assumption can be relaxed by replacing the eigenvalue decomposition A = P ΛP −1 with Jordan matrix decomposition A = SJS −1 . J is block diagonal, where the blocks J i ∈ C ni×ni are elementary Jordan blocks. Then A t−k = SJ t−k S −1 and J t−k is again block diagonal with the blocks J t−k i as upper triangular Toeplitz matrices. This structure of the J t−k i matrix blocks can be exploited to again collect the unknowns in B, C into a single matrix. Similarly, which of the basis functions are utilized in the linear regression can be seen through (potentially a subset of) the elements inC = CS.
The optimality results for the B, C, D matrices presented in this section are also true for frequency-domain identification when a cost function of the following form is utilized: The cost function Eq. is typical for many estimation approaches including the subspace methods. Here z k = e jTsω k is a point on the unit disk and G k is a frequency response measurement corresponding to the k th frequency point ω k rad/s. It follows that since also this cost function only depends on the input-output properties of the statespace system the same argumentation as made for the time-domain case also holds here.

DISCRETE-TIME MODELS WITH TIME DOMAIN DATA
The following experiment was performed in MATLAB Release 2018a (2018): (1) Generate 200 discrete-time systems randomly with drss. The systems are of order 7 with 4 inputs and 4 outputs. The feedthrough matrix D was set to zero. (2) Simulate each system with a Gaussian white input with unit variance matrix. No attempt was made to tailor the input to the system properties. The noisefree data are denoted by z. Add white Gaussian noise with unit covariance matrix to the output, giving the estimation data sequence zn.
Compute the error in that model compared to the true system by executing the commands: en=pe(z,mn); mnn=norm(en.y'*en.y/N) (4) Readjust the B, C parameters to data: mpBC=nn; mpBC.Structure.A.Free=zeros(7,7); mpBC=ssest(zn,mpBC); (5) Compute the error of that model: en=pe(z,mpBC); mnp=norm(en.y'*en.y/N) (6) Readjust the A, B, C parameters of the model by further Maximum Likelihood iterations on mn: mp=ssest(zn,mn); en=pe(z,mp); mnp=norm(en.y'*en.y/N) The results are summarized in the boxplot in Fig. 1. The medians of the error norms are The error for mn was larger than mpBC in 73.5 % of the cases. The error in mnBC was larger than mp in 80.5 % of the cases and mp outperformed mn for 83 % of the data sets.
A scatter plot for the error norms for the 200 tested systems is given in Fig. 2. It is in accordance with the observations in Ljung (2003). (1) Generate 200 seventh order continuous-time systems with 4 input and 4 outputs using G = rss(7,4,4) possibly with feedthrough terms D. Compute their frequency response functions Gf = idfrd(G,FG) at 410 linearly spaced frequencies FG. No attempts were made to adjust FG to the dynamics of G. Add 20 % random, multiplicative noise to the response to form the frequency domain data Gfn.
(2) Estimate 7th order state-space models from the noisy frequency response functions, using n4sid and the B, C adjusted models mnBC and A, B, C adjusted models mp and their errors compared to the true system as in the previous section.
The medians of the errors were found to be: The results are summarised in the boxplots in Fig. 3.
The error for mn was larger than mpBC in 83.5 % of the cases. The error in mnBC was larger than mp in 74.5 % of the cases and mp outperformed mn for 87 % of the data sets.
A scatter plot for the error norms for the 200 tested systems is given in Fig. 4.

REAL DATA FROM FLEXIBLE STRUCTURES
Mechanical vibration testing was performed on a Volvo XC90 (2015) rear subframe structure. During the testing the subframe was equipped with 26 uniaxial and 10 triaxial accelerometers yielding a total of 56 measurement channels. A shaker was used to provide excitation at two different locations. Frequency data was obtained by employing a single input multiple output (SIMO) stepped sine testing procedure for each shaker location. This testing directly yields measurements of the frequency response function at the excited frequencies. For each of the shaker locations, a total of 2998 different frequencies were excited in the frequency range between 60 and 500 Hz with a frequency spacing derived according to the method described in Vakilzadeh et al. (2015). For further details of the experiment, refer to Gibanica and Abrahamsson (2017). Using the frequency domain subspace method described in McKelvey et al. (1996), an initial continuous-time model was derived of order 40 with 56 outputs and 2 inputs. Subsequently the C and D matrices are iteratively reestimated followed by an estimate of the B and D matrices. After the first B, C, D iteration the least-squares loss function is reduced to 58.5% of the initial value. After four additional B, C, D iterations the loss function is slightly decreased to 56.9%. In Fig. 5 the magnitude of the final frequency response function corresponding to input 1 and output 7 is plotted together with the error magnitude, that is, the magnitude of the difference between the data and model. The error magnitude of the initial model is also plotted in the graph. For this rather large-scale example it is clear that it is advisable to employ at least a first B, C, D iteration after the initial estimate delivered by the subspace method.  et al. (1996), is used as the initial model for all approaches.
The cost function value obtained in the subsequent steps for all methods are normalized by the value attained by this initial step. For this specific data and the initial model delivered by the subspace methods, the block coordinate descent approach B, C, D iterations approximately converge in one step to a (local) minimum. This method was also the fastest. Gauss-Netwon approach for optimizing B, C, D matrices jointly approximately converged in two steps. Optimizing over A, B, C, D jointly can potentially find a better (local) minimum compared to the previous two methods since the A matrix is also optimized, but the progress was slower. Indeed once any optimization approach for B, C, D iterations converges to a minimum, this could be followed with an optimization over the A matrix.

CONCLUSIONS
This contribution shows that it is useful to do further iterations on the B, C, and D state-space matrix estimates obtained by the subspace identification when the system has multiple outputs. This is tested with both simulated and real data. It is also recommended to adjust the A estimate using maximum likelihood iterations. It leads to improvements in a majority of cases, but not always. One reason is that subspace identification algoritms contain several design variables (prediction horizons and prefilters) that can be difficult to choose in an optimal way. In a sense, maximum likelihood iterations (such as in ssest) perform such choices automatically.