Geometry-Informed Optimal Parameter Selection for Phase Space Warping with Application to Early Fatigue Damage Detection in 3D-Printed Beams

Phase space warping (PSW) methodology reconstructs a non-stationary hidden process from quasi-stationary observ- 8 able dynamics, where these two coupled dynamical processes have disparate time scales. PSW has been applied to 9 multivariate damage identiﬁcation and tracking, bio-mechanics, and manifold characterization in nonlinear dynami- 10 cal systems. However, its theory is not clearly connected to its practice. Furthermore, there is no associated sampling 11 theory or guidelines for optimal parameter selection to estimate the hidden dynamics reliably. This paper focuses 12 on a geometrical interpretation of PSW that coherently bridges its theory and practice by providing the needed 13 theoretical insights and explaining practical constraints. The corresponding algorithm’s parameter space is explored 14 to provide reliable and accurate estimates of the PSW function guided by the obtained geometrical properties and 15 insights. Numerical examples of a nonlinear hierarchical dynamical system with various hidden processes and ob- 16 servable dynamics are used to guide the parameter selection for the PSW algorithm. Parameter selection guidelines 17 are obtained through global sensitivity analysis to the estimation accuracy of the simulation results. The established 18 guidelines are used to extract fatigue damage evolution in 3D-printed beams from experimentally obtained vibration 19 data. The obtained results show how the PSW-based fatigue tracking can be used for early fatigue damage detection. 20 as indicators for early damage detection


23
• the inference of the reference fast-time flow as if the slow-time subsystem is static, and 83 • estimation of the difference between the observed fast-time flow and the inferred reference flow using the phase 84 space warping function (PSWF). 85 In the following section, we emphasize the continuous-time formulation of the PSW, the assumptions made, and the 86 approximations required to obtain reliable PSW estimation. Later, practical considerations of the PSW theory and its 87 algorithmic development through sampling are discussed. Aside from the theoretical developments, a comprehensive 88 graphical illustration of the PSW theory's geometry and its algorithmic considerations is provided. 89 PSW is applicable to the following HDS 91ẋ = f (x, t; µ(φ)), φ = ǫg(φ; η(x)), where x ∈ X ⊂ R m is the state variables that evolves in fast-time and it is directly observable (e.g., the structural 92 dynamics); t ∈ T ⊂ R is the time variable; φ ∈ D ⊂ R n is the slow-time dynamic variable which is unobservable or 93 partially observable (e.g., the damage evolution); µ(φ) ∈ R k is the fast-time parameter vector affected by the slow-94 time dynamics; η(x) ∈ R q is the slow-time loading function; 0 < ǫ ≪ 1 constant reflects the difference between the 95 fast and slow dynamics (for the rest of the paper, slow/fast system, slow-time/fast-time system, and damage/system 96 response are used interchangeably); and s is the observation of the fast dynamic variable (e.g., measurement function 97 of structural dynamics). 98 PSW quantifies the deformation (i.e., warping) of the fast-time phase space S X × T described by the vector field 99 f (·) (please, refer to Fig. 1 for the following discussion in this section). This deformation is assumed to be caused 100 by the change in the fast-time system parameters as a function of the slow-time variable, µ(φ). More explicitly, 101 we assume that for all φ ≥ 0, there exists an invariant, attractive, and finite-dimensional manifold M φ ⊂ S, that 102 smoothly depends on the change in φ. To track this deformation, the key operation in PSW is to estimate how the 103 fast system would have evolved for a healthy reference system (i.e., whenφ = 0 and µ = µ(φ R ) µ R ) and we call 104 this the reference dynamics. If one knows exactly how the fast dynamics evolves for some fixed reference damage 105 state (φ = φ R ), there is no need to estimate this reference dynamics. However, in practice, we do not expect to 106 have enough knowledge of the reference dynamics. As a result, one has to model this reference dynamics from a 107 designated reference portion of the acquired field measurements s. 108 To quantify the PSW, we first denote a fast-time flow in phase space S staring from initial conditions x(t 0 ) = x t0 109 and φ(t 0 ) = φ t0 , with t 0 an arbitrary initial time, as 110 x(t 0 + t, x t0 ; φ t0 ) = F t (x t0 ; µ(φ t0 )) .
To see the utility of this function in tracking the changes in slow-time variable ∆φ φ t0 − φ R , we expand the flow 112 about the reference value of damage variable Provided that the matrix exists and O ∆φ 2 is negligible, the PSWF is an affine transformation of the change in slow-time variable Therefore, given the reference damage variable φ R and if m ≥ n, we can write where P † R (x t0 ) is the pseudo-inverse of P R (x t0 ). If m < n, we cannot obtain P † R (x t0 ) and thus we do not have full 117 observabilily of φ t0 if we use only one initial condition x t0 . However, we may recover this observability by considering 118 a set of initial points in the fast-time phase space and estimating P † R in the least-squares sense over all those points.
x ∈ X ⊂ R m x(t 0 + t; µ R ) (x t0 , φ t0 , t 0 ) x(t0 + t; µ(φt 0 )) x(t0 + t; µR) + ∆φ x and φ axes are meant to represent m-and n-dimensional coordinate axes, respectively. The two instants of fast-time phase space S φ R and S φt are shown as orange sections through the reference and current values of slow-time variable φ R and φ t , respectively. Two different flows are originated from their initial conditions, where the blue curve depicts the flow in the reference fast phase space starting at (x t0 , φ R , t 0 ) and the red curve shows the flow in the current fast phase space starting at (x t0 , φ t0 , t 0 ), where t 0 is an arbitrary initial point of time. The image of the reference flow embedded into the current phase space S φt is shown by the dashed blue curve. The phase space warping function e t (x t0 , φ t0 ) is defined as difference between the fast flows starting from the same point but for the reference and current values of slow variable φ R and φ t , respectively. In practice, we will not have access to the reference flow starting exactly from (x t0 , φ R , t 0 ), and we have to use other available reference fast trajectories (shown by dashed thin blue lines on the light blue reference manifold M φ R ) to estimate it. Please note that M φ R from the actual system will have some small thickness along the φ axis, as damage evolves ad interim, which is less than O(ǫT R ).

