Application of Anti-Diagonal Averaging in Response Reconstruction

Abstract

Response reconstruction is used to obtain accurate replication of vehicle structural responses of field recorded measurements in a laboratory environment, a crucial step in the process of Accelerated Destructive Testing (ADA). Response Reconstruction is cast as an inverse problem whereby an input signal is inferred to generate the desired outputs of a system. By casting the problem as an inverse problem we veer away from the familiarity of symmetry in physical systems since multiple inputs may generate the same output. We differ in our approach from standard force reconstruction problems in that the optimisation goal is the recreated output of the system. This alleviates the need for highly accurate inputs. We focus on offline non-causal linear regression methods to obtain input signals. A new windowing method called AntiDiagonal Averaging (ADA) is proposed to improve the regression techniques’ performance. ADA introduces overlaps within the predicted time signal windows and averages them. The newly proposed method is tested on a numerical quarter car model and shown to accurately reproduce the system’s outputs, which outperform related Finite Impulse Response (FIR) methods. In the nonlinear configuration of the numerical quarter car, ADA achieved a recreated output Mean Fit Function Error (MFFE) score of $0.40\%$ compared to the next best performing FIR method, which generated a score of $4.89\%$. Similar performance was shown for the linear case.

1. Introduction

In the drive for continual improvement in vehicle engineering design, optimised structures and components with lower safety margins and greater reliability are sought [1]. Advances in computational design such as finite element analysis and dynamic modelling combined with fatigue prediction have furthered this goal tremendously during the design phase. Nevertheless, there is still a need to dynamically test physical prototypes or existing designs in a controlled laboratory environment. For the analysis to be worthwhile the excitation of the structure in the laboratory environment must induce responses in the structure as though it were being tested under real-world operating conditions. The end goal is to enable Accelerated Destructive Testing (ADT) of the structure. In ADT, a vehicle’s chassis is mounted with its suspension system on a set of hydraulic actuators. The hydraulic actuators then excite the system vertically. Laterally acting forces are simulated with additional actuators. An example of an ADT set-up is shown in Figure 1.

The structure’s excitation is then carried out for extended periods, allowing for the degradation of the structure to be measured in a controlled environment [1]. The structure is not typically excited until catastrophic failure, but rather until the degradation measured as vibration or noise has met a specified threshold [2]. This indicates possible failure points of the system and a means of predicting the component’s healthy lifespan. Other insights can be gained from dynamic testing, such as a better understanding of the system dynamics, vibration isolation [1] and vibration severities for passenger ride comfort [3].

The biggest hurdle with ADT is that the inputs to the system, such as the displacements or the forces acting on the vehicle’s tyres, are difficult or impossible to measure directly in the field. This means the problem must be cast as an inverse modelling or response reconstruction problem [4]. In inverse problems, the outputs of the system Z are used in conjunction with model parameters $\beta$ to determine the inputs U, i.e.,

$$U=f(Z,\beta ).$$

(1)

There are two possible choices for creating a model of the system. A mapping of the system can be constructed so that the system’s inputs are used to predict the outputs of the system. This is referred to as the forward problem. The forward problem is then inverted. If the model used to map the problem is nonlinear, an iterative optimisation scheme is employed to invert the system. However, the optimisation scheme may be prone to local minima. The second approach is to create a direct inverse of the model whereby the system’s outputs are used to predict the inputs of the system. The inverse method has an inherent stability check since the solution will only be obtained if the direct inverse model is stable [1]. However, we quickly find that most inverse problems are ill-posed. For a problem to be well-posed it needs to meet the following criteria: the solution is unique, the global solution exists for all data and the solution to the problem is continuously dependent on the given data [5]. The first criterion is normally the offending culprit since it is easy to construct a forward problem where two different inputs result in the same output. Therefore, the inverse solution is typically not unique. This introduces an asymmetry into the problem whereby the assumption of an one-to-one mapping is broken. If the problem is ill-posed we may use regularisation techniques to cast the problem as a more well-behaved problem. Regularisation techniques include: cross validation, SVD, iterative methods, data filtering and Tikhonov regularisation [4].

Most common reconstruction techniques are implemented in the frequency domain [6], whereby the discrete Fourier response is multiplied by an inverse or pseudoinverse frequency response function [7]. Raath’s Ph.D. thesis [1] highlighted the then known issues of using frequency response techniques in accelerated fatigue testing. It was shown that the frequency response was inaccurate for several reasons, that include:

• Assuming that the input and output signals are periodic when often they were not. These include sharp impulses from random impacts.
• Being unable to model nonlinear models since frequency response analyses assume a linear model.
• Requiring long time signals of the order of hours as opposed to minutes or seconds needed for the time domain. This ties in with the issue that low frequency information is easily lost due to spectral leakage where the energy in the lower frequencies is spread over to higher frequencies.
• Failing to capture the sequence or causal effects which play an important role in crack propagation.

Various time-domain techniques have been developed to overcome this. However, they have been shown to be slow or inaccurate [8]. The vehicle structures of interest typically contain many nonlinear components such as springs and pneumatic dampers. Typical control systems will overcome this issue by linearising the system around the operation point. It is then assumed that the system will experience small perturbations around this point. However, it is expected that the system will experience impact loadings and large displacements, which will force the system out of its linear region [1]. Another issue associated with response reconstruction is that of model mismatch, whereby the identified system does not truly represent the physical test rig. In response reconstruction the misrepresentation occurs when the physical test rig is taken from the real-world and recreated and simulated in the laboratory environment. Typically the degrees of freedom are not fully represented in the laboratory or the test rig parameters, such as mass, may vary. In this laboratory environment, the process of system identification occurs; therefore, we will have mapped a domain that differs from the real-world domain. When the real-world outputs need to be recreated we may find that the mapped inverse model may be forced to extrapolate into regions in the mapped domain to find a solution. In other words, the inverse model has over-fitted to the laboratory domain and generalises poorly with regard to the real-world domain. Regularisation can be employed to minimise this error [9]. A related field to response reconstruction is force identification whereby the inputs of the system are of interest. However, the inputs of the system for a given output are not unique [10]. Force identification tackles this problem by enforcing some prior knowledge of the system dynamics to constrain the inputs to reasonable solutions. Bayesian methods have become prevalent in force identification literature since they allow the experimenter to systematically incorporate prior knowledge [11]. Another benefit of incorporating Bayesian methods is that it provides for confidence intervals on the input predictions and model parameters [12]. A noticeable distinction in force identification literature is that a known finite element model of the structure is typically assumed, i.e., a known forward model.

