Identification of unknown , time-varying forces / moments in dynamics and vibration problems using a new approach to deconvolution

In this paper an alternative approach to the classical deconvolution idea is used to obtain a new and practical method for real-time identification of unknown, time-varying forces/moments in a general class of linear (linearized) dynamics and vibration problems with multiple-inputs and multiple-measurements. This new method for force/moment identification is unique in the respect that the uncertainty in the force/moment time-variations is not characterized by random-process methods, but rather by a generalized splinemodel with totally unknown weighting coefficients and completely known basis-functions. The basis-functions are custom chosen in each application to reflect, qualitatively, the known characteristics of the force/moment time-variations to be identified. The method does not involve explicit identification of the unknown weighting coefficients. Generalpurpose identification algorithms for both continuous-time and discrete-time measurements are developed, and a worked example including computer simulation results is presented.


Introduction
The problem of identifying or estimating the forces/ moments that acted on a dynamic system to produce an observed system response is of interest in many areas of dynamics and vibrations and has become an active research topic in recent years [5,27,29].Conceptually, this problem can be neatly formulated as a classical deconvolution problem but the solu-tion, using conventional deconvolution ideas, is fraught with numerous computational difficulties [2,25,28].In this paper, a recently published [18] alternative approach to the formulation and solution of deconvolution problems is used to obtain a practical, general method (general-purpose algorithm) for processing dynamic system response measurements to obtain realtime identification (estimation) of the unknown, timevarying forces/moments that produced the response.Our results are applicable to a broad class of practical dynamics and vibration problems involving multiple-(unknown) inputs, multiple-measurements and linear (linearized) equations of motion.Identification algorithms for both continuous-time and discrete-time measurements are presented, and a numerical example is worked.

The idea of "waveform-structured" variations in uncertain, time-varying forces/moments
There is a strong tendency in science and engineering to label as "random" any variations that are not accurately known a priori, and to model those variations in terms of probabilistic characterizations [30].However, the unknown force/moment variations of interest in "input-identification" problems in dynamics and vibrations are typically not totally random or "arbitrary" functions of time but rather, in each application, belong to some restricted class of time-functions that are related to an underlying physical process and have distinguishable patterns of characteristic waveform behavior, at least over short intervals of time.For instance, the unknown force (pressure) variations f (t) on missile launchers, due to launch "back-blasts", or on submerged structures due to underwater explosion "shocks", are not random variations but, in fact, Shock and Vibration 5 (1998) [181][182][183][184][185][186][187][188][189][190][191][192][193][194][195][196][197] ISSN 1070-9622 / $8.00 © 1998, IOS Press.All rights reserved are known to have a characteristic "pulse" pattern of time behavior that can be represented qualitatively by a mathematical expression of the form [6] (1) where the value of the positive integer k is determined by the application specifics, the set of weightingcoefficients {C 1 , . . ., C k } are totally unknown "constants", that can abruptly jump in value from time-totime, and the characteristic decay parameter α > 0 is essentially a known constant whose value is application specific.Unknown "recoil-forces" f (t) associated with weapon firings are also closely modeled by Eq. (1).Likewise, unknown force variations whose characteristic time-behavior is known to consist of uncertain, weighted combinations of step, ramp (linearin-time), and acceleration (∼ t 2 ) "modes" of behavior can be qualitatively represented by an expression of the form where values of the "constant" weighting coefficients {C 1 , C 2 , C 3 } are totally unknown and may occasionally jump in value in an uncertain manner.Expressions of the type Eqs (1), (2) serve to characterize, qualitatively, distinctive features of the waveform shape of the otherwise uncertain time-variations in f (t) and are examples of what we will call variations with "waveform-structure".Some additional examples of waveform-structured uncertain variations f (t) are: steady-state, damped or growing sinusoidal variations with known frequency and unknown amplitude and phase, growing or decaying exponentials, and arbitrary linear combinations of all the above.A-priori knowledge of the type Eqs (1), (2), regarding the class of unknown functions f (t) one is dealing with, can lead to important simplifications in force/moment identification procedures using deconvolution ideas.To demonstrate, generally, how this simplification occurs we will consider a broad class of unknown, waveform-structured force/moment timevariations f (t) that can be represented by an expression of the form where the functions {φ 1 (t), φ 2 (t), . . ., φ M (t)}, are chosen by the user to reflect, qualitatively, the known waveform-characteristics of the f (t) associated with each application and are therefore completely known.
The "constant" weighting coefficients {C 1 , . . ., C M } in Eq. ( 3) are totally unknown and may occasionally jump in value; the latter sparse-in-time jump-behavior of the C i is hereafter referred to as "stepwise-constant" and will be discussed in more detail in the Remarks below Eq. (12).Thus the set of functions {φ i } M 1 in Eq. ( 3) constitutes a known (finite) basis set for the unknown force/moment time-variations f (t) and the values of the set {C i } M 1 of unknown weighting coefficients determines just how the individual basis functions φ i (t) are weighted and linearly combined, at each moment of time, to produce the actual unknown function f (t).
Mathematical expressions of the form Eq. ( 3) are generalized "spline-models" for unknown functions f (t) and have a 30-year record of successful applications in the design of control systems which automatically reject or "accommodate" unknown, unmeasurable disturbance inputs [7,19].Hereafter Eq. ( 3) will be called a waveform model for unknown, waveformstructured input-variations f (t).Clearly, the φ i (t) in Eq. ( 3) must be linearly independent to form a basis for the variations f (t); however, in the novel deconvolution procedure we will present here, it is not necessary to explicitly identify the C i in Eq. (3) and consequently there is no particular advantage to choosing the φ i (t) as orthogonal functions.
In each practical application the particular basis functions one should use in Eq. (3) can usually be determined by consideration of the qualitative, characteristic features of experimental data, and/or the dynamic nature of the underlying physical process that produces f (t).Alternatively, if f (t) is unfamiliar, or has no distinctive characteristic time-behavior, or has basis functions φ i (t) with parametric uncertainties (i.e., uncertain frequencies, time-constants, etc.), one can proceed as in Johnson [7] and choose the set {φ i } M 1 as the polynomial basis set in which case Eq. (3) then becomes the (M − 1)-th degree polynomial spline model It is well-known [1] that polynomial spline-models of the type Eq. ( 5) can accurately represent a broad class of unknown, meandering functions f (t), even for relatively small values of M in the range 3 M 5.The numerical example presented later in this paper illustrates the effectiveness of the case M = 3 of Eq. ( 5); see Eq. ( 29) and Fig. 1.In fact, the use of polynomial basis-functions Eq. ( 4), with M = 3, 4 to augment an existing basis set {φ i } that involves parametric uncertainties is an effective way to model (account-for) the effects of those parameter uncertainties on the timebehavior of f (t), and in the new input-identification algorithms presented in this paper.
Although the basis functions {φ i (t)} in Eq. ( 3) can be chosen as rather complicated time-functions, it turns-out that, from the deconvolution-realization point-of-view, it is advantageous to choose the φ i (t) from among the class of functions that (individually) satisfy homogeneous, constant coefficient, linear differential equations.In that case, over each timeinterval where the C i in Eq. ( 3) are all constant, the unknown force/moment variations f (t) modeled by Eq. ( 3) also satisfy some known (knowable) ρ-th order, homogeneous, constant coefficient, linear differential equation which can be written as where the constant coefficients (β 1 , . . ., β ρ ) in Eq. ( 6) are independent of the (constant) values of the C i in Eq. ( 3) and are completely determined by the chosen basis set {φ i } M 1 .Expression (3) will be referred to as a linear dynamic waveform model in the case Eq. ( 6).The unknown force/moment variations f (t) that typically arise in practical applications can usually be represented by a linear dynamic waveform model Eqs (3), (6).For instance, Eqs (1), (2), (5) are examples of such models.Moreover, in that case the coefficients β i in Eq. ( 6) can be easily computed by Laplace transform techniques applied to the set {φ i (t)} [13].
In summary, the class of unknown force/moment variations f (t) we will consider in this paper consists of the set of functions f (t) which have waveformstructure and can be effectively represented (modeled) by a spline-type, linear dynamic waveform model of the form Eqs (3), (6) where the {φ i (t)} are custom chosen in each application to reflect qualitatively the characteristic modes of waveform behavior that the unknown time-variations f (t) can exhibit, and the stepwise-constant weighting coefficients C i are totally unknown.For simplicity we will refer to this class of unknown force/moment variations f (t) as "linear dynamic variations".The novel way, Eqs (3), ( 6), of mod-eling the uncertainty in unknown force/moment timevariations f (t) neatly avoids the subtle pitfalls of conventional probabilistic characterizations of uncertainty, as discussed by Kalman [24].In particular, we will show that, by virtue of Eq. ( 6), explicit identification of the unknown weighting coefficients C i in Eq. ( 3) is not required to effectively identify f (t).This feature of our deconvolution method results in a significant reduction in computational complexity.
In the case where two or more unknown inputs f (t) act on a dynamic system simultaneously, it will be assumed that each independent f i (t) has a known, linear dynamic waveform model Eqs (3), (6).Also, in all cases it will be necessary to assume that certain "information-theoretic" technical conditions are satisfied, to guarantee that the f i (t) are indeed "identifiable" from the information in the y(t) measurements; see Eqs (13), (19) in the sequel.

