Non-recursive sequential input deconvolution

Mechanical Systems and Signal Processing 82 (2017) 296–306

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Non-recursive sequential input deconvolution Dionisio Bernal Northeastern University, Civil and Environmental Engineering Department, Center for Digital Signal Processing, Northeastern University, Boston, MA 02115, United States

a r t i c l e in f o

abstract

Article history: Received 24 April 2015 Received in revised form 14 March 2016 Accepted 17 May 2016 Available online 8 June 2016

A scheme for sequential deconvolution of inputs from measured outputs is presented. The key feature in the formulation is elimination of the initial state from the input–output relations by projecting the output in the left null space of the observability block. Removal of the initial state allows the sequential format of the deconvolution, essential for computational reasons, to be implemented non-recursively, assuring unconditional stability. Identifiability is realized when the input-output arrangement does not have transmission zeros, and observability and controllability are shown immaterial. Comparison of results from the scheme with those from Dynamic Programming highlights the benefits of eliminating the initial state. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Non-recursive sequential deconvolution Input reconstruction Regularization Dynamic programming

1. Introduction Input reconstruction algorithms are of interest in the estimation of interface forces between rails and train wheels for condition monitoring and maintenance scheduling [1,2]; in acoustics in the estimation of forces at points causing structureborne noise [3,4], in the design of equivalent load regimes for the analysis of structures subjected to wind or blast pressures [5,6], and in estimation of vehicle weights, in weight in motion applications [7,8], among others. Another setting where input estimation arises is where the forces identified are not loads acting on structures but are corrective actions that bring the response of a model in line with observations, carrying information on discrepancies between the model and the structure [9]. Accepting replacement of the structure by a finite dimensional model and restricting the scope to linear time invariant (LTI) systems, the input reconstruction problem in the time domain boils down to a linear system of equations relating the unknown inputs, plus the initial condition, to the measurements. As is always the case, one is interested on whether the solution is unique and on its sensitivity to noise and model errors. The identifiability question translates to whether the outputs over some interval, T, constraints the inputs over a shorter interval (T
d), for some d. As pointed out in [10], the answer is encoded in the structure of the kernel of the matrix that implements the convolution in discrete time, namely, inputs are identifiable over the part of the time axis (if any) where the kernel is essentially zero; and we use the term “essentially”, instead of “identically”, because, as shown in [10], the finite dimensional model cannot capture dead time exactly. Adjustments to remove rank deficiency by time shifts have been proposed but the resulting schemes are inevitably approximate and the performance has been generally poor [11,12]. Conditioning difficulties in the input estimation problem are restricted to the frequency band where the magnitude of the inputs needed to match the noise are substantial compared to the inputs to be identified. The standard approach to E-mail address: [email protected] http://dx.doi.org/10.1016/j.ymssp.2016.05.024 0888-3270/& 2016 Elsevier Ltd. All rights reserved.

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

297