The approach taken in this paper of solving the issue of non-uniqueness of the input is mitigated by

• not focusing on the reconstructed input accuracies.
• using cross validation of the system’s reconstructed outputs to determine whether a given inverse model of the system is satisfactory i.e., using a forward pass through the physical system in each cross validation step to determine the model accuracy.

A potential drawback to this approach is that, if implemented naively, the cross validation can induce undue stress on the system before any ADT occurs. An overview of the response reconstruction methodology used in this paper is given in Figure 2.

This paper focuses on linear regression methods for mapping the relationship between the outputs X and inputs Y for response reconstruction, i.e., $X\beta =Y$. The core contribution of this paper is the proposed method of extending the capabilities of said linear regression methods by introducing overlaps and merging them using averaging in a process called AntiDiagonal Averaging (ADA), encapsulated in Equation (12). We show that ADA is closely related to FIR methods. We benchmark ADA in terms of its response reconstruction ability as well as its performance against the related FIR methods. We focus on Tikhonov regularisation with cross validation through the use of Ridge Regression (RR) to regularise the inversion of the system. Any suitable linear regression method can be employed with ADA; however, RR is needed for the FIR methods we cover.

We first give a brief overview of RR and how it enforces regularisation. The theory behind ADA is then introduced and compared against related FIR methods. The design of the investigation is then given with an overview of the numerical quarter car model, with which the reconstruction methods are benchmarked. The results of the benchmarks are then discussed. Finally, an illustrative comparison of the different regression methods is conducted showing the performance of the regression method on a challenging response reconstruction problem.

1.1. Ridge Regression

As opposed to discretely truncating the singular values, RR instead smoothly decays the singular values through the use of a regularisation matrix $\mathsf{\Gamma }$ which results in the solution

$$\begin{array}{c}\hfill \beta ={(X^{\prime }X+{\mathsf{\Gamma }}^{\prime }\mathsf{\Gamma })}^{-1}{X}^{\prime }Y,\end{array}$$

(2)

where $\mathsf{\Gamma }$ is typically chosen as a scaling of the identity matrix though the use of the regularisation constant $\alpha$, i.e., $\mathsf{\Gamma }=\alpha I$. RR has the solution in terms of the SVD of X [13]

$$\widehat{Y}=X\left(V_{x}D{U}_x^{\prime }\right)Y,$$

(3)

where the entries of the diagonal matrix D are given by

$$D_{ii}=\frac{s_i}{s_i^2+{\alpha }^2}.$$

(4)

$U_x$ and $V_x$ are the left and right singular vectors of X, respectively, with the corresponding singular values $s_i$.

1.2. ADA

Windowing methods are needed to represent the responses Z as the predictor matrix $X\in {\mathbb{R}}^{n\times p}$ and the inputs U as the target matrix $Y\in {\mathbb{R}}^{n\times r}$ in any suitable linear regression method. This is achieved by windowing said signals and treating each window as an observation. The original input and response measurements are given by $U\in {\mathbb{R}}^{m\times q}$ and $Z\in {\mathbb{R}}^{m\times o}$ where m is the original sequence length in samples and q and o are the number of actuator and sensor channels, respectively.

In ADA we introduce overlap between these observations. The overlap sample length $s_{\gamma }$ is defined by a proportion $\gamma$ of the proposed window sample length $s_w$, i.e.,

$$s_{\gamma }=⌊\gamma s_w⌋.$$

(5)

where $s_w$ is the window sample length given by the desired window length in seconds $T_w$ multiplied by the sampling frequency $f_s$

$$s_w=⌊f_s{T}_w⌋.$$

(6)

The stride of the window, $s_{\tau }$, is then given by

$$s_{\tau }=s_w-s_{\gamma }.$$

(7)

This occurs for each time sequence for either an actuator or sensor signal, being appended column-wise, resulting in

$$\begin{array}{c}\hfill p=s_w\times o,\end{array}$$

(8)

$$\begin{array}{c}\hfill r=s_w\times q\end{array}$$

(9)

and the number of rows or observations, n, equal to

$$n=⌊\frac{m-s_{\gamma }}{s_{\tau }}⌋.$$

(10)

We set the amount of overlap to the extreme such that the stride is one sample, i.e., $s_{\tau }=1$. This results in the following windowed target matrix Y

(11)

The windowed predictor matrix X takes on a similar form (not shown). We can simply average over the anti-diagonals of the windowed data matrix $\widehat{Y}$ to reconstruct the approximated input $\widehat{U}$. To compute the average response $\widehat{u}\left(k\right)$ we average all the anti-diagonal terms of ${\widehat{Y}}_{i,j}$, such that

$$\widehat{u}\left(k\right)=\frac{1}{n_{\mathrm{diag}}}\sum {\widehat{Y}}_{i,j},$$

(12)

for which $i+j=k+1$ and $n_{\mathrm{diag}}$ is the number of elements in the anti-diagonal. This process is known as Hankelization which is the same process followed in Singular Spectral Analysis (SSA) [14]. The corresponding windowed matrix is referred to as the trajectory matrix. This is known as the embedding step in SSA. The windowed matrix is then decomposed using SVD. In this case, we are merely borrowing the ADA concept from SSA for the regression problem, whereas SSA typically uses this process for autoregressive models.

An example signal with $m=7$ samples and window length $s_w=3$ and windowed with ADA results in the following equation