The force/moment identification problem for linear dynamic systems with linear dynamic input variations
The physical systems of interest here are multipleinput/multiple-measurement (MIMM) dynamic systems assumed to be modeled by a given, finite-dimensional, linear (linearized) vector-matrix, "equation of motion" in the standard state-variable format where {A, F , C, G} are known, constant matrices, x = (x 1 , . . ., x n ) is the "state" of the dynamical system (typically, x = an ordered n-tuple of generalized coordinates and associated momenta), ) is the vector of unknown, unmeasurable variations of the forces/moments that act on the system, and y = (y 1 , . . ., y m ) is the vector of system response measurements.The class of force/moment identification problems considered here, can now be precisely stated as follows.
Problem statement.Given: (i) the dynamic system's equation of motion Eq. ( 7), (ii) a linear dynamic waveform model Eqs (3), (6), for each independent element f i (t) of the unknown, unmeasurable force/moment input vector f (t), and (iii) a record of the system response vector y(t) = (y 1 (t), . . ., y m (t)), measured over a positive interval of time t 0 t T f , determine the corresponding time-variation of the input vector f (t), t 0 t T f , assuming the initial value x(t 0 ) in Eq. ( 7) is unknown and x(t) is not measurable.It will be initially assumed that any noise effects associated with the measurements {(y 1 (t), y 2 (t), . . ., y m (t)} are sufficiently low-level to be negligible.The accommodation of non-negligible measurement noise is discussed later in the paper.
It is remarked that (linearized) equations of motion for dynamic systems, involving "input-derivative" terms, df (t)/dt, d 2 f (t)/dt 2 , etc., can always be put into the standard format Eq. ( 7) by introducing special definitions for the state-variables x i (t), as explained in Zadeh and Desoer [32, pp.231-232].