mitigate noise related high frequency components in the reconstructed input is Tikhonov Regularization [13,14]. In the most commonly used Tikhonov regularization the solution is obtained by minimizing a weighted sum of the norms of the residual and the solution, and the weighting constant is referred to as the regularization parameter. Although some methods to select the regularization parameter have gained wide acceptance, the optimization of this selection remains an open problem [15]. In this paper the need for regularization is reduced by low-pass filtering the noisy outputs, prior to processing, and low pass filtering the estimated inputs, as a post-processing step. It is evident that the batch solution of the input estimation problem boils down to solving a system of coupled linear equations on the inputs and the initial state. A batch solution, however, is not computationally feasible when the inputs are long, and is not an option in applications where the time between the appearance of an input, and its estimation, is important. These limitations can be circumvented by a solution on a window whose size and rate of advance respects causality and ensures stability, as is done in the Sequential Deconvolution Reconstruction (SDR) scheme, presented in [10]. The objective of this paper is to introduce a modified version of SDR that is unconditionally stable. A version that is designated here as the Non-Recursive SDR or NR-SDR, with the acronym highlighting the fact that results in any window are independent of prior estimates. The key feature of NR-SDR is elimination of the initial state from the input output relations by projection of the output on the left null space of the observability block. A reasonable question to ask is whether the NR-SDR approach offers any advantage compared to a Dynamic Programming (DynProg) solution [16–19]. The DynProg approach is a computationally efficient scheme to obtain the regularized least square solution, as the dimension of the matrices that have to be inverted do not depend on the input duration, and duration affects only the storage requirements. Inspection of the theory suggests that if the initial state is adequately estimated by the DynProg algorithm, the performance of the two algorithms should be similar, and this is what simulations show. The estimation of the initial state, however, is a poorly conditioned problem for long inputs and there are, at the time of writing, no guidelines for selecting an optimal truncated length for its estimation in DynProg. Therefore, when the anticipated magnitude of the initial state is significant, NR-SDR appears to be the preferred approach. The fact that NR-SDR can be applied online while DynProg is restricted to batch applications is also worthy of explicit note. Before outlining the paper we make a parenthesis for two observations that, while not essential for the developments, are important for clarity. The first is that techniques that attempt to reconstruct inputs with one time step lag, in the socalled non-collocated case, cannot work in general because they violate causality constraints [20–24]. The second observation is a contention that the reason why this limitation hasn't been more widely recognized is because finite dimensional models have no dead time [10]. Needless to say, there are applications where the finite dimensional discrete time model is not formulated to approximate an infinitely dimensional one, and in these cases a one-time step lag prediction of causes from effects can be reasonable. The paper is organized as follows: Section 2 contains a brief description of the problem considered and sets some of the notation. Section 3 examines the issue of identifiability and focuses on the fact that identifiability implies the absence of loading conditions for which outputs are identically zero [25–27]. It is noted that observability, which places limits in system identification applications is immaterial in the input-identification problem and that the same holds true for controllability. Section 4 presents the essential contribution, the NR-SDR scheme. Regularization is examined in Section 5 and some observations on the DynProg solution are presented in Section 6. The one time step lag deterministic estimator is presented in Section 7, for completeness. Section 8 summarizes the essential results as a lead into Section 9 which contains the numerical results. The paper closes in Section 10 with a brief critical review of the content.

2. Problem considered The structure is assumed linear and time invariant and subjected to loads that are fixed at known positions. If inputs with known histories exist, it is assumed that their contribution has been subtracted from the output so there is no need to define them explicitly. There are instances where it may be appropriate to decompose the unknown inputs into two sets: one that is to be estimated plus a set that is to be treated as “nuisance loads” but we do not consider this case. Although measurement noise and model error are included in the numerical section, the derivations are carried out in a deterministic setting. In standard notation, the finite dimensional state-space representation in discrete time that replaces the structure and the corresponding state to output equations write

xk + 1 = Ad x x + Bduk

(1)

yk = Cdxk + Dd uk

(2)

and

where Ad ∈ R2nx2n, Bd ∈ R2nxr , Cd ∈ Rmx2n and Dd ∈ Rmxr with n¼number of degrees of freedom (DOF), m ¼number of outputs and r ¼number of (unknown) inputs. The matrices Bd and Dd depend on the operating assumption made on the intersample behavior of the input and, in the numerical section, we use linear interpolation between points, i.e. the so called “first-order-hold” [28]. No constraints regarding observability or controllability are imposed but we require, without loss of generality, that Bd and Cd be full column rank, and full row rank, respectively.

298

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

3. Identifiability of input histories Determining the number of independent inputs from observations of the output is a problem that can be solved in a data-driven fashion and the location and direction of the loads ascertained without resorting to combinatorics [29,30]. Here we assume that these questions have been answered and focus on extraction of the input histories. In the general case, where the initial state is not known, inputs are identifiable, if and only if, the set of Output Zeroing Loading Conditions (OZLC) defined by Eq. (3) is empty.

OZLC: =

{ (x , u ) 0

k

}

yk = 0 ∀ k

(3)

3.1. Output zeroing loading conditions (OZLC) OZLC can be computed in continuous or discrete time. Here we develop the more directly relevant to the practical implementation, discrete time version. Consider an input whose samples are

uk = u0z 0k

(4)

where u0 is a distribution, z0 is a complex scalar and the inter-sample behavior is the one assumed to pass from continuous time to discrete time (C2D) when formulating Eqs. (1) and (2). One can confirm that for this input the state recurrence is satisfied if the state trajectory evolves according to

xk = x 0z 0k

(5)

Substituting Eqs. (4) and (5) in Eqs. (1) and (2) and equating the output to zero writes

⎡ A − I⋅z B ⎤⎧ x ⎫ ⎧ 0⎫ 0 d ⎢ d ⎥⎨ 0 ⎬ = ⎨ ⎬ Cd Dd ⎦⎩ u0 ⎭ ⎩ 0⎭ ⎣

(6)