$$\left[\begin{array}{ccc}\widehat{u}\left(1\right)& \widehat{u}\left(2\right)& \widehat{u}\left(3\right)\\ \widehat{u}\left(2\right)& \widehat{u}\left(3\right)& \widehat{u}\left(4\right)\\ \widehat{u}\left(3\right)& \widehat{u}\left(4\right)& \widehat{u}\left(5\right)\\ \widehat{u}\left(4\right)& \widehat{u}\left(5\right)& \widehat{u}\left(6\right)\\ \widehat{u}\left(5\right)& \widehat{u}\left(6\right)& \widehat{u}\left(7\right)\end{array}\right]=\left[\begin{array}{ccc}z\left(1\right)& z\left(2\right)& z\left(3\right)\\ z\left(2\right)& z\left(3\right)& z\left(4\right)\\ z\left(3\right)& z\left(4\right)& z\left(5\right)\\ z\left(4\right)& z\left(5\right)& z\left(6\right)\\ z\left(5\right)& z\left(6\right)& z\left(7\right)\end{array}\right]\left[\begin{array}{ccc}{\beta }_{1,1}& {\beta }_{1,2}& {\beta }_{1,3}\\ {\beta }_{2,1}& {\beta }_{2,2}& {\beta }_{2,3}\\ {\beta }_{3,1}& {\beta }_{3,2}& {\beta }_{3,3}\end{array}\right].$$

(13)

Here z is the response signal used to predict the inputs u. The linear coefficients $\beta$ are computed using any suitable linear regression method. To gain insight into the workings of ADA we can write out the set of equations that infer $\widehat{u}\left(3\right)$, i.e.,

$$\begin{array}{ccc}\hfill \widehat{u}{\left(3\right)}_1& =& {\beta }_{1,3}z\left(1\right)+{\beta }_{2,3}z\left(2\right)+{\beta }_{3,3}z\left(3\right),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill \widehat{u}{\left(3\right)}_2& =& {\beta }_{1,2}z\left(2\right)+{\beta }_{2,2}z\left(3\right)+{\beta }_{3,2}z\left(4\right),\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \widehat{u}{\left(3\right)}_3& =& {\beta }_{1,1}z\left(3\right)+{\beta }_{2,1}z\left(4\right)+{\beta }_{3,1}z\left(5\right).\hfill \end{array}$$

(16)

We can then average over all the $\widehat{u}\left(3\right)$ predictions to obtain the final prediction of $\widehat{u}\left(3\right)$

$$\begin{array}{ccc}\hfill \widehat{u}\left(3\right)& =& \frac{1}{3}\left(\widehat{u}{\left(3\right)}_1+\widehat{u}{\left(3\right)}_2+\widehat{u}{\left(3\right)}_3\right)\hfill \\ & =& z\left(1\right)\left(\frac{{\beta }_{1,3}}{3}\right)+z\left(2\right)\left(\frac{{\beta }_{2,3}+{\beta }_{1,2}}{3}\right)+z\left(3\right)\left(\frac{{\beta }_{3,3}+{\beta }_{2,2}+{\beta }_{1,1}}{3}\right)\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& & +z\left(4\right)\left(\frac{{\beta }_{3,2}+{\beta }_{2,1}}{3}\right)+z\left(5\right)\left(\frac{{\beta }_{3,1}}{3}\right).\hfill \end{array}$$

(18)

If we rewrite the average of the $\beta$ multiplying with a particular z term as a new constant, e.g., ${\beta }_2=\frac{{\beta }_{2,3}+{\beta }_{1,2}}{2}$, we obtain

$$\widehat{u}\left(3\right)={\beta }_1\frac{1}{3}z\left(1\right)+{\beta }_2\frac{2}{3}z\left(2\right)+{\beta }_{3}1z\left(3\right)+{\beta }_4\frac{2}{3}z\left(4\right)+{\beta }_5\frac{1}{3}z\left(5\right).$$

(19)

Here we note that the ADA emphasises the middle term with decreasing emphasis placed on proceeding and preceding terms. It in effect creates a triangular windowing function. If we add a corresponding weight term w, e.g., $w_2=\frac{2}{3}$ we can rewrite the equation generally as

$$\widehat{u}\left(k\right)={\beta }_1{w}_{1}z(k-s_w)+\cdots +{\beta }_{s_w}{w}_{s_w}z\left(k\right)+\cdots +{\beta }_{k+s_w}{w}_{k+s_w}z(k+s_w).$$

(20)

This result demonstrates that ADA is an indirect method of creating a weighted moving average filter. In system identification this is known as a Finite Impulse Response (FIR) model. More specifically this an example of a non-causal weighted FIR model. The weights can be arbitrary and are a prior design choice. If we forgo the ADA method and use the weighted FIR model, we can be more creative with the weighting.

1.3. FIR Models

In FIR models the current output of the system is a function of past inputs such that

$$z\left(k\right)=f(u(k-1),\dots ,u(k-s_w)).$$

(21)

This is in contrast to other models such as Autoregressive eXogenous (ARX) which includes output feedback as well, i.e.,

$$z\left(k\right)=f(u(k-1),\dots ,u(k-s_w),z(k-1),\dots ,z(k-s_w)).$$

(22)

This paper focuses on non-causal inverse implementations of FIR models, where the current input is a function of both past and future outputs, written as

$$u\left(k\right)=f(z(k-s_w/2),\dots ,z(k+s_w/2)).$$

(23)

By using the FIR model, the predictor matrix X takes on the form

(24)

with the corresponding target matrix Y written as

$$Y=\left[\begin{array}{cccc}{u}_1(s_w/2)& \cdots & u_{q-1}(s_w/2)& u_q(s_w/2)\\ u_1(s_w/2+1)& \cdots & u_{q-1}(s_w/2+1)& u_q(s_w/2+1)\\ ⋮& \ddots & ⋮& ⋮\\ u_1(m-s_w/2-1)& \cdots & u_{q-1}(m-s_w/2-1)& u_q(m-s_w/2-1)\\ u_1(m-s_w/2)& \cdots & u_{q-1}(m-s_w/2)& u_q(m-s_w/2)\end{array}\right].$$

(25)

It is worth noting that we lose the first and last $s_w/2$ samples of the target matrix Y since we shifted the inputs to make the system non-causal.

The lack of feedback means that FIR methods are inherently stable. This is suitable and sometimes sought after if the system under consideration is stable. However, if the system is unstable, it will only approximate the instability for a short period before diverging [15]. FIR models come with the cost of needing significantly more terms than what output feedback models need to map the same system [15]. A similar approach to ADA can be achieved through the use of FIR models combined with Tikhonov regularisation. Using Tikhonov regularisation, the $\beta$ coefficients can be penalised and thus shaped by choice of the $\mathsf{\Gamma }$ matrix in Equation (2). To this end three options for the $\mathsf{\Gamma }$ matrix are implemented in this paper, namely: Finite Impulse Response with TriangularWeighting (FIR-T), Finite Impulse Response with Difference Smoothing and Triangular Weighting (FIR-DT) and Finite Impulse Response with Ridge Regression (FIR-RR).