122
The measurement s in Eq. (1) almost always has a lower dimensionality compared to the fast state vector x. Assuming 123 that s ∈ R and using ideas from the embedology [17], we can obtain a reconstructed state variable y topologically 124 equivalent to x using the delay coordinate embedding y = y(s) = y(h(x)) = y(x) or 125 systems is topologically equivalent. Now, we denote the corresponding flows of the dynamical subsystems as 132 y(t 0 + t, y t0 ; µ(φ t0 )) = F t (y t0 ; µ(φ t0 )) , where y t0 = y(t 0 ) and φ t0 = φ(t 0 ). By definition, Thus, the corresponding slow-time dynamics over a short time interval ∆t can be approximated as If the same map is applied k − 1 more times to the above slow-time state variable, the following is deduced (please 135 refer to Appendix A for detailed derivations), Therefore, the change in the damage state variable ||φ(t 0 + t, φ t0 ) − φ t0 || ∼ O(ǫt).

137
Due to the inherent time-scale separation, PSW focuses on reconstructing averaged slow-time damage variable where T represents the window size chosen for averaging over which fast-time dynamics can be treated as stationary.

139
Using the result from above, we can see that for any ν ∈ [t, t + T ], we have In this framework, we can specify the reference damage state of the slow dynamics as where T R ≥ T as we expect the rate of damage evolution to be the slowest initially, we may need more data to build 142 a suitable reference model. In general, as illustrated in Fig. 2, the total change in the averaged slow-time variable is 143 of the following order and it indicates that the change in the slow variable is detectable only if the observation time in the fast variable is 145 above that order. Fig. 2 Yt0 . PSWF e t (dashed red vector) is estimated by approximating F t (y Q ; µ(φ R )) (dashed blue trajectory) using the tangent space (straight blue line) of the reference trajectories starting from the nearest neighboring reference points y (i) (t 0 ; φ R ) r i=1 of y Q = y(t 0 ; φ t0 ) (only one of those trajectories is shown as solid thin blue curve). The rest of the nearest neighbor trajectories (NNTs) are encircled by the blue-shaded areas which is mapped into an ellipse over a short period of time, t. When the short-time assumption breaks, the NNTs are affected by the globally nonlinear dynamics (please refer to the distorted s-shaped ellipse), making the prediction unreliable. Base on those reference trajectories we build a LLM to estimate y(t 0 + t; φ R ), this introduces the prediction error in the estimate ε k = y(t 0 + t; φ R ) − y(t 0 + t; φ R ). If this error is small, then the estimated PSWFê k is close to the true value e k =ê k + ε k .