The values of z0 that satisfy Eq. (6), if any, are known as invariant zeros, x0 are the initial state zero directions and u0 are the zero direction load distributions [26]. Provided z0 is not a pole of the transfer matrix, so there is no pole-zero cancellation, the triplets { z0, x0, u0} define the OZLC. Since the dimension of the polynomial matrix on the left side of Eq. (6) is (2n + m) × (2n + r ) it follows that if r 4m the full complex plane satisfies the equation and, therefore, a necessary condition for input identifiability is that m ≥ r , namely, that the number of outputs be no less than the number of inputs to be identified. If the initial conditions are known their effect can be subtracted from the output and one concludes that lack of identifiability requires that Eq. (6) be satisfied for some z0 with x0 ¼0, which implies that the second column partition in Eq. (6) has to be rank deficient. This, however, is prevented by the fact that we've required that Bd be full column rank so it follows that when initial conditions are known the inputs are identifiable if m ≥ r . We close by noting that asymptotic stability of the system does not imply that errors in the estimated inputs from initial condition discrepancies vanish. An example is that where only acceleration measurements are available. In this instance one can check, by substitution in Eq. (6), that the triplet { z0, x0, u0} = 1 {K−1u0 , 0}T u0 is a solution, where K is the system stiffness matrix. The physical interpretation being that static load components are unobservable from acceleration measurements if the initial condition is unknown.

{

}

4. Non-recursive sequential deconvolution reconstruction (NR-SDR) The well-known relation between the initial state, inputs, and outputs, in a LTI system is k

yk = CAdk x 0 +

∑ Yjuk − j

(7)

j=0

where the so-called Markov parameters are

Y0 = Dd

Y j = CdAdj − 1Bd

(8a,b)

Stacking the inputs and outputs in columns writes

⎡Y ⎧ y0 ⎫ ⎡ Cd ⎤ 0 ⎥ ⎢ 0 ⎪ ⎪ ⎢ ⎪ y1 ⎪ ⎢ CdAd ⎥ ⎢ Y1 Y0 ⎨ ⎬−⎢ ⎥x 0 = ⎢ ⋮ ⎪ ⋮⎪ ⎢ ⋮ ⎥ ⎢⋮ ⎪ ⎢⎣ Yℓ Yℓ− 1 ⎩ yℓ ⎪ ⎭ ⎢⎣ CdAdℓ ⎥⎦

⋯ 0 ⎤⎧ u0 ⎫ ⎥ ⎪ u1 ⎪ ⎪ ⋯ 0 ⎥⎪ ⎨ ⎬ ⎥ ⋱ 0 ⎥⎪ ⋮ ⎪ ⎪ ⎪ ⋯ Y0 ⎥⎦⎩ uℓ ⎭

(9)

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

299

or, with obvious notation

y[0, ℓ] − Obℓ ⋅x 0 = Hu[0, ℓ]

(10)

(ℓ1m)x(ℓ1r )

where H ∈ R , Obℓ is the observability block of order ℓ , ℓ=total number of time steps and ℓ1 = ℓ + 1 is the total number of time stations. Let ℓ1 > 2n/m so the extended observability block is guaranteed to have a left null space, namely, a matrix Z ∈ R(mℓ1− g )xmℓ1 where g ¼number of observable modes, such that

Z⋅Obℓ = 0

(11)

Pre-multiplying both sides of Eq. (10) by Z writes

yz, ℓ = Hzur ⋅ℓ1

(12)

2

where

yz, ℓ = Z⋅y⎡⎣ 0, ℓ⎤⎦ and Hz = Z⋅H

(13a,b)

2

(mℓ1− g )xℓ1r

with Hz ∈ R and where the subscripts without brackets are the length of the associated vectors, namely, vector yz is ℓ2 = mℓ1 − g and the input vector is r⋅ℓ1. The general solution of Eq. (12) is

ur ⋅ℓ1 = Hz− *yz, ℓ + Ker (Hz )⋅h

(14)

2

where




* stands for pseudo-inversion and h is an arbitrary vector of appropriate dimension. Defining

Q p = ⎡⎣ Ip ⋅ r 0p ⋅ rx(ℓ1⋅ r − p) ⎤⎦

(15)

and multiplying Eq. (14) by Eq. (15) gives

ur ⋅ p = Q pHz− *yz, ℓ

(16)

2

provided p is sufficiently small so that (ℓ − p) is larger than any delay and the inputs are identifiable. We note that the matrix Qp need not be explicitly formed as the product on the right hand side of Eq. (16) can be expressed in terms of the economy SVD of Hz. Namely, let