In FIR-T, the coefficients relating to the outputs further away from the required input (both forwards and backwards in time) are penalised. This is achieved by setting

$${\mathsf{\Gamma }}^{\prime }\mathsf{\Gamma }=\alpha W$$

(26)

where W is an inverted triangular set of penalty weights, given as

$$W=diag\left(\left[\begin{array}{ccccccccc}{s}_w/2& s_w/2-1& \cdots & 2& 1& 2& \cdots & s_w/2-1& s_w/2\end{array}\right]\right),$$

(27)

and $\alpha$ scales the amount of regularisation we wish to impose. This should ideally mimic the weighting function achieved by ADA in Equation (20). FIR-DT further modifies the triangular weighting matrix through the use of a first difference matrix A, given as

$$A=\left[\begin{array}{ccccc}1& 0& \cdots & 0& 0\\ -1& 1& \cdots & 0& 0\\ 0& -1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& 0& \cdots & -1& 1\end{array}\right].$$

(28)

The first difference matrix ensures that the difference between each successive $\beta$ coefficient is small [16]. The difference matrix is then combined with the weighting matrix W to obtain the final form of the regularisation matrix such that

$${\mathsf{\Gamma }}^{\prime }\mathsf{\Gamma }=\alpha A^{\prime }WA.$$

(29)

This weighting scheme was initially implemented and developed for a causal FIR system where the penalty weights increased linearly further back in time [16]. Finally, the last choice of penalty matrix $\mathsf{\Gamma }$ is that of FIR-RR, i.e.,

$$\mathsf{\Gamma }=\alpha I$$

(30)

where we only limit the magnitude of the weights to act as a reference. This enables us to determine whether the regularisation of $\beta$ contributes to the accuracy of the response reconstruction.

2. Method

This section describes the general experimental design procedure for the numerical investigations.

2.1. Numerical Quarter Car Model

A simple two-degree-of-freedom nonlinear mass–spring–damper system, representing a quarter car model is used to investigate the methods explored in this paper. The numerical model employed is shown schematically in Figure 3. The sprung mass $M_A$ and unsprung mass $M_R$ represent the mass of the vehicle’s body and the suspension–tyre system, respectively. These bodies are connected by springs and dampers, which represent the dynamics of the suspension system. The unsprung mass is then connected to the road via a spring characterising the tyre stiffness. The system is excited by a road profile $u_{\mathrm{road}}$.

The system behaves according to the following equations of motion:

$$\begin{array}{ccc}\hfill {\ddot{z}}_A& =& -\frac{b_A}{M_A}({\dot{z}}_A-{\dot{z}}_R)-\frac{k_{NL}}{M_A}{(z_A-z_R)}^3-\frac{k_A}{M_A}(z_A-z_R),\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill {\ddot{z}}_R& =& +\frac{b_A}{M_R}({\dot{z}}_A-{\dot{z}}_R)+\frac{k_{NL}}{M_R}{(z_A-z_R)}^3+\frac{k_A}{M_R}(z_A-z_R)-\frac{k_R}{M_R}(z_R-u_{road}),\hfill \end{array}$$

(32)

where the k and b terms are the stiffness and damping coefficients, respectively. The nonlinearity is introduced by having cubic stiffening of the sprung mass spring controlled by the $k_{NL}$ term. The sprung mass spring force, $f_A$, is given by

$$\begin{array}{c}\hfill f_A\left({\mathsf{\Delta }}_z\right)=-k_A{\mathsf{\Delta }}_z-k_{NL}{\mathsf{\Delta }}_z^3,\end{array}$$

(33)

where we define a new state of the system representing the deflection of the spring, ${\mathsf{\Delta }}_z$, such that

$${\mathsf{\Delta }}_z\equiv z_A-z_R.$$

(34)

The $k_{NL}$ term can be varied to change the severity of the system’s nonlinearity or switch it completely off for linear behaviour. A hardening spring is modelled by choosing $k_{NL}>0$. This results in a spring that becomes stiffer as it undergoes compression or tension. Likewise, a softening spring can be implemented by choosing $k_{NL}<0$. In this study the linear component will always be restorative such that $k_A>0$. The default parameters chosen for the numerical quarter car are given in

Table 1. Numerical quarter car model’s default parameters.

${\mathit{M}}_{\mathit{A}}\left(\mathbf{k}\mathbf{g}\right)$	${\mathit{M}}_{\mathit{R}}\left(\mathbf{k}\mathbf{g}\right)$	${\mathit{k}}_{\mathit{A}}\left(\mathbf{N}{\mathbf{m}}^{-1}\right)$	${\mathit{k}}_{\mathit{R}}\left(\mathbf{N}{\mathbf{m}}^{-1}\right)$	${\mathit{k}}_{\mathit{NL}}\left(\mathbf{N}{\mathbf{m}}^{-3}\right)$	${\mathit{b}}_{\mathit{A}}\left(\mathbf{N}\mathbf{s}{\mathbf{m}}^{-1}\right)$
70	12	$1.6\times {10}^3$	$80\times {10}^3$	$12.8\times {10}^6$	500

.

2.2. Choice of Excitation Signals

Before we can begin building a direct inverse model of our plant we need informative data since the excitation signal’s quality places an upper bound on the accuracy of any subsequent model that we wish to build [15]. For response reconstruction, we can design the signals on which we want to train. There are two possible methods of designing excitation signals: model-free and model-based methods. In model-based methods subsequent excitation signals are chosen to improve the accuracy of the model [17]. In model-free methods we design an excitation signal that offers the best distributed coverage of the operating condition. Initially we have little prior knowledge of the system and of the real world input signals; therefore, we need to employ model-free methods. We assume we have some prior knowledge of the range of the operating condition. A suitable choice is the Amplitude Modulated Pseudo Random Binary Signal (APRBS).

2.2.1. Amplitude Modulated Pseudo Random Binary Signal

Since we are working with nonlinear structural systems, we know that the system responses are functions of input frequencies and the amplitude at which we excite the system. Therefore a signal that covers the necessary frequencies and the expected amplitude range of operations is required. The APRBS attempts to cover the amplitude operating conditions with a series of step responses that are fairly well distributed over the input range. An example of an APRBS is shown in Figure 4.