A state model for unknown, linear dynamic force/moment variations
When the uncertain force/moment time-variations f i (t) in Eq. ( 7) each have waveform structure and can be described qualitatively by a linear dynamic waveform model Eq.(3), the associated set of linear differential equations Eq. ( 6) can be used to associate a statevector z (i) = (z i1 , z i2 , . . ., z iρi ) with each unknown input f i (t).In particular, one can choose the individual state-variables z ij as for each independent f i (t), i = 1, 2, . . ., r.Thus, over those intervals of time in which all the associated unknown C ij in Eq. ( 3) are indeed constant in value, the equations of motion for the state-variables z ij (t) in Eq. ( 8) are To mathematically account for the uncertain, once-ina-while jumps that can occur in the C ij in Eq. ( 3), a symbolic representation σ ij (t) for a totally unknown, time-sparse sequence of random-like Dirac impulses (having totally unknown (sparse) arrival times and intensities) can be added to each żij expression in the state equations (9a) to obtain It is remarked that our deconvolution method does not involve "identification" of the σ ij (t).
If the set of state vectors {z (1) , z (2) , . . ., z (r) } associated with the vector f = (f 1 , f 2 , . . ., f r ) of unknown inputs in Eq. ( 7) are "stacked" to form one large combined state-vector z = (z (1) |z (2) | • • • |z (r) ) the associated state equations Eq. (9b) can be combined in a similar manner to obtain one vector-matrix "equation of motion" for the combined input state z(t) in the form where and (D, H) have the block structure and where each block D i is a known (knowable) ρ i ×ρ i companion matrix, and each H i a known ρ i -dim.row vector, given by , In summary, the vector f (t) of unknown waveformstructured inputs f i (t) in Eq. ( 7), with a known linear dynamic waveform model Eqs (3), (6) for each independent component f i (t) can be represented equivalently in the state-model format Eq. ( 10) where (D, H) are known, the vector z(t) is the "state" of f (t) and each component σ i (t) of σ(t) is a symbolic representation for a totally unknown, time-sparse sequence of random-like Dirac impulses which "cause" the sparsein-time jumps that naturally occur in the "constants" C ij in Eq. ( 3), associated with each f i (t).

Representation of coupling-effects in the f i (t).
It is remarked that in some practical applications involving multiple-input forces/moments (f i (t), . . ., f r (t)) one or more of the inputs f i (t) may be dynamicallycoupled with other of the inputs f k (t).In such cases the general representation Eq. (10a) still applies [assuming the "coupling" is linear (linearizable)], but the matrices (D, H) then do not necessarily have the block-diagonal structure indicated in Eqs (10b), (10c).Moreover in some applications, particularly those involving bodies moving in fluids, gases, etc., the unknown dynamic variations in the f i (t) may be coupled-with (influenced by) the kinetic or kinematic state-variables x j (t) of the dynamic system Eq.(7).If this latter coupling is linear (linearizable), a representation of the type Eq.(10a) can still be used to model the uncertain behavior of f (t) by generalizing Eq. (10a) to read: ż = Dz+Mx+σ(t), f (t) = Hz(t) + Lx(t), where the coupling-matrices (L, M) presumably can be determined by appropriate modeling procedures.At the expense of some additional algebraic complexity, this latter generalization of Eq. (10a) can be easily incorporated into the forceidentification algorithms presented here by following the procedures outlined in Johnson [13, pp. 424, 432, 481].For simplicity in the present paper, it is hereafter assumed that f (t) is modeled by Eq. (10a).

Endogenization of the unknown input f (t)
The introduction of a state-vector z(t) and associated state-model Eq. (10) for the unknown, lineardynamic waveform-structured input f (t) in Eq. ( 7) is the first of two key ideas that make the identification of f (t) by our deconvolution technique physically realizable.The second idea is that f (t) in Eq. ( 10) can be incorporated into the model Eq. ( 7) in such a way that f (t) is made to "appear" as an internal-variable of the system Eq.( 7) (called "endogenizing" the input f (t)).This conversion is achieved by first substituting Eq. (10a) into Eq.( 7) to obtain the coupled-models and then introducing the new (n + R)-dimensional "composite state" x = (x|z) so that Eq. ( 11) can be re-written as the single, enlarged, "sparse-impulse forced" linear dynamic system model where Ā, C are completely known and σ(t) is totally unknown.Thus, the general class of force/moment identification problems we are considering here can now be restated in terms of Eq. ( 12) as follows.
Problem restatement.Given the enlarged dynamic system model Eq. ( 12) where Ā, C are known and σ(t) is a vector of totally unknown, time-sparse sequences of Dirac impulses, and given a record of the (vector) response measurement y(t) = Cx(t) + Gf (t) over a positive time-interval t 0 t T f , determine (identify, deconvolve) the corresponding input vector f (t) = Hz(t), t 0 t T f , where x(t 0 ), x(t), and z(t) are unknown and not directly measurable.

Remarks.
It is important to note that, even though f (t) is unknown, the modeling procedure of Eqs ( 3), ( 6), ( 10) leads to a completely known composite model Eq.(12), where the only "unknown" is the impulse sequence σ(t) which models the occasional jumps in the C i and is immaterial to the estimation of f (t) -as long as the σ(t) impulses are sparse in time.A more precise technical explanation of the term "sparse", as used here, is given below.
It is clear from the totally unknown nature of the presumed sparse, random-like jumps in the C i in Eq. ( 3) that it is physically impossible for any deconvolution algorithm to generate accurate estimates f(t) of f (t) at the (sparse, isolated) moments of time t k when the C i in Eq. ( 3) abruptly jump in value.Moreover, the speed at which f (t) → f (t) after each such jump will be determined by the algorithm's convergence rate which, for the new algorithms we will propose here, can be controlled by algorithm parameter-values chosen by the user, subject to implicit constraints related to measurement "noises" as discussed later in this paper.However, in practice that convergence rate will always be finite and therefore to achieve f(t) ≈ f (t) "most" of the time it is essential to invoke the tacit assumption that the physical process which produces f (t) is such that the associated, unknown C i in Eq. ( 3) do not jump in value too frequently (i.e., those jumps turn-out to be sufficiently sparse-in-time, relative to the algorithm's designed settling-time for f(t) → f (t)).Otherwise, the convergence f (t) → f (t) will be thwarted by the relentless jumping of the C i and, apparently, no deconvolution procedure can accurately identify f (t) in such cases; however, see Johnson [13, pp. 419, 420] for some ways to mitigate the relentless jumping of the C i , when it occurs, by using alternative basis-functions φ i (t) in Eq. ( 4).
It is remarked that the inability to find a waveformmodel Eq. ( 3) that can adequately represent the behavior of an uncertain function f (t) -with C i that do not rapidly jump in value -is a practical way to conclude that f (t) should be modeled instead as a conventional random process.It is recalled that in random-process methodologies the individual timebehavior of any one random function f (t) (i.e., any one "sample-function") cannot be effectively characterized; rather, it is only possible to characterize certain long-term, time-averaged features of a family (ensemble) of "similar" f (t), such as mean-value, variance, power-spectral density, etc., [30].The latter methodologies are often based on the tacit presumption that the ensemble of sample-functions {f (t)} obeys the ergodic hypothesis -an elusive assumption that can be difficult to verify in practical applications.In contrast, when the waveform-model Eq. ( 3) and associated identification methodology as used in this paper is applicable, it applies equally well to families of unknown f (t) that are non-ergodic and allows qualitative characterization of the time-behavior of each sample-function f (t) in a family.
We will now present two physically realizable algorithms that can accomplish the deconvolution process required to identify unknown, linear-dynamic input variations f (t) in real time, from real-time measurements of the response vector y(t), t 0 t T f , in (7).