ℓ1 ≥

2n m−r

(17)

so, Hz is a tall matrix, and let its economy SVD be

Hz = UzSzVzT where Uz ∈ R

(mℓ1− g )x(r ⋅ℓ1)

(18) , Sz ∈ R

r . ℓ1xr ⋅ℓ1

and Vz ∈ R

(r ⋅ℓ1)x(r ⋅ℓ1)

. One can then confirm that

→ Q pHz− * = VzSz− *UzT

(19)

→ where Vz ∈ R(pxr )x(ℓ1r ) is the matrix obtained by retaining the first (pxr) rows of Vz. Taking the number of non-zero singular → ↔ values in Sz (the rank of Hz) as rz and defining the matrices Uz ∈ R(mℓ1− g )x(rz ) and Vz ∈ R(pxr )xrz where rz are the first rz columns, it follows that

↔→−1→T ur ⋅ p = ( Vz S z Uz )yz, ℓ

(20)

2

→ where Sz ∈ Rrzxrz . Note that Eq. (20) is non-recursive and independent of the initial state, implementation to solve for inputs of arbitrary length is evident. 4.1. Identifiability of the input at the first station of each window The first input of any window is often poorly conditioned because it can be replaced by “nearly equivalent” initial conditions. To illustrate consider the situation

x 0 = 0, uk = u0 k = 0; uk = 0 k = 1, 2...

(21)

for which

x 0 = 0, x1 = Bdu0 , xk = x1e Ac (k − 1)Δt

(k = 1, 2... ), y0 = Dd u0 , yk = Cxk

and the situation where the loading is zero but the initial condition is x0 =

(22a-e) Ad−1Bdu0 ;

in this instance

300

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

x 0 = Ad−1Bdu0 , x1 = Bdu0 , xk = x1e Ac (k − 1)Δt

(k = 1, 2... ), y0 = CAd−1Bdu0 , yk = Cxk

(23a-e)

comparing Eqs. (22) and Eqs. (23) shows that the outputs in the two cases are identical everywhere except at k ¼0 where the difference depends on ε , where

ε = ⎡⎣ Dd − CAd−1Bd⎤⎦

(24)

The matrix in Eq. (24) is identically zero when expressed in continuous time if accelerations or velocity measurements are used and this suggests that the norm of ε is likely small in discrete time when these types of sensors are used.

5. Regularization Inspection shows that Eq. (20) can be written as

⎛ rz ↔→T ⎞ vj u j ⎟ ⎜ ur ⋅ p = ⎜ ∑ ⋅y ⎜ j = 1 sj ⎟⎟ z, ℓ2 ⎝ ⎠

(25)

where the lower case symbols with the j subscript are the jth columns of the associated matrices and sj is the jth diagonal → entry in Sz . The reconstructed input is, therefore, the sum of the projection of the transformed output on each of the rankone matrices in the parenthesis. A well-known difficulty arises when there are very small singular values since in this case the noise in the transformed output can, for some values of j, have large contributions compared to that of the real output. The standard approach to prevent large noise contributions in linear problems is Tikhonov regularization [13], which, in the common case where the weighting matrix is the identity is readily implemented in Eq. (25) by changing sj to sj, m using

s j, m =

s 2j + α sj

(26)

where α is the regularization parameter. Selection of the regularization parameter to minimize some norm of the error is an open problem. Several recommended criteria were explored during the course of this work but none was found to provide results that consistently matched the optimal regularization determined by trial and error (which is feasible, of course, only in the academic setting where the solution is known). In any case, we settled on the use of the Generalize Cross Validation (GCV) approach [31] but with a reduction of two orders of magnitude. The fact that all methods suggested higher regularization than the optimal value is consistent with findings previously reported in [18]. It should be clarified, however, that the solutions with the higher regularization values are not far from the results for the regularization selected. For implementation of GCV we used the Matlab Toolbox in [32]. 5.1. Frequency dependency Conditioning difficulties arise because the high frequency noise components in the output demand large amplitude high frequency inputs to match them. This suggests that the amount of regularization can be reduced by low pass filtering the output prior to its use in Eq. (25). Selection of the cutoff frequency requires that one have some idea about the frequency content of the signals to be estimated but this is usually not a problem. Note that conditioning is also promoted by selecting a sampling rate that is no larger than necessary to describe the inputs -in the sense of the C2D transfer – as the sampling rate fixes the bandwidth of the noise.