To specify the profile, a set of N design points $d_n$ are chosen to define the step’s amplitude. The design points are sampled from the desired range $[u_{\min },u_{\max }]$ using Latin Hypercube Sampling (LHS). LHS splits the design space into N intervals with one design point placed randomly in each interval. LHS then iteratively optimises the design points such that each design point is the maximum distance away from its neighbours. This provides a random but equally spread set of design points. Since no physical system can achieve an instantaneous change in displacements required for a true step input, the step is instead approximated by a ramp function. The slope of the ramp is determined by the maximum allowed velocity $v_{\max }$ that can safely or accurately be performed by the actuator. The ramp’s slope affects the frequency content of the signals with higher velocities resulting in higher frequencies being excited [18]. The length of the step is then specified by the hold time $T_h$. Since the testing time is limited, the maximum number of steps that best cover the input space in the shortest time is sought. The hold time $T_h$ must be small enough to fit as many steps in but must be long enough that the steps actively excite the system at that point. The hold time $T_h$ is typically set to be at least the length of the largest time constant $T_{c,\max }$ of the system [15]. This can be determined with a simple step test of the system if no prior knowledge is known. The parameters used for the investigations are given in

Table 2. APRBS parameters used to generate the training and validation signals used in the non-overlapping windows numerical investigation.

${\mathit{T}}_{\mathit{s}}\left(\mathit{s}\right)$	${\mathit{T}}_{\mathit{h}}\left(\mathit{s}\right)$	${\mathit{v}}_{\mathbf{max}}\left(\mathit{m}{\mathit{s}}^{-1}\right)$	${\mathit{u}}_{\mathbf{max}}\left(\mathit{m}\right)$	${\mathit{u}}_{\mathbf{max}}\left(\mathit{m}\right)$
$0.001$	$0.2$	10	$-0.1$	$0.1$

.

2.2.2. Road Profile

The ISO 8608 standard [19] for specifying road profiles is used to generate a separate test set to determine how well the direct inverse model performs on unseen data. The ISO 8608 standard defines inputs that are distinct from APRBS while still being representative of real-world operating conditions. The profiles are characterised by the standard in the frequency domain where the spectral density $S_z$ is given by

$$S_z\left(\varphi \right)=A{\left(\varphi \right)}^{-n},$$

(35)

for the given spatial frequency $\varphi$ with units $m^{-1}$. The A term represents the road’s roughness coefficient, whereas n represents the road index of the profile. The A coefficient controls how large the amplitudes are at each frequency whereas n controls how quickly the amplitudes decay as frequency functions. Varying types of profiles such as ploughed agricultural land to smooth gravel highways can be produced by altering these two coefficients. The spatial frequencies $\varphi$ are limited between $0.5$ and 10 $m^{-1}$. The former represents the broad changes in the landscape which have negligible effects on vehicle dynamics. In contrast, the upper limit on the frequency represents small variations which are filtered out by the tyre [20]. When generating the profiles only the amplitude information is given by the ISO 8608 standard; therefore, in order to generate time signals, a uniformly random signal is generated for the phase signal with spatial frequencies sampled at discrete intervals. This generates a displacement signal as a function of distance. The vehicle’s velocity must be chosen to generate a displacement signal as a function of time. The parameters of the road profile used are given in

Table 3. Road profile parameters used to generate the test signal used in the non-overlapping windows numerical investigation.

n	A	${\mathit{\varphi }}_{\mathit{min}}\left({\mathit{m}}^{-1}\right)$	${\mathit{\varphi }}_{\mathit{max}}\left({\mathit{m}}^{-1}\right)$	${\mathit{\varphi }}_{\mathit{int}}\left({\mathit{m}}^{-1}\right)$	$\mathit{v}\left(\mathit{m}{\mathit{s}}^{-1}\right)$
10	$6.5\times {10}^{-4}$	$0.5$	10	$3.5\times {10}^{-4}$	5

.

2.2.3. Preprocessing

The windowing techniques covered in this paper will truncate some of the testing and training set samples. To ensure a fair comparison between the different data sets, a dead time is appended and prepended. The dead times will be excluded when calculating the cost function during cross validation and reporting the final accuracy of the predictions. The constant initial and final conditions also allow for different signals to be concatenated without introducing unwanted jumps.

3. Scaling

The windowed inputs Y and windowed outputs X of the system record different types of signals which will have different variances across them. We may also find that constant biases need to be accounted for from the sensors. Therefore, the inputs and outputs are z-scored normalised to scale the rows to have a mean of zero and a variance of one [13].

4. Cross Validation

To determine the optimal regularisation constant $\alpha$ for RR, cross validation is used. However, cross validation can be misleading if it is implemented without considering the correlation between observations. Suppose the validation set is removed once the data have already been windowed with overlaps. In that case, the validated set will be correlated to the training set due to the overlaps introduced. If the validation set is first removed from the middle portion of the dataset and then windowed, then care must be taken when splitting and merging the training set to ensure that no unintended overlap is introduced between the separated training segments. A simpler solution to this problem is implemented by removing a single validation set from either the beginning or end of the dataset before windowing. In this work, a validation set was created independently of the training set.

4.1. Choice of Cost Function

We have the choice of either using the errors of the approximated inputs or the approximated outputs as the cost function of the optimisation scheme. In response reconstruction, we are interested in producing an accurate output response since a unique input may not exist. The downside of this is that, to obtain the output error, the approximated input needs to be passed through the test rig. This needs to occur for every loop in the cross validation step. The numerical model is computationally efficient to compute. However, this would result in significant fatigue of the experimental rig in the real world and would take considerable time to run. Therefore, it is necessary to limit the number of forward evaluations in the cross validation step. In evaluating these methods for response reconstruction, the output error is used during cross validation. Since we need to measure and compare response and input reconstruction accuracies across different types of signals, we need a normalised measure of error. The Mean Fit Function Error (MFFE) [21] is used to report the final test accuracies of the reconstructed input and output signals. MFFE is defined as

$$\mathrm{MFFE}=100\times \frac{{\sum }_{m=1}^M\left|e_0\right|}{{\sum }_{m=1}^M\left|z_0\right|}[\%],$$

(36)