147
If a model of the fast-time dynamics for the initial damage state φ(0) is not available, we can use the measured 148 reference fast-time dynamics to approximate it over short time interval t ≪ T and for any t 0 ∈ [0, F t (y; µ R ) can be viewed as a model for a fast-time subsystem for the reference value of φ R ≈ φ(0).

150
To infer how the fast-time flow would evolve for the reference slow variable, we work in the reconstructed phase 151 space as shown in Fig. 3. In practical applications, we use local linear models (LLMs) to estimate PSW-based 152 that the tangent ∂F t ∂y is a sufficient approximation to the actual reference flow, and the second assumption ensures 155 slow dynamic variable does not significantly affect the fast-time flow during the prediction. Specifically, in the 156 reconstructed phase space, we consider some arbitrary state vector on the current flow, y Q = y(t 0 ; φ t0 ), and look for 157 a set of nn state vectors on the reference flow, which are similar to y Q , we find them using y (i) (t 0 ; φ R ) = y (i) (σ(y Q )) 158 some similarity metric σ(·). Then, we approximate the globally nonlinear autonomous map with an LLM-based on 159 a set of known state vectors tuple that describes the linearized map. These model parameters can be estimated in the least-squares sense using 162 all points in the neighborhood {y (i) } nn i=1 . Therefore, the inferred flow can be expressed as an estimate based on this 163 LLM, The PSWF for the current fast-time phase space point y Q = y(t 0 ; φ t0 ) and the current slow-time state φ t0 with 165 respect to the reference value φ R is given by where e t (·) is the short-time reference model prediction error or the damage tracking function in damage estimation 167 applications; F t (y Q ; µ(φ t0 )) is a reconstructed fast-time phase space point that is t-time ahead of y Q point, which is 168 directly available from measurements, F t (y Q ; µ(φ t0 )) = y(t 0 + t; φ t0 ). However, the corresponding fast-time evolution 169 that would have occurred for the reference/healthy state of the system F t (y Q ; µ R ) needs to be estimated. Using whereA t (y Q ) and b t (y Q ) are estimated from the nearest neighbors y (i) (t 0 ; φ R ) nn i=1 found in the reference fast-time 172 data record and their images t-times later y (i) (t 0 + t; φ R ) nn i=1 . In particular, we find the best fit parameters in the 173 least squares sense to the following map for all points in the neighborhood of y Q . The solution to this in the least squares sense is: where (·) † is the pseudo-inverse of (·),

177
In practice, the estimation quality of the PSW is affected by the sampling and the subsequent data analysis method.

178
The acquired fast-time series are first preprocessed, including data visualization, data truncation, and filtering.
damage extraction, which is referred to as the PSW algorithm as shown in Fig. 4.  and (20), the estimated PSWF is the difference between the estimated reference state y i+k (φ R ) and the current state 206 y i+k (φ i ), which can be expressed in the discrete time as, This is called the training stage of the reference LLM, based on the reference reconstructed phase space points. As 216 a result, Eq. (20) can be expressed as the following estimate,  Estimating PSWF e k i for a reconstructed phase space point y i (φ i ) entails (1) finding the nearest neighbors of y i (φ i ) in the reference phase space for φ = φ R contained inside a blue ball NN i and their images k-samples later contained in the blue ellipsoid NN i+k ; and (2) finding the optimal LLM between these neighborhoods to estimate the the same as the error in the reference LLM prediction.
As long as the dimensionality of e k is greater than or equal to the dimensionality of φ, its observability condition is 223 valid. Reference LLMs are obtained through ordinary least-squares using the reference data inside the two ellipsoids Fast-time Dynamics Averaged Slow-time Dynamics where c i (φ j ) ∈ R d×1 and N i is the number of current points in the hypercuboid B i . One should realize that Eq. (29) 233 involves averaging over the sampled time in each of the hypercuboids. the averaged PSWF values in the empty boxes for the current step will be replaced by their predecessors (i.e., 243 last observation carried forward method in [20]). To help extract the slow-time dynamics from the obtained local 244 averaged PSWFs, they are reshaped into a each windowed data records j = 1, · · · , nr: and then concatenated into a tracking matrix, which can also be written as where each column E i , i = 1, 2, · · · , N d represents the projection of the variation of the locally averaged PSWF onto 247 each independent local embedding dimension in Fig. 6 (c).

248
Reference [10] claims this phase space sectioning suppresses spurious fluctuations that are not related to the slow-time 249 dynamics without giving any explanations. Here, we address this by providing a supporting explanation as to why 250 the tracking accuracy favors partitioning. Consider a system whose reconstructed phase space is two-dimensional, in more empty subspaces, which will cause rank-deficiency in the error matrix E.

264
After the error matrix is obtained, multivariate data analysis is applied to obtain a desirable representative of the derivative DE, usually estimated using a finite difference operator to get: where U and V are unitary matrices, D is some finite difference operator in time, S and C are diagonal, and X is 270 a full rank matrix composed of smooth modes. Then, the corresponding smooth orthogonal coordinates (SOCs) are 271 obtained using: where Φ ∈ R nr×N d is used to estimate the hidden slow dynamics trajectory by considering only the smoothest The PSW algorithm is a non-parametric estimator for the slow-time dynamics in the HDS that requires numerous 277 input parameters that determine the algorithm's performance. However, the parameter selection was not explicitly No partition 4 partition Uniform Expansion Local Warping Global Warping the variance of the output. In particular, the GSA helps one to prioritize (or skip) the selection of the parameters 288 (i.e., the factor prioritization). Given a square integrable function f (·)-in this context, the PSW algorithm-over 289 a m-dimensional parameter space (usually represented by a m-dimensional hypercuboid), where the random variable X represents the input variables to the PSW algorithm, and for each of the components 291 X i , the domain is designated by some type of probability distributions. Sobol proposed the higher-order model 292 representation which decompose the function f (·) into, in which each individual term is also square integrable over the domain of the functions. Based on the higher-order 294 model representation, the model can be treated as the function described in Eq. (36). To characterize the sensitivity, 295 one first square integrates every term in Eq. (36) and by implementing the notation of the variances, one can obtain 296 the ANOVA-HDMR (analysis of variance-higher dimensional model representation) decomposition, where V (·) now is the variance estimator, is the variance to the bivariate function f ij , and so do the rest of the components.
represents the sensitivity of the first order (or the univariate part of the function decom- where U (·, ·) represents the uniform distribution whose lower bound and upper bound are separated by a comma.

352
All the samples drawn from these distributions are of integer values.  The ramp and linear-exponential damage function can be expressed explicitly as, where φ e is the critical damage value when the simulations stop; Θ(·) is the Heaviside function; • indicates a 368 Hadamard product; T e is the total time to failure; r φ represents the damage initiation time as a ratio to the total 369 time to failure; t is the time variable; and a i , i = 1, · · · , 3 are the parameters that control the proportionality of the and in the case of linear-exponential damage accumulation, one can set the damage initiation ratio to infinity (i.e., 375 r φ ∼ ∞) and obtain, To examine the parameter selection under various measurable fast-time dynamics, a Duffing oscillator with deterio-377 rating stiffness is considered,  Table 1: Parameter used to produce system responses with various complexity averaging parameters rl = 2 6 and nr ensuring 0% overlapping. The first dominant SOC from each PSW estimation, 400 obtained from Eq. (34) is treated as the estimate of the damage-time history and is scaled to its original range.

401
Eventually, the prediction quality is quantified through the mean absolute error between the estimated damage 402 and the actual damage,   Fig. 10. Plots depicting the responses due to prescribed damage evolution. Upper: damage-time histories of the ramp type with the damage initiation at 50% of the total life span; lower: system responses with the prescribed damage. Left: system response is of period-1 type; middle: system response is of period-3 type; right: system response is of chaotic type.
significantly drops when the initial crack growth rate is low, see the decrements in the sensitivity of nn in the first 422 row of Fig. 13. Among all the parameters considered in the sensitivity study, k or step has the least contribution to 423 the changes in the prediction error variance, under the considered range of k. Since the prediction error is insensitive 424 to k, this parameter may be ruled out from further consideration. In addition, we are free to choose any value of k, 425 within the considered range, without significantly affecting the prediction quality under the considered distribution.   Table 4 are embedology-parameter-normalized ; temporally, nm is normalized by the delayed  to be included in the reference model, i.e., a larger number of nm is preferred. Some cases may seem 449 like outliers, e.g., the P-1 cases in the LE damage group; however, they have a large standard deviation, 450 indicating the tracking quality is less sensitive to the selection of nm, which does not contradict the general 451 conclusion made herein.

452
(b) For the normalized statistics, the dependence of the optimal nb on the complexity of the response is 453 nullified when compared to the original statistics in Table 4. Generally, the normalized nb does not 454 correlate with the complexity of the response, and one should use less nb to reduce computation time, if 455 possible. However, when considering non-normalized nb, a significant separation in parameter sensitivity 456 between the chaotic case and the period-1 response is expected; this supports the hypothesis made in 457 section 3.2, see the difference between P-3 cases and chaotic cases in Table 4 in the Appendix.

458
(c) The nn is less sensitive to the responses' complexity than to the other two parameters. However, it favors 459 a higher number of nearest neighbors as the complexity of the response increases for accurate reference 460 modeling. Moreover, if the normalized statistics are used, the nn has a negative correlation to the response 461 complexity but spans a limited range of selection, see the rows of nn in Table 2.  Fig. 14. An example of visual inspection of the reduced parameter space from which the sampled parameters result in high prediction accuracy with the prediction errors within 5-percentile (5PCTL) of the total variance of the prediction error among 10 4 trials from the quasi-Monte-Carlo simulations. The upper-left figure is the parameter triplets' scatter plot, leading to very low prediction error for the prescribed damage. Larger and darker dots indicate the triplets that yield a lower error. The rest of the subplots show the 2-dimensional marginal PDF of the parameters, which gives the 5-percentile error. exponential growth rate increases. This observation is also consistent with the other two cases; for the period-3 response, its higher bound decreases, whereas for the chaotic one, its lower bound fades into lower 471 values of nb. As the growth rate increases, the damaged manifold diverges from the reference manifold 472 drastically, which causes empty local hypercuboids in the warped phase space if large nb is used. Therefore, 473 a conservative number of nb is always preferred.

474
(c) From the non-normalized statistics for both ramp and LE damage dynamics in Table 4, one can conclude 475 that a higher number of nn is preferred for the chaotic responses as the quasi-stationarity breaks down. 476 For example, as the exponential growth rate decreases from a 3 = 30 to a 3 = 13, the fast dynamics is less 477 stationary, more nn is preferred by the chaotic cases. And the opposite is true for the periodic cases.

478
Based on the preceding statistical inference, the general guidelines for input parameter selection for the PSW algo-479 rithm can be summarized as follows:     (a) illustrates the three distinct raster angles used during the specimen preparation where x is the reference direction, and r is the direction along which the raster is deposited; therefore, the three raster angles are named according to the angle α, the angle between the raster direction and the reference direction. (b1) illustrates the setup of the structural health monitoring system based on a Moon beam setup, where the specimen, C, is cantilevered to the frame, B, which provides harmonic base motion to the beam; two magnets, F, provides nonlinear restoring force; two laser displacement sensors are used to record the beam motion and the base motion, respectively. (b2) illustrates the boundary condition and the initial damage created about 1 inch below the fixture. The effective load at the damage site (initial notch in the red circle) is tension-compression of variable amplitude. The analog signal is low-pass filtered and properly sampled for its use for the PSW algorithm.
normalized parameters, nn = 5 to 6 can be used for the first trial. To reliably obtain the PSW estimate, the authors advocate using the normalized values of nm, nb, and nn, listed 506 in Table 2 during the parameter selection. Since the spatially related parameters are less affected by the type of The PSW is applied to three experimentally obtained displacement-time histories to validate the parameter selection 513 guidelines based on the empirically determined "trust-region" for the PSW algorithm. These displacement-time 514 histories are recorded from a magneto-elastic beam system with polymeric beam specimens whose cross-sections 515 undergo area reduction due to an initial notch and fatigue loading. The beams are 3D-printed with three distinct 516 raster orientations, and their responses to the harmonic forcing are all chaotic, see the graphical abstract of the 517 experiment setup in Fig. 15. This system can be perceived as nonlinear with (1) nonlinear stiffness and (2) time-518 varying linear stiffness whose sub-structure uses anisotropic materials. Since obtaining the damage evolution in such 519 non-conductive material is difficult, PSW using non-destructive health monitoring signals is more convenient for such 520 applications. The selection of PSW parameters is guided by the delay-coordinate embedding parameters (τ and d).

538
The 90-degree case shares a similar damage propagation rate with the 45-degree case up until the first dashed line 539 and diverges with a higher rate than the 45-degree case, which deteriorates with the same rate until the very late The authors declare that they have no known competing financial interests or personal relationships that could have 576 appeared to influence the work reported in this paper. The data sets generated during and/or analysed during the current study are not publicly available but can be made 579 available from the corresponding author on reasonable request.
As a result, one can obtain the iteration of any integer multiples k ∈ Z + and its subsequent iteration k + 1; by 585 mathematical induction, the iteration of the flow in the continuous time with a constant time frame ∆t can be 586 obtained as,

587
G k∆t (φ t0 ; η(y t0 )) = φ t0 + ǫk∆tg(φ t0 ; η(y t0 )) + O (ǫ∆t) 2 .  Fig. 17. Histograms to the parameters which yields first-dencile error, the slow dynamics is of the ramp type. From top to bottom, each row specifies ramp damage dynamics with various length of damage-free responses; and from left to right, each column represents one parameter considered in the reduced parameter space.  Fig. 18. Histograms to the parameters which yields first-decile error, the slow dynamics is of the linear-exponential type. From top to bottom, each row specifies ramp damage dynamics with various length of damage-free responses; and from left to right, each column represents one parameter considered in the reduced parameter space.