6. A few observations on the dynamic programming solution The contribution from the first time station to the cost function that is minimized in a DynProg solution is

J1 = x 0T R1x 0 + x 0T s1 + q1

(27)

where R1, s1 and q1 are results obtained from the backward pass and x0 is the (unknown) initial state. Since the objective is to minimize the total cost, which is the sum over all time steps, the DynProg estimate of the initial state is

x0 = −

1 −1 R1 s1 2

(28)

Eq. (28) often fails, however, because R1 diverges. One suspects on physical grounds that the issue may be resolved by restricting the length of the output used to estimate x0 and an exploratory investigation showed this to be the case. Research to identify an approach to select the optimal length appears justified but was deemed too far from the central theme of this

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

301

paper. In the numerical section we simply illustrate that durations for which R1 does not diverge, and which are otherwise reasonable, can nonetheless result in significant error in the initial state, and as a consequence, in the early part of the input estimation. We close this brief section with two observations: the first is that while R1 cannot attain full rank if there are unobservable states this is not an issue in input estimation because replacing the inverse with the minimum norm solution gives the initial state projection that influences the reconstruction. The second is that while DynProg violates causality, since it offers input estimates for the same length of time as the output sequence, the bogus values at the (generally few) time stations not constrained by the output are typically of little consequence in the overall solution.

7. Deterministic one lag input identification filter Although the one lag input predictor can only work only under very restricted conditions, we present it for completeness, and to allow reference to it in the numerical section. For simplicity we present the case where the feedthrough matrix is zero, although the case where this matrix is non-zero but rank deficient can be easily fitted in the derivation. The derivation shown is inspired in the decomposition introduced by Ho and Müller [33] and subsequently used in [21]; it begins by factoring the matrix Bd as

⎡ R⎤ Bd = N ⎢ ⎥P T ⎣ 0⎦

(29) T

where N and P are orthogonal and R is diagonal, substituting Eq. (29) into Eq. (1), pre-multiplying by N and defining

x¯ k = NT xk

u¯ k = P T uk

(30a,b)

gives

⎡ R⎤ x¯ k + 1 = A¯ x¯ k + ⎢ ⎥u¯ k ⎣ 0⎦

(31)

which shows that the bottom partition of the transformed state does not contain the input explicitly. Eq. (31) can be written as

⎧ x¯1, k + 1⎫ ⎡ A¯11 A¯12 ⎤⎧ x¯1, k ⎫ ⎡ R ⎤ ⎥⎨ ⎬ + ⎢ ⎥u¯ k ⎨ ⎬=⎢ ⎩ x¯2, k + 1⎭ ⎢⎣ A¯ 21 A¯ 22 ⎥⎦⎩ x¯2, k ⎭ ⎣ 0 ⎦ ⎪

⎪

⎪

⎪

(32)

where the dimensions of the various partitions are apparent. Consider now the output equation, namely

yk = C xk

(33)

substituting Eq. (30a) into Eq. (33) gives

yk = C N x¯ k

(34)

which, with obvious definitions, can be written as

⎧ x¯1, k ⎫ ⎬ yk = ⎡⎣ C¯1 C¯2 ⎤⎦⎨ ⎩ x¯2, k ⎭ ⎪

⎪

(35)

where, recognizing that R is square and equal to the number of inputs, shows that the number of columns in C¯1 is also equal to the number of inputs. Requiring that the number of outputs be no less than the number of inputs ensures that C¯1 is full column rank and one can thus write, without approximation − x¯1, k = C¯1 * yk − C¯2x¯2, k

(

)

(36)

Substituting Eq. (36) in the second partition of Eq. (32) gives, after some simple re-arranging

⌢ ⌢ x¯2, k + 1 = A22 x¯2, k + A21yk

(37)

where

⌢ − A22 = A¯ 22 − A¯ 21C¯1 *C¯2

⌢ − A21 = A¯ 21C¯1 *

(38)

The first partition of Eq. (32) reads

x¯1, k + 1 = A¯11x¯1, k + A¯12 x¯2, k + Ru¯ k so it follows from Eq. (39) and (30b) that

(39)

302

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

uk = PR−1 x¯1, k + 1 − A¯11x¯1, k − A¯12 x¯2, k

(

)

(40)

given a known initial condition one can solve Eq. (37) to get x¯2. The results are subsequently used to solve for x¯1 in Eq. (36) and the unknown inputs are obtained from Eq. (40). 7.1. Stability of recursive solutions The numerical stability of the recursive version of sequential deconvolution, as given in [10], writes