where $e_0$ is the error between the true output $z_0$ and the approximate output ${\widehat{z}}_0$, i.e.,

$$e_0=z_0-{\widehat{z}}_0.$$

(37)

The signals under consideration have been mean centred such that

$$\begin{array}{ccc}\hfill z_0& =& z-{\mu }_z,\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill {\widehat{z}}_0& =& \widehat{z}-{\mu }_{\widehat{z}}.\hfill \end{array}$$

(39)

4.2. Training Procedure

The cross validation algorithm consists of two sub-routines: an outer routine that incorporates the windowing parameter grid search for the optimal window length $T_w$ and an inner subroutine which optimises the regularisation constant $\alpha$.

4.2.1. Window Loop

A graphical overview of the training process is shown in Figure 5 with focus on the window parameter search. The window optimisation loops over the window length, $T_{w,i}$, where i represents the ith iteration of the loop. The training set $U_{\mathrm{train}}$ and $Z_{\mathrm{train}}$ as well as validation output $Z_{\mathrm{val}}$ are then windowed accordingly. The z-score parameters, ${\sigma }_i$ and ${\mu }_i$, are then calculated using only the training dataset and applied to both the training and validation set. The training set is then decomposed using SVD according to the regression method specified, in this case RR. The decomposed SVD is then passed to the regularisation optimisation loop.

4.2.2. Latent Variable Loop

Figure 6 depicts a graphical overview of the regularisation constant optimisation. The regression coefficients ${\beta }_k$ are then calculated and weighted with ${\alpha }_k$, where k is the $k^{th}$ iteration of the loop. The approximated windowed validation inputs ${\widehat{Y}}_{\mathrm{val}}$ are then predicted using the windowed validation outputs ${\widehat{X}}_{\mathrm{val}}$. The approximated windowed validation inputs are rescaled and then merged using the specified windowing methods to obtain the approximated input ${\widehat{U}}_{\mathrm{val}}$. The merged inputs are then passed through the test rig to obtain the approximated output ${\widehat{Z}}_{\mathrm{val}}$. The MFFE is then calculated between the true output $Z_{\mathrm{val}}$ and the approximated output ${\widehat{Z}}_{\mathrm{val}}$. The optimised regularisation constant ${\alpha }_k$ and the corresponding minimum MFFE are then returned from this loop to the windowing loop as seen in Figure 5. This minimum MFFE result is then used in the window loop to find the corresponding optimal window length $T_{w,\min }$.

4.2.3. Final Training Step

In the final training step, the training set is concatenated with the validation set. This newly combined set is then windowed with the optimised window parameters $T_{w,\min }$. The new z-score parameters $[\sigma ,\mu ]$ are then calculated. The combined set is decomposed and used in the regression step with the optimised regularisation constant ${\alpha }_{\min }$ to determine the final regression coefficients ${\beta }_{\mathrm{final}}$.

4.3. Prediction

A graphical overview of the prediction step and the approximation of the output is shown in Figure 7. Once the training step is complete, it is relatively straightforward to use the optimised parameters to make further predictions. The test output signal $Z_{\mathrm{test}}$ needs to be preprocessed first before predictions can be made. To obtain the predictor matrix, $X_{\mathrm{test}}$, the test signal is windowed and z-scored normalised using the parameters determined during the training phase. The prediction step then occurs using the regression coefficients ${\beta }_{\mathrm{final}}$ obtained during training to obtain the approximate target matrix, ${\widehat{Y}}_{\mathrm{test}}$. The windowing and z-score normalisation are then reversed before passing the approximated input ${\widehat{U}}_{\mathrm{test}}$ into the test rig to obtain the approximated output, ${\widehat{Z}}_{\mathrm{test}}$.

5. Comparison against Finite Impulse Response (FIR) Models

This section aims to benchmark ADA against FIR in terms of response reconstruction since ADA can be seen as a subset of FIR. The idea behind this benchmark is to ensure that ADA is not an indirect method of achieving an FIR implementation. If so, it needs to be determined whether ADA offers any substantial benefits over using FIR.

5.1. Finite Impulse Response (FIR) Comparison Procedure

The three different regularised FIR implementations will be compared against ADA combined with RR for two different test cases, linear and nonlinear. The experiment will be performed with a system configuration more representative of a typical test rig. In this case the sprung mass acceleration and spring displacement (i.e., the delta between the sprung and unsprung mass) of the quarter car will be used. The first being the linear system and the second being the default nonlinear system. The inputs and responses will be sampled at 250 $\mathrm{Hz}$ and 350 $\mathrm{Hz}$ for the linear and nonlinear systems, respectively. The window lengths will be determined via grid search cross validation with the window lengths being sampled from $T_w\in [0.1,12]$ $\mathrm{s}$ with a grid of 50 equally spaced intervals. The potential $\alpha$ values used to regularise RR will be spaced equally on a log scale within the range $\alpha \in [s_{\min }\times {10}^{-5},s_{\max }]$, where s are the singular values. Thirty equally spaced divisions will be used. An overview of the numerical experiment parameters is given in

Table 4. Experimental design benchmarking ADA against different forms of FIR models. (Variables of interest shown first).

Variable	Details
Reg. method	FIR-T, FIR-DT, FIR-RR and ADA-RR
${\mathbf{k}}_{\mathbf{NL}}$	Linear: 0, Nonlinear: $1.28\times {10}^7$ $\mathrm{N}/{\mathrm{m}}^3$
Sensor config.	Sprung mass acceleration + spring displacement
Window length ${\mathbf{T}}_{\mathbf{w}}$	$\in [0.1,12]$ $\mathrm{s}$ with a grid of 50 equally spaced intervals
Sampling frequency ${\mathbf{f}}_{\mathbf{s}}$	Linear: 250 $\mathrm{Hz}$, Nonlinear: 350 $\mathrm{Hz}$
Window proportional overlap $\gamma$	Maximum
Ridge regression regularisation constant$\alpha$	$\in [{10}^{-16},{10}^5]$ with 30 divisions spaced logarithmically
QC parameters	Default values; Table 1
Noise level ${\eta }_{\%}$	0 %
Train. Set	APRBS; Table 2
Val. Set	APRBS; Table 2
Test Set	Road profile; Table 3

.

5.2. FIR Comparison Numerical Results

