- Email: [email protected]
Asymptotic approximation method of force reconstruction: Proof of concept
Mechanical Systems and Signal Processing 92 (2017) 39–63
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Asymptotic approximation method of force reconstruction: Proof of concept J. Sanchez a,⇑, H. Benaroya b a b
Sanchez Advanced Engineering Concepts, LLC, South Plainfield, NJ 07080, United States Department of Mechanical and Aerospace Engineering, Rutgers University, 98 Brett Road, Piscataway, NJ 08854, United States
a r t i c l e
i n f o
Article history: Received 24 March 2016 Received in revised form 16 January 2017 Accepted 18 January 2017 Available online 3 February 2017 Keywords: Dynamic force estimation Force reconstruction Asymptotic Approximation
a b s t r a c t An important problem in engineering is the determination of the system input based on the system response. This type of problem is difficult to solve as it is often ill-defined, and produces inaccurate or non-unique results. Current reconstruction techniques typically involve the employment of optimization methods or additional constraints to regularize the problem, but these methods are not without their flaws as they may be suboptimally applied and produce inadequate results. An alternative approach is developed that draws upon concepts from control systems theory, the equilibrium analysis of linear dynamical systems with time-dependent inputs, and asymptotic approximation analysis. This paper presents the theoretical development of the proposed method. A simple application of the method is presented to demonstrate the procedure. A more complex application to a continuous system is performed to demonstrate the applicability of the method. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction In many engineering applications, knowledge of the input or applied force is critical to the analysis of a system or structure. Unfortunately, at times, this knowledge may be difficult to obtain as it may not be possible to measure the applied force directly, and we must rely on measurements of the system response to determine the input or force applied to a system. This presents the problem of developing a model that estimates the input based on the measured output. Various approaches to this problem are outlined in [19], where such reconstruction methods are divided into three categories: direct methods, regularization methods, and probabilistic/statistical methods. Direct methods involve the application of the governing equation explicitly to determine an estimate for the applied force. These methods were employed by Jacquelin and Hamelin [11] to a split Hopkinson pressure bar using strain measurements, Law et al. [10] to an Euler-Bernoulli beam with a moving force applied, and were more generally developed by Karlsson [14], and Avitabile et al. [1]. While these methods are simplest to apply, they generally are ill-conditioned and can produce inadequate results. To overcome this ill-conditioning, additional constraints are often applied to transform the system into a wellconditioned problem, which is known as regularization. Various methods have been developed to regularize a problem, for example, through modal truncation, and optimization. Methods of truncation involve the isolation and removal of characteristics deemed to be attributed to noise, such as mode selection techniques developed by Jiang and Hu [12,13], and the
⇑ Corresponding author. E-mail addresses: [email protected] (J. Sanchez), [email protected] (H. Benaroya). http://dx.doi.org/10.1016/j.ymssp.2017.01.022 0888-3270/Ó 2017 Elsevier Ltd. All rights reserved.
40
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
application of wavenumber filtering by Djamaa et al. [4]. Additionally, the application of optimization methods were presented by Huang et al. [9] that utilized and optimized additional terms to filter noise inherent in the system, and the conjugate gradient method was detailed by Huang [8] and Yen and Wu [22,23]. Another type of optimization solution involves the application of the Levenberg-Marquardt iterative regularization method to solve inverse problems, and were presented by Gunawan [7] and Gasparo et al. [5]. While the application of these methods are generally successful, they do have their drawbacks. Some regularization methods possess additional parameters for regularization that lack a clearly defined method of determination while others require an intense level of computation to produce an estimate of the applied force. In addition to regularization methods, probabilistic/statistical methods have been applied to the solution of the inverse problem. These methods generally incorporate uncertainties such as measurement or system noise into the analysis and allows for a more holistic perspective of the inverse problem. One method that is becoming prevalent is the application of the Kalman filter to the inverse problem, and has been applied in conjunction with a recursive least squares algorithm by Ma et al. [16], by Xu et al. [21], and in an augmented form by Lourens et al. [15]. While the application of these methods have proven successful, like regularization methods, probabilistic methods can be computationally intensive as they may require a significant amount of data and time for computation. Fundamentally, these methods of estimation involve the truncation of system characteristics attributed to noise or the development of an optimization problem to produce a well-defined formulation. The dismissal of these assumptions would allow more freedom in the estimation, and allow us to explore new avenues of analysis that have otherwise gone unexplored. Our objective in developing our dynamic force estimation method was to develop a procedure that is inherently stable independent of any optimization or assumptions on the characteristics of the temporal component of the applied force throughout the recovery process. The work presented in this paper should be seen as a proof-of-concept where the estimation of dynamic forces is viewed from an asymptotic perspective. From this view, we aim to develop a method for estimating forces that has the potential to asymptotically converge to the applied force. To achieve this goal, we will compare the excited system to an unexcited estimate system under the assumptions that the system parameters are known and that sampling is continuous in time. Under these assumptions, and concepts from the following subject areas:
linear dynamical systems theory, control theory, equilibrium analysis of linear dynamical systems, asymptotic approximation analysis,
we will derive a relationship between the excited system, the unexcited system, and the system input. As we are interested in the development and proof of the concepts, we begin the development of this work under ideal circumstances where measurement noise is neglected, and extend these concepts to a real model, which includes measurement noise. Additional details are provided in [18], and while the focus of our applications in this paper are limited to the recovery of forces on mechanical systems, it is noted that the method is applicable to any linear system that satisfies the conditions discussed in this method. In this paper, the preliminary concepts developed for this analysis are summarized in Section 2. The derivation and significant results of the estimation method for an ideal model are discussed in Section 3. A simple example is developed in Section 4, and the application of the method to the estimation of distributed loads on continuous systems by using modal analysis is shown in Section 5. The concepts of the proposed method are applied to a real model in Section 6. Additional discussion and future work is provided in Section 7. 2. Preliminaries The development of the proposed estimation method utilizes concepts from multiple disciplines, but while these methods are critical to the development of the method, linear dynamical systems analysis is the starting point for the development. Throughout this work, we use bold face variables to denote matrices and underlined variables denote vectors. 2.1. Equilibrium analysis of linear dynamical systems The development of the proposed force reconstruction method begins with an analysis of the equilibrium of systems of the form:
x_ ¼ Ax þ Q ðtÞ;
ð1Þ
where Q ðtÞ is a time-dependent vector input, x is an n 1 state vector, and A is an invertible, well-conditioned matrix. Our aim is to transform Eq. (1) to a simplified equation of the form:
z_ ¼ Az;
ð2Þ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
41
where z ¼ x xH , and xH is the equilibrium state. By reducing Eq. (1) to Eq. (2), we aim to derive an expression for the equilibrium state, xH , in terms of the input, Q ðtÞ . We are additionally assuming that Q ðtÞ is infinitely differentiable and varying with time on the interval t 2 ½0; 1Þ. To achieve this, we begin by finding where x_ ¼ 0 for t 2 ½0; 1Þ. This equilibrium xH 0 of the linear dynamical system is given by: 1 xH 0 ¼ A QðtÞ:
ð3Þ
To better see what is happening, we can transform the state vector to be centered around the equilibrium xH 0 by defining _ _ _H the following variables: x1 ¼ x xH 0 and x1 ¼ x x0 . The variables would thus be transformed from x to x1 by substituting the relations for x1 and x_ 1 into Eq. (1), which produces:
_ x_ 1 ¼ Ax1 þ A1 QðtÞ:
ð4Þ
From this procedure, we see something interesting happening. Through the change of coordinates to center Eq. (1) around 1 _ the equilibrium xH 0 , we obtain another system of the same form but Q ðtÞ is replaced by A Q ðtÞ. Thus, for Eq. (3) to hold, Eq. (4) must hold. Additionally, as Eq. (4) is not of the form given by Eq. (2), the simplication of Eq. (1) to the form given by Eq. (2) is incomplete. Comparing Eq. (1) to Eq. (4), this suggests a recursive process. Repeating this process, we can find the equilibrium of Eq. (4), which is given by x_ 1 ¼ 0. This is true when:
2 1 _ xH QðtÞ: 1 ¼ A
ð5Þ
Defining x2 ¼ x1 xH 1 , we can rewrite Eq. (4) in terms of x2 . This produces the following:
€ x_ 2 ¼ Ax2 þ A1 QðtÞ:
ð6Þ
As we can see, Eq. (6) is of the same form of Eqs. (1) and (4). Beginning with Eq. (6), we repeat the procedure of changing the variables n times, with the following relations resulting by induction:
xn ¼ x þ
n1 X mþ1 ððAÞ1 Þ Q ðmÞ
ð7aÞ
m¼0 n
x_ n ¼ Axn þ ðA1 Þ Q ðnÞ ;
ð7bÞ
where ðmÞ denotes the mth order derivative with respect to time and n is a positive integer. Additionally, x_ n can be directly related to x by substitution of Eq. (7a) into Eq. (7b), to produce:
x_ n ¼ A x þ
! n X mþ1 ðA1 Þ Q ðmÞ :
ð8Þ
m¼0
Taking the limit as n goes to infinity, Eqs. (7a) and (8) become:
lim xn ¼ x þ
n!1
1 X mþ1 ðA1 Þ Q ðmÞ
ð9aÞ
m¼0
lim x_ n ¼ A x þ A1
n!1
1 X
! m
ðA1 Þ Q ðmÞ :
ð9bÞ
m¼0
Defining x1 ¼ limn!1 xn , Eqs. (9a) and (9b) become:
x1 ¼ x þ
1 X mþ1 ðA1 Þ Q ðmÞ
ð10aÞ
m¼0
x_ 1 ¼ Ax1 ;
ð10bÞ
and thus we have transformed Eq. (1) into the form of Eq. (2), where z ¼ x1 . The equilibrium of Eq. (10b) is found by setting x_ 1 ¼ 0. Since A is invertible and well-conditioned, the only state in which this is true is x1 ¼ 0. Thus, the equilibrium state of Eq. (10b) is xH 1 ¼ 0. To determine what happens at equilibrium, we set Eq. (10a) equal to xH , which produces the following: 1 1 X mþ1 xH ¼ ðA1 Þ Q ðmÞ :
ð11Þ
m¼0
Thus, the equilibrium state of the linear dynamical system given in Eq. (1) is given by Eq. (11), and the result of this analysis is presented as Theorem 1.
42
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Theorem 1. For a linear dynamical system of the form of Eq. (1) with the matrix A invertible and well-conditioned, and an input Q ðtÞ that is infinitely differentiable, the dynamic equilibrium xH of the system is given by: 1 mþ1 X xH ¼ A1 Q ðmÞ ;
ð12Þ
m¼0 m
where Q ðmÞ ¼ dtd m Q ðtÞ. The series expansion of xH as stated by Theorem 1 is important to the development of the proposed estimation method. It enables us to see the relationship between the system parameters of A and the system input Q ðtÞ, and how they relate to the equilibrium of the dynamical system given by Eq. (1). 2.1.1. Discussion of Theorem 1 The form for the equilibrium of Eq. (1) given by Theorem 1 allows us to transform the dynamical system into the simpler form of Eq. (2). To see this, we begin from the discussion after Eq. (2):
z ¼ x xH 1 X mþ1 ðA1 Þ Q ðmÞ : ¼xþ
ð13aÞ ð13bÞ
m¼0
Substitution of Eq. (13a) into Eq. (1) produces:
x_ ¼ Ax þ Q ðtÞ 1 1 X X m mþ1 z_ ðA1 Þ Q ðmÞ ¼ A z ðA1 Þ Q ðmÞ
z_ z_ z_
m¼1 1 X
ð14bÞ
m
m
ð14cÞ
m
m¼0 1 X
m
ð14dÞ
m
m¼1 1 X
m
ð14eÞ
ðA1 Þ Q ðmÞ ¼ Az
m¼1 1 X
þ QðtÞ
m¼0 1 X
ðA1 Þ Q ðmÞ ¼ Az
m¼1 1 X
ð14aÞ
!
ðA1 Þ Q ðmÞ ¼ Az
m¼1
z_ ¼ Az:
ðA1 Þ Q ðmÞ þ QðtÞ ðA1 Þ Q ðmÞ Q ðtÞ þ Q ðtÞ ðA1 Þ Q ðmÞ
m¼1
ð14fÞ
As we can see in Eqs. (14), the series expansion of the equilibrium state given in Theorem 1 allows for the cancelation of Q ðtÞ on the right hand side of Eq. (1), and the higher order derivatives of the series balance out on each side of the equality. As a result, the system is reduced from an excited to an unexcited linear system. Additionally, the solution for Eq. (14f) is:
z ¼ expðAtÞz0 ;
ð15Þ
where z0 is the initial state of z. Since we know the relationship between z and x, we can rewrite Eq. (15) using Eq. (13a) as:
x xH ¼ expðAtÞz0 :
ð16Þ
From this, we can see that z0 is essentially the initial condition of the transient response for the system given in Eq. (1). This allows us to rewrite Eq. (16) as:
x ¼ expðAtÞz0 þ xH :
ð17Þ
Looking at the right hand side of Eq. (17), we can see in this form that the series expansion of the equilibrium state given by Theorem 1 is a alternative form of the particular solution for nonhomogeneous linear dynamical systems derived using the concept of equilibrium from dynamical systems theory. From these multiple forms, another important conclusion that can be drawn is that the response of the system is related to the input and its infinitely many derivatives, if they exist. It is not just the value of the input that determines the evolution of the system at steady state or equilibrium. 2.1.2. Existence of Theorem 1 While we have shown computationally that Theorem 1 is an expression of the equilibrium, we have yet to address the issue of the convergence of the series representation given by Theorem 1. In this section, we will provide the conditions in which the series from Theorem 1 converges. Since we already assume A in Eq. (1) is invertible and well-conditioned, we will assume, for simplicity, that it is therefore diagonalizable, and we can simplify this analysis by performing this analysis on a scalar system, given by:
43
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
x_ ¼ ax þ QðtÞ:
ð18Þ
Theorem 1 dictates the equilibrium of Eq. (18) is given by: 1 X Q ðmÞ xH ¼ : amþ1 m¼0
ð19Þ
Now we can apply the ratio test to determine the conditions for convergence of the series. This states that an infinite series converges absolutely if:
Q ðmÞ 1=m lim < 1; m!1amþ1
ð20Þ
where jj represents the absolute value. Eq. (20) simplifies to:
1=m lim Q ðmÞ < jaj:
m!1
ð21Þ
Eq. (21) is the general condition for the existence of Theorem 1. While Eq. (21) is a general result, we can produce more explicit results with more assumptions on the form of Q ðtÞ and its derivatives. First, we examine the case in which there exists a constant Bm > 0 such that:
ðmÞ Q 6 Bm
ð22Þ
for all t 2 ½0; 1Þ for each m 2 ½0; 1Þ. Combining Eqs. (22) and (21) produces:
1=m 1=m lim Q ðmÞ 6 lim Bm < jaj:
m!1
m!1
ð23Þ
Thus, for all Q ðtÞ that satisfies Eq. (22), the criterion for the existence of Theorem 1 is given by:
lim B1=m m!1 m
< jaj:
ð24Þ
We can relax the assumption on Bm to exponentially bounded functions. We define the followin:
Bm ¼ C m expðMm tÞ;
ð25Þ
meaning Q ðtÞ and its derivatives are bounded as:
ðmÞ Q 6 C m expðM m tÞ;
ð26Þ
where M m and C m are positive constants for each m 2 ½0; 1Þ. Substituting Eq. (26) into Eq. (24) produces the following: 1=m lim C m exp
m!1
Mm t < jaj: m
ð27Þ
For Eq. (27) to hold, we can conclude the following criteria must be true:
lim C 1=m m!1 m lim
m!1
< jaj
Mm ¼ 0: m
ð28aÞ ð28bÞ
While we have conducted this analysis on a scalar system, it should be noted these results hold for matrices as well. The only difference between the scalar system proof and the proof for higher order systems would be that we are using the eigenvalues of the matrix A throughout the discussion. Otherwise, the proof holds identically for scalar systems and higher order matrix systems. 2.2. Asymptotic approximation analysis Now that we have derived an expression for the equilibrium of Eq. (1), we want to address the question as to what would occur as the eigenvalues of A are varied. In particular, we are most interested in what occurs as the eigenvalues of A tend to infinity. We restrict our analysis to inputs with derivatives of exponential order. This means that for each Q ðmÞ , there exist positive constants, C m and M m , for each m 2 ½0; 1Þ such that:
ðmÞ Q i 6 C m expðM m tÞ for m 2 ½0; 1Þ8t 2 ½0; 1Þ;
ð29Þ
ðmÞ where Q i is representative of the absolute value of each term of the vector, Q ðmÞ . Additionally, we will assume that A is diagonalizable, and as a result, we can simplify this analysis by examining scalar systems, given by:
44
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
x_ ¼ ax þ Q ðtÞ:
ð30Þ
and inputs as we can perform this analysis term by term in the vectors Q ðmÞ and xH . First, we examine what happens to the magnitude of the equilibrium xH . By using Theorem 1, Eq. (29), and the triangle inequality, we find: 1 H X C m expðM m tÞ x 6 ; jajmþ1 m¼0
ð31Þ
where jxH j represents the absolute value of xH . Now, we can examine what happens to each term at a tends to infinity. For all t 2 ½0; 1Þ, each term of the inequality on the right hand side in Eq. (31) tends to zero as jaj goes to infinity. Thus, we can conclude from Eq. (31) that jxH j goes to zero as jaj tends to infinity. This result is summarized in Proposition 1. Proposition 1. For an input Q such that:
ðmÞ Q i 6 C m expðM m tÞ for m 2 ½0; 1Þ;
ð32Þ
xH goes to zero as the eigenvalues of A tend to infinity. Proposition 1 shows what occurs to xH , but we must determine what occurs to the term AxH as the eigenvalues of A go to infinity. To achieve this, we must rearrange Eq. (12) into the following form:
AxH Q ¼
1 X
A1
m
Q ðmÞ ;
ð33Þ
m¼1
and determine what occurs to the magnitude of AxH Q as the eigenvalues of A tend to infinity. Again, the proof is performed for a scalar system and input given by Eq. (30). Thus, Eq. (33) in scalar form reduces to:
axH Q ¼
1 X
a1
m
Q ðmÞ ;
ð34Þ
m¼1
Applying the Triangle Inequality to Eq. (34), we obtain: 1 H X C m expðMm tÞ ax Q 6 : jajm m¼1
ð35Þ
Again, we can see that for all t 2 ½0; 1Þ, each term of the inequality on the right hand side in Eq. (35) goes to zero as jaj goes to infinity. Thus, from Eq. (35), we can see that axH tends to Q as jaj tends to infinity. This result is summarized in Proposition 2. Proposition 2. For an input Q such that:
ðmÞ Q i 6 C m expðM m tÞ for m 2 ½0; 1Þ;
ð36Þ
AxH goes to Q as the eigenvalues of A tend to infinity. It should be noted that Propositions 1 and 2 also apply to inputs with impulsive derivatives, and further analysis of this is found in [18]. 3. Derivation of the proposed method – ideal case We are now able to pull together all the essential aspects of the prior discussion to formulate a method for input reconstruction for the ideal case. We begin with a linear dynamical system of the form:
x_ ¼ Ax þ Q ðtÞ
ð37aÞ
y ¼ Cx;
ð37bÞ
where Q ðtÞ is an n 1 time-dependent input vector, x is an n 1 state vector, A is an n n matrix, C is an m n matrix, and y is the m 1 measured output vector. To find the relationship between the measured output and the input, we define an estimation system with observer feedback, where we neglect the excitation term Q ðtÞ. This is done to draw a comparison between the excited and unexcited systems, and to produce a relationship between the state and the excitation Q ðtÞ, as the result of the following derivation will show. Thus, we define the estimation system as:
^x_ ¼ A^x þ LðC ^x yÞ;
ð38Þ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
45
where ^ x is the state estimate, and L is the n m observer gain matrix, which is a design parameter. With this form, we are able to generate a direct relationship between the states of the actual and estimation systems and the input of the actual system. Substitution of Eq. (37b) into Eq. (38) produces:
^x_ ¼ A^x þ LCð^x xÞ:
ð39Þ
x x, we have the following error system: Defining e ¼ ^
e_ ¼ ðA þ LC Þe QðtÞ:
ð40Þ
For an observable system, we can design the eigenvalues of A þ LC with an appropriately chosen L [17]. As a result, we can design the eigenvalues of A þ LC such that the matrix is well-conditioned and invertible. This would allow us to use Theorem 1 to determine the equilibrium state of the system defined by Eq. (40), which would give us:
eH ¼
1 mþ1 X ðA þ LC Þ1 Q ðmÞ ;
ð41Þ
m¼0
where eH is the equilibrium state of Eq. (40). The parameter L could be designed such that: the real part of the eigenvalues of A þ LC are negative so that the equilibrium state of the error system is asymptotically stable, and the magnitude of the eigenvalues of A þ LC are sufficiently large in magnitude such that: 1 m X ðA þ LC Þ1 Q ðmÞ 0:
ð42Þ
m¼1
We can apply Proposition 2, and as the magnitudes of the eigenvalues of A þ LC go to infinity, we have:
ðA þ LC ÞeH ! QðtÞ;
ð43Þ
^: and we could define our estimate for the input Q
b ¼ ðA þ LC ÞeH : Q
ð44Þ H
Additionally, from Proposition 1, we know that e goes to zero as the eigenvalues of A þ LC tend to infinity. Since the terms of the matrix L are of sufficiently large magnitude, we can also conclude in Eq. (44) that AeH is negligible in relation to LCeH , and we can define our estimate for the input as:
^ ¼ LCeH : Q
ð45Þ
This can be rewritten as:
^ ¼L y ^H yH : Q
ð46Þ
3.1. Force reconstruction procedure Now that we have an expression for the input estimate, we can discuss the process of estimating the input. To better visualize what is occurring in the method, the block diagram in Fig. 1 shows us that the method consists of three fundamental steps. First, the system response y must be measured. Once the measurement is made, this is substituted into the unexcited b to x. Then, the measurement y and the estimate state ^ x are substituted into Q estimate system, resulting in the estimate state ^ determine the estimate for the system input. 3.2. Interpretation The method developed here defines an estimate system with observer feedback, and in that sense, it is comparable to the Luenberger observer. The two methods differ in their assumptions and application. The Luenberger observer assumes knowledge of the input and measured output in order to determine an estimate for the state, but in the proposed method, we know the measured output in order to determine an estimate for the input. As a result, the proposed method is not a Luenberger observer. It may be considered to be a complement to the Luenberger observer, and an extension of the application of observer feedback. Additionally, the new application of observer feedback in the method includes new interpretations for the parameters. The matrix L consists of design parameters and the matrix C acts to filter terms of the error equilibrium, eH , that we cannot
46
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Fig. 1. Block diagram of force reconstruction method.
observe or measure. When functioning in tandem with each other in the matrix product LC, we can say that LC acts as a high gain observer filter to approximate the input. It achieves this in two ways: as a filter of the derivatives of the input, and as a weighting function to assign more value to the observed terms in relation to the unobserved terms. Its ability to act as a filter is apparent in Eq. (33). As the eigenvalues of A þ LC go to infinity, each term in the series with Q ðmÞ for m 2 ½1; 1Þ becomes increasingly negligible. As a result, the terms in the series with the input differentiated are filtered, leaving us with the original input Q ðtÞ. To see the second purpose of L, we examine what is occurring within the matrix A þ LC. The eigenvalues of A þ LC control the magnitude of eH , thus, for larger eigenvalues, eH becomes predictably smaller. And since the matrix L controls the eigenvalues of A þ LC, then the terms of L are large. As L works in tandem with C, and C is related to the measurements, then LC acts to assign more weight to the measurements. As the terms of L become large, the term AeH becomes negligible in the estimation of the input.
4. Example: Single degree–of–freedom system To detail the manner in which the proposed method would be applied, numerical simulations were conducted in [18] for a scalar system, a single degree–of–freedom system, and a 2 degree–of–freedom system. We present here the application of the proposed method to a single degree–of–freedom example for the system shown in Fig. 2. In this example, we assume particular inputs and use these to find the outputs. It is the output that are the ‘‘data” used in the proposed method to reconstruct the inputs.
47
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Fig. 2. Single degree-of-freedom dynamical system.
We select the particular single degree-of-freedom system:
€ þ 4y_ þ 4y ¼ FðtÞ; y
ð47Þ
with the observed displacement:
z ¼ y:
ð48Þ
_ resulting in: Our method requires that we first rewrite the system in the form of Eq. (37) by defining x1 ¼ y and x2 ¼ y,
A¼ Q¼
0
1
ð49aÞ
4 4
0
ð49bÞ
FðtÞ
C ¼ ½1 0
ð49cÞ
z ¼ Cx
x1 : x¼ x2
ð49dÞ ð49eÞ
We can define our estimation system by:
^x_ ¼ A^x þ LðC ^x zÞ;
ð50Þ
where ^ x is the estimation state. Since we are observing the displacement and z is a scalar quantity, L takes the form:
L¼
l1 l2
:
ð51Þ
We can define the error system as e ¼ ^ x x, producing:
e_ ¼ ðA þ LC Þe Q:
ð52Þ
According to our method for an observable system, we can design L such that the matrix A þ LC is well-conditioned so that A þ LC is invertible, and the eigenvalues of A þ LC sufficiently large such that:
Q LCeH ;
ð53Þ
^: and we can use Eq. (45) to define our force estimate Q
^ ¼ LCeH : Q
ð54Þ
Due to the size of the system, we are able to design the eigenvalues through direct computation, but it should be noted there are other methods that can be employed to design the eigenvalues which would more readily lend itself to systems of higher dimension, such as the Bass-Gura method or linear quadratic estimator method [20]. To design the eigenvalues of A þ LC, we begin with the eigenvalue problem for A þ LC:
det ½kI 22 ðA þ LCÞ ¼ 0;
ð55Þ
where k represents the eigenvalues of the matrix A þ LC and:
A þ LC ¼
l1
1
4 þ l2
4
:
ð56Þ
48
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
The eigenvalues of A þ LC are determined by finding the zeroes of Eq. (55), which when expanded produces the following:
2 1 1 2 k þ ð4 l1 Þ ð4 l1 Þ 4l1 þ 4 l2 ¼ 0; 2 4
ð57Þ
and are specified by assigning values to the parameters l1 and l2 . For simplicity, we choose the eigenvalues that produce the double-root k ¼ ð4 l1 Þ=2. Thus, we can define l2 in terms of l1 :
1 2 l2 ¼ ð4 l1 Þ 4l1 þ 4 4
ð58Þ
and we can assign values to the parameter l1 . Eqs. (47) and (50) were simulated with the input: FðtÞ ¼ 10 sinð2tÞ N, FðtÞ is a sawtooth wave with period 2p and amplitude of 10 N, and with the various values of l1 : l1 ¼ 10 s1, l1 ¼ 100 s1, l1 ¼ 1000 s1, to exhibit the convergence characteristics of the proposed estimation method. The frequency for the harmonic applied force was chosen to coincide with the natural frequency of the system to show that the proposed method is applicable at the natural frequency of the system. The sawtooth wave was chosen to highlight the applicability of the method to an applied force with rapid changes in amplitude. To obtain the displacement data, Eqs. (47) and (50) were simulated in Matlab with the prescribed input FðtÞ and system parameters defined in A and L. The total period of simulation was 10 s, and was divided into segments of .0001 s each. Using a fifth-order Dormand-Prince method given by the Matlab function ode5, the terms, x and ^ x, were solved at 51 points between time segments. Beginning at time t ¼ 0 with initial conditions set as zero and using the fifth-order DormandPrince method, we solved for final condition of the first segment, which ends at time t ¼ :0001. This is then the initial condition of the next segment, and this process is continued for the total period of 10 s. This provided us with the system state x, which consists of the displacement and velocity of the mass. To obtain the displacement data, we multiplied x by C to determine z. This was substituted into Eq. (54), and can be written as:
^ ¼ LC ð^x xÞ: Q
ð59Þ
x and x by C and using Eq. (51), Eq. (59) can be expressed as: Multiplying ^
^ ¼ l1 ð^x1 x1 Þ: Q l2
ð60Þ
Using:
" ^ ¼ Q
# b1 Q ; b2 Q
ð61Þ
^ approximates Q , our approximation for FðtÞ is given by: and since Q
b 2 ¼ l2 ð^x1 x1 Þ: Q
ð62Þ
The systems were simulated for each value of l1 , leading to a value of l2 via Eq. (58), and the results for each of these cases are shown in Figs. 3 and 4. It is readily observed that as the magnitude of l1 is increased, the estimate for the input converges to the actual input. It should be noted that in Fig. 4, the initial value of the force estimate is does not agree as well. While this is dependent on the initial conditions, if the magnitude of l1 is further increased, the initial difference will remain but will converge to the applied force rapidly. Also, we should note that the application of the proposed method is inherently asymptotic, meaning that as the magnitude of l1 goes to infinity Eq. (62) will converge to a particular set of values. And while the values of l1 were selected to highlight the asymptotic nature of the proposed method, l1 can be selected through the assignment of a convergence criteria. As b 2 will tend to zero, and once it reaches an assigned threshold value, this value the magnitude of l1 increases, the variation in Q of l1 will be the determined value for acceptable convergence. From these results, we are able to conclude that the proposed method is able to provide an estimate of the input based only on the observed parameters. This is an important result because it shows that we are able to produce a relationship
49
Applied Force (N)
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
l1= −10 10 Actual Estimate
5 0 −5 −10 0
1
2
3
4
5
6
7
8
9
10
Applied Force (N)
Time (s)
l1= −100 10 Actual Estimate
5 0 −5 −10 0
1
2
3
4
5
6
7
8
9
10
Applied Force (N)
Time (s)
l1= −1000 10 Actual Estimate
5 0 −5 −10 0
1
2
3
4
5
6
7
8
9
10
Time (s) Fig. 3. Single degree-of-freedom system sinusoidally forced.
Applied Force (N)
l1= −10 10 Actual Estimate
5 0 −5 −10
0
1
2
3
4
5
6
7
8
9
10
Time (s)
Applied Force (N)
l1= −100 10 Actual Estimate
5 0 −5 −10
0
1
2
3
4
5
6
7
8
9
10
Time (s)
Applied Force (N)
l1= −1000 10 Actual Estimate
5 0 −5 −10
0
1
2
3
4
5
6
7
Time (s) Fig. 4. Single degree-of-freedom system with a sawtooth wave force.
8
9
10
50
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
between the observed measurements and the input directly. Additionally, it shows that the estimate for the input can be independent of the system parameters, given by A. 5. Application to continuous systems using modal analysis So far, we have developed the theory and applied the proposed method to a single degree-of-freedom system. In this section, our aim is to apply the proposed method to more complex systems by estimating loads on continuous systems. While the scope of this section is strictly limited to modal analysis as the method of discretizing a continuous system, we are confident that the proposed estimation method can work in tandem with other commonly accepted discretization methods, such as Rayleigh-Ritz or finite element. This ability to work with discretization methods broadens the applicability of the proposed estimation method. It is not limited in scope to simple models, but is capable of being applied to more complex structures or applications. To better express the development of the proposed method with modal analysis, we apply the proposed method to the example of estimating the force applied to an Euler-Bernoulli beam, which is governed by the following equation:
qA
@2w @w @4w þc þ EI 4 ¼ pðx; tÞ; 2 @t @x @t
ð63Þ
where wðx; tÞ is the transverse displacement of the beam at location x along the beam and at time t; pðx; tÞ is the transverse load along the length of the beam, E is the Young’s modulus, I is the moment of inertia of the beam, q is the material density, A is the cross sectional area, and c is the viscous damping coefficient [3]. In the example presented in this section, we are assuming the full displacement field of the beam is known. While this assumption is not realistic, we use this assumption strictly in the context of demonstrating the application of the proposed method to continuous systems. As the proposed method is developed with the aim of estimating the temporal component of a dynamic force, our objective in this section is to show that the proposed method is capable of being applicable to continuous systems using discretization methods to estimate time-varying loads. In reality, we would only be able to take measurements at discrete points along the surface of the beam. In this context, the proposed method would have to be merged with a modal observer, such as that used in [6], but by assuming access to the continuous displacement field, we are able to develop the theory and determine the behavior of the proposed method when used with more complex structures. To develop the theory, the proposed method will be applied to estimate the loads as shown in Figs. 5 and 6, and given by the following:
x pðx; tÞ ¼ 100 CðtÞ l pðx; tÞ ¼ 100dðx aÞCðtÞ; where dðx aÞ is the Dirac delta function, l is the length of the beam, a is a point on the beam, and:
Fig. 5. Simply supported Euler-Bernoulli beam of length l with a linearly distributed load.
ð64aÞ ð64bÞ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
51
Fig. 6. Simply supported Euler-Bernoulli beam of length l with a point load.
CðtÞ ¼ sinð2tÞ:
ð65Þ
Loads characterized by Eqs. (64a) and (64b) were chosen specifically to highlight certain attributes of modal analysis. In Section 5.1, we will describe the theoretical analysis of estimating the load using modal analysis. In Section 5.2, we will derive an exact solution to the estimation of the temporal components of the distributed load, and in Section 5.3, we present results of the simulation of this method and discuss the important attributes of these simulations. 5.1. Modal asymptotic approximation: theory To apply the proposed method to more complex structures using modal analysis, we begin by expanding the displacement, wðx; tÞ, and the distributed load, pðx; tÞ in terms of orthogonal functions f n ðxÞ, where:
Z
l
f n ðxÞf m ðxÞdx ¼ dnm ;
ð66Þ
0
and f n ðxÞ satisfy the given boundary conditions. This property is critical to the application of the proposed method as it enables us to transform the governing equation, given by Eq. (63), into the necessary form. Using the orthogonality property of Eq. (66), the parameters wðx; tÞ and pðx; tÞ can be expanded in infinite series as: 1 X
wðx; tÞ ¼
yn ðtÞf n ðxÞ;
ð67aÞ
1 X Q n ðtÞf n ðxÞ;
ð67bÞ
m¼1
pðx; tÞ ¼
m¼1
where:
Z
l
yn ðtÞ ¼
wðx; tÞf n ðxÞdx
ð68aÞ
0
Z Q n ðtÞ ¼
0
l
pðx; tÞf n ðxÞdx:
ð68bÞ
By multiplying Eq. (63) by f n ðxÞ, integrating over the domain, and using Eqs. (67a) and (67b), we find:
qAy€n þ cy_ n þ EI
np4 l
yn ¼ Q n ðtÞ:
In this form, we are able to apply our reconstruction method to each mode yn . We need to put Eq. (69) into the same form as Eq. (37) by defining X 1n ¼ yn and X 2n ¼ y_ n , yielding:
ð69Þ
52
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
" An ¼ Xn ¼
0 EI np 4
qA
l
X 1n ðtÞ
#
1 qcA
ð70aÞ ð70bÞ
X 2n ðtÞ " # 0 Qn ¼ 1 : qA Q n ðtÞ
ð70cÞ
Since the displacement is the parameter being observed, it needs to be put in terms of the modal parameters. We define the measurement as z ¼ wðx; tÞ. To obtain the modal measurement, we need to multiply z by f n ðxÞ and integrate over the domain to produce zn ¼ yn . Rearranging zn to be in the same form as Eq. (37), we obtain:
zn ¼ CX n ;
ð71Þ
C ¼ ½ 1 0 ;
ð72Þ
where
and zn is the modal measurement. Now we can apply the reconstruction method for each mode, which begins by defining an unperturbed estimation system with observer feedback for each mode by:
^ ^ n þ Ln C X ^ n Xn ; X_ n ¼ An X
ð73Þ
where Ln is a design parameter. The nth mode error system is:
E_ n ¼ ðAn þ Ln C ÞEn Qn ;
ð74Þ
b n X n . For well-conditioned An þ Ln C with sufficiently large eigenvalues, we can estimate Qn by: where En ¼ X
^ n ¼ Ln CEn ; Q
ð75Þ
^ n is the estimate of Qn . To obtain an approximation for the input, we must multiply Eq. (75) by f ðxÞ and sum over where Q n the domain n from 1 to infinity. As a result, our estimate becomes:
^ tÞ ¼ Qðx;
1 X ^ n ðtÞf ðxÞ; Q n
ð76Þ
n¼1
^ is the estimate for: where Q
"
Qðx; tÞ ¼
0 1
qA pðx; tÞ
#
:
ð77Þ
We have approximated the input vector Q, but we still have not quite approximated the distributed load pðx; tÞ. To do so, we must acknowledge that Q consists of the term pðx; tÞ=ðqAÞ. Thus, to obtain the distributed load, we must define a vector P such that:
P ¼ qAQ;
ð78Þ
where P contains the distributed load, and our estimate for the distributed load becomes:
^ ^ ¼ qAQ; P
ð79Þ
^ is the estimate for the distributed load, pðx; tÞ. where P 5.2. Modal asymptotic approximation: exact solution Now that we have outlined the theory, we can demonstrate its usage by applying it to the linearly-varying load from Fig. 5. Our aim in this section is to develop an exact solution so that we may be able to better develop the analysis of the proposed method in this type of application. For this example, the Euler-Bernoulli beam is simply supported, which gives us the following boundary conditions:
@2w ð0; tÞ ¼ 0 @x2 2 @ w wðl; tÞ ¼ 2 ðl; tÞ ¼ 0: @x
wð0; tÞ ¼
ð80aÞ ð80bÞ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
53
The initial conditions for the undeformed beam are assumed to be:
wðx; 0Þ ¼
@w ðx; 0Þ ¼ 0: @x
ð81Þ
We are developing the solution for the applied load given by Eq. (64a), and are assuming the transverse displacement field wðx; tÞ is measured. First, we must determine the mode shape, f n ðxÞ. For the given boundary conditions, the mode shape is given by:
rffiffiffi np 2 sin f n ðxÞ ¼ x : l l
ð82Þ
Then, we must select the matrix LC. Due to the size of the system, we are able to design the eigenvalues through direct computation, but we should again note there are other methods which can be employed to design the eigenvalues which would more readily apply to systems of higher dimension, such as the Bass-Gura method or linear quadratic estimator method [20]. Since we are assuming that we can only observe the displacement field of the beam, Ln C takes the form:
Ln C ¼
ln1
0
ln2
0
:
ð83Þ
We need to design the terms ln1 and ln2 such that the eigenvalues of An þ Ln C are sufficiently large with a negative real part to ensure asymptotic stability. The process to determine the eigenvalues as a function of l1 is identical to that of the example for a single degree-of-freedom system in Section 4. Using Eq. (55), the eigenvalues for each mode are computed from the characteristic polynomial:
det ½kI 22 ðAn þ Ln CÞ ¼ k2 þ
c c EI np4 þ ln2 : ln1 k ln1 qA qA qA l
ð84Þ
We determine the eigenvalues by finding the zeros of Eq. (84), and since ln1 and ln2 are design parameters, we can parameterize these to be functions of a single variable, b, such that ln1 ðbÞ and ln2 ðbÞ. To achieve this, we first need the quadratic equation with the desired zeros k ¼ a and k ¼ b. This is given by:
k2 þ ða þ bÞk þ ab ¼ 0:
ð85Þ
From here, we equate the coefficients of Eqs. (84) and (85) as follows:
c aþb¼ l qA n1 c EI np4 ab ¼ þ ln2 : ln1 qA qA l
ð86aÞ ð86bÞ
To further simplify our eigenvalue selection, we define a ¼ 2b to produce:
3b ¼
c
ln1 c EI np4 2 2b ¼ þ ln2 : ln1 qA qA l
qA
ð87aÞ ð87bÞ
This gives us the values for ln1 and ln2 in terms of b:
ln1 ¼
c
qA
3b
ð88aÞ
2
ln2 ¼ 2b
c c EI np4 : 3b þ qA qA qA l
ð88bÞ
For simplicity, we define:
EI np4 l c c ^ 2 ¼ 2b2 þ x 3b : qA qA
x2n ¼
qA
ð89aÞ ð89bÞ
Using Eqs. (89), ln2 can be rewritten as:
^ 2: ln2 ¼ x2n x
ð90Þ
This allows us to use only a single parameter, b, to define the eigenvalues of An þ Ln C. Using Eqs. (88)–(90), Ln C becomes:
54
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
" Ln C ¼
qA 3b c
0
#
^2 0 x2n x
:
ð91Þ
From here, we can use Eq. (75) as the estimate for each modal input, and while it can be applied exactly as is, we can gain more insights by returning to the modal governing equation, given by Eq. (69), in combination with the state-space analysis to derive the solution explicitly in closed form. For this analysis, we will provide the solution to pðx; tÞ corresponding to Eq. (64a). By multiplying Eq. (64a) by f n ðxÞ and integrating from 0 to l, we obtain the modal input as follows:
Q n ðtÞ ¼ T n sinðxtÞ; where:
10 Tn ¼ qA
ð92Þ
rffiffiffi 2 cosðnpÞ : l np
ð93Þ
For the initial conditions given by Eq. (81), the solution of Eq. (69) for yn ðtÞ is given by:
yn ðtÞ ¼
c
x2n x2 2
x2n x2 þ
c qA
2 T n sinðxtÞ
x
qA x
x2n x2 2 þ
c qA
x
2 T n cosðxtÞ
ð94Þ
Using Eq. (92), Eq. (94) can be rewritten in terms of Q n ðtÞ as:
yn ðtÞ ¼
c
x2n x2 2
x2n x2 þ
c qA
2 Q n ðtÞ
x
qA x
x2n x2 2 þ
c qA
x
2 T n cosðxtÞ:
ð95Þ
Now, we must rewrite the estimate system, given by Eq. (73), in terms of yn ðtÞ. To do so, we must first rewrite Eq. (73) as individual equations. Using Eq. (91), this produces:
c b n2 þ 3b c X n1 b n1 þ X 3b X qA qA 2 c 2 2 b n1 b þ x ^ X ^ xn X n1 : ¼ x X qA n2
b_ n1 ¼ X b_ n2 X
ð96aÞ ð96bÞ
b n2 and X b n1 are estimates for X n2 and X n1 , respectively, and since X n1 ¼ y , we can conclude that X b n1 ¼ y ^n and Since X n _ € b ^n . This means we can express Eq. (96b) as: X n2 ¼ y
€^ þ y n
c
qA
2 ^ 2y ^ x2n yn ðtÞ: ^n ¼ x ^n þ x y
ð97Þ
^n is given If we assign the initial conditions of zero for the estimation system and use Eq. (95), the solution of Eq. (97) for y by:
2 3 2 2 2 c 2 2 2 2 2 2 ^ ^ ^ x x x x x x x x x 6 n n n 7 qA ^n ðtÞ ¼ 6 7Q ðtÞ y 2 4 2 c 2 5 n 2 c 2 2 2 2 ^ x Þ þ qA x ðx xn x þ qA x 3 2 2 2 c c 2 2 2 2 2 2 ^ ^ ^ x x x x x x x x x x þ 7 6 n n n qA qA 6 2 2 7 5T n cosðtÞ: 4 2 2 ^ 2 x2 Þ þ qcA x ðx x2n x2 þ qcA x
ð98Þ
Using Eq. (75), the estimate for each modal input is expressed as:
b n ¼ ln2 ðy ^n yn Þ: Q
ð99Þ
Substitution of Eqs. (88b), (95), and (98) into Eq. (99) and simplifying produces:
2 2 2 c ^2 ^ x2n x ^ x2 x qA x xn x b Q n ðtÞ ¼ 2 Q n ðtÞ þ 2 T n cosðxtÞ: ^ 2 x2 Þ2 þ qcA x ^ 2 x2 Þ2 þ qcA x ðx ðx
The distributed load, pðx; tÞ, is estimated by multiplying Eq. (100) by f n ðxÞ and summing from 1 to infinity:
ð100Þ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
2 3 2 2 2 c 2 1 2 2 ^ X x x x ^ ^ x x x x n qA 7 n ^ðx; tÞ ¼ qA 6 p 4 2 Q n ðtÞ þ 2 T n cosðxtÞ5f n ðxÞ: n¼1 ðx ^ 2 x2 Þ2 þ qcA x ^ 2 x2 Þ2 þ qcA x ðx
55
ð101Þ
From Eqs. (89) and (100), we can see several important attributes of the application of the proposed method with modal b n to Q as b goes to infinity. As b goes to infinity in Eq. analysis. First, this provides direct evidence of the convergence of Q n 2 2 b ^ ^ (89b), x goes to infinity, and in Eq. (100), as x goes to infinity, Q n converges to Q n . In the context of the continuous problem of an Euler-Bernoulli beam using modal analysis, we can see convergence occurs for the temporal component of each mode in the ideal case as b goes to infinity. Additionally, we can see the connection between the selection of a finite value for b and the modal natural frequency. Since we’re choosing a finite value for b, there can be a coupling between b and the number of modes in which we choose in the approximation of the load, pðx; tÞ. This is apparent in Eqs. (89) and (100). In Eq. (100), the coefficient of Q n is dependent ^ 2 is a function of b and x2n is a function of n. From Eq. (89a), we can see that x2n grows with ^ 2 x2n , where x on the term x each additional term used in the modal sum. For a large but finite value for b, each additional mode used in the approxima^ 2 x2n . From this, we can conclude that the number of tion of pðx; tÞ will result in a decrease in the magnitude of the term x modes used to approximate pðx; tÞ should be taken into account during the selection of the parameter b. In practice, the value bn of b can be selected through the assignment of a convergence criteria. As the magnitude of b increases, the variation in Q will diminish, and once it reaches an assigned threshold value, this value of b will be the determined value for acceptable convergence. This procedure should not need to be applied to every mode in the approximation of pðx; tÞ. It can be applied to the mode with the largest modal natural frequency, and once a suitable value for b is determined for this mode, it should be suitable for all other modes. 5.3. Modal asymptotic approximation: numerical simulations To exhibit the application of the proposed method in estimating the load applied to an Euler-Bernoulli beam, we applied the results from Section 5.2 to show the resulting estimation of the distributed load. This will be performed for the linearly varying load shown in Fig. 5, and for a point load as shown in Fig. 6. These two spatial distributions were chosen to highlight certain attributes of modal analysis. The linearly distributed load is an odd function, and the point load is not. As a result, the characteristics of the functions must be taken into account. To run the simulations, the following physical parameters are defined:
EI ¼ 107 Nm2, qA ¼ 0:1 kg/m, l ¼ 10 m, c ¼ 4 Ns/m2, x ¼ 2 radians/s.
Additionally, for the linearly varying load, we varied the value for b to show how the estimate varies with a fixed modal summation, and for the point load, we varied the number of modes summed to show how the estimate varies for fixed b. As we have already shown the convergence in the temporal domain, we are only showing the convergence of the estimate for pðx; tÞ in the spatial domain. b n ðtÞ. Eq. (82) was evalTo obtain data, we computed Eqs. (95) and (98), and substituted them into Eq. (99) to determine Q b n ðtÞ. This was performed for n 2 ½1; 20, and then uated at 101 points on the interval ½0; d, and we multiplied this by each Q summed to obtain the estimate for pðx; tÞ. This was performed for the values of b: b ¼ 10; 000 s1, b ¼ 100; 000 s1, b ¼ 1; 000; 000 s1. The results for each b are shown in Fig. 7, and an enlarged view of b ¼ 1; 000; 000 s1 is shown in Fig. 9. As we can see from Fig. 7, as we increase the value of b, the estimate for the distributed load converges actual distributed load. This actually agrees very well with the theory we derived in Section 5.2. If we define the coefficient for Q n in Eq. (100) as follows: (Fig. 10)
bn ¼
^ 2 x2n x ^ 2 x2 x 2
^ 2 x2 Þ þ ðx
c
qA x
2 ;
ð102Þ
we can use this relation as the metric for the quality of the estimation. To rewrite bn in terms of b, we substitute Eqs. (89) into Eq. (102):
56
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
b= 10,000
4
Estimated Distribution (N/m)
1.5
x 10
1 0.5 0 −0.5 −1 −1.5
0
1
2
3
4
5
6
7
8
9
10
6
7
8
9
10
6
7
8
9
10
Distance (m)
Estimated Distribution (N/m)
b= 100,000 200
100
0
−100
−200
0
1
2
3
4
5
Distance (m)
Estimated Distribution (N/m)
b= 1,000,000 150 100 50 0 −50 −100 −150
0
1
2
3
4
5
Distance (m) Fig. 7. Estimated linear load distribution under various values of b.
bn ¼
2
2b þ qcA
2 2 2 4 2 2b þ qcA qcA 3b 3b qEIA nlp x2 : 2 2 2 2 x2 þ qcA x 2b þ qcA qcA 3b
c qA
ð103Þ
From Eq. (103), we can see that as b goes to infinity, bn goes to 1. So, the closer to 1 bn is, the better the approximation we obtain. For n ¼ 20, Eq. (103) can be evaluated for each value of b. For b ¼ 10; 000 and b ¼ 100; 000; bn is on the order of 773 and 6:8 respectively, which indicates these values are severely underperforming. But for b ¼ 1; 000; 000; bn ¼ 0:922, indicating a high level of agreement for the 20th mode. This is readily seen in Fig. 7, as the first two graphs show poor results and the third graph exhibits a high level of convergence. For the case of b ¼ 1; 000; 000, the load distribution peaks at 112 N/m while the actual peak value is 100 N/m. This means there is only a 12% difference between the estimate peak and actual peak values, and higher accuracy could be achieved by increasing the value of b and the number of modes used in the estimation. Now that we have seen the behavior of the proposed method as b varies, we want to see the behavior of the proposed method as the number of modes varies. For this case, we will apply the proposed method to the estimation of a point load. As we have derived everything we need regarding the linearly distributed load in Section 5.2, we show here only the nec-
57
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Estimated Distribution (N/m)
10 Modes Estimation 1500 1000 500 0 −500 −1000 −1500 0
1
2
3
4
5
6
7
8
9
10
7
8
9
10
7
8
9
10
Distance (m)
Estimated Distribution (N/m)
15 Modes Estimation 1500 1000 500 0 −500 −1000 −1500 0
1
2
3
4
5
6
Distance (m)
1500 1000 500 0 −500 −1000 −1500 0
1
2
3
4
5
6
Distance (m) Fig. 8. Estimated point load under various modal sum conditions.
Estimated Distributed Load with b=1,000,000 150
Estimated Distribution (N/m)
Estimated Distribution (N/m)
20 Modes Estimation
100 50 0
−50 −100 −150 0
2
4
6
8
Distance (m) Fig. 9. Estimated load distribution with b ¼ 1; 000; 000.
10
58
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Estimated Point Load with b=1,000,000
Estimated Distribution (N/m)
2000
1000
0
−1000
−2000
0
2
4
6
8
10
Distance (m) Fig. 10. Estimated point load with b ¼ 1; 000; 000 and 20 modes summed.
essary results. We repeated the procedure performed in Section 5.2 for a point load at a distance a ¼ 5 m from the end of the b n ðtÞ are presented below: ^n ðtÞ, and Q beam, which is given by Eq. (64b). The results for Q n ðtÞ; yn ðtÞ; y
Q n ðtÞ ¼ T n sinðxtÞ rffiffiffi anp 10 2 sin Tn ¼ qA l l yn ðtÞ ¼
ð104aÞ ð104bÞ c
x2n x2
x2n x2 2 þ
c
2 T n sinðxtÞ
qA x
qA x
x2n x2 2 þ
c
qA x
2 T n cosðxtÞ
2 3 2 2 2 c 2 2 2 2 2 2 ^ ^ ^ x x x x x x x x x 6 n n n 7 qA ^n ðtÞ ¼ 6 7Q ðtÞ y 2 4 2 c 2 5 n 2 c 2 2 2 ^ 2 x Þ þ qA x ðx xn x þ qA x 3 2 2 2 c c 2 2 2 2 2 2 ^ ^ ^ 6 x x x xn qA x þ qA x x xn xn x 7 6 2 2 7 5T n cosðtÞ: 4 ^ 2 x2 Þ2 þ qcA x ðx x2n x2 2 þ qcA x 2 2 2 c 2 2 ^2 ^ ^ qA x xn x b n ðtÞ ¼ x xn x x Q ðtÞ þ Q 2 T n cosðxtÞ: n 2 ^ 2 x2 Þ2 þ qcA x ^ 2 x2 Þ2 þ qcA x ðx ðx
ð104cÞ
ð104dÞ
ð104eÞ
We repeated the procedure of the linearly distributed load to obtain data, but we computed Eq. (104c) to determine the modal displacements and (104d) for the modal estimate. In this case, we kept b constant at the value of 1; 000; 000 s1, but varied the number of modes used in the estimation of pðx; tÞ. The number of modes used are as follows: n ¼ 10, n ¼ 15, n ¼ 20. As we can see in Fig. 8, the number of modes greatly affects the estimate for the point load. While 20 modes were adequate for the linearly varying load, the point load has is not as fortunate. While we can see some sense of convergence, some concerns do remain. For starters, the width of the peak grows more narrow with an increase in the number of modes used, but it might require a much larger value for n before an acceptable estimate for the applied load is computed. In addition to the narrowing of the peak, a significant number of minor peaks surround the peak at the location of the point load. The difference in the convergence between the linearly varying load and the point load is a consequence of the agreement of the types of loads with the types of problem. A simply-supported beam is essentially an odd-functioned problem, and the best result will come from odd-functioned loads, such as the linearly varying load. The point load is not as fortunate in that it is not a purely odd function, and convergence to the true value can be limited in this respect. This is one of the limitations of applying the proposed method with modal analysis to estimate distributed loads. For the linearly-varying load, modal analysis is a perfectly acceptable method of discretizing the continuous problem, but a more appropriate method of discretization would be required to better estimate the point load.
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
59
6. Derivation of proposed method – real case Thus far, we have assumed perfect knowledge of the response in our measurement, but it is important to understand what would occur in the likely event that uncertainty in the measurements is present. To make such a determination, we must alter the original system given by Eq. (37) to include measurement noise as:
x_ ¼ Ax þ Q ðtÞ
ð105aÞ
y ¼ Cx þ g;
ð105bÞ
where g is an m 1 measurement noise vector. Using the estimate system given by Eq. (39) in conjunction with Eq. (105), the error system becomes:
e_ ¼ ðA þ LC Þe QðtÞ Lg:
ð106Þ
We apply Theorem 1 and rearrange to produce:
ðA þ LC ÞeH ¼ Q ðtÞ þ Lg þ ðA þ LC Þ1 Q_ ðtÞ þ ðA þ LC Þ1 Lg_ þ
1 h X m¼2
ðA þ LC Þ1
im dm h i m Q ðtÞ þ Lg : dt
ð107Þ
By applying Propositions 1 and 2 to Eq. (107), we produce the following approximation:
ðA þ LC ÞeH Q ðtÞ þ Lg þ ðA þ LC Þ1 Lg_ :
ð108Þ
The term ðA þ LC ÞeH can be further simplified. So far, the proposed method has been discussed only in terms of the error system, but it is just as applicable to the estimate system given by Eq. (39). When we apply Theorem 1 and Propositions 1 and 2 to Eq. (39), we obtain:
^xH ðA þ LC Þ1 Ly:
ð109Þ
Using Eqs. (105b) and (109), we can further simplify Eq. (108) to produce the approximation:
^ y Q ðtÞ AðA þ LC Þ1 Lg þ ðA þ LC Þ1 Lg_ : L y
ð110Þ
In the form given by Eq. (110), we can define our estimate by:
b g ðtÞ ¼ L y ^y ; Q
ð111Þ
^ ðtÞ is the estimate of the input with measurement noise. where Q g 6.1. Discussion The concepts developed in Section 2 provided useful insight into the inverse problem. While measurement noise remains an issue, the concepts of the proposed method allow us redefine the force reconstruction problem in a different context. The concepts developed allow the force reconstruction problem to be defined as a time series analysis problem. In this treatment of measurement noise, alternative signal processing techniques can be applied to Eq. (111) to reduce the effects of measurement noise and recover an estimate of the original input. In addition, the recovery of the forces is not limited to strictly forces in the time domain, but the concepts herein can be extended further. For example, if random forces are present, the objective would be to recover the spectral density of the random forces applied to the system. The current formulation readily lends itself to such problems, and additional work in the future will be performed to address random forces from a probabilistic perspective. 6.2. Example: Exact solution of a scalar system To show that the approximation given by Eq. (110) holds analytically, we apply the proposed method to the scalar system:
x_ ¼ Ax þ expðMtÞ y ¼ Cx þ g;
ð112aÞ ð112bÞ
where M is a non-zero constant, A < 0, and to simplify the analysis, we let C ¼ 1. Since our analysis does not make any assumptions on the form of g, we define g for the purposes of demonstration as:
g ¼ B expðixtÞ;
ð113Þ
60
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
pffiffiffiffiffiffiffi where B and x are constants and i ¼ 1. This definition of g will allow us to produce an exact analysis of the system using the proposed method, and as a result, it will allow us to better demonstrate the characteristics of the proposed method and the concepts developed. At equilibrium, Eq. (112b) becomes:
y¼
expðMtÞ þ B expðixtÞ: MA
ð114Þ
To apply the proposed method, the following estimate system must be defined:
^x_ ¼ A^x þ Lð^x yÞ:
ð115Þ
Using Eq. (114), we apply Theorem 1 to Eq. (115) to produce:
^xH ¼
1 X ðA þ LÞðmþ1Þ m¼0
L m Mm expðMtÞ þ BLðixÞ expðixtÞ : MA
ð116Þ
Using Eqs. (114) and (116), we can apply Eq. (111) to produce our estimate. This yields:
bg ¼ Q
L Lðix AÞ expðMtÞ þ B expðixtÞ: AþLM A þ L ix
ð117Þ
For L sufficiently large, it can be said that:
b g expðMtÞ þ ðix AÞB expðixtÞ: Q
ð118Þ
Now if we use the right-hand side of Eq. (110) with g defined as in Eq. (113), we proceed as follows:
AL Lix B expðixtÞ þ B expðixtÞ AþL AþL Lðix AÞ B expðixtÞ ¼ expðMtÞ þ AþL expðMtÞ þ ðix AÞB expðixtÞ:
QðtÞ AðA þ LÞ1 Lg þ ðA þ LÞ1 Lg_ ¼ expðMtÞ
ð119aÞ ð119bÞ ð119cÞ
Comparing Eqs. (118) and (119c), we can see the approximation given by Eq. (110) holds and that our approximation defined by Eq. (111) behaves as predicted. 6.3. Example: Single degree-of-freedom (continued) To see the application of the proposed method to force estimation with realistic measurement noise, we return to the single degree-of-freedom example from Section 4. While we still assume that we are measuring the displacement of the mass, we amend Eq. (49e) to include noise as follows:
z ¼ Cx þ g;
ð120Þ
where here g is uncorrelated zero-mean Gaussian white noise with variance 1, which was used as we wanted to have a large ratio of variance to maximum amplitude. The rest of the analysis from Section 4 remains unchanged, and this analysis is simulated with the applied force:
FðtÞ ¼ 10 sinð2tÞ N;
ð121Þ
and at values of l1 given by: l1 ¼ 10000 s1, l1 ¼ 100000 s1, l1 ¼ 1000000 s1. With the prescribed FðtÞ, the max amplitude of the displacement of the mass is 1.49 m. This gives a ratio of noise variance to max amplitude of 0.6727. To obtain the displacement data, Eqs. (47) and (50) were simulated in Matlab with the prescribed force FðtÞ and system parameters defined in A and L. The total period of simulation was 10 s, and was divided into segments of .0001 s each. Using a fifth-order Dormand-Prince method solved by the Matlab function ode5, the terms, x and ^ x, were solved at 51 points simultaneously between time segments. Beginning at time t ¼ 0 with initial conditions set as zero and using the fifth-order Dormand-Prince method, we solved for final condition of the first segment, which ends at time t ¼ :0001. This is then the initial condition of the next segment beginning at t ¼ :0001 and ending at t ¼ :0002. This process is continued for the total period of 10 s. Additionally, using a random number generator with a zero-mean Gaussian distribution with variance of 1, a static random sample was added to the observer feedback during each segment to represent the measurement noise. At this point, we have produced an estimate for the input with a significant amount of noise. To smooth the signal and reproduce a cleaner estimate for the input, we introduced a low-pass filter to this result. Due to its simplicity, we chose to
61
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
apply a moving average to the signal to smooth out some of the high-frequency noise in the signal [2]. In continuous time, the moving average can be characterized as:
gðtÞ ¼
1=T : 0:
jtj < T=2 jt j P T=2
ð122Þ
The frequency attributes of the moving average are exhibited taking the Fourier transform of gðtÞ. To determine the frequencies to be emphasized by the low-pass filter, we must find the Fourier transform of the moving average:
GðxÞ ¼
sin
xT 2
xT
:
ð123Þ
2
where GðxÞ is the Fourier transform of gðtÞ. The frequency range most emphasized by the moving average is defined by the area between x ¼ 0 and the first x such that GðxÞ ¼ 0. This is given by:
xT 2
¼ p:
ð124Þ
Thus, the cutoff frequency, xcutoff , can be determined by rearranging Eq. (124) as:
xcutoff ¼
2p : T
ð125Þ
From Eq. (125), we can see that the frequency range to be emphasized is dependent on T. In this example, we choose T ¼ 0:01 s since it would give xcutoff ¼ 628:3 radians per second, which is a wide enough frequency range while reducing the effects of high frequency noise. The results of these simulations will be represented with estimates given by Eq. (111), as presented in Fig. 13, and with a moving average applied to Eq. (111), which is presented in Fig. 14. The moving average is used as a simple example of the potential for post-processing the estimate using noise reduction techniques applied in conjunction with the proposed method. Due to the presence of measurement noise, the values of l1 are larger in magnitude than those used in Section 4. Thus, the top two graphs in Figs. 13 and 14 are unable to be properly scaled, and are on orders of magnitude between 105 and 108 . The influence due to measurement noise is too great, and once we reach the order of magnitude for l1 ¼ 1; 000; 000, we are able to obtain an adequate approximation for the applied force relationship given by Eq. (110). An enlarged view of the case for l1 ¼ 1; 000; 000 is shown in Figs. 11 and 12 for better visibility.
l1= −10000
Applied Force (N)
7
4
x 10
Actual Estimate
2 0 −2 −4 0
1
2
3
4
5
6
7
8
9
10
Time (s)
l1= −100000
Applied Force (N)
8
4
x 10
Actual Estimate
2 0 −2 −4 0
1
2
3
4
5
6
7
8
9
10
Time (s)
Applied Force (N)
l1= −1000000 20 Actual Estimate
10 0 −10 −20 0
1
2
3
4
5
6
7
Time (s) Fig. 11. Single degree-of-freedom system with measurement noise – unfiltered.
8
9
10
62
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
l1= −10000
Applied Force (N)
5
5
x 10
Actual Averaged
0
−5 0
1
2
3
4
5
6
7
8
9
10
Time (s)
l1= −100000
Applied Force (N)
6
4
x 10
Actual Averaged
2 0 −2 −4 0
1
2
3
4
5
6
7
8
9
10
l1= −1000000 20 Actual Averaged
10 0 −10 −20 1
2
3
4
5
6
7
8
Time (s) Fig. 12. Single degree-of-freedom system with measurement noise – moving average.
l1= −1000000 20
Applied Force (N)
0
Actual Estimate
10
0
−10
−20
0
2
4
6
8
10
Time (s) Fig. 13. Single degree-of-freedom system with measurement noise for l1 ¼ 1000000 – unfiltered.
l1= −1000000 15
Actual Averaged
10
Applied Force (N)
Applied Force (N)
Time (s)
5 0 −5 −10 −15 0
2
4
6
8
10
Time (s) Fig. 14. Single degree-of-freedom system with measurement noise for l1 ¼ 1000000 – moving average.
9
10
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
63
7. Discussion and conclusion In this paper, we present the development of concepts in application to the problem of force reconstruction. It involves the usage of observer feedback, and in the present analysis, the parameter, given by L, is asymptotic in nature, and allows for the design of the eigenvalues of A þ LC such that approximations in terms of the applied force can be determined. The concepts were developed for an ideal model, and extended to a real model. The concepts developed present a number of benefits. The equilibrium analysis allowed us to develop an alternative representation of the relationship between the applied force and the system response. The solution is generally given by the convolution integral, and this analysis has allowed us to produce a differential series representation of the equilibrium, which proved critical to the development of the asymptotic approximation analysis. Additionally, the asymptotic approximation analysis allows for a simplified analysis of the inverse problem. The concepts developed simplify the analysis by allowing the inverse problem to be analyzed algebraically. This allows for the ease of computation and analysis, and effectively, circumventing the act of inversion in the force reconstruction problem. In application to the force reconstruction problem, the concepts developed herein alleviate the measurement noise issue by redefining the force reconstruction problem in the context of a time series analysis. This allows for the application of postprocessing to reduce the influence of noise and estimate the applied force. Additionally, probabilistic methods can potentially be applied to develop additional methods of modeling the measurement noise in relation to the concepts developed. The current analysis described herein is developed under the assumption of continuous sampling. Moving forward, we aim to relax this assumption, and reconcile the concepts developed with discrete sampling theory. Additionally, we aim to further develop these concepts in relation to probabilistic applied forces, beginning with stationary forces and the generalization of these concepts for non-stationary forces. In these circumstances, the recovery of a time series is not as useful as the recovery of the probabilistic aspects of this series, such as the autocorrelation function or the spectral density. The concepts developed herein can be extended to the estimation of the spectral density of a random force based on the spectral density of the measurements. Acknowledgements We would like to thank the reviewers for their detailed review of this manuscript. Their feedback was incredibly helpful as it allowed us to greatly improve the clarity of the presentation and the depth of the discussion of our work. References [1] Peter Avitabile, Fabio Piergentili, Ken Lown, Identifying dynamic loadings from measured responses, Sound Vib. 33 (1999) 24–28. [2] Christopher Terman, George Verghese, Hari Balakrishnan, Bits, signals, and packets: an introduction to digital communications and networks, From MIT open courseware.
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Asymptotic approximation method of force reconstruction: Proof of concept J. Sanchez a,⇑, H. Benaroya b a b
Sanchez Advanced Engineering Concepts, LLC, South Plainfield, NJ 07080, United States Department of Mechanical and Aerospace Engineering, Rutgers University, 98 Brett Road, Piscataway, NJ 08854, United States
a r t i c l e
i n f o
Article history: Received 24 March 2016 Received in revised form 16 January 2017 Accepted 18 January 2017 Available online 3 February 2017 Keywords: Dynamic force estimation Force reconstruction Asymptotic Approximation
a b s t r a c t An important problem in engineering is the determination of the system input based on the system response. This type of problem is difficult to solve as it is often ill-defined, and produces inaccurate or non-unique results. Current reconstruction techniques typically involve the employment of optimization methods or additional constraints to regularize the problem, but these methods are not without their flaws as they may be suboptimally applied and produce inadequate results. An alternative approach is developed that draws upon concepts from control systems theory, the equilibrium analysis of linear dynamical systems with time-dependent inputs, and asymptotic approximation analysis. This paper presents the theoretical development of the proposed method. A simple application of the method is presented to demonstrate the procedure. A more complex application to a continuous system is performed to demonstrate the applicability of the method. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction In many engineering applications, knowledge of the input or applied force is critical to the analysis of a system or structure. Unfortunately, at times, this knowledge may be difficult to obtain as it may not be possible to measure the applied force directly, and we must rely on measurements of the system response to determine the input or force applied to a system. This presents the problem of developing a model that estimates the input based on the measured output. Various approaches to this problem are outlined in [19], where such reconstruction methods are divided into three categories: direct methods, regularization methods, and probabilistic/statistical methods. Direct methods involve the application of the governing equation explicitly to determine an estimate for the applied force. These methods were employed by Jacquelin and Hamelin [11] to a split Hopkinson pressure bar using strain measurements, Law et al. [10] to an Euler-Bernoulli beam with a moving force applied, and were more generally developed by Karlsson [14], and Avitabile et al. [1]. While these methods are simplest to apply, they generally are ill-conditioned and can produce inadequate results. To overcome this ill-conditioning, additional constraints are often applied to transform the system into a wellconditioned problem, which is known as regularization. Various methods have been developed to regularize a problem, for example, through modal truncation, and optimization. Methods of truncation involve the isolation and removal of characteristics deemed to be attributed to noise, such as mode selection techniques developed by Jiang and Hu [12,13], and the
⇑ Corresponding author. E-mail addresses: [email protected] (J. Sanchez), [email protected] (H. Benaroya). http://dx.doi.org/10.1016/j.ymssp.2017.01.022 0888-3270/Ó 2017 Elsevier Ltd. All rights reserved.
40
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
application of wavenumber filtering by Djamaa et al. [4]. Additionally, the application of optimization methods were presented by Huang et al. [9] that utilized and optimized additional terms to filter noise inherent in the system, and the conjugate gradient method was detailed by Huang [8] and Yen and Wu [22,23]. Another type of optimization solution involves the application of the Levenberg-Marquardt iterative regularization method to solve inverse problems, and were presented by Gunawan [7] and Gasparo et al. [5]. While the application of these methods are generally successful, they do have their drawbacks. Some regularization methods possess additional parameters for regularization that lack a clearly defined method of determination while others require an intense level of computation to produce an estimate of the applied force. In addition to regularization methods, probabilistic/statistical methods have been applied to the solution of the inverse problem. These methods generally incorporate uncertainties such as measurement or system noise into the analysis and allows for a more holistic perspective of the inverse problem. One method that is becoming prevalent is the application of the Kalman filter to the inverse problem, and has been applied in conjunction with a recursive least squares algorithm by Ma et al. [16], by Xu et al. [21], and in an augmented form by Lourens et al. [15]. While the application of these methods have proven successful, like regularization methods, probabilistic methods can be computationally intensive as they may require a significant amount of data and time for computation. Fundamentally, these methods of estimation involve the truncation of system characteristics attributed to noise or the development of an optimization problem to produce a well-defined formulation. The dismissal of these assumptions would allow more freedom in the estimation, and allow us to explore new avenues of analysis that have otherwise gone unexplored. Our objective in developing our dynamic force estimation method was to develop a procedure that is inherently stable independent of any optimization or assumptions on the characteristics of the temporal component of the applied force throughout the recovery process. The work presented in this paper should be seen as a proof-of-concept where the estimation of dynamic forces is viewed from an asymptotic perspective. From this view, we aim to develop a method for estimating forces that has the potential to asymptotically converge to the applied force. To achieve this goal, we will compare the excited system to an unexcited estimate system under the assumptions that the system parameters are known and that sampling is continuous in time. Under these assumptions, and concepts from the following subject areas:
linear dynamical systems theory, control theory, equilibrium analysis of linear dynamical systems, asymptotic approximation analysis,
we will derive a relationship between the excited system, the unexcited system, and the system input. As we are interested in the development and proof of the concepts, we begin the development of this work under ideal circumstances where measurement noise is neglected, and extend these concepts to a real model, which includes measurement noise. Additional details are provided in [18], and while the focus of our applications in this paper are limited to the recovery of forces on mechanical systems, it is noted that the method is applicable to any linear system that satisfies the conditions discussed in this method. In this paper, the preliminary concepts developed for this analysis are summarized in Section 2. The derivation and significant results of the estimation method for an ideal model are discussed in Section 3. A simple example is developed in Section 4, and the application of the method to the estimation of distributed loads on continuous systems by using modal analysis is shown in Section 5. The concepts of the proposed method are applied to a real model in Section 6. Additional discussion and future work is provided in Section 7. 2. Preliminaries The development of the proposed estimation method utilizes concepts from multiple disciplines, but while these methods are critical to the development of the method, linear dynamical systems analysis is the starting point for the development. Throughout this work, we use bold face variables to denote matrices and underlined variables denote vectors. 2.1. Equilibrium analysis of linear dynamical systems The development of the proposed force reconstruction method begins with an analysis of the equilibrium of systems of the form:
x_ ¼ Ax þ Q ðtÞ;
ð1Þ
where Q ðtÞ is a time-dependent vector input, x is an n 1 state vector, and A is an invertible, well-conditioned matrix. Our aim is to transform Eq. (1) to a simplified equation of the form:
z_ ¼ Az;
ð2Þ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
41
where z ¼ x xH , and xH is the equilibrium state. By reducing Eq. (1) to Eq. (2), we aim to derive an expression for the equilibrium state, xH , in terms of the input, Q ðtÞ . We are additionally assuming that Q ðtÞ is infinitely differentiable and varying with time on the interval t 2 ½0; 1Þ. To achieve this, we begin by finding where x_ ¼ 0 for t 2 ½0; 1Þ. This equilibrium xH 0 of the linear dynamical system is given by: 1 xH 0 ¼ A QðtÞ:
ð3Þ
To better see what is happening, we can transform the state vector to be centered around the equilibrium xH 0 by defining _ _ _H the following variables: x1 ¼ x xH 0 and x1 ¼ x x0 . The variables would thus be transformed from x to x1 by substituting the relations for x1 and x_ 1 into Eq. (1), which produces:
_ x_ 1 ¼ Ax1 þ A1 QðtÞ:
ð4Þ
From this procedure, we see something interesting happening. Through the change of coordinates to center Eq. (1) around 1 _ the equilibrium xH 0 , we obtain another system of the same form but Q ðtÞ is replaced by A Q ðtÞ. Thus, for Eq. (3) to hold, Eq. (4) must hold. Additionally, as Eq. (4) is not of the form given by Eq. (2), the simplication of Eq. (1) to the form given by Eq. (2) is incomplete. Comparing Eq. (1) to Eq. (4), this suggests a recursive process. Repeating this process, we can find the equilibrium of Eq. (4), which is given by x_ 1 ¼ 0. This is true when:
2 1 _ xH QðtÞ: 1 ¼ A
ð5Þ
Defining x2 ¼ x1 xH 1 , we can rewrite Eq. (4) in terms of x2 . This produces the following:
€ x_ 2 ¼ Ax2 þ A1 QðtÞ:
ð6Þ
As we can see, Eq. (6) is of the same form of Eqs. (1) and (4). Beginning with Eq. (6), we repeat the procedure of changing the variables n times, with the following relations resulting by induction:
xn ¼ x þ
n1 X mþ1 ððAÞ1 Þ Q ðmÞ
ð7aÞ
m¼0 n
x_ n ¼ Axn þ ðA1 Þ Q ðnÞ ;
ð7bÞ
where ðmÞ denotes the mth order derivative with respect to time and n is a positive integer. Additionally, x_ n can be directly related to x by substitution of Eq. (7a) into Eq. (7b), to produce:
x_ n ¼ A x þ
! n X mþ1 ðA1 Þ Q ðmÞ :
ð8Þ
m¼0
Taking the limit as n goes to infinity, Eqs. (7a) and (8) become:
lim xn ¼ x þ
n!1
1 X mþ1 ðA1 Þ Q ðmÞ
ð9aÞ
m¼0
lim x_ n ¼ A x þ A1
n!1
1 X
! m
ðA1 Þ Q ðmÞ :
ð9bÞ
m¼0
Defining x1 ¼ limn!1 xn , Eqs. (9a) and (9b) become:
x1 ¼ x þ
1 X mþ1 ðA1 Þ Q ðmÞ
ð10aÞ
m¼0
x_ 1 ¼ Ax1 ;
ð10bÞ
and thus we have transformed Eq. (1) into the form of Eq. (2), where z ¼ x1 . The equilibrium of Eq. (10b) is found by setting x_ 1 ¼ 0. Since A is invertible and well-conditioned, the only state in which this is true is x1 ¼ 0. Thus, the equilibrium state of Eq. (10b) is xH 1 ¼ 0. To determine what happens at equilibrium, we set Eq. (10a) equal to xH , which produces the following: 1 1 X mþ1 xH ¼ ðA1 Þ Q ðmÞ :
ð11Þ
m¼0
Thus, the equilibrium state of the linear dynamical system given in Eq. (1) is given by Eq. (11), and the result of this analysis is presented as Theorem 1.
42
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Theorem 1. For a linear dynamical system of the form of Eq. (1) with the matrix A invertible and well-conditioned, and an input Q ðtÞ that is infinitely differentiable, the dynamic equilibrium xH of the system is given by: 1 mþ1 X xH ¼ A1 Q ðmÞ ;
ð12Þ
m¼0 m
where Q ðmÞ ¼ dtd m Q ðtÞ. The series expansion of xH as stated by Theorem 1 is important to the development of the proposed estimation method. It enables us to see the relationship between the system parameters of A and the system input Q ðtÞ, and how they relate to the equilibrium of the dynamical system given by Eq. (1). 2.1.1. Discussion of Theorem 1 The form for the equilibrium of Eq. (1) given by Theorem 1 allows us to transform the dynamical system into the simpler form of Eq. (2). To see this, we begin from the discussion after Eq. (2):
z ¼ x xH 1 X mþ1 ðA1 Þ Q ðmÞ : ¼xþ
ð13aÞ ð13bÞ
m¼0
Substitution of Eq. (13a) into Eq. (1) produces:
x_ ¼ Ax þ Q ðtÞ 1 1 X X m mþ1 z_ ðA1 Þ Q ðmÞ ¼ A z ðA1 Þ Q ðmÞ
z_ z_ z_
m¼1 1 X
ð14bÞ
m
m
ð14cÞ
m
m¼0 1 X
m
ð14dÞ
m
m¼1 1 X
m
ð14eÞ
ðA1 Þ Q ðmÞ ¼ Az
m¼1 1 X
þ QðtÞ
m¼0 1 X
ðA1 Þ Q ðmÞ ¼ Az
m¼1 1 X
ð14aÞ
!
ðA1 Þ Q ðmÞ ¼ Az
m¼1
z_ ¼ Az:
ðA1 Þ Q ðmÞ þ QðtÞ ðA1 Þ Q ðmÞ Q ðtÞ þ Q ðtÞ ðA1 Þ Q ðmÞ
m¼1
ð14fÞ
As we can see in Eqs. (14), the series expansion of the equilibrium state given in Theorem 1 allows for the cancelation of Q ðtÞ on the right hand side of Eq. (1), and the higher order derivatives of the series balance out on each side of the equality. As a result, the system is reduced from an excited to an unexcited linear system. Additionally, the solution for Eq. (14f) is:
z ¼ expðAtÞz0 ;
ð15Þ
where z0 is the initial state of z. Since we know the relationship between z and x, we can rewrite Eq. (15) using Eq. (13a) as:
x xH ¼ expðAtÞz0 :
ð16Þ
From this, we can see that z0 is essentially the initial condition of the transient response for the system given in Eq. (1). This allows us to rewrite Eq. (16) as:
x ¼ expðAtÞz0 þ xH :
ð17Þ
Looking at the right hand side of Eq. (17), we can see in this form that the series expansion of the equilibrium state given by Theorem 1 is a alternative form of the particular solution for nonhomogeneous linear dynamical systems derived using the concept of equilibrium from dynamical systems theory. From these multiple forms, another important conclusion that can be drawn is that the response of the system is related to the input and its infinitely many derivatives, if they exist. It is not just the value of the input that determines the evolution of the system at steady state or equilibrium. 2.1.2. Existence of Theorem 1 While we have shown computationally that Theorem 1 is an expression of the equilibrium, we have yet to address the issue of the convergence of the series representation given by Theorem 1. In this section, we will provide the conditions in which the series from Theorem 1 converges. Since we already assume A in Eq. (1) is invertible and well-conditioned, we will assume, for simplicity, that it is therefore diagonalizable, and we can simplify this analysis by performing this analysis on a scalar system, given by:
43
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
x_ ¼ ax þ QðtÞ:
ð18Þ
Theorem 1 dictates the equilibrium of Eq. (18) is given by: 1 X Q ðmÞ xH ¼ : amþ1 m¼0
ð19Þ
Now we can apply the ratio test to determine the conditions for convergence of the series. This states that an infinite series converges absolutely if:
Q ðmÞ 1=m lim < 1; m!1amþ1
ð20Þ
where jj represents the absolute value. Eq. (20) simplifies to:
1=m lim Q ðmÞ < jaj:
m!1
ð21Þ
Eq. (21) is the general condition for the existence of Theorem 1. While Eq. (21) is a general result, we can produce more explicit results with more assumptions on the form of Q ðtÞ and its derivatives. First, we examine the case in which there exists a constant Bm > 0 such that:
ðmÞ Q 6 Bm
ð22Þ
for all t 2 ½0; 1Þ for each m 2 ½0; 1Þ. Combining Eqs. (22) and (21) produces:
1=m 1=m lim Q ðmÞ 6 lim Bm < jaj:
m!1
m!1
ð23Þ
Thus, for all Q ðtÞ that satisfies Eq. (22), the criterion for the existence of Theorem 1 is given by:
lim B1=m m!1 m
< jaj:
ð24Þ
We can relax the assumption on Bm to exponentially bounded functions. We define the followin:
Bm ¼ C m expðMm tÞ;
ð25Þ
meaning Q ðtÞ and its derivatives are bounded as:
ðmÞ Q 6 C m expðM m tÞ;
ð26Þ
where M m and C m are positive constants for each m 2 ½0; 1Þ. Substituting Eq. (26) into Eq. (24) produces the following: 1=m lim C m exp
m!1
Mm t < jaj: m
ð27Þ
For Eq. (27) to hold, we can conclude the following criteria must be true:
lim C 1=m m!1 m lim
m!1
< jaj
Mm ¼ 0: m
ð28aÞ ð28bÞ
While we have conducted this analysis on a scalar system, it should be noted these results hold for matrices as well. The only difference between the scalar system proof and the proof for higher order systems would be that we are using the eigenvalues of the matrix A throughout the discussion. Otherwise, the proof holds identically for scalar systems and higher order matrix systems. 2.2. Asymptotic approximation analysis Now that we have derived an expression for the equilibrium of Eq. (1), we want to address the question as to what would occur as the eigenvalues of A are varied. In particular, we are most interested in what occurs as the eigenvalues of A tend to infinity. We restrict our analysis to inputs with derivatives of exponential order. This means that for each Q ðmÞ , there exist positive constants, C m and M m , for each m 2 ½0; 1Þ such that:
ðmÞ Q i 6 C m expðM m tÞ for m 2 ½0; 1Þ8t 2 ½0; 1Þ;
ð29Þ
ðmÞ where Q i is representative of the absolute value of each term of the vector, Q ðmÞ . Additionally, we will assume that A is diagonalizable, and as a result, we can simplify this analysis by examining scalar systems, given by:
44
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
x_ ¼ ax þ Q ðtÞ:
ð30Þ
and inputs as we can perform this analysis term by term in the vectors Q ðmÞ and xH . First, we examine what happens to the magnitude of the equilibrium xH . By using Theorem 1, Eq. (29), and the triangle inequality, we find: 1 H X C m expðM m tÞ x 6 ; jajmþ1 m¼0
ð31Þ
where jxH j represents the absolute value of xH . Now, we can examine what happens to each term at a tends to infinity. For all t 2 ½0; 1Þ, each term of the inequality on the right hand side in Eq. (31) tends to zero as jaj goes to infinity. Thus, we can conclude from Eq. (31) that jxH j goes to zero as jaj tends to infinity. This result is summarized in Proposition 1. Proposition 1. For an input Q such that:
ðmÞ Q i 6 C m expðM m tÞ for m 2 ½0; 1Þ;
ð32Þ
xH goes to zero as the eigenvalues of A tend to infinity. Proposition 1 shows what occurs to xH , but we must determine what occurs to the term AxH as the eigenvalues of A go to infinity. To achieve this, we must rearrange Eq. (12) into the following form:
AxH Q ¼
1 X
A1
m
Q ðmÞ ;
ð33Þ
m¼1
and determine what occurs to the magnitude of AxH Q as the eigenvalues of A tend to infinity. Again, the proof is performed for a scalar system and input given by Eq. (30). Thus, Eq. (33) in scalar form reduces to:
axH Q ¼
1 X
a1
m
Q ðmÞ ;
ð34Þ
m¼1
Applying the Triangle Inequality to Eq. (34), we obtain: 1 H X C m expðMm tÞ ax Q 6 : jajm m¼1
ð35Þ
Again, we can see that for all t 2 ½0; 1Þ, each term of the inequality on the right hand side in Eq. (35) goes to zero as jaj goes to infinity. Thus, from Eq. (35), we can see that axH tends to Q as jaj tends to infinity. This result is summarized in Proposition 2. Proposition 2. For an input Q such that:
ðmÞ Q i 6 C m expðM m tÞ for m 2 ½0; 1Þ;
ð36Þ
AxH goes to Q as the eigenvalues of A tend to infinity. It should be noted that Propositions 1 and 2 also apply to inputs with impulsive derivatives, and further analysis of this is found in [18]. 3. Derivation of the proposed method – ideal case We are now able to pull together all the essential aspects of the prior discussion to formulate a method for input reconstruction for the ideal case. We begin with a linear dynamical system of the form:
x_ ¼ Ax þ Q ðtÞ
ð37aÞ
y ¼ Cx;
ð37bÞ
where Q ðtÞ is an n 1 time-dependent input vector, x is an n 1 state vector, A is an n n matrix, C is an m n matrix, and y is the m 1 measured output vector. To find the relationship between the measured output and the input, we define an estimation system with observer feedback, where we neglect the excitation term Q ðtÞ. This is done to draw a comparison between the excited and unexcited systems, and to produce a relationship between the state and the excitation Q ðtÞ, as the result of the following derivation will show. Thus, we define the estimation system as:
^x_ ¼ A^x þ LðC ^x yÞ;
ð38Þ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
45
where ^ x is the state estimate, and L is the n m observer gain matrix, which is a design parameter. With this form, we are able to generate a direct relationship between the states of the actual and estimation systems and the input of the actual system. Substitution of Eq. (37b) into Eq. (38) produces:
^x_ ¼ A^x þ LCð^x xÞ:
ð39Þ
x x, we have the following error system: Defining e ¼ ^
e_ ¼ ðA þ LC Þe QðtÞ:
ð40Þ
For an observable system, we can design the eigenvalues of A þ LC with an appropriately chosen L [17]. As a result, we can design the eigenvalues of A þ LC such that the matrix is well-conditioned and invertible. This would allow us to use Theorem 1 to determine the equilibrium state of the system defined by Eq. (40), which would give us:
eH ¼
1 mþ1 X ðA þ LC Þ1 Q ðmÞ ;
ð41Þ
m¼0
where eH is the equilibrium state of Eq. (40). The parameter L could be designed such that: the real part of the eigenvalues of A þ LC are negative so that the equilibrium state of the error system is asymptotically stable, and the magnitude of the eigenvalues of A þ LC are sufficiently large in magnitude such that: 1 m X ðA þ LC Þ1 Q ðmÞ 0:
ð42Þ
m¼1
We can apply Proposition 2, and as the magnitudes of the eigenvalues of A þ LC go to infinity, we have:
ðA þ LC ÞeH ! QðtÞ;
ð43Þ
^: and we could define our estimate for the input Q
b ¼ ðA þ LC ÞeH : Q
ð44Þ H
Additionally, from Proposition 1, we know that e goes to zero as the eigenvalues of A þ LC tend to infinity. Since the terms of the matrix L are of sufficiently large magnitude, we can also conclude in Eq. (44) that AeH is negligible in relation to LCeH , and we can define our estimate for the input as:
^ ¼ LCeH : Q
ð45Þ
This can be rewritten as:
^ ¼L y ^H yH : Q
ð46Þ
3.1. Force reconstruction procedure Now that we have an expression for the input estimate, we can discuss the process of estimating the input. To better visualize what is occurring in the method, the block diagram in Fig. 1 shows us that the method consists of three fundamental steps. First, the system response y must be measured. Once the measurement is made, this is substituted into the unexcited b to x. Then, the measurement y and the estimate state ^ x are substituted into Q estimate system, resulting in the estimate state ^ determine the estimate for the system input. 3.2. Interpretation The method developed here defines an estimate system with observer feedback, and in that sense, it is comparable to the Luenberger observer. The two methods differ in their assumptions and application. The Luenberger observer assumes knowledge of the input and measured output in order to determine an estimate for the state, but in the proposed method, we know the measured output in order to determine an estimate for the input. As a result, the proposed method is not a Luenberger observer. It may be considered to be a complement to the Luenberger observer, and an extension of the application of observer feedback. Additionally, the new application of observer feedback in the method includes new interpretations for the parameters. The matrix L consists of design parameters and the matrix C acts to filter terms of the error equilibrium, eH , that we cannot
46
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Fig. 1. Block diagram of force reconstruction method.
observe or measure. When functioning in tandem with each other in the matrix product LC, we can say that LC acts as a high gain observer filter to approximate the input. It achieves this in two ways: as a filter of the derivatives of the input, and as a weighting function to assign more value to the observed terms in relation to the unobserved terms. Its ability to act as a filter is apparent in Eq. (33). As the eigenvalues of A þ LC go to infinity, each term in the series with Q ðmÞ for m 2 ½1; 1Þ becomes increasingly negligible. As a result, the terms in the series with the input differentiated are filtered, leaving us with the original input Q ðtÞ. To see the second purpose of L, we examine what is occurring within the matrix A þ LC. The eigenvalues of A þ LC control the magnitude of eH , thus, for larger eigenvalues, eH becomes predictably smaller. And since the matrix L controls the eigenvalues of A þ LC, then the terms of L are large. As L works in tandem with C, and C is related to the measurements, then LC acts to assign more weight to the measurements. As the terms of L become large, the term AeH becomes negligible in the estimation of the input.
4. Example: Single degree–of–freedom system To detail the manner in which the proposed method would be applied, numerical simulations were conducted in [18] for a scalar system, a single degree–of–freedom system, and a 2 degree–of–freedom system. We present here the application of the proposed method to a single degree–of–freedom example for the system shown in Fig. 2. In this example, we assume particular inputs and use these to find the outputs. It is the output that are the ‘‘data” used in the proposed method to reconstruct the inputs.
47
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Fig. 2. Single degree-of-freedom dynamical system.
We select the particular single degree-of-freedom system:
€ þ 4y_ þ 4y ¼ FðtÞ; y
ð47Þ
with the observed displacement:
z ¼ y:
ð48Þ
_ resulting in: Our method requires that we first rewrite the system in the form of Eq. (37) by defining x1 ¼ y and x2 ¼ y,
A¼ Q¼
0
1
ð49aÞ
4 4
0
ð49bÞ
FðtÞ
C ¼ ½1 0
ð49cÞ
z ¼ Cx
x1 : x¼ x2
ð49dÞ ð49eÞ
We can define our estimation system by:
^x_ ¼ A^x þ LðC ^x zÞ;
ð50Þ
where ^ x is the estimation state. Since we are observing the displacement and z is a scalar quantity, L takes the form:
L¼
l1 l2
:
ð51Þ
We can define the error system as e ¼ ^ x x, producing:
e_ ¼ ðA þ LC Þe Q:
ð52Þ
According to our method for an observable system, we can design L such that the matrix A þ LC is well-conditioned so that A þ LC is invertible, and the eigenvalues of A þ LC sufficiently large such that:
Q LCeH ;
ð53Þ
^: and we can use Eq. (45) to define our force estimate Q
^ ¼ LCeH : Q
ð54Þ
Due to the size of the system, we are able to design the eigenvalues through direct computation, but it should be noted there are other methods that can be employed to design the eigenvalues which would more readily lend itself to systems of higher dimension, such as the Bass-Gura method or linear quadratic estimator method [20]. To design the eigenvalues of A þ LC, we begin with the eigenvalue problem for A þ LC:
det ½kI 22 ðA þ LCÞ ¼ 0;
ð55Þ
where k represents the eigenvalues of the matrix A þ LC and:
A þ LC ¼
l1
1
4 þ l2
4
:
ð56Þ
48
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
The eigenvalues of A þ LC are determined by finding the zeroes of Eq. (55), which when expanded produces the following:
2 1 1 2 k þ ð4 l1 Þ ð4 l1 Þ 4l1 þ 4 l2 ¼ 0; 2 4
ð57Þ
and are specified by assigning values to the parameters l1 and l2 . For simplicity, we choose the eigenvalues that produce the double-root k ¼ ð4 l1 Þ=2. Thus, we can define l2 in terms of l1 :
1 2 l2 ¼ ð4 l1 Þ 4l1 þ 4 4
ð58Þ
and we can assign values to the parameter l1 . Eqs. (47) and (50) were simulated with the input: FðtÞ ¼ 10 sinð2tÞ N, FðtÞ is a sawtooth wave with period 2p and amplitude of 10 N, and with the various values of l1 : l1 ¼ 10 s1, l1 ¼ 100 s1, l1 ¼ 1000 s1, to exhibit the convergence characteristics of the proposed estimation method. The frequency for the harmonic applied force was chosen to coincide with the natural frequency of the system to show that the proposed method is applicable at the natural frequency of the system. The sawtooth wave was chosen to highlight the applicability of the method to an applied force with rapid changes in amplitude. To obtain the displacement data, Eqs. (47) and (50) were simulated in Matlab with the prescribed input FðtÞ and system parameters defined in A and L. The total period of simulation was 10 s, and was divided into segments of .0001 s each. Using a fifth-order Dormand-Prince method given by the Matlab function ode5, the terms, x and ^ x, were solved at 51 points between time segments. Beginning at time t ¼ 0 with initial conditions set as zero and using the fifth-order DormandPrince method, we solved for final condition of the first segment, which ends at time t ¼ :0001. This is then the initial condition of the next segment, and this process is continued for the total period of 10 s. This provided us with the system state x, which consists of the displacement and velocity of the mass. To obtain the displacement data, we multiplied x by C to determine z. This was substituted into Eq. (54), and can be written as:
^ ¼ LC ð^x xÞ: Q
ð59Þ
x and x by C and using Eq. (51), Eq. (59) can be expressed as: Multiplying ^
^ ¼ l1 ð^x1 x1 Þ: Q l2
ð60Þ
Using:
" ^ ¼ Q
# b1 Q ; b2 Q
ð61Þ
^ approximates Q , our approximation for FðtÞ is given by: and since Q
b 2 ¼ l2 ð^x1 x1 Þ: Q
ð62Þ
The systems were simulated for each value of l1 , leading to a value of l2 via Eq. (58), and the results for each of these cases are shown in Figs. 3 and 4. It is readily observed that as the magnitude of l1 is increased, the estimate for the input converges to the actual input. It should be noted that in Fig. 4, the initial value of the force estimate is does not agree as well. While this is dependent on the initial conditions, if the magnitude of l1 is further increased, the initial difference will remain but will converge to the applied force rapidly. Also, we should note that the application of the proposed method is inherently asymptotic, meaning that as the magnitude of l1 goes to infinity Eq. (62) will converge to a particular set of values. And while the values of l1 were selected to highlight the asymptotic nature of the proposed method, l1 can be selected through the assignment of a convergence criteria. As b 2 will tend to zero, and once it reaches an assigned threshold value, this value the magnitude of l1 increases, the variation in Q of l1 will be the determined value for acceptable convergence. From these results, we are able to conclude that the proposed method is able to provide an estimate of the input based only on the observed parameters. This is an important result because it shows that we are able to produce a relationship
49
Applied Force (N)
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
l1= −10 10 Actual Estimate
5 0 −5 −10 0
1
2
3
4
5
6
7
8
9
10
Applied Force (N)
Time (s)
l1= −100 10 Actual Estimate
5 0 −5 −10 0
1
2
3
4
5
6
7
8
9
10
Applied Force (N)
Time (s)
l1= −1000 10 Actual Estimate
5 0 −5 −10 0
1
2
3
4
5
6
7
8
9
10
Time (s) Fig. 3. Single degree-of-freedom system sinusoidally forced.
Applied Force (N)
l1= −10 10 Actual Estimate
5 0 −5 −10
0
1
2
3
4
5
6
7
8
9
10
Time (s)
Applied Force (N)
l1= −100 10 Actual Estimate
5 0 −5 −10
0
1
2
3
4
5
6
7
8
9
10
Time (s)
Applied Force (N)
l1= −1000 10 Actual Estimate
5 0 −5 −10
0
1
2
3
4
5
6
7
Time (s) Fig. 4. Single degree-of-freedom system with a sawtooth wave force.
8
9
10
50
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
between the observed measurements and the input directly. Additionally, it shows that the estimate for the input can be independent of the system parameters, given by A. 5. Application to continuous systems using modal analysis So far, we have developed the theory and applied the proposed method to a single degree-of-freedom system. In this section, our aim is to apply the proposed method to more complex systems by estimating loads on continuous systems. While the scope of this section is strictly limited to modal analysis as the method of discretizing a continuous system, we are confident that the proposed estimation method can work in tandem with other commonly accepted discretization methods, such as Rayleigh-Ritz or finite element. This ability to work with discretization methods broadens the applicability of the proposed estimation method. It is not limited in scope to simple models, but is capable of being applied to more complex structures or applications. To better express the development of the proposed method with modal analysis, we apply the proposed method to the example of estimating the force applied to an Euler-Bernoulli beam, which is governed by the following equation:
qA
@2w @w @4w þc þ EI 4 ¼ pðx; tÞ; 2 @t @x @t
ð63Þ
where wðx; tÞ is the transverse displacement of the beam at location x along the beam and at time t; pðx; tÞ is the transverse load along the length of the beam, E is the Young’s modulus, I is the moment of inertia of the beam, q is the material density, A is the cross sectional area, and c is the viscous damping coefficient [3]. In the example presented in this section, we are assuming the full displacement field of the beam is known. While this assumption is not realistic, we use this assumption strictly in the context of demonstrating the application of the proposed method to continuous systems. As the proposed method is developed with the aim of estimating the temporal component of a dynamic force, our objective in this section is to show that the proposed method is capable of being applicable to continuous systems using discretization methods to estimate time-varying loads. In reality, we would only be able to take measurements at discrete points along the surface of the beam. In this context, the proposed method would have to be merged with a modal observer, such as that used in [6], but by assuming access to the continuous displacement field, we are able to develop the theory and determine the behavior of the proposed method when used with more complex structures. To develop the theory, the proposed method will be applied to estimate the loads as shown in Figs. 5 and 6, and given by the following:
x pðx; tÞ ¼ 100 CðtÞ l pðx; tÞ ¼ 100dðx aÞCðtÞ; where dðx aÞ is the Dirac delta function, l is the length of the beam, a is a point on the beam, and:
Fig. 5. Simply supported Euler-Bernoulli beam of length l with a linearly distributed load.
ð64aÞ ð64bÞ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
51
Fig. 6. Simply supported Euler-Bernoulli beam of length l with a point load.
CðtÞ ¼ sinð2tÞ:
ð65Þ
Loads characterized by Eqs. (64a) and (64b) were chosen specifically to highlight certain attributes of modal analysis. In Section 5.1, we will describe the theoretical analysis of estimating the load using modal analysis. In Section 5.2, we will derive an exact solution to the estimation of the temporal components of the distributed load, and in Section 5.3, we present results of the simulation of this method and discuss the important attributes of these simulations. 5.1. Modal asymptotic approximation: theory To apply the proposed method to more complex structures using modal analysis, we begin by expanding the displacement, wðx; tÞ, and the distributed load, pðx; tÞ in terms of orthogonal functions f n ðxÞ, where:
Z
l
f n ðxÞf m ðxÞdx ¼ dnm ;
ð66Þ
0
and f n ðxÞ satisfy the given boundary conditions. This property is critical to the application of the proposed method as it enables us to transform the governing equation, given by Eq. (63), into the necessary form. Using the orthogonality property of Eq. (66), the parameters wðx; tÞ and pðx; tÞ can be expanded in infinite series as: 1 X
wðx; tÞ ¼
yn ðtÞf n ðxÞ;
ð67aÞ
1 X Q n ðtÞf n ðxÞ;
ð67bÞ
m¼1
pðx; tÞ ¼
m¼1
where:
Z
l
yn ðtÞ ¼
wðx; tÞf n ðxÞdx
ð68aÞ
0
Z Q n ðtÞ ¼
0
l
pðx; tÞf n ðxÞdx:
ð68bÞ
By multiplying Eq. (63) by f n ðxÞ, integrating over the domain, and using Eqs. (67a) and (67b), we find:
qAy€n þ cy_ n þ EI
np4 l
yn ¼ Q n ðtÞ:
In this form, we are able to apply our reconstruction method to each mode yn . We need to put Eq. (69) into the same form as Eq. (37) by defining X 1n ¼ yn and X 2n ¼ y_ n , yielding:
ð69Þ
52
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
" An ¼ Xn ¼
0 EI np 4
qA
l
X 1n ðtÞ
#
1 qcA
ð70aÞ ð70bÞ
X 2n ðtÞ " # 0 Qn ¼ 1 : qA Q n ðtÞ
ð70cÞ
Since the displacement is the parameter being observed, it needs to be put in terms of the modal parameters. We define the measurement as z ¼ wðx; tÞ. To obtain the modal measurement, we need to multiply z by f n ðxÞ and integrate over the domain to produce zn ¼ yn . Rearranging zn to be in the same form as Eq. (37), we obtain:
zn ¼ CX n ;
ð71Þ
C ¼ ½ 1 0 ;
ð72Þ
where
and zn is the modal measurement. Now we can apply the reconstruction method for each mode, which begins by defining an unperturbed estimation system with observer feedback for each mode by:
^ ^ n þ Ln C X ^ n Xn ; X_ n ¼ An X
ð73Þ
where Ln is a design parameter. The nth mode error system is:
E_ n ¼ ðAn þ Ln C ÞEn Qn ;
ð74Þ
b n X n . For well-conditioned An þ Ln C with sufficiently large eigenvalues, we can estimate Qn by: where En ¼ X
^ n ¼ Ln CEn ; Q
ð75Þ
^ n is the estimate of Qn . To obtain an approximation for the input, we must multiply Eq. (75) by f ðxÞ and sum over where Q n the domain n from 1 to infinity. As a result, our estimate becomes:
^ tÞ ¼ Qðx;
1 X ^ n ðtÞf ðxÞ; Q n
ð76Þ
n¼1
^ is the estimate for: where Q
"
Qðx; tÞ ¼
0 1
qA pðx; tÞ
#
:
ð77Þ
We have approximated the input vector Q, but we still have not quite approximated the distributed load pðx; tÞ. To do so, we must acknowledge that Q consists of the term pðx; tÞ=ðqAÞ. Thus, to obtain the distributed load, we must define a vector P such that:
P ¼ qAQ;
ð78Þ
where P contains the distributed load, and our estimate for the distributed load becomes:
^ ^ ¼ qAQ; P
ð79Þ
^ is the estimate for the distributed load, pðx; tÞ. where P 5.2. Modal asymptotic approximation: exact solution Now that we have outlined the theory, we can demonstrate its usage by applying it to the linearly-varying load from Fig. 5. Our aim in this section is to develop an exact solution so that we may be able to better develop the analysis of the proposed method in this type of application. For this example, the Euler-Bernoulli beam is simply supported, which gives us the following boundary conditions:
@2w ð0; tÞ ¼ 0 @x2 2 @ w wðl; tÞ ¼ 2 ðl; tÞ ¼ 0: @x
wð0; tÞ ¼
ð80aÞ ð80bÞ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
53
The initial conditions for the undeformed beam are assumed to be:
wðx; 0Þ ¼
@w ðx; 0Þ ¼ 0: @x
ð81Þ
We are developing the solution for the applied load given by Eq. (64a), and are assuming the transverse displacement field wðx; tÞ is measured. First, we must determine the mode shape, f n ðxÞ. For the given boundary conditions, the mode shape is given by:
rffiffiffi np 2 sin f n ðxÞ ¼ x : l l
ð82Þ
Then, we must select the matrix LC. Due to the size of the system, we are able to design the eigenvalues through direct computation, but we should again note there are other methods which can be employed to design the eigenvalues which would more readily apply to systems of higher dimension, such as the Bass-Gura method or linear quadratic estimator method [20]. Since we are assuming that we can only observe the displacement field of the beam, Ln C takes the form:
Ln C ¼
ln1
0
ln2
0
:
ð83Þ
We need to design the terms ln1 and ln2 such that the eigenvalues of An þ Ln C are sufficiently large with a negative real part to ensure asymptotic stability. The process to determine the eigenvalues as a function of l1 is identical to that of the example for a single degree-of-freedom system in Section 4. Using Eq. (55), the eigenvalues for each mode are computed from the characteristic polynomial:
det ½kI 22 ðAn þ Ln CÞ ¼ k2 þ
c c EI np4 þ ln2 : ln1 k ln1 qA qA qA l
ð84Þ
We determine the eigenvalues by finding the zeros of Eq. (84), and since ln1 and ln2 are design parameters, we can parameterize these to be functions of a single variable, b, such that ln1 ðbÞ and ln2 ðbÞ. To achieve this, we first need the quadratic equation with the desired zeros k ¼ a and k ¼ b. This is given by:
k2 þ ða þ bÞk þ ab ¼ 0:
ð85Þ
From here, we equate the coefficients of Eqs. (84) and (85) as follows:
c aþb¼ l qA n1 c EI np4 ab ¼ þ ln2 : ln1 qA qA l
ð86aÞ ð86bÞ
To further simplify our eigenvalue selection, we define a ¼ 2b to produce:
3b ¼
c
ln1 c EI np4 2 2b ¼ þ ln2 : ln1 qA qA l
qA
ð87aÞ ð87bÞ
This gives us the values for ln1 and ln2 in terms of b:
ln1 ¼
c
qA
3b
ð88aÞ
2
ln2 ¼ 2b
c c EI np4 : 3b þ qA qA qA l
ð88bÞ
For simplicity, we define:
EI np4 l c c ^ 2 ¼ 2b2 þ x 3b : qA qA
x2n ¼
qA
ð89aÞ ð89bÞ
Using Eqs. (89), ln2 can be rewritten as:
^ 2: ln2 ¼ x2n x
ð90Þ
This allows us to use only a single parameter, b, to define the eigenvalues of An þ Ln C. Using Eqs. (88)–(90), Ln C becomes:
54
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
" Ln C ¼
qA 3b c
0
#
^2 0 x2n x
:
ð91Þ
From here, we can use Eq. (75) as the estimate for each modal input, and while it can be applied exactly as is, we can gain more insights by returning to the modal governing equation, given by Eq. (69), in combination with the state-space analysis to derive the solution explicitly in closed form. For this analysis, we will provide the solution to pðx; tÞ corresponding to Eq. (64a). By multiplying Eq. (64a) by f n ðxÞ and integrating from 0 to l, we obtain the modal input as follows:
Q n ðtÞ ¼ T n sinðxtÞ; where:
10 Tn ¼ qA
ð92Þ
rffiffiffi 2 cosðnpÞ : l np
ð93Þ
For the initial conditions given by Eq. (81), the solution of Eq. (69) for yn ðtÞ is given by:
yn ðtÞ ¼
c
x2n x2 2
x2n x2 þ
c qA
2 T n sinðxtÞ
x
qA x
x2n x2 2 þ
c qA
x
2 T n cosðxtÞ
ð94Þ
Using Eq. (92), Eq. (94) can be rewritten in terms of Q n ðtÞ as:
yn ðtÞ ¼
c
x2n x2 2
x2n x2 þ
c qA
2 Q n ðtÞ
x
qA x
x2n x2 2 þ
c qA
x
2 T n cosðxtÞ:
ð95Þ
Now, we must rewrite the estimate system, given by Eq. (73), in terms of yn ðtÞ. To do so, we must first rewrite Eq. (73) as individual equations. Using Eq. (91), this produces:
c b n2 þ 3b c X n1 b n1 þ X 3b X qA qA 2 c 2 2 b n1 b þ x ^ X ^ xn X n1 : ¼ x X qA n2
b_ n1 ¼ X b_ n2 X
ð96aÞ ð96bÞ
b n2 and X b n1 are estimates for X n2 and X n1 , respectively, and since X n1 ¼ y , we can conclude that X b n1 ¼ y ^n and Since X n _ € b ^n . This means we can express Eq. (96b) as: X n2 ¼ y
€^ þ y n
c
qA
2 ^ 2y ^ x2n yn ðtÞ: ^n ¼ x ^n þ x y
ð97Þ
^n is given If we assign the initial conditions of zero for the estimation system and use Eq. (95), the solution of Eq. (97) for y by:
2 3 2 2 2 c 2 2 2 2 2 2 ^ ^ ^ x x x x x x x x x 6 n n n 7 qA ^n ðtÞ ¼ 6 7Q ðtÞ y 2 4 2 c 2 5 n 2 c 2 2 2 2 ^ x Þ þ qA x ðx xn x þ qA x 3 2 2 2 c c 2 2 2 2 2 2 ^ ^ ^ x x x x x x x x x x þ 7 6 n n n qA qA 6 2 2 7 5T n cosðtÞ: 4 2 2 ^ 2 x2 Þ þ qcA x ðx x2n x2 þ qcA x
ð98Þ
Using Eq. (75), the estimate for each modal input is expressed as:
b n ¼ ln2 ðy ^n yn Þ: Q
ð99Þ
Substitution of Eqs. (88b), (95), and (98) into Eq. (99) and simplifying produces:
2 2 2 c ^2 ^ x2n x ^ x2 x qA x xn x b Q n ðtÞ ¼ 2 Q n ðtÞ þ 2 T n cosðxtÞ: ^ 2 x2 Þ2 þ qcA x ^ 2 x2 Þ2 þ qcA x ðx ðx
The distributed load, pðx; tÞ, is estimated by multiplying Eq. (100) by f n ðxÞ and summing from 1 to infinity:
ð100Þ
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
2 3 2 2 2 c 2 1 2 2 ^ X x x x ^ ^ x x x x n qA 7 n ^ðx; tÞ ¼ qA 6 p 4 2 Q n ðtÞ þ 2 T n cosðxtÞ5f n ðxÞ: n¼1 ðx ^ 2 x2 Þ2 þ qcA x ^ 2 x2 Þ2 þ qcA x ðx
55
ð101Þ
From Eqs. (89) and (100), we can see several important attributes of the application of the proposed method with modal b n to Q as b goes to infinity. As b goes to infinity in Eq. analysis. First, this provides direct evidence of the convergence of Q n 2 2 b ^ ^ (89b), x goes to infinity, and in Eq. (100), as x goes to infinity, Q n converges to Q n . In the context of the continuous problem of an Euler-Bernoulli beam using modal analysis, we can see convergence occurs for the temporal component of each mode in the ideal case as b goes to infinity. Additionally, we can see the connection between the selection of a finite value for b and the modal natural frequency. Since we’re choosing a finite value for b, there can be a coupling between b and the number of modes in which we choose in the approximation of the load, pðx; tÞ. This is apparent in Eqs. (89) and (100). In Eq. (100), the coefficient of Q n is dependent ^ 2 is a function of b and x2n is a function of n. From Eq. (89a), we can see that x2n grows with ^ 2 x2n , where x on the term x each additional term used in the modal sum. For a large but finite value for b, each additional mode used in the approxima^ 2 x2n . From this, we can conclude that the number of tion of pðx; tÞ will result in a decrease in the magnitude of the term x modes used to approximate pðx; tÞ should be taken into account during the selection of the parameter b. In practice, the value bn of b can be selected through the assignment of a convergence criteria. As the magnitude of b increases, the variation in Q will diminish, and once it reaches an assigned threshold value, this value of b will be the determined value for acceptable convergence. This procedure should not need to be applied to every mode in the approximation of pðx; tÞ. It can be applied to the mode with the largest modal natural frequency, and once a suitable value for b is determined for this mode, it should be suitable for all other modes. 5.3. Modal asymptotic approximation: numerical simulations To exhibit the application of the proposed method in estimating the load applied to an Euler-Bernoulli beam, we applied the results from Section 5.2 to show the resulting estimation of the distributed load. This will be performed for the linearly varying load shown in Fig. 5, and for a point load as shown in Fig. 6. These two spatial distributions were chosen to highlight certain attributes of modal analysis. The linearly distributed load is an odd function, and the point load is not. As a result, the characteristics of the functions must be taken into account. To run the simulations, the following physical parameters are defined:
EI ¼ 107 Nm2, qA ¼ 0:1 kg/m, l ¼ 10 m, c ¼ 4 Ns/m2, x ¼ 2 radians/s.
Additionally, for the linearly varying load, we varied the value for b to show how the estimate varies with a fixed modal summation, and for the point load, we varied the number of modes summed to show how the estimate varies for fixed b. As we have already shown the convergence in the temporal domain, we are only showing the convergence of the estimate for pðx; tÞ in the spatial domain. b n ðtÞ. Eq. (82) was evalTo obtain data, we computed Eqs. (95) and (98), and substituted them into Eq. (99) to determine Q b n ðtÞ. This was performed for n 2 ½1; 20, and then uated at 101 points on the interval ½0; d, and we multiplied this by each Q summed to obtain the estimate for pðx; tÞ. This was performed for the values of b: b ¼ 10; 000 s1, b ¼ 100; 000 s1, b ¼ 1; 000; 000 s1. The results for each b are shown in Fig. 7, and an enlarged view of b ¼ 1; 000; 000 s1 is shown in Fig. 9. As we can see from Fig. 7, as we increase the value of b, the estimate for the distributed load converges actual distributed load. This actually agrees very well with the theory we derived in Section 5.2. If we define the coefficient for Q n in Eq. (100) as follows: (Fig. 10)
bn ¼
^ 2 x2n x ^ 2 x2 x 2
^ 2 x2 Þ þ ðx
c
qA x
2 ;
ð102Þ
we can use this relation as the metric for the quality of the estimation. To rewrite bn in terms of b, we substitute Eqs. (89) into Eq. (102):
56
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
b= 10,000
4
Estimated Distribution (N/m)
1.5
x 10
1 0.5 0 −0.5 −1 −1.5
0
1
2
3
4
5
6
7
8
9
10
6
7
8
9
10
6
7
8
9
10
Distance (m)
Estimated Distribution (N/m)
b= 100,000 200
100
0
−100
−200
0
1
2
3
4
5
Distance (m)
Estimated Distribution (N/m)
b= 1,000,000 150 100 50 0 −50 −100 −150
0
1
2
3
4
5
Distance (m) Fig. 7. Estimated linear load distribution under various values of b.
bn ¼
2
2b þ qcA
2 2 2 4 2 2b þ qcA qcA 3b 3b qEIA nlp x2 : 2 2 2 2 x2 þ qcA x 2b þ qcA qcA 3b
c qA
ð103Þ
From Eq. (103), we can see that as b goes to infinity, bn goes to 1. So, the closer to 1 bn is, the better the approximation we obtain. For n ¼ 20, Eq. (103) can be evaluated for each value of b. For b ¼ 10; 000 and b ¼ 100; 000; bn is on the order of 773 and 6:8 respectively, which indicates these values are severely underperforming. But for b ¼ 1; 000; 000; bn ¼ 0:922, indicating a high level of agreement for the 20th mode. This is readily seen in Fig. 7, as the first two graphs show poor results and the third graph exhibits a high level of convergence. For the case of b ¼ 1; 000; 000, the load distribution peaks at 112 N/m while the actual peak value is 100 N/m. This means there is only a 12% difference between the estimate peak and actual peak values, and higher accuracy could be achieved by increasing the value of b and the number of modes used in the estimation. Now that we have seen the behavior of the proposed method as b varies, we want to see the behavior of the proposed method as the number of modes varies. For this case, we will apply the proposed method to the estimation of a point load. As we have derived everything we need regarding the linearly distributed load in Section 5.2, we show here only the nec-
57
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Estimated Distribution (N/m)
10 Modes Estimation 1500 1000 500 0 −500 −1000 −1500 0
1
2
3
4
5
6
7
8
9
10
7
8
9
10
7
8
9
10
Distance (m)
Estimated Distribution (N/m)
15 Modes Estimation 1500 1000 500 0 −500 −1000 −1500 0
1
2
3
4
5
6
Distance (m)
1500 1000 500 0 −500 −1000 −1500 0
1
2
3
4
5
6
Distance (m) Fig. 8. Estimated point load under various modal sum conditions.
Estimated Distributed Load with b=1,000,000 150
Estimated Distribution (N/m)
Estimated Distribution (N/m)
20 Modes Estimation
100 50 0
−50 −100 −150 0
2
4
6
8
Distance (m) Fig. 9. Estimated load distribution with b ¼ 1; 000; 000.
10
58
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
Estimated Point Load with b=1,000,000
Estimated Distribution (N/m)
2000
1000
0
−1000
−2000
0
2
4
6
8
10
Distance (m) Fig. 10. Estimated point load with b ¼ 1; 000; 000 and 20 modes summed.
essary results. We repeated the procedure performed in Section 5.2 for a point load at a distance a ¼ 5 m from the end of the b n ðtÞ are presented below: ^n ðtÞ, and Q beam, which is given by Eq. (64b). The results for Q n ðtÞ; yn ðtÞ; y
Q n ðtÞ ¼ T n sinðxtÞ rffiffiffi anp 10 2 sin Tn ¼ qA l l yn ðtÞ ¼
ð104aÞ ð104bÞ c
x2n x2
x2n x2 2 þ
c
2 T n sinðxtÞ
qA x
qA x
x2n x2 2 þ
c
qA x
2 T n cosðxtÞ
2 3 2 2 2 c 2 2 2 2 2 2 ^ ^ ^ x x x x x x x x x 6 n n n 7 qA ^n ðtÞ ¼ 6 7Q ðtÞ y 2 4 2 c 2 5 n 2 c 2 2 2 ^ 2 x Þ þ qA x ðx xn x þ qA x 3 2 2 2 c c 2 2 2 2 2 2 ^ ^ ^ 6 x x x xn qA x þ qA x x xn xn x 7 6 2 2 7 5T n cosðtÞ: 4 ^ 2 x2 Þ2 þ qcA x ðx x2n x2 2 þ qcA x 2 2 2 c 2 2 ^2 ^ ^ qA x xn x b n ðtÞ ¼ x xn x x Q ðtÞ þ Q 2 T n cosðxtÞ: n 2 ^ 2 x2 Þ2 þ qcA x ^ 2 x2 Þ2 þ qcA x ðx ðx
ð104cÞ
ð104dÞ
ð104eÞ
We repeated the procedure of the linearly distributed load to obtain data, but we computed Eq. (104c) to determine the modal displacements and (104d) for the modal estimate. In this case, we kept b constant at the value of 1; 000; 000 s1, but varied the number of modes used in the estimation of pðx; tÞ. The number of modes used are as follows: n ¼ 10, n ¼ 15, n ¼ 20. As we can see in Fig. 8, the number of modes greatly affects the estimate for the point load. While 20 modes were adequate for the linearly varying load, the point load has is not as fortunate. While we can see some sense of convergence, some concerns do remain. For starters, the width of the peak grows more narrow with an increase in the number of modes used, but it might require a much larger value for n before an acceptable estimate for the applied load is computed. In addition to the narrowing of the peak, a significant number of minor peaks surround the peak at the location of the point load. The difference in the convergence between the linearly varying load and the point load is a consequence of the agreement of the types of loads with the types of problem. A simply-supported beam is essentially an odd-functioned problem, and the best result will come from odd-functioned loads, such as the linearly varying load. The point load is not as fortunate in that it is not a purely odd function, and convergence to the true value can be limited in this respect. This is one of the limitations of applying the proposed method with modal analysis to estimate distributed loads. For the linearly-varying load, modal analysis is a perfectly acceptable method of discretizing the continuous problem, but a more appropriate method of discretization would be required to better estimate the point load.
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
59
6. Derivation of proposed method – real case Thus far, we have assumed perfect knowledge of the response in our measurement, but it is important to understand what would occur in the likely event that uncertainty in the measurements is present. To make such a determination, we must alter the original system given by Eq. (37) to include measurement noise as:
x_ ¼ Ax þ Q ðtÞ
ð105aÞ
y ¼ Cx þ g;
ð105bÞ
where g is an m 1 measurement noise vector. Using the estimate system given by Eq. (39) in conjunction with Eq. (105), the error system becomes:
e_ ¼ ðA þ LC Þe QðtÞ Lg:
ð106Þ
We apply Theorem 1 and rearrange to produce:
ðA þ LC ÞeH ¼ Q ðtÞ þ Lg þ ðA þ LC Þ1 Q_ ðtÞ þ ðA þ LC Þ1 Lg_ þ
1 h X m¼2
ðA þ LC Þ1
im dm h i m Q ðtÞ þ Lg : dt
ð107Þ
By applying Propositions 1 and 2 to Eq. (107), we produce the following approximation:
ðA þ LC ÞeH Q ðtÞ þ Lg þ ðA þ LC Þ1 Lg_ :
ð108Þ
The term ðA þ LC ÞeH can be further simplified. So far, the proposed method has been discussed only in terms of the error system, but it is just as applicable to the estimate system given by Eq. (39). When we apply Theorem 1 and Propositions 1 and 2 to Eq. (39), we obtain:
^xH ðA þ LC Þ1 Ly:
ð109Þ
Using Eqs. (105b) and (109), we can further simplify Eq. (108) to produce the approximation:
^ y Q ðtÞ AðA þ LC Þ1 Lg þ ðA þ LC Þ1 Lg_ : L y
ð110Þ
In the form given by Eq. (110), we can define our estimate by:
b g ðtÞ ¼ L y ^y ; Q
ð111Þ
^ ðtÞ is the estimate of the input with measurement noise. where Q g 6.1. Discussion The concepts developed in Section 2 provided useful insight into the inverse problem. While measurement noise remains an issue, the concepts of the proposed method allow us redefine the force reconstruction problem in a different context. The concepts developed allow the force reconstruction problem to be defined as a time series analysis problem. In this treatment of measurement noise, alternative signal processing techniques can be applied to Eq. (111) to reduce the effects of measurement noise and recover an estimate of the original input. In addition, the recovery of the forces is not limited to strictly forces in the time domain, but the concepts herein can be extended further. For example, if random forces are present, the objective would be to recover the spectral density of the random forces applied to the system. The current formulation readily lends itself to such problems, and additional work in the future will be performed to address random forces from a probabilistic perspective. 6.2. Example: Exact solution of a scalar system To show that the approximation given by Eq. (110) holds analytically, we apply the proposed method to the scalar system:
x_ ¼ Ax þ expðMtÞ y ¼ Cx þ g;
ð112aÞ ð112bÞ
where M is a non-zero constant, A < 0, and to simplify the analysis, we let C ¼ 1. Since our analysis does not make any assumptions on the form of g, we define g for the purposes of demonstration as:
g ¼ B expðixtÞ;
ð113Þ
60
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
pffiffiffiffiffiffiffi where B and x are constants and i ¼ 1. This definition of g will allow us to produce an exact analysis of the system using the proposed method, and as a result, it will allow us to better demonstrate the characteristics of the proposed method and the concepts developed. At equilibrium, Eq. (112b) becomes:
y¼
expðMtÞ þ B expðixtÞ: MA
ð114Þ
To apply the proposed method, the following estimate system must be defined:
^x_ ¼ A^x þ Lð^x yÞ:
ð115Þ
Using Eq. (114), we apply Theorem 1 to Eq. (115) to produce:
^xH ¼
1 X ðA þ LÞðmþ1Þ m¼0
L m Mm expðMtÞ þ BLðixÞ expðixtÞ : MA
ð116Þ
Using Eqs. (114) and (116), we can apply Eq. (111) to produce our estimate. This yields:
bg ¼ Q
L Lðix AÞ expðMtÞ þ B expðixtÞ: AþLM A þ L ix
ð117Þ
For L sufficiently large, it can be said that:
b g expðMtÞ þ ðix AÞB expðixtÞ: Q
ð118Þ
Now if we use the right-hand side of Eq. (110) with g defined as in Eq. (113), we proceed as follows:
AL Lix B expðixtÞ þ B expðixtÞ AþL AþL Lðix AÞ B expðixtÞ ¼ expðMtÞ þ AþL expðMtÞ þ ðix AÞB expðixtÞ:
QðtÞ AðA þ LÞ1 Lg þ ðA þ LÞ1 Lg_ ¼ expðMtÞ
ð119aÞ ð119bÞ ð119cÞ
Comparing Eqs. (118) and (119c), we can see the approximation given by Eq. (110) holds and that our approximation defined by Eq. (111) behaves as predicted. 6.3. Example: Single degree-of-freedom (continued) To see the application of the proposed method to force estimation with realistic measurement noise, we return to the single degree-of-freedom example from Section 4. While we still assume that we are measuring the displacement of the mass, we amend Eq. (49e) to include noise as follows:
z ¼ Cx þ g;
ð120Þ
where here g is uncorrelated zero-mean Gaussian white noise with variance 1, which was used as we wanted to have a large ratio of variance to maximum amplitude. The rest of the analysis from Section 4 remains unchanged, and this analysis is simulated with the applied force:
FðtÞ ¼ 10 sinð2tÞ N;
ð121Þ
and at values of l1 given by: l1 ¼ 10000 s1, l1 ¼ 100000 s1, l1 ¼ 1000000 s1. With the prescribed FðtÞ, the max amplitude of the displacement of the mass is 1.49 m. This gives a ratio of noise variance to max amplitude of 0.6727. To obtain the displacement data, Eqs. (47) and (50) were simulated in Matlab with the prescribed force FðtÞ and system parameters defined in A and L. The total period of simulation was 10 s, and was divided into segments of .0001 s each. Using a fifth-order Dormand-Prince method solved by the Matlab function ode5, the terms, x and ^ x, were solved at 51 points simultaneously between time segments. Beginning at time t ¼ 0 with initial conditions set as zero and using the fifth-order Dormand-Prince method, we solved for final condition of the first segment, which ends at time t ¼ :0001. This is then the initial condition of the next segment beginning at t ¼ :0001 and ending at t ¼ :0002. This process is continued for the total period of 10 s. Additionally, using a random number generator with a zero-mean Gaussian distribution with variance of 1, a static random sample was added to the observer feedback during each segment to represent the measurement noise. At this point, we have produced an estimate for the input with a significant amount of noise. To smooth the signal and reproduce a cleaner estimate for the input, we introduced a low-pass filter to this result. Due to its simplicity, we chose to
61
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
apply a moving average to the signal to smooth out some of the high-frequency noise in the signal [2]. In continuous time, the moving average can be characterized as:
gðtÞ ¼
1=T : 0:
jtj < T=2 jt j P T=2
ð122Þ
The frequency attributes of the moving average are exhibited taking the Fourier transform of gðtÞ. To determine the frequencies to be emphasized by the low-pass filter, we must find the Fourier transform of the moving average:
GðxÞ ¼
sin
xT 2
xT
:
ð123Þ
2
where GðxÞ is the Fourier transform of gðtÞ. The frequency range most emphasized by the moving average is defined by the area between x ¼ 0 and the first x such that GðxÞ ¼ 0. This is given by:
xT 2
¼ p:
ð124Þ
Thus, the cutoff frequency, xcutoff , can be determined by rearranging Eq. (124) as:
xcutoff ¼
2p : T
ð125Þ
From Eq. (125), we can see that the frequency range to be emphasized is dependent on T. In this example, we choose T ¼ 0:01 s since it would give xcutoff ¼ 628:3 radians per second, which is a wide enough frequency range while reducing the effects of high frequency noise. The results of these simulations will be represented with estimates given by Eq. (111), as presented in Fig. 13, and with a moving average applied to Eq. (111), which is presented in Fig. 14. The moving average is used as a simple example of the potential for post-processing the estimate using noise reduction techniques applied in conjunction with the proposed method. Due to the presence of measurement noise, the values of l1 are larger in magnitude than those used in Section 4. Thus, the top two graphs in Figs. 13 and 14 are unable to be properly scaled, and are on orders of magnitude between 105 and 108 . The influence due to measurement noise is too great, and once we reach the order of magnitude for l1 ¼ 1; 000; 000, we are able to obtain an adequate approximation for the applied force relationship given by Eq. (110). An enlarged view of the case for l1 ¼ 1; 000; 000 is shown in Figs. 11 and 12 for better visibility.
l1= −10000
Applied Force (N)
7
4
x 10
Actual Estimate
2 0 −2 −4 0
1
2
3
4
5
6
7
8
9
10
Time (s)
l1= −100000
Applied Force (N)
8
4
x 10
Actual Estimate
2 0 −2 −4 0
1
2
3
4
5
6
7
8
9
10
Time (s)
Applied Force (N)
l1= −1000000 20 Actual Estimate
10 0 −10 −20 0
1
2
3
4
5
6
7
Time (s) Fig. 11. Single degree-of-freedom system with measurement noise – unfiltered.
8
9
10
62
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
l1= −10000
Applied Force (N)
5
5
x 10
Actual Averaged
0
−5 0
1
2
3
4
5
6
7
8
9
10
Time (s)
l1= −100000
Applied Force (N)
6
4
x 10
Actual Averaged
2 0 −2 −4 0
1
2
3
4
5
6
7
8
9
10
l1= −1000000 20 Actual Averaged
10 0 −10 −20 1
2
3
4
5
6
7
8
Time (s) Fig. 12. Single degree-of-freedom system with measurement noise – moving average.
l1= −1000000 20
Applied Force (N)
0
Actual Estimate
10
0
−10
−20
0
2
4
6
8
10
Time (s) Fig. 13. Single degree-of-freedom system with measurement noise for l1 ¼ 1000000 – unfiltered.
l1= −1000000 15
Actual Averaged
10
Applied Force (N)
Applied Force (N)
Time (s)
5 0 −5 −10 −15 0
2
4
6
8
10
Time (s) Fig. 14. Single degree-of-freedom system with measurement noise for l1 ¼ 1000000 – moving average.
9
10
J. Sanchez, H. Benaroya / Mechanical Systems and Signal Processing 92 (2017) 39–63
63
7. Discussion and conclusion In this paper, we present the development of concepts in application to the problem of force reconstruction. It involves the usage of observer feedback, and in the present analysis, the parameter, given by L, is asymptotic in nature, and allows for the design of the eigenvalues of A þ LC such that approximations in terms of the applied force can be determined. The concepts were developed for an ideal model, and extended to a real model. The concepts developed present a number of benefits. The equilibrium analysis allowed us to develop an alternative representation of the relationship between the applied force and the system response. The solution is generally given by the convolution integral, and this analysis has allowed us to produce a differential series representation of the equilibrium, which proved critical to the development of the asymptotic approximation analysis. Additionally, the asymptotic approximation analysis allows for a simplified analysis of the inverse problem. The concepts developed simplify the analysis by allowing the inverse problem to be analyzed algebraically. This allows for the ease of computation and analysis, and effectively, circumventing the act of inversion in the force reconstruction problem. In application to the force reconstruction problem, the concepts developed herein alleviate the measurement noise issue by redefining the force reconstruction problem in the context of a time series analysis. This allows for the application of postprocessing to reduce the influence of noise and estimate the applied force. Additionally, probabilistic methods can potentially be applied to develop additional methods of modeling the measurement noise in relation to the concepts developed. The current analysis described herein is developed under the assumption of continuous sampling. Moving forward, we aim to relax this assumption, and reconcile the concepts developed with discrete sampling theory. Additionally, we aim to further develop these concepts in relation to probabilistic applied forces, beginning with stationary forces and the generalization of these concepts for non-stationary forces. In these circumstances, the recovery of a time series is not as useful as the recovery of the probabilistic aspects of this series, such as the autocorrelation function or the spectral density. The concepts developed herein can be extended to the estimation of the spectral density of a random force based on the spectral density of the measurements. Acknowledgements We would like to thank the reviewers for their detailed review of this manuscript. Their feedback was incredibly helpful as it allowed us to greatly improve the clarity of the presentation and the depth of the discussion of our work. References [1] Peter Avitabile, Fabio Piergentili, Ken Lown, Identifying dynamic loadings from measured responses, Sound Vib. 33 (1999) 24–28. [2] Christopher Terman, George Verghese, Hari Balakrishnan, Bits, signals, and packets: an introduction to digital communications and networks, From MIT open courseware.