ρ = ηj

max

≤1

(41)

where η are the eigenvalues of

Pp = Adp − Rcp − 1Q pH − *Obℓ

(42)

with Rc as the reversed controllability of order p
1 given by

Rcp − 1 = ⎡⎣ Adp − 1Bd Adp − 2 Bd ⋯ Bd ⎤⎦

(43)

One suspects that the stability of the one lag estimator can be evaluated by taking {ℓ, p} ¼{1,1} and numerical results show that this is in fact the case.

8. Summary We summarize the main results and contentions as a lead into the numerical section.

 The NR-SDR algorithm is encapsulated by Eq. (20).  NR-SDR is unconditionally stable and accounts for arbitrary initial conditions without the need to estimate them explicitly.

 Inputs are strictly identifiable if there are no transmission zeros.  Observability and controllability are immaterial in the reconstruction of inputs.  Pre-filtering the output and post-filtering the estimated inputs reduces the level of regularization needed and generally enhances accuracy. From a theoretical perspective one expects that NR-SDR and DynProg will perform similarly if the initial state is adequately estimated in the DynProg solution. The fundamental advantage of NR-SDR being that the initial condition need not be estimated. The fact that NR-SDR operates with a small prediction lags while DynProg is essentially an off-line scheme may be relevant in some applications.

9. Numerical illustrations 9.1. Example#1 Consider an 8-DOF uniform chain system with a fundamental frequency of 1 Hz and Rayleigh damping that gives 1% of critical in the first and the last mode. The sampling rate is 50 Hz. The initial state is selected randomly, with a magnitude that is of the order of the state amplitude during the response, and is assumed unknown. In all instances noise with 2% of the signal standard deviation is added to the measurements and discrepancy between the identification model and the model used to generate the output data is included. Model error is selected so masses are off by the addition of random perturbations with limits of 5% and stiffness values are perturbed with limits of plus or minus 10%. An input signal obtained by filtering an earthquake record is applied at coordinate #8 while outputs are: accelerometers at {1,3,5} and displacements at {3,5}. The system is such that the first DOF is connected by first spring to the only existing support. Solution: Examination shows that there are no OZLC so the input is strictly identifiable. The regularization parameter from GCV, reduced by two orders of magnitude, is 4.88e
7. In NR-SDR we choose { ℓ, p} = { 60, 50}, namely a 1.2 s window that advances 1.0 s at a time. The initial condition for the DynProg solution is estimated using two seconds of output. The results in Fig. 1 show, as expected, that both techniques work well, with the behavior discussed in Section 4 exemplified in Fig. 1a. The spectral radius of Eq. (41) for {ℓ, p} = {1, 1} is 7.7, indicating that the one-lag estimator would prove unstable and this is what is found in a simulation.

force (N)

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

20

20

10

10

0

0

303

Actual Input

NR-SDR -10

-20

-10

Actual Input 0

5

10

-20

DynProg

5

0

time(sec)

10

time(sec)

Fig. 1. Reconstructed input acting @8 vs actual signal.

20

20

NR-SDR

force(N)

10

0

0

-10

-10

0

NR-SDR

Input @8

Actual Input -20

Actual Input

10

5

time(sec)

Input @4

10

-20

0

5

10

time(sec)

Fig. 2. Reconstructed inputs vs actual signal using NR-SDR (left: input @8, right: input @4).

9.2. Example #2 Exactly the same situation as in example#1 is considered, except that the time reversed version of the input at coordinate#8 is applied as a second input at coordinate #4. Solution: The inputs prove identifiable and the regularization is α = 2.805e−6 . The solution with NR-SDR is contrasted with the exact signals in Fig. 2. The deterioration of accuracy relative to the results in Fig. 1 is evident and is rationalized by noting that the number of unknowns is double and the constraints are the same. Solutions obtain using DynProg for the first input when the response segment used to estimate the initial condition is 0.50, 1.0 and 1.50 s are depicted in Fig. 3. As can be seen, the solution for the first duration is good but for the second and third windows the results are poor in the first three to four seconds. Results for the one-time step lag predictor are not illustrated as the approach proved unstable. 9.3. Example #3 This example exemplifies a scenario when OZLC are feasible. The situation is the same as in example #1 except that the displacement measurements are removed. Examination shows that for this arrangement the system has 2 transmission zeros, both real, located at s¼{
0.0192, 1}. The first zero is close to the origin and is thus associated with loads at very small times (from the start of each window) but the zero at s ¼1 allows for arbitrary biases. The zero direction initial state for the zero at s ¼1 is x0,dis ¼{0.08,0.17,0.25,0.33,0.41,0.50,0.58,0.66} x0,vel ¼0 and the loading (for the selected scaling) is u0 ¼10 N; the initial state being, of course, the static position when a load of 10 N acts at coordinate #8. In the simulations we add the