The reconstructed inputs and outputs for the linear and nonlinear systems are shown in Figure 8 and, Figure 9 respectively. The response reconstruction results for the linear and nonlinear systems are shown in

Table 5. MFFE scores for the approximated input and output signals using different FIR methods. Best performing results shown in bold.

	Training			Validation			Test
	${\mathbf{u}}_{\mathbf{road}}$	${\ddot{\mathbf{z}}}_{\mathbf{A}}$	${\mathsf{\Delta }}_{\mathbf{z}}$	${\mathbf{u}}_{\mathbf{road}}$	${\ddot{\mathbf{z}}}_{\mathbf{A}}$	${\mathsf{\Delta }}_{\mathbf{z}}$	${\mathbf{u}}_{\mathbf{road}}$	${\ddot{\mathbf{z}}}_{\mathbf{A}}$	${\mathsf{\Delta }}_{\mathbf{z}}$
Linear Case
FIR-RR	12.40	1.26	0.26	12.71	1.63	0.35	$\mathbf{3}.\mathbf{50}$	0.91	1.13
FIR-T	12.33	1.38	0.34	12.41	0.46	0.25	3.80	0.90	1.15
FIR-DT	$\mathbf{10}.\mathbf{57}$	4.39	2.01	17.98	1.48	1.81	6.98	3.08	3.16
ADA-RR	11.14	$\mathbf{0}.\mathbf{71}$	$\mathbf{0}.\mathbf{16}$	$\mathbf{10}.\mathbf{48}$	$\mathbf{0}.\mathbf{38}$	$\mathbf{0}.\mathbf{23}$	5.49	$\mathbf{0}.\mathbf{16}$	$\mathbf{0}.\mathbf{35}$
Nonlinear Case
FIR-RR	25.66	11.33	1.99	43.72	5.03	2.82	20.76	8.13	7.25
FIR-T	38.86	5.92	4.32	59.83	2.34	1.60	38.90	17.56	13.44
FIR-DT	$\mathbf{5}.\mathbf{22}$	8.54	4.70	$\mathbf{35}.\mathbf{94}$	4.71	3.60	$\mathbf{7}.\mathbf{50}$	5.41	4.89
ADA-RR	5.91	$\mathbf{0}.\mathbf{55}$	$\mathbf{0}.\mathbf{13}$	36.60	$\mathbf{0}.\mathbf{21}$	$\mathbf{0}.\mathbf{16}$	12.55	$\mathbf{0}.\mathbf{39}$	$\mathbf{0}.\mathbf{40}$

. We treat FIR-RR as the bare minimum regression method since it does not impose a prior choice on the shape or smoothness of the $\beta$ parameters.

For the linear case it appears that the imposed smoothness offered by FIR-DT does not contribute any significant improvement and actually hinders the reconstruction performance. If we refer to the optimised hyper-parameters for the numerical experiment in

Table 6. Optimised hyper-parameter results for the numerical demonstration using different FIR methods and ADA-RR.

	$\mathit{\alpha }$	${\mathbf{T}}_{\mathbf{w}}\left(\mathbf{s}\right)$	$\mathit{\alpha }$	${\mathbf{T}}_{\mathbf{w}}\left(\mathbf{s}\right)$
	Linear Case		Nonlinear Case
FIR-RR	11.80	12.00	4.27	3.50
FIR-T	$4.74\times {10}^{-2}$	11.51	$4.80\times {10}^{-1}$	4.47
FIR-DT	$1.04\times {10}^{-9}$	8.11	$4.79\times {10}^{-2}$	7.63
ADA-RR	$8.52\times {10}^{-3}$	10.54	$4.66\times {10}^{-9}$	6.65

, we see that FIR-DT used a small amount of regularisation which further indicates the poor suitability of the methodology to the problem. We note that the triangular weighting offered by FIR-T performs similarly to FIR-RR, which suggests that the shape of the $\beta$ parameters, is not as important for the linear case. However, ADA still performs an order of magnitude better in terms of the recreated output MFFE scores. This suggests, at least for the linear case, that the ADA performance is not necessarily due to the shape factor or due to imposed smoothing of each successive $\beta$ parameter.

For the nonlinear case, we note that FIR-T, obtains the worst recreated output score with the default regression method, FIR-RR, performing significantly better. This suggests that the introduction of the triangular weighting is ill-suited for the nonlinear case. The introduction of the difference smoothing in the form of FIR-DT is an improvement over FIR-RR, which suggests that the smoothing of the $\beta$ parameters is an improvement in the FIR regression methods’ performance. If we refer to the optimised hyper-parameters in

Table 7. Optimised hyper-parameter results for the numerical demonstration using different FIR methods for an illustrative use case.

	$\mathit{\alpha }$	${\mathit{T}}_{\mathit{w}}\left(\mathit{s}\right)$
FIR-RR	0.64	8.03
FIR-T	7.23	2.08
FIR-DT	469.42	9.52
ADA-RR	$6.81\times {10}^{-5}$	7.04

, we note that ADA implemented a low amount of regularisation for the nonlinear case. This indicates that averaging used in ADA adds an extra form of regularisation since it performs an order of magnitude better than the other regression methods for the nonlinear case, without relying on a large regularisation constant. This is corroborated by the fact that ADA outperforms the other regression methods when either smoothing is better suited (nonlinear case) or imposing a shape is better suited (linear case), which suggests that the averaging inherent to ADA is the key factor for its performance regarding the problem at hand.

In general, we note that the MFFE results for the recreated outputs of the system (for both the linear and nonlinear case) are lower than their associated input MFFE results. This indicates the non-uniqueness of the inputs for the given response since a seemingly poor input can result in an accurate output. This justifies the need to incorporate the forward pass through the system to determine the suitability of the input by judging it by its associated recreated output.

6. Illustrative Use Case

In this section, we create a scenario whereby all the challenges to response reconstruction are introduced. These are noise, model mismatch and nonlinearity. In this experiment we focus on a narrower scope of model mismatch whereby the model parameters of the system are simply scaled from the real-world environment to that of the laboratory environment. A broader scope of model mismatch would be to add new dynamics going from one environment to the other. One such example would be to add or remove a discontinuity, i.e., a tyre separating from the road surface. This paper focuses on this narrower view of model mismatch. The default parameters given in Table 1 are modified such that