Algorithms for identification of f (t) by an asymptotic deconvolution process
In a little-known, limited circulation, but highly significant 1963 document [4], the mathematician/physicist R.W. Bass developed a novel algorithm (filter) for generating approximations of higher-order timederivatives of measured time-varying signals s(t).Bass called his algorithm "an asymptotic polynomial differentiator" [4, p. 19] and envisioned it would be useful in aerospace applications.According to Kalman [23], that algorithm was later rediscovered, refined, and generalized (by Luenberger [26] and other researchers) and has come to be called a state-observer, or state-estimator, in today's terminology.The deconvolution algorithms we are about to describe, for solving the force/moment identification problem, are based on those state-observer ideas as they have been adapted to solve the "disturbance-observer" problem in Disturbance-Accommodating Control Theory [7,9,13,14], and to solve a broad class of "inverse-system" realization problems [18].
An (n + R)-th order deconvolution algorithm for force/moment identification.The deconvolution algorithm presented in this section is analog (i.e., continuous-time) in nature and is intended for processing realtime analog measurement signals {y 1 (t), y 2 (t), . . ., y m (t)} in Eq. ( 7) to generate real-time analog estimates { f 1 (t), f 2 (t), . . ., f r (t)} of the unknown, unmeasurable force/moment time-variations {f 1 (t), . . ., f r (t)} in Eq. ( 7).The algorithm is entirely linear in structure and is composed of combinations of (analog) integration with respect to time, multiplication by constants, and summations, and is compactly defined in the form of a set of linear differential and algebraic equations of the type Eq. ( 12).The elegance of this algorithm makes it fitting to present it as a theorem.An all-digital (digital-signal-processing) version of this algorithm is presented later in this paper; see Theorem 3 below.Theorem 1.The vector of unknown force/moment variations f (t) in Eq. ( 7) is identifiable (i.e., f (t) can be identified or estimated), between successive, sparse jumps of the associated C ij in Eq. (3), from the vector of response measurements y(t) if, in Eq. (12), where (•) T denotes transpose.Moreover, in that case the estimate f (t) of f (t) can be generated in realtime by processing the vector y(t) of response measurements in the (analog) deconvolution algorithm defined by the set of linear algebraic and differential equations [compare to Eq. (12)] where x(0) is to be chosen, see remarks below Eq. (17), and K in Eq. (14b) denotes a constant (n + R) × m "gain" matrix to be chosen (designed) to make the real-part of each eigenvalue λ i of ( Ā + K C) have a sufficiently "large" negative value, in accordance with the required rate of convergence of f (t) → f (t) as determined by the application requirements.
Proof.The real-time identification error ε f (t) associated with the estimate f (t) can be expressed as.
Applying Eqs ( 12), (14), between successive arrivals of the (time-sparse) impulses in σ(t), it follows that the error variable ε(t) is governed by the linear homogeneous differential equation It is clear from Eq. ( 16) that ε(t) → 0 [and hence ε f (t) → 0, also] arbitrarily fast, between successive impulses in σ(t), only if the matrix K can be chosen such that the real-part of each eigenvalue λ i of [ Ā + K C] can be made an arbitrarily "large" negative value.It is well-known [31] that, for a given ( Ā, C), Eq. ( 13) is the necessary and sufficient condition for the existence of such a K.