304

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

20

20

Actual Input

force(N)

10

DynProg

20

DynProg Actual Input

10

0

0

-10

-10

-10

IC = 1.0 sec

IC = 0.5 sec 0

Actual Input

10

0

-20

DynProg

5

10

-20

0

5

IC = 1.5 sec 10

-20

0

5

10

time(sec) Fig. 3. Dynamic Programming estimation of the input@8 for 3 different response durations used to estimate the initial condition.

20

20

Actual Input

Actual Input 10

force (N)

force (N)

10

0

-10

0

-10

NR-SDR -20 0

DynProg

5

-20

10

0

5

time (sec)

10

time (sec)

Fig. 4. Same as Fig. 1 but with displacement measurements removed: in this case the “true loading” contains a contribution from an OZLC.

noted initial condition to the arbitrary previously selected one and a load of 10 N to the input and perform the identification. Fig. 4 illustrates the expected results. 9.4. Example #4 This last example illustrates input identification in an arrangement that is not observable or controllable. Inspection of the system depicted in Fig. 5 shows that the 2nd, 4th and 6th modes are uncontrollable and that the 4th mode is unobservable. Analysis shows, however, that there are no transmission zeros so the input is identifiable. The input history is Measurements Accelerometers @ ={2,6} Displacement @ 4

k

k m (1)

1% Rayleigh damping – 1 and last mode m = 1.0 kg , k = 1037 N/m

k

k m

m

m u(t)

(2)

k

k

k m

f = 0.5 Hz, 50Hz sampling

Fig. 5. System of example #3.

m

k m (7)

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

305

20

force(N)

10

NR-SDR 0

Actual Input

-10

-20

0

1

2

3

4

5

6

7

8

9

10

time(sec) Fig. 6. Comparison of NR-SDR solution and actual input for the system in Fig. 5.

the same used in example #1 and the initial state is an arbitrary state (with non-zero projection in the 4th mode). The result from NR-SDR when 2% noise is added to the output is depicted in Fig. 6. The reconstruction, as expected, is accurate.

10. Conclusions NR-SDR simplifies sequential deconvolution by eliminating the need to select window parameters that satisfy numerical stability and considers arbitrary initial conditions automatically. Numerical results confirm that the accuracy attained is of the same order of magnitude as that realized using DynProg, provided that in the DynProg solution the initial conditions are known. In practice knowledge of the initial conditions is usually tantamount to stating that it is reasonable to take them as zero - when this is not the case NR-SDR is the preferred alternative. The price paid for non-recursiveness and unconditional stability is (primarily) computation of the left null space of the observability block. For problems with DOF in the 100's this evaluation is not an issue but when the DOFs go into the thousands or 10's of thousands a direct assessment is impractical and a reduced order model becomes necessary. The paper brings attention to the fact that reconstruction algorithms formulated with one-time step lag do not respect causality and thus cannot work in general.

