## 2 Brief Review of MIMO-NARX Modeling
### Model Order Reduction  
Although MIMO-NARX is suitable for multiple outputs, it is computationally expensive to establish the NARX model for each degree of freedom. The fact that the first few modes dominate the dynamic response of a multi-degree-of-freedom structure can be utilized for eliminating the burden due to a large number of degrees of freedom. The MOR method [15; 16; 17] thus can be utilized to lower the degrees of freedom of the MDOF nonlinear dynamic system in this study. Using a small number of lower structural modes corresponding to the number of DOF, an accurate approximation of the response history can be acquired. The response vector in the reduced space can be obtained through the coordinate transformation matrix and the response vector in the original space:  
\[\mathbf{q}(t)=\mathbf{\Phi}_{n_{r}}^{\mathrm{T}}\mathbf{x}(t) \tag{7}\]  
where \(\mathbf{q}(t)\) is the response in the reduced \(n_{r}\)-dimensional space; \(\mathbf{x}(t)\) is the response vector in the original \(n\)-dimensional space; \(\mathbf{\Phi}_{n_{r}}^{\mathrm{T}}\) is the \(n_{r}\times n\) coordinate transformation matrix with \(n_{r}\ll n_{r}\). The transformation matrix can be obtained by singular value decomposition (SVD) [24; 25]. For the \(n_{s}\) samples of the simulation results for the MDOF dynamic system under uncertainties, a snapshot matrix can be formulated as follows:  
\[\mathbf{X}=[\mathbf{x}_{1}(t_{1}),\dots,\mathbf{x}_{1}(t_{n_{1}}),\dots,\mathbf{x}_{n_{s}}(t_ {1}),\dots,\mathbf{x}_{n_{s}}(t_{n_{r}})] \tag{8}\]  
where \(\mathbf{X}\in\mathbb{R}^{n\times(n_{r}n_{s})}\) contains \(n_{t}\) time step responses of the \(n\)-DOF dynamic system. The SVD can then be used to decompose the snapshot matrix \(\mathbf{X}\) as:  
\[\mathbf{X}=\mathbf{U}\mathbf{\Lambda}\mathbf{V}^{\mathrm{T}} \tag{9}\]  
where \(\mathbf{U}\) is a \(n\times n\) orthonormal matrix that contains the left singular vectors of \(\mathbf{X}\); \(\mathbf{V}\) is the \((n_{t}n_{s})\times(n_{t}n_{s})\) orthonormal matrix that contains the right singular vectors; \(\mathbf{\Lambda}\) is a \(n\times(n_{t}n_{s})\) pseudo-diagonal matrix with the singular values \(\mathbf{\Lambda}(j,j)=\lambda_{j}\), \(j=1,\dots,n\). A computationally effective economy-size SVD is used in view of the extremely large size of \(\mathbf{V}\) with \((n_{t}n_{s})\gg n\)[14]. The left singular vectors \(\mathbf{U}\) are the proper orthogonal decomposition (POD) modes.  
Then the coordinate transformation matrix is defined as the first \(n_{r}\) POD modes as:  
\[\mathbf{\Phi}_{n_{r}}=[\mathbf{U}_{1},\mathbf{U}_{2},\dots,\mathbf{U}_{n_{r}}] \tag{10}\]  
where \([\mathbf{U}_{1},\mathbf{U}_{2},\dots,\mathbf{U}_{n_{r}}]\) is the first \(n_{r}\) columns of \(\mathbf{U}\); \(n_{r}\) can be determined by the energy criterion which ensures the energy captured in the truncated representation of \(x\) reaches \(\eta\) of the total energy [17]:  
\[\sum_{j=1}^{n_{r}}\lambda_{j}\geq\eta\sum_{j=1}^{n}\lambda_{j} \tag{11}\]  
where \(\eta\) is typically assumed to be close to 1, for example, 0.99. A properly selected \(\eta\) can bring significant dimensional reduction to the system, that is, \(n_{r}\ll n\). The transformation error of MOR can be measured by the comparison between the original response \(\mathbf{x}(t)\) from Equation (7) and the reverse-transformation response \(\mathbf{x}_{-}(t)\), which can be calculated as:  
\[\mathbf{x}_{-}(t)=\mathbf{\Phi}_{n_{r}}\mathbf{q}(t) \tag{12}\]