Remarks.
As indicated, Eq. ( 13) is a simple, practical sufficient condition for existence of an algorithm which can identify f (t) from the y(t) data.The more precise necessary and sufficient conditions for the "identifiability" of f (t) are somewhat more complicated than Eq. ( 13) and are derived in Johnson [20].Note that our deconvolution algorithm Eq. ( 14) generates the estimate f (t) = H z(t) without explicitly identifying the unknown, stepwise-constant weighting coefficients C ij in Eq. ( 3).The problem of designing the (n + R) × m matrix K in Eq. ( 16) to "place" the eigenvalues λ i of [ Ā + K C] at locations sufficiently "deep" in the left-half of the complex plane can be approached by either of two ways.If (n+R) is relatively small, say (n+R) 7, it may be easier to employ the brute-force method of first picking the desired numerical values for the λ i , then setting and then solving Eq. ( 17) for the elements k ij of K. Unless m = 1, the matrix K that satisfies Eq. ( 17) is highly non-unique, in general.On the other hand, if (n + R) is relatively large it may be easier to first transpose [ Ā+K C], then transform ( ĀT , CT ) to a generalized block companion canonical form, such as de-scribed in Johnson [10], and then employ Eq. ( 17) with [ Ā+K C] in Eq. ( 17) replaced by its transformed transpose [ ĀT + CT K T ]; see Johnson ([10], Section 5) and recall that a matrix and its transpose share the same eigenvalues.With the numerical choice of K so determined (and transformed back), the deconvolution algorithm Eq. ( 14) can be physically implemented using computer-based numerical-integration programs or by op-amp type, linear analog circuitry.This latter circuit will employ (n + R) op-amps configured as "integrators".Assuming K is properly designed, the estimate f(t) produced by the algorithm Eq. ( 14) will converge f(t) → f (t), for any initial-condition value of ˆ x(0) = ( x(0)| z(0)) in Eq. (14b).Thus, it is attractive and convenient to simply set ˆ x(0) = (0|0).However, given a non-zero value of the response initial-condition y(0), the physically-knowable "best" choice for ˆ x(0) in Eq. ( 14b) is a certain non-zero value that depends linearly on y(0), as shown in Johnson [17].The matrix C in Eq. ( 12) reflects the physical positioning (placement) of the response-measurement sensors that produce the data y(t), and that sensor placement clearly can affect the critical rank condition, Eq. ( 13).A mathematical theory for optimal sensor placement (optimal choice of C) to "maximize" both the rank condition, Eq. ( 13), and the robustness of that rank against parameter-variations, is developed in Johnson [8].
An alternative form of (analog) algorithm Eq. ( 14), that requires only (n+R− m) integrators, m = rank C, will now be described.
An (n + R − m)-th order deconvolution algorithm for force/moment identification.The algorithm Eq. ( 14) requires (n + R) integration operations because it necessarily estimates each component x i (t) of the n-dimensional state x(t) of the dynamic system Eq.( 7) and each component z i (t) of the Rdimensional state z(t) of the unknown input f (t).Since the m-vector of response measurements y(t) represents a directly measurable "part" of the composite system state-vector x(t) in Eq. ( 12), it is redundant to "estimate" that measurable part of x(t) in the deconvolution process.Moreover, if rank C = m < m, (m − m) of the measurements {y 1 (t), y 2 (t), . . ., y m (t)} in Eq. ( 12) can be expressed as weighted linear combinations of the remaining m measurements.Consequently it is always possible to express the original measurement set (y 1 , . . ., y m ) in terms of an m-dimensional linearly-independent subset of measurements (y 1 , . . ., y m) which can in-turn be written in the form Eq. (12a) where then dim.C = m × (n + R) and rank C = m.For this reason, no generality is lost by assuming at the outset that rank C = maximum = m in the representation Eqs (11), (12).In that case the measurement vector y(t) constitutes an m-dimensional subset (sub-state) of the (n + R)-dimensional composite system state x(t) in Eq. (12).Then, as a consequence of Luenberger's reduced-order observer theory [26] applied to Eq. ( 12), there exists, in principle, an alternative deconvolution algorithm of the type Eq. ( 14) that requires only (n + R − m) integration operations, and therefore has less "computational inertia", enabling faster convergence.The general mathematical description of that reduced-order algorithm involves somewhat more matrix-algebra than Eq. ( 14), but is equally elegant and worthy of being stated as a theorem, as we shall now do.An all-digital version of this reduced-order deconvolution algorithm is presented later in Theorem 4 of this paper.Theorem 2. Suppose rank C = m = dim.y in Eq. (12a).Then, between successive jumps of the associated C ij in Eq. (3), the vector of unknown force/moment variations f (t) in Eq. ( 7) is identifiable, from the vector of response measurements y(t), by the (n + R − m)-th order deconvolution algorithm where and where T 12 , T 22 are, respectively, any n × provided the (n + R − m) × m matrix Σ in Eqs (18a), (18b), (18d) can be designed to make all eigenvalues λ i of (D+ΣH) have sufficiently "large" negative realparts.The latter is possible if Proof.As a consequence of Eqs (18e), (18g) the linear transformation associated with Eqs ( 11), ( 12) is non-singular for all choices of Σ.In fact, the analytical inverse of Eq. ( 20) is given explicitly by Let ξ(t) be any solution of Eq. (18b) and, using Eq. ( 20), define the associated estimates x(t), z(t), f(t) of x(t), z(t), f (t) = Hz(t) in Eq. ( 11) as It follows that where The equation of motion governing ε(t) is obtained by computing ε in Eq. (23b), using Eqs ( 21), ( 11), (18b)-(18g), to obtain the linear differential equation Thus, between successive jumps of the C ij in Eq. (3) [i.e., on the time-intervals where σ(t) ≡ 0 in Eq. ( 24)], the estimate f (t) given by Eq. (18a) approaches the true value f (t) = Hz(t) in accordance with the expression f (t) − f (t) = HT 22 ε(t) where, by Eq. ( 24), According to Eq. ( 19), the matrix Σ can be chosen by the user to make all eigenvalues of [D + ΣH] have "large" negative real-parts.In that case, it follows from Eq. ( 25) that ε(t) → 0 promptly, thus assuring that f (t) → f (t).
Remarks.As in the case of Theorem 1 and Eq. ( 13), the condition Eq. ( 19) is a simple, practical sufficient condition for identifiability of f (t) by the reducedorder algorithm Eq. ( 18).The more complicated necessary and sufficient condition for identifiability of f (t) is developed in Johnson [20].The design of Σ in Eq. ( 25) and the implementation of Eq. ( 18) is accomplished by using the same techniques described below the proof of Theorem 1.In the absence of any reliable information about the true value of f (0) in Eq. ( 7), the "best" choice for the initial-condition value ξ(0) in Eq. ( 18b) is zero [17, p. 864].