References [1] C.P. Ward, R.M. Goodall, R. Dixon, Contact force estimation in the railway vehicle wheel-rail interface, in: Proceedings of the 18th IFAC World Congress, Milano, Italy, 2011. [2] T. Zhu, S. Xiao, G. Yang, M. Wang, Estimation of wheel/rail contact forces based on an inverse technique, in: Proceedings of the 13th International Conference on Fracture, Beijing, China, 2013. [3] M.H.A. Janssens, J.W. Verheij, D.J. Thompson, The use of an equivalent forces method for the experimental quantification of structural sound transmission, J. Sound Vib. 226 (1999) 305–328. [4] A.N. Thite, Inverse Determination of Structure-borne Sound Sources (Ph.D. Thesis), University of Southampton, Highfield, UK, 2003. [5] Y. Niu, C.-P. Fritzen, Y.Q., Ni, Online wind load reconstruction study for Canton Tower, in: Proceedings of the 6th International Conference on Structural Health Monitoring of Intelligent Infrastructure, Hong Kong, 2013. [6] S. Xu, X. Deng, V. Tiwari, M.A. Sutton, W.L. Fourney, D. Bretall, An inverse approach for pressure load identification, Int. J. Impact Eng. 37 (2010) 865–877. [7] S.S. Law, Y.L. Fang, Moving force identification: optimal state estimation approach, J. Sound Vib. 239 (2) (2001) 233–254. [8] L. Yu, T.H.T. Chan, Moving force identification from bridge dynamic responses, Struct. Eng. Mech. 21 (3) (2005) 369–374. [9] R.J. Patton, P.M. Frank, R.N. Clark (Eds.), Issues of Fault Diagnosis for Dynamical Systems, Springer-Verlag, London, 2000. [10] D. Bernal, A. Usssia, Sequential deconvolution input reconstruction, Mech. Syst. Signal Process. 50 (2015) 41–55. [11] L.J.L. Nordstrom, T.P. Nordberg, A time delay method to solve non-collocated input estimation problems, Mech. Syst. Signal Process. 18 (2004) 1469–1483. [12] M.S. Allen, T.G. Carne, Delayed, multi-step inverse structural filter for robust force identification, Mech. Syst. Signal Process. 22 (2008) 1036–1054. [13] A.N. Tikhonov, A.V. Goncharsky, V.V. Stepanov, A.G. Yagola, Numerical Methods for the Solution of Ill-Posed Problems, Kluwer Academic Publishers, Berlin, 1995. [14] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 2003. [15] J.L. Muller, S. Siltanen, Linear and nonlinear inverse problems with practical applications, Comput. Sci. Eng. (2012). SIAM. [16] D.M. Trujillo, H.R. Busby, Practical Inverse Analysis in Engineering, CRC Press LLC, New York, 1997. [17] D.M. Trujillo, Application of dynamic programming to the general inverse problem, Int. J. Numer. Methods Eng. 12 (1978) 613–624. [18] L.J.L. Nordstrom, A dynamic programming algorithm for input estimation on linear time-variant systems, Comput. Methods Appl. Mech. Eng. 195 (2006) 6407–6427. [19] R. Adams, F. Doyle, Multiple force identification for complex structures, Exp. Mech. 42 (1) (2002) 25–36. [20] S. Gillijns, B. De Moor, Unbiased minimum-variance input and state estimation for linear discrete-time system with direct feedthrough, Automatica 43 (2007) 934–937. [21] D. Marquin, B. Gddouna, J. Ragot, Estimation of unknown inputs in linear systems, Am. Control Conf. 1 (1994). [22] E.M. Lourens, G. Lombaert, C. Papadimitriou, G. De Roeck, Joint Estimation of States and Input in Linear Structural Dynamics, 3rd International Conference on Computational Methods in Structural Dynamics and Earthquake Engineering (COMPDYN 2011), Corfu, Greece, 2011, pp. 25–28. [23] M.K. Sain, J.L. Massey, Invertibility of linear time-invariant dynamical systems, IEEE Trans. Autom. Control AC-14 (2) (1969) 141–149. [24] L.M. Silverman, Inversion of multivariable linear systems, IEEE Trans. Autom. Control. AC-14 (3) (1969) 270–276. [25] J. Tokarzewski, A general solution to the output zeroing problem for MIMOLTI systems, Int. J. Appl. Math. Comput. Sci. 12 (2) (2002) 161–171. [26] K. Zhou, Essentials of Robust Control, Prentice Hall, Upper Saddle River, New Jersey, 1988. [27] N. Karcanias, B. Kouvaritakis, The output-zeroing problem and its relationship to the invariant zero structure, Int. J. Control 30 (3) (1979) 395–415.

306

[28] [29] [30] [31] [32] [33]

D. Bernal / Mechanical Systems and Signal Processing 82 (2017) 296–306

D. Bernal, Optimal discrete to continuous transfer for band limited inputs, J. Eng. Mech. ASCE 133 (12) (2007) 1370–1377. D. Bernal, A. Ussia, O. Bursi, On Line and Off-line Identification of Inputs, Eurodyn, Porto, Portugal (2014), pp. , 2014, 3053–3056. D. Bernal, On the Identification of Inputs, ICEDyn, Faro, Portugal, 2015. G. Golub, M. Heath, G. Wahba, Generalized cross-validation as a method for choosing a good ridge parameter, Econometrics (1979). C. Hansen, Regularization Tools: A Matlab Package for Analysis and Solution of Discrete Ill-posed Problems, 2008. M. Hou, P.C. Müller, Design of observers for linear systems with unknown inputs, IEEE Trans. Autom. Control 37 (6) (1992) 871–875.