$$\begin{array}{ccc}\hfill M_{A,\mathrm{mis}}& =& M_A\left(1-\frac{m_{\%}}{100}\right)\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill M_{R,\mathrm{mis}}& =& M_R\left(1+\frac{m_{\%}}{100}\right)\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill b_{A,\mathrm{mis}}& =& b_A\left(1-\frac{m_{\%}}{100}\right)\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill k_{A,\mathrm{mis}}& =& k_A\left(1+\frac{m_{\%}}{100}\right)\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill k_{R,\mathrm{mis}}& =& k_R\left(1-\frac{m_{\%}}{100}\right).\hfill \end{array}$$

(44)

The investigation is not exhaustive but rather proposed to give an illustrative sense of the regression methods’ performance on a challenging response reconstruction problem. To this end, the numerical experiment will be performed with the FIR and ADA regression methods with noise, model-mismatch and nonlinearity implemented. This investigation’s level of noise is defined in percentage terms, ${\eta }_{\%}$, of the standard deviation for each channel o of the outputs z. The noise is assumed to be Gaussian with zero mean, resulting in

$${\mathbf{z}}_{o,\mathrm{noisy}}={\mathbf{z}}_o+\mathcal{N}\left(0,{\eta }_{\%}^2{\sigma }_{z_o}^2\right).$$

(45)

Noise ${\eta }_{\%}$ will be set to 5 %, model mismatch $m_{\%}$ set to $10\%$ and the non-linearity term $k_{NL}$ set to $1.28\times {10}^7$ $\mathrm{N}/{\mathrm{m}}^3$. In the case of model mismatch, the validation response set will come from a field recording instead of a laboratory recording. The idea behind this is to force the cross validation to only retain latent variables that allow the laboratory environment to recreate dynamics that are common to both the real-world and the lab environment. An overview of the numerical procedure is given in

Table 8. Experimental design comparing ADA against different FIR models for an illustrative use case. (Variables of interest shown first).

Variable	Details
Regression method	FIR-T, FIR-DT, FIR-RR and ADA-RR
Nonlinearity constant${\mathbf{k}}_{\mathbf{NL}}$	Nonlinear: $1.28\times {10}^7$ $\mathrm{N}/{\mathrm{m}}^3$
Model mismatch	Perturbed as per Equations (40) to (44) with $m_{\%}=10$
Noise level ${\eta }_{\%}$	5%
Sensor configuration	Sprung mass acceleration + spring displacement
Window lengths ${\mathbf{T}}_{\mathbf{w}}$	$\in [0.1,12]$ $\mathrm{s}$ with a grid of 50 equally spaced intervals
Sampling frequency ${\mathbf{f}}_{\mathbf{s}}$	350 $\mathrm{Hz}$
Window proportional overlap $\gamma$	Maximum
Ridge regression regularisation constant $\alpha$	$\in [{10}^{-16},{10}^5]$ with 30 divisions spaced logarithmically
QC parameters	Default values; Table 1
Training Set	APRBS; Table 2
Validation Set	APRBS; Table 2
Test Set	Road profile; Table 3

.

Illustrative Use Case Numerical Results

The response reconstruction results are shown in

Table 9. MFFE $(\%)$ scores for the approximated input and output signals using different FIR methods for an illustrative use case. Best performing results shown in bold.

	Training			Validation			Test
	${\mathbf{u}}_{\mathbf{road}}$	${\ddot{\mathbf{z}}}_{\mathbf{A}}$	${\mathsf{\Delta }}_{\mathbf{z}}$	${\mathbf{u}}_{\mathbf{road}}$	${\ddot{\mathbf{z}}}_{\mathbf{A}}$	${\mathsf{\Delta }}_{\mathbf{z}}$	${\mathbf{u}}_{\mathbf{road}}$	${\ddot{\mathbf{z}}}_{\mathbf{A}}$	${\mathsf{\Delta }}_{\mathbf{z}}$
FIR-RR	15.37	8.78	6.26	46.52	16.26	12.66	24.45	31.56	25.17
FIR-T	56.34	16.12	6.71	62.47	$\mathbf{13}.\mathbf{79}$	11.62	71.46	46.49	35.09
FIR-DT	25.72	$\mathbf{6}.\mathbf{59}$	11.84	46.65	14.91	$\mathbf{11}.\mathbf{58}$	$\mathbf{22}.\mathbf{20}$	15.53	14.41
ADA-RR	$\mathbf{18}.\mathbf{82}$	6.64	$\mathbf{4}.\mathbf{68}$	$\mathbf{29}.\mathbf{90}$	16.01	13.85	24.69	$\mathbf{14}.\mathbf{75}$	$\mathbf{14}.\mathbf{23}$

with the corresponding reconstructed inputs and outputs shown in Figure 10. By referring to the results in Table 9, we see that ADA and FIR-DT perform similarly well for the reconstructed test results. These results are achieved within a close enough margin to each other that it probably falls within the uncertainty introduced by noise. We see that FIR-T performs poorly for the problem at hand. This follows the general trend of reconstruction performance as found in the previous nonlinear benchmark. The optimised hyper-parameters for this numerical experiment are shown in Table 7. Here we note that the different regression methods use similar window lengths, save for FIR-T, which used a significantly shorter window length. Here we also note that the regularisation constants $\alpha$ are larger than those found in Table 6, which is to be expected since more regularisation is needed due the introduced model mismatch as well as the added output noise.

7. Conclusions

By introducing the overlapping windows inherent to the ADA implementation as well as focusing on the recreated outputs of the system, we overcome the asymmetry introduced by the inverse nature of response reconstruction. In summary, it is shown that ADA combined with an appropriate linear regression is a suitable black-box method of reconstructing responses in dynamic systems. It has wide application in response reconstruction in that it can be readily applied to practical sensor configurations as well as non-linear systems. We compared the performance of ADA to related FIR regression methods. Although the experiments were not exhaustive, the results indicate that ADA outperforms the related FIR methods in response reconstruction accuracy. By repeating the experiment with challenges that require better regularisation, insights into how ADA may be performing regularisation was gained. The current ADA implementation can be seen as a post-processing smoothing step that occurs after a linear regression prediction. An exciting avenue to explore would be to replace the linear regression with a non-linear regression method such as a neural network.