Solution and simulation of a numerical example
To illustrate the effectiveness of the unknown force/ moment identification (deconvolution) technique presented in this paper, consider the generic one degree of freedom, damped spring-mass system (rectilinear or rotational) modeled by where (k, c, M ) are known parameters, f (t) denotes unknown, unmeasurable force/moment variations and y = x 1 denotes the single, continuous-time response measurement.An n = 2, R = 1, m = 1 state-model Eq. ( 7) for this dynamic system can be obtained by defining the additional state-variable x 2 = ẋ1 so that Thus, in Eq. ( 7) It is assumed that the qualitative characteristics of the unknown time-variations of f (t) in Eq. ( 26) are known to consist of (or to be closely approximated by) weighted linear combinations (superpositions) of "constant", "ramp" and "acceleration" modes of timebehavior.In other words, the characteristic waveform behavior of f (t) can be modeled by the quadratic polynomial-spline expression where {C 1 , C 2 , C 3 } are totally unknown "constants" that may abruptly jump in value in a time-sparse, oncein-a-while (stepwise-constant) manner.It is clear from Eq. ( 29) that, between successive jumps in the C i , the unknown input f (t) satisfies the linear differential equation (compare to Eq. ( 6)) Thus R = ρ = 3 in Eq. ( 6) and f (t) can be modeled by a state-type equation of motion Eq. ( 10) where z 1 = f , z 2 = ḟ , z 3 = f and D, H are the known matrices For this example, an (n + R − m) = (2 + 3 − 1) = 4-th order continuous-time deconvolution algorithm of the reduced-order type Eq. ( 18) can be designed as follows.Using Eqs (12b), ( 28), (31), the matrices T 11 , T 21 defined by Eq. (18e) are computed to be and the (non-unique) matrices T 12 , T 22 defined by Eq. (18g) can be chosen as Consequently, Eq. (18f) yields The matrices (D, H) can now be computed from Eq. (18c), using Eqs (12b), ( 33), (34), to obtain It is easily checked that Eq. ( 19) is satisfied for the (D, H) in Eq. ( 35) and thus f (t) can indeed be identified by the deconvolution algorithm Eq. ( 18), from the "information" encoded in the response y(t).
For this example, the gain-matrix Σ in Eq. ( 18) turns-out to be a column 4-vector so that the error dynamics Eq. ( 25) associated with estimation of f (t) is governed by the matrix The characteristic polynomial P(λ) of Eq. ( 36) is calculated to be Although not required by Theorem 2, it is performanceeffective, and computationally convenient, to design the elements of Σ to "place" (locate) all 4 eigenvalues {λ i } 4 1 of Eq. ( 37) at a common location on the negative real-axis, i.e., where h = real-value > 0 is chosen sufficiently large.For that particular choice, the associated "desired" characteristic polynomial P d (λ) is given by The required values of the elements Σ i in Eqs (36), (37) are then determined by equating Eq. (38) to Eq. ( 37) to obtain the gain-design equations where one should choose the numerical value of h > 0 sufficiently large to achieve a satisfactory rate of convergence f(t) → f (t), as determined by Eqs ( 23)-( 25) and the application requirements.Finally, computing Ψ by Eq. ( 18d), the main algorithm equations Eqs (18a), (18b) for this example are calculated to be where, in the absence of any reliable knowledge of the value f (0) in Eq. ( 26), the initial-condition values ξ i (0) in Eq. (40b) should all be set to zero.This completes the design of the force/moment identification algorithm for this example.Some representative results of a simulation exercise of the system Eq.( 26) with a simulated unknown input f (t), and using a simulation of the deconvolution algorithm Eq. ( 40), with non-dimensional parameter values M = 1, c = 0.2, k = 100, and Σ chosen to "place" all four eigenvalues of (D + ΣH) at λ 1 = λ 2 = λ 3 = λ 4 = −5, are shown in Fig. 1.It is clear from Fig. 1 that, even though the actual function f (t) is not a quadratic polynomial in time, the estimate f(t) based on the quadratic-spline model Eq. ( 29) and the associated deconvolution algorithm Eq. ( 40) is an accurate, real-time identification of the unknown input f (t), except in a small right-neighborhood of t = 0 and at the other sparse, isolated times t k where the "constants" C i in Eq. ( 29) jump in value, thereby introducing a sequence of real-time "initial-conditions" on ε(t) in Eq. ( 24) which thereafter rapidly decay toward zero.Those brief, transient errors in f(t), occurring in the right neighborhoods of the t k , can be further reduced, in principle, by re-locating the λ i of (D + ΣH) to positions "deeper" into the left-half of the complex-plane (i.e., Re (λ i ) < −5).This latter step will generally result in an increase in Σ .Owing to the way Σ multiplies y(t) in Eqs (18a), (18b), such progressive increases in Σ can eventually become intolerable if the response measurement y(t) in Eq. ( 27) is "noisy", thus leading to "optimal" compromises and trade-offs in the practical design of Σ (i.e., choice of the λ i of D+ΣH); see Section 9 entitled "Other generalizations".
The simulation results shown in Fig. 1 clearly demonstrate that the new force-identification algorithms proposed in this paper can accurately identify unknown-inputs f (t) even when f (t) is loosely modeled by a quadratic polynomial spline Eq. ( 29) and when the system's response measurement-data y(t) = x 1 (t) in Eq. ( 26) contains a strong component of structural-oscillations. Similar accurate estimates f (t) of unknown input-variations can be obtained in more complex examples/applications where the system's response measurement-data y(t) may contain many different transient and/or steady-state modes of structuraloscillations (or other response modes) -provided the essential characteristics (frequencies, time-constants, etc.) of each such response mode are properly represented in the state-vector x and A matrix of the dynamic system model, Eq. (7).For instance, complex, unknown periodic inputs f (t), that produce resonant structural-oscillations in the response measurement y(t), can be effectively identified by the algorithms de-veloped in this paper, provided all relevant modal frequencies are represented in the system model Eq. ( 7).
Approximate digital-versions of Eqs ( 14), ( 18) can be developed by replacing the derivatives ˙ x and ˙ ξ by suitable finite-difference approximations to obtain a difference-equation algorithm that operates on a sequence of periodically sampled-values y(kT ), k = 0, 1, 2, 3, . . .; T = sample-period > 0, of the system response measurement y(t) in Eq. ( 7) to generate a corresponding sequence of discrete-time estimates f (kT ) of the unknown input f (t).Due to the derivative approximations involved, the accuracy of such algorithms tends to suffer unless the value of T > 0 is chosen sufficiently small.On the other hand, "exact" digital-versions of Eqs ( 14), (18), which can produce theoretically exact estimates f(kT ) for all values of the sample-period T > 0, can be derived by starting with exact discretetime versions of the system dynamics model Eq. ( 7) and unknown-input model Eq.(10).The final result is that one then obtains exact difference equation counterparts of Eq. ( 14) and of Eq. ( 18).Procedures for deriving such digital-versions of the algorithms Eqs ( 14), (18) are discussed in Johnson [14, pp. 242, 245], and the essential results are summarized in the following subsections.

The exact digital version of algorithm (14)
The exact digital-version of the force-identification algorithm Eq. ( 14) is based on the exact equation counterparts of: (i) the underlying physical system equation of motion Eq. ( 7) and (ii) the state model Eq. ( 10) of the unknown input f (t), assuming the vector y(t) of system response measurements is periodically sampled every T units of time to obtain the sequence of "sampled-data" measurement vectors {y(0), y(T ), y(2T ), y(3T ), . .., y(kT ), . ..},T = constant > 0, hereafter denoted by y(kT ), k = 0, 1, 2, . . . .Those exact difference equation counterparts of Eqs ( 7), (10), which model the "discrete-time" evolution of x(t), z(t), for t = kT , k = 0, 1, 2, . . ., are derived in Johnson [14, p. 237] and turn-out to be x((k + 1)T ) = Ax(kT ) + F Hz(kT ) + γ, (41a) where and The difference equation models Eqs (41a), (41c) describe how the exact "next" values x((k + 1)T ), z((k + 1)T ) are related to the current (real-time) values of x(kT ) and z(kT ) and to the values of the terms γ, σ.However, at the current time t = kT , the values of γ, σ cannot be determined because, according to Eq. (41e), (41f), those values depend on the (unknown, unpredictable) "future behavior" of the sparse sequence of Dirac impulses σ(t), over the interval kT t (k + 1)T .Thus, since γ, σ in Eq. ( 41) are never "knowable" at t = kT , the only useable part of Eqs (41a)-(41c) is Note that Eq. (42) becomes "exact" at each t = kT where ( γ, σ) vanish; i.e., where the sparse impulse sequence σ(t) turns-out to be zero over the corresponding sample-interval kT t (k + 1)T .This latter fact is an important consideration in the users choice of numerical value for the data sample-period T > 0, relative to the sparseness of the σ(t) impulses.That is, to obtain reliable estimates f(t) of f (t) using Eq.(42), the value of T should be chosen sufficiently smaller than the average time-spacing between the σ(t) impulses so that Eq. (42) becomes exact for "most" of the k values, k = 0, 1, 2, . . ., associated with the total intervalof-time during which measurements are taken in each particular application.
The exact digital-version of the force-identification algorithm Eq. ( 14) is based on Eq. (42) and can be concisely described as follows.
Theorem 3. Assuming the y(t) data sample-period T is significantly smaller than the average time-spacing between arrivals of the sparse impulses of σ(t) in Eq. (10), the vector f (t) of unknown force/moment variations in Eq. ( 7) is identifiable in discrete-time t = kT , k = 0, 1, 2, . . ., from the vector of sampleddata response measurements y(kT where Moreover, in that case the discrete-time estimate f (kT ) of f (kT ) can be generated in real-time by processing the sampled-measurements y(kT ) in the digital deconvolution algorithm defined by the set of linear algebraic and difference equations where K d in Eq. (44b) denotes an (n+R)×m constant "gain" matrix to be chosen (designed) to make each eigenvalue λ i of ( Ād + K d C) lie inside the unit-radius circle in the complex-plane, and sufficiently close to the origin, in accordance with the required rate of convergence f (kT ) → f (kT ).
Then, the (discrete-time) real-time identification error ε f (kT ) = f (kT ) − f (kT ) can be expressed as ε f (kT ) = [O|H] ε(kT ) where ε = x − ˆ x, just as in Eq. (15).Using Eqs (42), (44), the difference equation governing ε(kT ) is determined to be It follows from the stability theory for linear homogeneous, constant-coefficient difference equations of the form Eq. (45) that ε(kT ) → 0 [and hence ε f (kT ) → 0, also] arbitrarily fast (in discrete-time) only if the matrix K d can be chosen such that the modulus | λ i | of each eigenvalue λ i of ( Ād + K d C) can be made arbitrarily small.It is well-known [31] that, for arbitrarily given ( Ād , C), Eq. (43a) is the necessary and sufficient condition for the existence of such a K d .Like Eq. ( 13), the condition in Eq. ( 43a) is a simple, practical sufficient condition for identifiability of f (t) by the algorithm Eq. (44); see Johnson [22] for the more complicated necessary and sufficient conditions.It is important to note that the matrix Ād in Eq. ( 43) is an explicit function of the users chosen value for the sample-period T ; see Eq. (41d).Thus satisfaction of the sufficient condition Eq. (43a) will depend, in part, on the chosen value of T .The same is true for the more-complicated necessary and sufficient conditions developed in Johnson [22], and is a natural con-sequence of information-losses that can occur due to "sampling" of dynamic data at certain inappropriate frequencies (1/T ).
Remarks.The design of the (n + R) × m matrix K d to place all eigenvalues λ i of ( Ād + K d C) sufficiently close to zero, and the design of the "best" choice for the initial-condition of Eq. (44b), can be accomplished by the same methods previously described in the remarks below the proof of Theorem 1.The design of K d to place all the λ i at λ i = 0 is a particularly attractive option because then, for any initial-condition ε(0), the solution ε(kT ) of Eq. ( 45) [and hence ε f (kT ) also] reaches zero, and stays there, in a finite number N of sample-periods T ; i.e., then ε(NT ) = 0 for any ε(0).This desirably rapid convergence of f(kT ) → f (kT ) in the finite time-interval [0, NT ] is called "deadbeat response" of the identification algorithm Eq. (44).Linear (vector) differenceequation algorithms of the type Eq. ( 44) can be easily programmed into a general-purpose digital computer or special-purpose digital signal processing (DSP) circuit/chip to achieve essentially real-time generation of the discrete-time force/moment estimate f(kT ), k = 0, 1, 2, . . . .In practice, the sequence of discrete-time estimates f(kT ) generated by Eq. (44) produces a "stairstep"type time plot which may not be smooth enough for some application requirements.For such cases, the algorithm equations Eq. (44) can be augmented by an additional equation which eliminates that stairstep characteristic by producing a smooth, natural-waveform estimation interpolation f (τ ), kT τ < (k +1)T , between each successive pair of sample-times.As shown in Johnson [18, pp. 516-517], that interpolation equation, based on the homogeneous part of the model Eq. ( 10), is where the symbol I in Eq. ( 46) denotes the R×R identity matrix.As mentioned below Eq. (42), impulses of σ(t) in Eq. ( 10) that happen to arrive between a particular pair of (successive) sample-times will introduce unavoidable errors in the corresponding "next" estimate f(kT ) generated by Eq. (44), and the same comment applies also to Eq. ( 46).However, such errors do not propagate beyond the next sample-interval over which σ(t) remains "quiet".Thus, if σ(t) is sparse and T is chosen sufficiently small, such quiet-intervals will occur frequently and consequently Eq. (46) will then produce reliably accurate interpolation estimates f(τ ) "most" of the time.The exact digital version of Eq. ( 18) will now be presented.

The exact digital version of algorithm (18)
The exact digital-version of the (n + R − m)-th order continuous-time force-identification algorithm Eq. ( 18) is based on the "exact" system (and unknowninput) difference-equation models Eq. ( 42), with A, F H, D defined as in Eq. (41d).Under the assumptions stated above Eq.( 41) and below Eq. ( 42) that digitalversion of Eq. ( 18) can be concisely stated as follows; see Johnson [14, pp. 244-245] for details.

Theorem 4. Under the same hypothesis associated with Theorem 3, and satisfaction of the (sufficient) condition rank H
where where the matrices Ād , {T 11 , T 12 , T 21 , T 22 , T12 , T22 } are defined as in Eqs (43b), (18e)-(18g) and the matrix Σ in Eq. (48) is designed to make each eigenvalue λ i of ( D + Σ H) lie inside the unit-radius circle in the complex-plane, and sufficiently close to the origin, in accordance with the required rate of convergence f (kT ) → f (kT ).
Proof.As in the proof of Theorem 2, introduce the non-singular linear transformation Eq. ( 20) with the explicit inverse Eq. ( 21).Let ξ(kT ) be any solution of Eq. (48b) and, using Eq. ( 20), define the associated discrete-time estimates { x(kT ), z(kT ), f(kT )} of {x(kT ), z(kT ), f (kT ) = Hz(kT )} in Eq. (42) as It follows that where where ξ is defined in Eq. ( 21).The difference equation governing the motion of ε(kT ) in Eq. ( 51) is found by forward-shifting Eqs ( 21) and (51) once, using Eqs (42), (48), to obtain Therefore ε(kT ) → 0, and thus f(kT ) → f (kT ) also, if, and only if, the (n + R − m) × m design matrix Σ is chosen to place all eigenvalues of λ i of ( D + Σ H) inside the unit-circle of the complex plane.The condition Eq. ( 47) is necessary and sufficient for the existence of a Σ that will place all those λ i precisely at the origin, thus assuring ε(kT ) → 0 as fast as possible, in the discrete-time sense.

Other generalizations
In some applications, the system dynamics model Eq. ( 7) includes a controllable input u(t), as well as the unknown, uncontrollable input f (t).Moreover, one or more of the matrices {A, F , C, G} may contain elements that can vary with time in some known manner.In addition, the time-behavior of f (t) as modeled by Eq. ( 10) may involve known time-varying matrices (D, H).Procedures for generalizing Eqs ( 14), (18) to accommodate such cases are presented in Johnson [12,13].
In this paper the unknown forces/moments f (t) have been presumed to vary with respect to the temporal independent variable (time) t.However, the methods and identification procedures used here are equally applicable to problems of identification of unknown spatialvariations in static load-distributions from measured static deflections of beams, plates, shells, structures, etc., where in the latter class of problems the counterpart of the system equations of motion Eq. ( 7) are the well-known "static deflection" differential equations, for beams, plates, shells, etc., and the unknown force/moment static load-distributions to be identified vary with respect to some spatial independent variable.Application of the methodologies described in this paper to the identification of unknown, unmeasurable static load-distributions on beams, from measurements of beam static deflections, is treated in Johnson [21].Generalizations of the identification methods presented here, to accommodate more general distributed-parameter (continua) problems involving unknown force/moment variations and distributions that depend on both temporal and spatial independent variables, i.e., f = f (t, x), are presented in Balas and Johnson [3].

Summary and conclusions
In this paper it has been shown that a broad class of multi-input/multi-measurement, force/moment identification (deconvolution) problems in (linear) dynamics and vibration can be effectively solved by use of a new form of physically-realizable, real-time deconvolution algorithm, based on linear dynamic waveformmodels Eqs (3), ( 6) for time-variations of the vector of unknown inputs f (t) = (f 1 (t), . . ., f r (t)).A worked example, with simulation results, has been presented to illustrate the effectiveness of the proposed method of force-identification. Various generalizations of the class of problems, and deconvolution algorithms, considered here have also been presented, including "exact" digital-based (DSP) versions of the algorithms.
The force/moment identification problem addressed in this paper arises in a variety of practical applications where the unknown, unmeasurable inputvariations f (t) that "caused" an observed system response y(t) must be determined.The capability of the force/moment identification algorithms presented here to generate accurate estimates f(t) in real-time should be useful in contemporary dynamics and vibration applications where reliable estimates f(t) are needed to make real-time control decisions in "loadadaptive", "smart" and "reconfigurable" structures, active vibration-control systems, active isolators, motionabsorbers, etc., as explained and illustrated in the 1973 paper [11].abama, for his technical advice and support of this research effort.Finally, the author is deeply indebted to Prof. H.H.E.Leipholz, University of Waterloo, Ontario, who, 25 years ago, kindly provided the author a professional forum to advocate, and demonstrate, the (then) relatively new idea of applying modern control techniques to active-vibration and structural-control (smart-structures) problems in a keynote presentation at the University of Waterloo 1973 Symposium on Stochastic Problems in Mechanics; see Johnson [11].

)
the discrete-time estimate f (kT ) of f (kT ) can be generated in real-time by processing the sampledmeasurements y(kT ) in the (n + R − m)-th order digital deconvolution algorithm defined by the set of linear algebraic and difference equations f(kT ) = H (T 21 − T 22 Σ)y(kT ) + T 22 ξ(kT ) , (48a)