# Online Structural Health Monitoring and Parameter Estimation for Vibrating Active Cantilever Beams Using Low-Priced Microcontrollers

Authors: Gergely Takács; Ján Vachálek; Boris Rohal’-Ilkiv
Journal: Shock and Vibration
 (2015)

---

## Abstract

This paper presents a structural health monitoring and parameter estimation system for vibrating active cantilever beams using low-cost embedded computing hardware. The actuator input and the measured position are used in an augmented nonlinear model to observe the dynamic states and parameters of the beam by the continuous-discrete extended Kalman filter (EKF). The presence of undesirable structural change is detected by variations of the first resonance estimate computed from the observed equivalent mass, stiffness, damping, and voltage-force conversion coefficients. A fault signal is generated upon its departure from a predetermined nominal tolerance band. The algorithm is implemented using automatically generated and deployed machine code on an electronics prototyping platform, featuring an economically feasible 8-bit microcontroller unit (MCU). The validation experiments demonstrate the viability of the proposed system to detect sudden or gradual mechanical changes in real-time, while the functionality on low-cost miniaturized hardware suggests a strong potential for mass-production and structural integration. The modest computing power of the microcontroller and automated code generation designates the proposed system only for very flexible structures, with a first dominant resonant frequency under 4 Hz; however, a code-optimized version certainly allows much stiffer structures or more complicated models on the same hardware.


---

## Body

## 1. Introduction

Recent engineering trends foreshadow a tendency, where mechanical structures are becoming more vibration prone. The design approach, in which lighter materials and more slender shapes are used, may simply have an aesthetic justification, but the savings in material costs are an important factor as well. Aircraft wing surfaces with progressive design may benefit from reduced fuel and maintenance costs or increased maneuverability, while robotic arms and manipulators are faced with increasing demands on working speeds. Often, these improvements are plagued by the curse of lowered stiffness and a susceptibility to demonstrate undesired vibration response.Smart materials and the availability of cheap computing technology enable designers to balance or even completely eliminate the undesirable vibration response of modern structures by active vibration control (AVC) measures [1, 2]. Microcontroller units (MCU) are essentially miniaturized computing platforms that are capable of independent operation and pose the built in peripherals necessary to read sensors or drive actuators. As both the price and the size of MCUs are decreasing [3], it is now possible to integrate vibration control and diagnostics systems into mechanical structures on a mass produced scale, even with the possibility to use several detached systems in a decentralized manner.In addition to the control functionality of structures with AVC, safety-critical applications demand highly integrated online diagnostics functionality [4]. Most engineering structures may experience changes in properties and response dynamics due to unexpected outside influences, such as sensor failures or degradation, material buildup, fatigue, cracks, icing, and a plethora of other factors. In addition to these, process variables such as changing load or working conditions may induce substantial changes in the structure, providing a good reason to use structural health monitoring and diagnostic systems.The requirements for diagnostic systems can range from simple fault indicators to complex parameter estimation procedures. Online parameter estimation allows for a more detailed assessment of the properties of the monitored structures, representing unexpected changes as variations in physically interpretable parameters. Moreover, the estimated parameters can be used to adjust the controllers to counteract structural changes, thereby creating an adaptive vibration control system [5].Several exciting fault and damage detection strategies have appeared for cantilever beams and other structures in the recent literature. Localization of the damaged section is possible by employing progressive optimization-based algorithms, such as genetic algorithms (GA) [6, 7], fuzzy-genetic algorithms [8] particle swarm optimization [9], wavelet transform [10, 11], and artificial neural networks [12]. Crack location and depth were detected by Vakil Baghmisheh et al. [13] using a hybrid particle swarm-Nelder-Mead algorithm and modal testing. Although most of these methods offer much more than a simple status signal on the presence or absence of a fault in the structure, they require a large computational effort, suggest but do not demonstrate online operation, or are simply not suitable for online application due to design and technical limitations [7–13].This paper introduces a prototype of an online diagnostic system intended to detect and monitor changes in the mechanical properties of actively controlled structures resembling the dynamic characteristics of thin cantilever beams. In addition to monitoring the occurrence of changes caused by unexpected outside influences, the system based on a nonlinear joint state and parameter observer is capable of providing real-time estimates that can be used to create adaptive mechanical structures. The proposed estimation scheme is implemented on a low-cost electronic prototyping platform using automatized code generation tools.The approach presented here assumes that the dynamics of active mechanical structures resembling the properties thin cantilever beams can be sufficiently approximated using a single degree of freedom (SDOF) system [14]. This assumption holds well in the case of flexible structures with single digit first—fundamental—resonant frequencies, dominating the overall response in comparison with the contribution of the higher modes. The state-space model of the driven mass-spring-damper is augmented by the monitored physical parameters, namely, the unknown equivalent stiffness, damping, mass, and force conversion coefficients.The introduction of unknown parameters in the linear SDOF dynamic system renders the augmented model nonlinear. This paper uses the nonlinear expansion of the acknowledged Kalman filter [15–17]—the well-known extended Kalman filter (EKF) [18, 19]—to observe the states and parameters of this augmented system. Though it is possible to directly compute and output the equivalent physical parameters for the needs of adaptive control or detailed monitoring, here the parameters estimated by the EKF are used to assess the first damped resonant frequency of the active structure. It is assumed that this remains bounded if there is no change in the structure but will diverge from the nominal levels as soon as a mechanical failure or an undesired change of physical configuration occurs [20].The proposed structural monitoring system is tested in a laboratory setting. The test structure uses the classical example of the clamped cantilever (see Figure1) that is often employed to represent the general dynamic behavior of real-life structures resembling flexible thin beams [14], such as fixed and rotary aerodynamic surfaces [21, 22], antennae [23], manipulators [24], and others [25]. The thin aluminum cantilever is equipped with two large double-layer piezoceramic transducers, mounted on opposing sides close to the clamping mechanism. The transducers are driven by a power amplifier, receiving a pseudorandom-binary (PRBS) signal. The motion of the beam is measured at a single location near the end of the structure. The actuator inputs and the measurements are fed to a low-cost microcontroller and electronics prototyping platform featuring a 8-bit MCU (see Figure 2).Figure 1
Simplified hardware connection scheme of the experimental system.Figure 2
Embedded microcontroller unit and the electronics prototyping board.The experimental validation procedure of the proposed parameter monitoring and fault detection scheme involves supplying the beam with the PRBS signal, emulating the control input to the actuator. In the first stage of the experiments, the beam is in its nominal configuration without any changes and the parameters converge to nominal levels. At the beginning of the second phase of the validation experiments, different weights are added to the end of the beam in order to emulate the occurrence of a mechanical fault. The algorithm running on the MCU estimates the parameters and calculates the first resonant frequency of the beam in real-time, based on the measured actuator input and position. As soon as the frequency leaves its user determined tolerance band, the fault is detected and a warning signal is sent. The experiments featured in this paper demonstrate the estimates of the dynamic states and parameters for the different levels of mechanical faults and indicate the time when the error is detected.The use of the extended Kalman filter for the estimation of vibration dynamics is, with a few exceptions, limited to offline a posteriori system identification procedures. Real-time application of EKF for diagnostics and monitoring is presented by Lourens et al. [26] for disturbance force estimation and by Mu et al. [27] for the monitoring of structures under seismic loads, stiffness, and damping coefficient estimation in automotive applications [28] and for the adaptive control of torsional vibration in a two-mass rotating drive system [29, 30]. While Williams suggested the use of unscented EKF to assess damage and damage location based on accelerometer measurements [31], the work is limited to offline numerical studies. Stress and strain fields are observed experimentally using a Kalman filter-based methods by Erazo and Hernandez [32], but the work is again limited to an a posteriori test. The EKF was used in previous works to provide state and parameter estimates for an adaptive vibration control system for thin beams [5, 33, 34]; however, this was implemented on an expensive rapid control prototyping platform offering extensive computational power. The main contribution of this paper is therefore not only the proposition to use the EKF for the estimation of vibration dynamics, but also the online and real-time implementation and validation of an EKF-based scheme on a low-cost platform with a potential for mass-production and close structural integration.Having stated the contributions of the paper, it is equally important to realize what is not the ambition of this work. The proposed monitoring system is capable of providing real-time estimates of the parameters; however, this paper does not attempt to solve the question of the usage of this online information. The observed parameters may be employed to generate simple warning signals, drive fail-safe mechanisms, or serve as a basis for adaptive vibration control systems. Likewise, the paper does not attempt to make recommendations on the mechanical integration of the system into structures, since it implements the algorithm on a prototyping platform; nevertheless, the miniaturized MCU used on the board suggests that this is a likely possibility.It is also important to note that structural failure detection in this paper is only solved on the level of its occurrence. Localization or quantification of the structural failure is not within the scope of this work, and the proposed method will only provide information on the presence or absence of a certain change.


## 2. Materials and Methods

### 2.1. Theoretical Foundations

#### 2.1.1. Dynamic Model of the Beam

Due to the flexible nature of the structure, the dynamics of the actively controlled thin beam is dominated by its first resonant frequency [14, 35]. Even though from the perspective of mechanics the SDOF assumption may seem simplistic, it is sufficient for our purposes and is often used in combination with EKF to estimate vibration motion [36–38].Naturally, a first order model and the resulting EKF may cover only structures with a single and significantly dominant resonant mode. Nonetheless, a better representation of a structure yielding models that are still feasible for computational implementation can be obtained from larger models by reduced order modelling techniques or, for example, by the use of the proper order decomposition method. A model covering approximately 3 modes of vibration and augmented by unknown parameters is plausible for EKF and its implementation on a low-cost MCU, however, beyond that one may encounter difficulties with algorithm tuning and computational feasibility.Let us assume now that, from the viewpoint of the proposed fault detection algorithm, the cantilever beam can be approximated by the second order differential equation describing the motion of a mass-spring-damper system with an external driving force [39] (1)mq¨t+bq˙t+kqt=Ft,where the driving actuator force (2)Ft=cutis a linear multiple of the input voltage to the actuators u(t) (V), while c (N/V) is the voltage-force conversion coefficient. The mass of the structure m (kg) is a sum of the mass of the beam and the actuators with an additional weight attached at its end. The damping coefficient is given by b (Ns/m), while k (N/m) is the stiffness coefficient. The true position of the beam at the free end is denoted by q(t) (m), while its velocity is q˙(t) (ms−1) and acceleration is q¨(t) (ms−2).The second order differential equation in (1) may be deconstructed into two first order differential equations to create a state-space representation of the beam dynamics [40]. The choice of position x1(t)=q(t) and velocity x2(t)=q˙(t) as dynamic states in the representation will yield the familiar continuous state-space equations for the mass-spring-damper [41] (3)x˙t=Fxt+Gutzt=Hxt,where z(t) (m) is the position measured at a single location close to the end of the active beam. The position output can be obtained from the state vector x(t)=x1(t)x2(t)T by using H=10 as the R2×1 output matrix. Furthermore, F∈R2×2 is the dynamic matrix given by (4)F=01-km-bm,and G∈R2×1 is the input matrix given by (5)G=0cm.


#### 2.1.2. Augmented Dynamic Model

The dynamic states of the structure—position and velocity as measured at the free end—are contained in the state variablex(t). This state variable can be augmented by a vector, incorporating the unknown and potentially changing parameters to obtain (6)xat=xtptT.The system dynamics given by (3) is clearly linear; however, by including the unknown parameters in the augmented state formulation this new estimation problem will be nonlinear and will require methods designed for nonlinear systems [42–44]. In the case of the active beam, none of the parameters are known beforehand and all of them are likely to vary during the operation of the structure, so the parameter vector will be (7)pt=ktbtct1mtT.The process noise and the measurement noise are assumed to be sequentially uncorrelated and have a Gaussian distribution with zero mean. This white noise disturbance assumption comes from the theoretical formulation of the EKF [44] but is also a common premise in operational modal analysis (OMA). If the process noise for the motion dynamics components is denoted by ws(t) and the process noise for the parameters is wp(t), the nonlinear continuous function expressing the dynamics of the augmented model with noise will be (8)x˙t=f~cxt,pt,ut+wstp˙t=0+wptfor the motion and parameter dynamics which are assumed to remain constant between the individual samples. Combining the motion dynamics and parameter noise into a single term w(t)=ws(t)wp(t)T enables us to express the augmented dynamics compactly by (9)x˙at=fcxat,ut+wtzt=hcxat,ut+vt,where the continuous function fc defines the augmented nonlinear system dynamics, hc is a continuous measurement function, and v(t) is the measurement noise. The properties of the process noise are defined using the covariance matrix [42] (10)wt~N0,Qt,where the discretized analogy of the process noise is obtained by the relation (11)Q=QtTs,where Ts is the measurement sampling period. Similar to the process noise covariance, the properties of the measurement noise are given by [42] (12)vt~N0,R.This work assumes the covariance matricesQt and R to remain constant during the operation of the fault detection system. This implies that the statistical characteristics of the process noise of the dynamic states and parameters and the measurement noise are assumed to remain unchanged. While it is possible to conceive a modification to the EKF with robust, rescalable and possibly adaptive process, and measurement noise covariance matrices [45–47], this is beyond the scope of this paper.Modern control systems operate on a discrete basis. The nonlinear estimation procedure introduced in the following system will be also implemented on an embedded microcontroller unit, performing measurements at discrete intervalst=T(k), k=1,2,3,(…). The continuous augmented dynamics in (9) will be propagated using the Euler’s method by the system (13)xak+1=fdxak,uk+wkzk=hdxak,uk+vk.


#### 2.1.3. Joint State and Parameter Estimation

The extended Kalman filter is proposed in this work to estimate the nonlinear augmented states and parameters of the active structure [43, 44]. The algorithm considered here is often referred to as the continuous-discrete or hybrid formulation in literature [42–44], since the EKF takes its measurements and yields new state estimates at discrete sample times, while the system model is evaluated by numerical integration. The observer algorithm is initialized by a first state estimate at zero time (14)x^a0+=Exa0that is determined by the user as an initial guess of the states and parameters. Similarly, the initial covariance matrix of the error of this state estimate is given by the user as (15)P0+=Exa0-x^a0+xa0-x^a0+T.The primary estimate—ora priori estimate, as it is often denoted in literature—of the augmented state is obtained using simulation. The dynamics of the driven point-mass-damper with the unknown parameter vector is described by the continuous augmented function f(x^a(k-1)+,u(k)). Moreover, another assumption is that the parameters p(t)=k(t)b(t)c(t)1/m(t)T remain constant between the samples. The continuous state is then propagated using (16)x^ak-=fcx^ak-1+,uk,between samples t=(k-1)Ts to t=(k)Ts to obtain a priori estimates, denoted here by the - subscript. Though several methods exist to solve ordinary differential equations numerically, this work utilizes the first order Euler method. The continuous augmented function fc(x^a(k-1)+,u(k)) can be safely considered as (17)q˙t=x˙1t=x2tq¨t=x˙2t=-x1tx3tx6t-x2tx4tx6t+x5tx6tutk˙t=x˙3t=0b˙t=x˙4t=0c˙t=x˙5t=01m˙t=x˙6t=0for flexible active structures with dominant first resonant frequencies.In addition to the computation of the a priori state estimate, it is also necessary to evaluate the a priori update of the covariance matrix of the state estimates(18)P˙t=Zx^atPt+PtZTx^at+Qt,where Z(x^a(t)) is the Jacobian of (17) with respect to the augmented state (19)Zx^at=∂fcx^at∂x^atxat=x^at,which for the formerly described augmented beam model it is given by (20)Zx^at=010000-x6tx3t-x6tx4t-x6tx1t-x6tx2t-x6tut-x1tx3t--x2tx4t--utx5t000000000000000000000000.Similar to the numerical propagation of the augmented state, the covariance matrix is obtained by simulating (18) from sample t=(k-1)Ts to t=(k)Ts.The discrete part of the EKF uses the current measurement to update the a priori state to get a reviseda posteriori state estimate, which is denoted by the + subscript. The first step is to compute the Kalman gain K(k) using (21)Kk=Pk-LTLPk-LT+MRMT-1.The a priori state estimate is then revised by the measurement at time(k) with its relative importance weighted to the simulation estimate by the Kalman gain K(k), which is computed at the next step. The a posteriori update of the augmented state is (22)x^ak+=x^ak-+Kkzk-hx^ak-,uk.The state error covariance matrix is updated by(23)Pk+=I-KkLPk-I-KkLT+KkMRMTKkT,where the expression L is the Jacobian of the measurement equation hc(xa(t)) with respect to the state xa(t)(24)L=∂hcxat∂xatxat=x^at-=00000x6t000001000001000001000001,and matrix M is the Jacobian of hc(xa(t)) with respect to the measurement noise term va(t)(25)M=∂hcxat∂vatxat=x^at-=1.The a posteriori augmented state estimate x^a(k)+ and covariance matrix estimate P(k)+ are passed onto the next computation step, where they become the basis for the next a priori estimates.


#### 2.1.4. Resonance Estimates

The first dominant resonant frequency of the beam is approximated on the basis of the estimated parametersp^(k) outside the observer algorithm. Given the parameter estimates k^(k),b^(k),m^(k) at t=kTs and assuming that the damped natural frequency of the dynamically analogous SDOF system in (1) is equal to the dominant resonance frequency of the structure [14], the frequency is (26)f^k=4k^km^k-b^k24πm^k.The algorithm monitors this frequency estimate throughout the operation of the system. As it has been noted by Jassim et al. [20] and others, monitoring the change of natural frequency is a feasible tool to indicate damage occurrence. In case that the estimate is outside the tolerance band (27)ek=0iffn-ft≤f^k≥fn+ft1iffn-ft≥f^k1iffn+ft≤f^k,a fault detection signal e(k) is produced and retained unless cleared. The resonance of the nominal structure is denoted by fn and ft is a symmetric tolerance band.


#### 2.1.5. Algorithm Summary

The resulting real-time fault detection strategy with joint state and parameter estimation can be summarized in the following algorithm.Algorithm 1.
Perform the following operations at each sampling instant(k).(1)
Propagate (16) and (18) by numerical simulation using the Euler method. Obtain the a priori estimates of the augmented state x^a(k)- and covariance matrix P(k)-.(2)
Sample the actual deflectionz(k) at the end of the cantilever beam.(3)
Compute the Kalman gainKf(k) using (21).(4)
In order to compute the a posteriori state estimatex^a(k)+, use the new measurement sample z(k) and Kalman gain update Kf(k) to revise the a priori state estimate through (22).(5)
Compute the a posteriori update of the covariance matrixP(k)+ in (23).(6)
Retain the a posteriori state and covariance estimatesx^a(k)+ and P(k)+ in the memory for the next estimation step.(7)
Approximate the damped natural frequencyf^(k) of the beam according to (26).(8)
If the frequency estimatef^(k) is outside the tolerance band in (27), enable the fault detection signal e(k).(9)
At the next sample restart the procedure from 1, use the estimates from Step 6.


### 2.2. Experimental Hardware and Software

#### 2.2.1. Experimental Hardware

The proposed fault detection system was tested on a thin flexible cantilever beam with piezoceramic actuation. The test hardware composition will be described only briefly as the test bed was previously introduced in detail in previous work, including [5, 35, 48]. The schematic setup of the experimental hardware, including the active thin structure and the prototyping system, is illustrated in Figure 1.The cantilever is a 550 × 40 × 3 mm beam, fixed at one end and freely moving at the other. As it is indicated in the side view of the beam in Figure1, it undergoes flexural vibration in the direction of its thickness. The beam is made of commercially pure aluminium, having mechanical properties which render the structure more flexible than it would be possible with regularly available aluminium alloys.To further decrease its first resonant frequency, a heavy nut, bolt, and a set of four washers are mounted at the free end of the beam. The nominal weight installed at the end is necessary to shift its original resonance (8.17 Hz) to the area, where is possible to run the proposed fault detection system in real-time on the hardware with limited computing power and automatically deployed code. The weight of the aluminum cantilever beam wasmb=0.182 kg, and the total weight of the nut-washer-bolt assembly is mw=0.305 kg, making the total physical weight of m=0.487 kg. This base configuration assumed here therefore consists of the beam with a weight, with a first nominal resonance decreased to 2.26 Hz. This nominal configuration is taken as the initial configuration in each experiment, which is then further modified by the experiments.An industrial laser triangulation device is mounted near the free end of the beam, to provide position measurements for the fault detection algorithm. The Keyence LK-G 32 laser sensor head is pointed at a single location at the beam end, while its signal is preprocessed with a low pass filter in the Keyence LK-G3001V processing unit. The Keyence LK-G3001V unit provides an analog signal scaling the ±5 mm maximal measurement range into a 0–5 V signal, compatible with the analog measurement range of the prototyping board. The laser processing unit is directly connected to the AI1 analog input port of MCU. The measurements are rescaled and shifted to readings in meters within the initial portion of the estimation algorithm.The actuation of the active thin structure is realized with a set of two double layer piezoceramic transducers. The MIDÉ QuickPack QP45n double layer actuators are mounted on both fixed sides of the beam and receive the same electric signal with opposing polarity. The high voltage driving signal is provided through a MIDÉ 1225 capacitive precision amplifier with a 20× V/V gain on the signal input. The active control signal is emulated by a pseudorandom binary voltage signal, having its low level at 0 V and high level at 4 V, resulting in 80 V impulses to the transducers. The PRBS signal has a 10 bit/s base sample (in order to excite the first resonance of the beam) and is generated by an Agilent 3351A signal generator unit. In addition to the measured position inputz(k), the extended Kalman filter requires the input signal u(k) as well; therefore the output from the signal generator is directly connected to the AI1 analog input port of the prototyping board and then rescaled to the amplified levels within the EKF algorithm.The low-cost electronics prototyping board employed in this work is an Arduino Mega 2560 Revision 3 device, featuring a simple 8-bit Atmel ATmega 2560 microcontroller unit (Figure2). The MCU used in this board has a clock speed of 16 MHz and a 256 kB large flash memory for storing user algorithms [49], enabling its online external use with the Simulink product range. The detected fault signal e(k) is passed on to the DO9 digital output of the board with an LED connected to it for signal indication. Although this is only a binary information relating to the occurrence of a mechanical fault, the parameter estimates may also be output using a pulse-width modulated signal or encoded into several digital outputs. The board is powered through the USB port, which also serves as the means for programming and online data transfer for logging and storage. Programming of the board and data logging is realized using a laptop computer with an Intel Core i3 CPU running at 2.4 GHz and 4 GB operating memory. The computing power of the design computer has no effect on the speed of the proposed algorithm, as it is running in a stand alone mode on the Arduino Mega 2560 prototyping board.


#### 2.2.2. Experimental Software

The proposed estimation and fault detection algorithm was developed within the Mathworks Matlab and Simulink framework. The software was implemented in Matlab and Simulink 8.0 (R2012b), running under the Microsoft Windows 7 × 64 operating system. In addition to the programming role of this software implementation setup, the prototyping board was also providing the estimates and other signals through the USB connection in external mode. There was an approximately 4 s communication delay between the prototyping board and the Simulink data logging and visualizing interface. Even though this portion of the experiments is not shown on figures here, the proposed algorithm was still fully functional on the MCU during this initialization period.The simplified Simulink scheme of the software implementation is featured in Figure3. Communication between the analog input and digital output ports of the Arduino Mega 2560 board and the real-time logging capability was ensured using the Simulink Target for use with Arduino Hardware (v 2.0) suite, providing access to the peripherals of the MCU and a possibility to deploy Simulink schemes onto the hardware. The 10-bit in-processor A/D converter [49] provides a decimal integer reading of the voltages, which is converted using type conversion blocks. The fault signal is transferred into the digital port using the corresponding output block. Internal state estimates are displayed and saved using a scope block, with its logging capability enabled.Figure 3
Block scheme of the algorithm running on the microcontroller unit.In addition to the readily available standardized Simulink blocks, the core of the algorithm featured previously was developed in the Matlab M-script language and then automatically transcribed into a C/C++ version for the embedded MCU. The custom algorithm block noted as EKF in Figure3 contains means to manipulate voltage inputs to physical units, the extended Kalman filter providing state and parameter estimates, and the simple algorithm generating the fault signal based on the calculated damped first resonant frequency. This M-script code is automatically recompiled into embedded C/C++ code and downloaded to the target hardware using the Embedded Coder (v 6.3) and Simulink Coder (v 8.3). The compilation itself is performed by the Microsoft Visual C++ 2010 compiler.


### 2.3. Experiment Design and Settings

#### 2.3.1. Experiment Design and Procedures

The experiments featured in this paper were designed to demonstrate the functionality of the proposed estimation and fault detection algorithm and the viability of its real-time implementation on low-cost microcontroller hardware. Experiments were lasting 80 s with the pseudorandom binary driving signal supplied to the actuators during the entire duration of the test. After initializing the experiments with completely identical tuning and settings, the estimated dynamic states and the parameters converged to the values in agreement with the mechanical properties of the nominal beam. The majority of this starting convergence is not indicated in this paper, due to the previously mentioned communication issue between the MCU and the laptop computer, which is performing data logging and visualization.In addition to a nominal beam configuration with the beam itself and the heavy bolt-washer-nut combination, 6 types of different weights ranging from 16 g to 345 g were added to the end of the beam approximately at the ~40 s time mark (see Figure4). The added weights emulated the effects of mechanical failures and changes with different severity. As these weights were placed manually, there may be a ±1-2 s variance in this time between the individual tests. The weights and individual experiment types are summarized in Table 1. The added weights changed the overall dynamic behavior of the beam and yielded parameters which diverged from the nominal values. If the approximate of the first damped resonance f^(k) left the ft=±0.1 Hz tolerance band of the nominal first resonance at fn=2.27 Hz, the alarm signal e(k) was switched on and this fact was logged along with the rest of the parameters. The first resonant frequency of the thin structure in its nominal configuration was determined by supplying the piezoceramic transducers with a PRBS signal and performing a position measurement. This time domain measurement was then analysed using the fast Fourier transform (FFT) to extract the first fundamental resonant frequency.Table 1
Summary of experimental scenarios.
Number
Weight
Change
Nominal
Color
1
0 g
0%
✓
Blue
2
16 g
3%
✗
Red
3
50 g
10%
✗
Green
4
99 g
20%
✗
Black
5
195 g
40%
✗
Magenta
6
246 g
50%
✗
Orange
7
345 g
70%
✗
GreyFigure 4
Estimated dynamic states for the experiment with the 246 g added weight.


#### 2.3.2. Experimental Parameters and Settings

All experiments were initialized with the initial state estimate of(28)x^a0+=0010011E-410.5T,denoted in the experiments with an upside down black triangle located at zero time. Likewise, the initial covariance estimates remained the same throughout the experiments with (29)P0+=1E-50000001E-40000001E+50000001E-20000001E-50000001E+0.The values in the P0+ matrix provide a compromise between a realistic guess of the variance in the initial state estimate error and an incremental tuning of the EKF algorithm.The covariance of the single measurement signal was estimated based on the square of the maximal amplitude resolution of the A/D input of the prototyping board in meters, yielding ~2.5 V at 10 bits [49] or in the terms of the measurement covariance (30)R=R=1E-5.The greatest freedom in tuning the behavior of the algorithm can be achieved by manipulating the continuous-time process noise covariance matrix. The process noise matrix was tuned for best performance based on incremental testing and the expected range of changes in the physical variables and then fixed for all experiments featured in this paper as (31)Qt=5E-3000001E+1000001E-1000001E-12000003E-1,with its discrete-time analogy recomputed by Q=Qt/Ts.Global sampling times were set asTs = 0.1 s, while the continuous model was numerically integrated using the Euler method with a dt=0.02 s time step.


## 2.1. Theoretical Foundations

### 2.1.1. Dynamic Model of the Beam

Due to the flexible nature of the structure, the dynamics of the actively controlled thin beam is dominated by its first resonant frequency [14, 35]. Even though from the perspective of mechanics the SDOF assumption may seem simplistic, it is sufficient for our purposes and is often used in combination with EKF to estimate vibration motion [36–38].Naturally, a first order model and the resulting EKF may cover only structures with a single and significantly dominant resonant mode. Nonetheless, a better representation of a structure yielding models that are still feasible for computational implementation can be obtained from larger models by reduced order modelling techniques or, for example, by the use of the proper order decomposition method. A model covering approximately 3 modes of vibration and augmented by unknown parameters is plausible for EKF and its implementation on a low-cost MCU, however, beyond that one may encounter difficulties with algorithm tuning and computational feasibility.Let us assume now that, from the viewpoint of the proposed fault detection algorithm, the cantilever beam can be approximated by the second order differential equation describing the motion of a mass-spring-damper system with an external driving force [39] (1)mq¨t+bq˙t+kqt=Ft,where the driving actuator force (2)Ft=cutis a linear multiple of the input voltage to the actuators u(t) (V), while c (N/V) is the voltage-force conversion coefficient. The mass of the structure m (kg) is a sum of the mass of the beam and the actuators with an additional weight attached at its end. The damping coefficient is given by b (Ns/m), while k (N/m) is the stiffness coefficient. The true position of the beam at the free end is denoted by q(t) (m), while its velocity is q˙(t) (ms−1) and acceleration is q¨(t) (ms−2).The second order differential equation in (1) may be deconstructed into two first order differential equations to create a state-space representation of the beam dynamics [40]. The choice of position x1(t)=q(t) and velocity x2(t)=q˙(t) as dynamic states in the representation will yield the familiar continuous state-space equations for the mass-spring-damper [41] (3)x˙t=Fxt+Gutzt=Hxt,where z(t) (m) is the position measured at a single location close to the end of the active beam. The position output can be obtained from the state vector x(t)=x1(t)x2(t)T by using H=10 as the R2×1 output matrix. Furthermore, F∈R2×2 is the dynamic matrix given by (4)F=01-km-bm,and G∈R2×1 is the input matrix given by (5)G=0cm.


### 2.1.2. Augmented Dynamic Model

The dynamic states of the structure—position and velocity as measured at the free end—are contained in the state variablex(t). This state variable can be augmented by a vector, incorporating the unknown and potentially changing parameters to obtain (6)xat=xtptT.The system dynamics given by (3) is clearly linear; however, by including the unknown parameters in the augmented state formulation this new estimation problem will be nonlinear and will require methods designed for nonlinear systems [42–44]. In the case of the active beam, none of the parameters are known beforehand and all of them are likely to vary during the operation of the structure, so the parameter vector will be (7)pt=ktbtct1mtT.The process noise and the measurement noise are assumed to be sequentially uncorrelated and have a Gaussian distribution with zero mean. This white noise disturbance assumption comes from the theoretical formulation of the EKF [44] but is also a common premise in operational modal analysis (OMA). If the process noise for the motion dynamics components is denoted by ws(t) and the process noise for the parameters is wp(t), the nonlinear continuous function expressing the dynamics of the augmented model with noise will be (8)x˙t=f~cxt,pt,ut+wstp˙t=0+wptfor the motion and parameter dynamics which are assumed to remain constant between the individual samples. Combining the motion dynamics and parameter noise into a single term w(t)=ws(t)wp(t)T enables us to express the augmented dynamics compactly by (9)x˙at=fcxat,ut+wtzt=hcxat,ut+vt,where the continuous function fc defines the augmented nonlinear system dynamics, hc is a continuous measurement function, and v(t) is the measurement noise. The properties of the process noise are defined using the covariance matrix [42] (10)wt~N0,Qt,where the discretized analogy of the process noise is obtained by the relation (11)Q=QtTs,where Ts is the measurement sampling period. Similar to the process noise covariance, the properties of the measurement noise are given by [42] (12)vt~N0,R.This work assumes the covariance matricesQt and R to remain constant during the operation of the fault detection system. This implies that the statistical characteristics of the process noise of the dynamic states and parameters and the measurement noise are assumed to remain unchanged. While it is possible to conceive a modification to the EKF with robust, rescalable and possibly adaptive process, and measurement noise covariance matrices [45–47], this is beyond the scope of this paper.Modern control systems operate on a discrete basis. The nonlinear estimation procedure introduced in the following system will be also implemented on an embedded microcontroller unit, performing measurements at discrete intervalst=T(k), k=1,2,3,(…). The continuous augmented dynamics in (9) will be propagated using the Euler’s method by the system (13)xak+1=fdxak,uk+wkzk=hdxak,uk+vk.


### 2.1.3. Joint State and Parameter Estimation

The extended Kalman filter is proposed in this work to estimate the nonlinear augmented states and parameters of the active structure [43, 44]. The algorithm considered here is often referred to as the continuous-discrete or hybrid formulation in literature [42–44], since the EKF takes its measurements and yields new state estimates at discrete sample times, while the system model is evaluated by numerical integration. The observer algorithm is initialized by a first state estimate at zero time (14)x^a0+=Exa0that is determined by the user as an initial guess of the states and parameters. Similarly, the initial covariance matrix of the error of this state estimate is given by the user as (15)P0+=Exa0-x^a0+xa0-x^a0+T.The primary estimate—ora priori estimate, as it is often denoted in literature—of the augmented state is obtained using simulation. The dynamics of the driven point-mass-damper with the unknown parameter vector is described by the continuous augmented function f(x^a(k-1)+,u(k)). Moreover, another assumption is that the parameters p(t)=k(t)b(t)c(t)1/m(t)T remain constant between the samples. The continuous state is then propagated using (16)x^ak-=fcx^ak-1+,uk,between samples t=(k-1)Ts to t=(k)Ts to obtain a priori estimates, denoted here by the - subscript. Though several methods exist to solve ordinary differential equations numerically, this work utilizes the first order Euler method. The continuous augmented function fc(x^a(k-1)+,u(k)) can be safely considered as (17)q˙t=x˙1t=x2tq¨t=x˙2t=-x1tx3tx6t-x2tx4tx6t+x5tx6tutk˙t=x˙3t=0b˙t=x˙4t=0c˙t=x˙5t=01m˙t=x˙6t=0for flexible active structures with dominant first resonant frequencies.In addition to the computation of the a priori state estimate, it is also necessary to evaluate the a priori update of the covariance matrix of the state estimates(18)P˙t=Zx^atPt+PtZTx^at+Qt,where Z(x^a(t)) is the Jacobian of (17) with respect to the augmented state (19)Zx^at=∂fcx^at∂x^atxat=x^at,which for the formerly described augmented beam model it is given by (20)Zx^at=010000-x6tx3t-x6tx4t-x6tx1t-x6tx2t-x6tut-x1tx3t--x2tx4t--utx5t000000000000000000000000.Similar to the numerical propagation of the augmented state, the covariance matrix is obtained by simulating (18) from sample t=(k-1)Ts to t=(k)Ts.The discrete part of the EKF uses the current measurement to update the a priori state to get a reviseda posteriori state estimate, which is denoted by the + subscript. The first step is to compute the Kalman gain K(k) using (21)Kk=Pk-LTLPk-LT+MRMT-1.The a priori state estimate is then revised by the measurement at time(k) with its relative importance weighted to the simulation estimate by the Kalman gain K(k), which is computed at the next step. The a posteriori update of the augmented state is (22)x^ak+=x^ak-+Kkzk-hx^ak-,uk.The state error covariance matrix is updated by(23)Pk+=I-KkLPk-I-KkLT+KkMRMTKkT,where the expression L is the Jacobian of the measurement equation hc(xa(t)) with respect to the state xa(t)(24)L=∂hcxat∂xatxat=x^at-=00000x6t000001000001000001000001,and matrix M is the Jacobian of hc(xa(t)) with respect to the measurement noise term va(t)(25)M=∂hcxat∂vatxat=x^at-=1.The a posteriori augmented state estimate x^a(k)+ and covariance matrix estimate P(k)+ are passed onto the next computation step, where they become the basis for the next a priori estimates.


### 2.1.4. Resonance Estimates

The first dominant resonant frequency of the beam is approximated on the basis of the estimated parametersp^(k) outside the observer algorithm. Given the parameter estimates k^(k),b^(k),m^(k) at t=kTs and assuming that the damped natural frequency of the dynamically analogous SDOF system in (1) is equal to the dominant resonance frequency of the structure [14], the frequency is (26)f^k=4k^km^k-b^k24πm^k.The algorithm monitors this frequency estimate throughout the operation of the system. As it has been noted by Jassim et al. [20] and others, monitoring the change of natural frequency is a feasible tool to indicate damage occurrence. In case that the estimate is outside the tolerance band (27)ek=0iffn-ft≤f^k≥fn+ft1iffn-ft≥f^k1iffn+ft≤f^k,a fault detection signal e(k) is produced and retained unless cleared. The resonance of the nominal structure is denoted by fn and ft is a symmetric tolerance band.


### 2.1.5. Algorithm Summary

The resulting real-time fault detection strategy with joint state and parameter estimation can be summarized in the following algorithm.Algorithm 1.
Perform the following operations at each sampling instant(k).(1)
Propagate (16) and (18) by numerical simulation using the Euler method. Obtain the a priori estimates of the augmented state x^a(k)- and covariance matrix P(k)-.(2)
Sample the actual deflectionz(k) at the end of the cantilever beam.(3)
Compute the Kalman gainKf(k) using (21).(4)
In order to compute the a posteriori state estimatex^a(k)+, use the new measurement sample z(k) and Kalman gain update Kf(k) to revise the a priori state estimate through (22).(5)
Compute the a posteriori update of the covariance matrixP(k)+ in (23).(6)
Retain the a posteriori state and covariance estimatesx^a(k)+ and P(k)+ in the memory for the next estimation step.(7)
Approximate the damped natural frequencyf^(k) of the beam according to (26).(8)
If the frequency estimatef^(k) is outside the tolerance band in (27), enable the fault detection signal e(k).(9)
At the next sample restart the procedure from 1, use the estimates from Step 6.


## 2.1.1. Dynamic Model of the Beam

Due to the flexible nature of the structure, the dynamics of the actively controlled thin beam is dominated by its first resonant frequency [14, 35]. Even though from the perspective of mechanics the SDOF assumption may seem simplistic, it is sufficient for our purposes and is often used in combination with EKF to estimate vibration motion [36–38].Naturally, a first order model and the resulting EKF may cover only structures with a single and significantly dominant resonant mode. Nonetheless, a better representation of a structure yielding models that are still feasible for computational implementation can be obtained from larger models by reduced order modelling techniques or, for example, by the use of the proper order decomposition method. A model covering approximately 3 modes of vibration and augmented by unknown parameters is plausible for EKF and its implementation on a low-cost MCU, however, beyond that one may encounter difficulties with algorithm tuning and computational feasibility.Let us assume now that, from the viewpoint of the proposed fault detection algorithm, the cantilever beam can be approximated by the second order differential equation describing the motion of a mass-spring-damper system with an external driving force [39] (1)mq¨t+bq˙t+kqt=Ft,where the driving actuator force (2)Ft=cutis a linear multiple of the input voltage to the actuators u(t) (V), while c (N/V) is the voltage-force conversion coefficient. The mass of the structure m (kg) is a sum of the mass of the beam and the actuators with an additional weight attached at its end. The damping coefficient is given by b (Ns/m), while k (N/m) is the stiffness coefficient. The true position of the beam at the free end is denoted by q(t) (m), while its velocity is q˙(t) (ms−1) and acceleration is q¨(t) (ms−2).The second order differential equation in (1) may be deconstructed into two first order differential equations to create a state-space representation of the beam dynamics [40]. The choice of position x1(t)=q(t) and velocity x2(t)=q˙(t) as dynamic states in the representation will yield the familiar continuous state-space equations for the mass-spring-damper [41] (3)x˙t=Fxt+Gutzt=Hxt,where z(t) (m) is the position measured at a single location close to the end of the active beam. The position output can be obtained from the state vector x(t)=x1(t)x2(t)T by using H=10 as the R2×1 output matrix. Furthermore, F∈R2×2 is the dynamic matrix given by (4)F=01-km-bm,and G∈R2×1 is the input matrix given by (5)G=0cm.


## 2.1.2. Augmented Dynamic Model

The dynamic states of the structure—position and velocity as measured at the free end—are contained in the state variablex(t). This state variable can be augmented by a vector, incorporating the unknown and potentially changing parameters to obtain (6)xat=xtptT.The system dynamics given by (3) is clearly linear; however, by including the unknown parameters in the augmented state formulation this new estimation problem will be nonlinear and will require methods designed for nonlinear systems [42–44]. In the case of the active beam, none of the parameters are known beforehand and all of them are likely to vary during the operation of the structure, so the parameter vector will be (7)pt=ktbtct1mtT.The process noise and the measurement noise are assumed to be sequentially uncorrelated and have a Gaussian distribution with zero mean. This white noise disturbance assumption comes from the theoretical formulation of the EKF [44] but is also a common premise in operational modal analysis (OMA). If the process noise for the motion dynamics components is denoted by ws(t) and the process noise for the parameters is wp(t), the nonlinear continuous function expressing the dynamics of the augmented model with noise will be (8)x˙t=f~cxt,pt,ut+wstp˙t=0+wptfor the motion and parameter dynamics which are assumed to remain constant between the individual samples. Combining the motion dynamics and parameter noise into a single term w(t)=ws(t)wp(t)T enables us to express the augmented dynamics compactly by (9)x˙at=fcxat,ut+wtzt=hcxat,ut+vt,where the continuous function fc defines the augmented nonlinear system dynamics, hc is a continuous measurement function, and v(t) is the measurement noise. The properties of the process noise are defined using the covariance matrix [42] (10)wt~N0,Qt,where the discretized analogy of the process noise is obtained by the relation (11)Q=QtTs,where Ts is the measurement sampling period. Similar to the process noise covariance, the properties of the measurement noise are given by [42] (12)vt~N0,R.This work assumes the covariance matricesQt and R to remain constant during the operation of the fault detection system. This implies that the statistical characteristics of the process noise of the dynamic states and parameters and the measurement noise are assumed to remain unchanged. While it is possible to conceive a modification to the EKF with robust, rescalable and possibly adaptive process, and measurement noise covariance matrices [45–47], this is beyond the scope of this paper.Modern control systems operate on a discrete basis. The nonlinear estimation procedure introduced in the following system will be also implemented on an embedded microcontroller unit, performing measurements at discrete intervalst=T(k), k=1,2,3,(…). The continuous augmented dynamics in (9) will be propagated using the Euler’s method by the system (13)xak+1=fdxak,uk+wkzk=hdxak,uk+vk.


## 2.1.3. Joint State and Parameter Estimation

The extended Kalman filter is proposed in this work to estimate the nonlinear augmented states and parameters of the active structure [43, 44]. The algorithm considered here is often referred to as the continuous-discrete or hybrid formulation in literature [42–44], since the EKF takes its measurements and yields new state estimates at discrete sample times, while the system model is evaluated by numerical integration. The observer algorithm is initialized by a first state estimate at zero time (14)x^a0+=Exa0that is determined by the user as an initial guess of the states and parameters. Similarly, the initial covariance matrix of the error of this state estimate is given by the user as (15)P0+=Exa0-x^a0+xa0-x^a0+T.The primary estimate—ora priori estimate, as it is often denoted in literature—of the augmented state is obtained using simulation. The dynamics of the driven point-mass-damper with the unknown parameter vector is described by the continuous augmented function f(x^a(k-1)+,u(k)). Moreover, another assumption is that the parameters p(t)=k(t)b(t)c(t)1/m(t)T remain constant between the samples. The continuous state is then propagated using (16)x^ak-=fcx^ak-1+,uk,between samples t=(k-1)Ts to t=(k)Ts to obtain a priori estimates, denoted here by the - subscript. Though several methods exist to solve ordinary differential equations numerically, this work utilizes the first order Euler method. The continuous augmented function fc(x^a(k-1)+,u(k)) can be safely considered as (17)q˙t=x˙1t=x2tq¨t=x˙2t=-x1tx3tx6t-x2tx4tx6t+x5tx6tutk˙t=x˙3t=0b˙t=x˙4t=0c˙t=x˙5t=01m˙t=x˙6t=0for flexible active structures with dominant first resonant frequencies.In addition to the computation of the a priori state estimate, it is also necessary to evaluate the a priori update of the covariance matrix of the state estimates(18)P˙t=Zx^atPt+PtZTx^at+Qt,where Z(x^a(t)) is the Jacobian of (17) with respect to the augmented state (19)Zx^at=∂fcx^at∂x^atxat=x^at,which for the formerly described augmented beam model it is given by (20)Zx^at=010000-x6tx3t-x6tx4t-x6tx1t-x6tx2t-x6tut-x1tx3t--x2tx4t--utx5t000000000000000000000000.Similar to the numerical propagation of the augmented state, the covariance matrix is obtained by simulating (18) from sample t=(k-1)Ts to t=(k)Ts.The discrete part of the EKF uses the current measurement to update the a priori state to get a reviseda posteriori state estimate, which is denoted by the + subscript. The first step is to compute the Kalman gain K(k) using (21)Kk=Pk-LTLPk-LT+MRMT-1.The a priori state estimate is then revised by the measurement at time(k) with its relative importance weighted to the simulation estimate by the Kalman gain K(k), which is computed at the next step. The a posteriori update of the augmented state is (22)x^ak+=x^ak-+Kkzk-hx^ak-,uk.The state error covariance matrix is updated by(23)Pk+=I-KkLPk-I-KkLT+KkMRMTKkT,where the expression L is the Jacobian of the measurement equation hc(xa(t)) with respect to the state xa(t)(24)L=∂hcxat∂xatxat=x^at-=00000x6t000001000001000001000001,and matrix M is the Jacobian of hc(xa(t)) with respect to the measurement noise term va(t)(25)M=∂hcxat∂vatxat=x^at-=1.The a posteriori augmented state estimate x^a(k)+ and covariance matrix estimate P(k)+ are passed onto the next computation step, where they become the basis for the next a priori estimates.


## 2.1.4. Resonance Estimates

The first dominant resonant frequency of the beam is approximated on the basis of the estimated parametersp^(k) outside the observer algorithm. Given the parameter estimates k^(k),b^(k),m^(k) at t=kTs and assuming that the damped natural frequency of the dynamically analogous SDOF system in (1) is equal to the dominant resonance frequency of the structure [14], the frequency is (26)f^k=4k^km^k-b^k24πm^k.The algorithm monitors this frequency estimate throughout the operation of the system. As it has been noted by Jassim et al. [20] and others, monitoring the change of natural frequency is a feasible tool to indicate damage occurrence. In case that the estimate is outside the tolerance band (27)ek=0iffn-ft≤f^k≥fn+ft1iffn-ft≥f^k1iffn+ft≤f^k,a fault detection signal e(k) is produced and retained unless cleared. The resonance of the nominal structure is denoted by fn and ft is a symmetric tolerance band.


## 2.1.5. Algorithm Summary

The resulting real-time fault detection strategy with joint state and parameter estimation can be summarized in the following algorithm.Algorithm 1.
Perform the following operations at each sampling instant(k).(1)
Propagate (16) and (18) by numerical simulation using the Euler method. Obtain the a priori estimates of the augmented state x^a(k)- and covariance matrix P(k)-.(2)
Sample the actual deflectionz(k) at the end of the cantilever beam.(3)
Compute the Kalman gainKf(k) using (21).(4)
In order to compute the a posteriori state estimatex^a(k)+, use the new measurement sample z(k) and Kalman gain update Kf(k) to revise the a priori state estimate through (22).(5)
Compute the a posteriori update of the covariance matrixP(k)+ in (23).(6)
Retain the a posteriori state and covariance estimatesx^a(k)+ and P(k)+ in the memory for the next estimation step.(7)
Approximate the damped natural frequencyf^(k) of the beam according to (26).(8)
If the frequency estimatef^(k) is outside the tolerance band in (27), enable the fault detection signal e(k).(9)
At the next sample restart the procedure from 1, use the estimates from Step 6.


## 2.2. Experimental Hardware and Software

### 2.2.1. Experimental Hardware

The proposed fault detection system was tested on a thin flexible cantilever beam with piezoceramic actuation. The test hardware composition will be described only briefly as the test bed was previously introduced in detail in previous work, including [5, 35, 48]. The schematic setup of the experimental hardware, including the active thin structure and the prototyping system, is illustrated in Figure 1.The cantilever is a 550 × 40 × 3 mm beam, fixed at one end and freely moving at the other. As it is indicated in the side view of the beam in Figure1, it undergoes flexural vibration in the direction of its thickness. The beam is made of commercially pure aluminium, having mechanical properties which render the structure more flexible than it would be possible with regularly available aluminium alloys.To further decrease its first resonant frequency, a heavy nut, bolt, and a set of four washers are mounted at the free end of the beam. The nominal weight installed at the end is necessary to shift its original resonance (8.17 Hz) to the area, where is possible to run the proposed fault detection system in real-time on the hardware with limited computing power and automatically deployed code. The weight of the aluminum cantilever beam wasmb=0.182 kg, and the total weight of the nut-washer-bolt assembly is mw=0.305 kg, making the total physical weight of m=0.487 kg. This base configuration assumed here therefore consists of the beam with a weight, with a first nominal resonance decreased to 2.26 Hz. This nominal configuration is taken as the initial configuration in each experiment, which is then further modified by the experiments.An industrial laser triangulation device is mounted near the free end of the beam, to provide position measurements for the fault detection algorithm. The Keyence LK-G 32 laser sensor head is pointed at a single location at the beam end, while its signal is preprocessed with a low pass filter in the Keyence LK-G3001V processing unit. The Keyence LK-G3001V unit provides an analog signal scaling the ±5 mm maximal measurement range into a 0–5 V signal, compatible with the analog measurement range of the prototyping board. The laser processing unit is directly connected to the AI1 analog input port of MCU. The measurements are rescaled and shifted to readings in meters within the initial portion of the estimation algorithm.The actuation of the active thin structure is realized with a set of two double layer piezoceramic transducers. The MIDÉ QuickPack QP45n double layer actuators are mounted on both fixed sides of the beam and receive the same electric signal with opposing polarity. The high voltage driving signal is provided through a MIDÉ 1225 capacitive precision amplifier with a 20× V/V gain on the signal input. The active control signal is emulated by a pseudorandom binary voltage signal, having its low level at 0 V and high level at 4 V, resulting in 80 V impulses to the transducers. The PRBS signal has a 10 bit/s base sample (in order to excite the first resonance of the beam) and is generated by an Agilent 3351A signal generator unit. In addition to the measured position inputz(k), the extended Kalman filter requires the input signal u(k) as well; therefore the output from the signal generator is directly connected to the AI1 analog input port of the prototyping board and then rescaled to the amplified levels within the EKF algorithm.The low-cost electronics prototyping board employed in this work is an Arduino Mega 2560 Revision 3 device, featuring a simple 8-bit Atmel ATmega 2560 microcontroller unit (Figure2). The MCU used in this board has a clock speed of 16 MHz and a 256 kB large flash memory for storing user algorithms [49], enabling its online external use with the Simulink product range. The detected fault signal e(k) is passed on to the DO9 digital output of the board with an LED connected to it for signal indication. Although this is only a binary information relating to the occurrence of a mechanical fault, the parameter estimates may also be output using a pulse-width modulated signal or encoded into several digital outputs. The board is powered through the USB port, which also serves as the means for programming and online data transfer for logging and storage. Programming of the board and data logging is realized using a laptop computer with an Intel Core i3 CPU running at 2.4 GHz and 4 GB operating memory. The computing power of the design computer has no effect on the speed of the proposed algorithm, as it is running in a stand alone mode on the Arduino Mega 2560 prototyping board.


### 2.2.2. Experimental Software

The proposed estimation and fault detection algorithm was developed within the Mathworks Matlab and Simulink framework. The software was implemented in Matlab and Simulink 8.0 (R2012b), running under the Microsoft Windows 7 × 64 operating system. In addition to the programming role of this software implementation setup, the prototyping board was also providing the estimates and other signals through the USB connection in external mode. There was an approximately 4 s communication delay between the prototyping board and the Simulink data logging and visualizing interface. Even though this portion of the experiments is not shown on figures here, the proposed algorithm was still fully functional on the MCU during this initialization period.The simplified Simulink scheme of the software implementation is featured in Figure3. Communication between the analog input and digital output ports of the Arduino Mega 2560 board and the real-time logging capability was ensured using the Simulink Target for use with Arduino Hardware (v 2.0) suite, providing access to the peripherals of the MCU and a possibility to deploy Simulink schemes onto the hardware. The 10-bit in-processor A/D converter [49] provides a decimal integer reading of the voltages, which is converted using type conversion blocks. The fault signal is transferred into the digital port using the corresponding output block. Internal state estimates are displayed and saved using a scope block, with its logging capability enabled.Figure 3
Block scheme of the algorithm running on the microcontroller unit.In addition to the readily available standardized Simulink blocks, the core of the algorithm featured previously was developed in the Matlab M-script language and then automatically transcribed into a C/C++ version for the embedded MCU. The custom algorithm block noted as EKF in Figure3 contains means to manipulate voltage inputs to physical units, the extended Kalman filter providing state and parameter estimates, and the simple algorithm generating the fault signal based on the calculated damped first resonant frequency. This M-script code is automatically recompiled into embedded C/C++ code and downloaded to the target hardware using the Embedded Coder (v 6.3) and Simulink Coder (v 8.3). The compilation itself is performed by the Microsoft Visual C++ 2010 compiler.


## 2.2.1. Experimental Hardware

The proposed fault detection system was tested on a thin flexible cantilever beam with piezoceramic actuation. The test hardware composition will be described only briefly as the test bed was previously introduced in detail in previous work, including [5, 35, 48]. The schematic setup of the experimental hardware, including the active thin structure and the prototyping system, is illustrated in Figure 1.The cantilever is a 550 × 40 × 3 mm beam, fixed at one end and freely moving at the other. As it is indicated in the side view of the beam in Figure1, it undergoes flexural vibration in the direction of its thickness. The beam is made of commercially pure aluminium, having mechanical properties which render the structure more flexible than it would be possible with regularly available aluminium alloys.To further decrease its first resonant frequency, a heavy nut, bolt, and a set of four washers are mounted at the free end of the beam. The nominal weight installed at the end is necessary to shift its original resonance (8.17 Hz) to the area, where is possible to run the proposed fault detection system in real-time on the hardware with limited computing power and automatically deployed code. The weight of the aluminum cantilever beam wasmb=0.182 kg, and the total weight of the nut-washer-bolt assembly is mw=0.305 kg, making the total physical weight of m=0.487 kg. This base configuration assumed here therefore consists of the beam with a weight, with a first nominal resonance decreased to 2.26 Hz. This nominal configuration is taken as the initial configuration in each experiment, which is then further modified by the experiments.An industrial laser triangulation device is mounted near the free end of the beam, to provide position measurements for the fault detection algorithm. The Keyence LK-G 32 laser sensor head is pointed at a single location at the beam end, while its signal is preprocessed with a low pass filter in the Keyence LK-G3001V processing unit. The Keyence LK-G3001V unit provides an analog signal scaling the ±5 mm maximal measurement range into a 0–5 V signal, compatible with the analog measurement range of the prototyping board. The laser processing unit is directly connected to the AI1 analog input port of MCU. The measurements are rescaled and shifted to readings in meters within the initial portion of the estimation algorithm.The actuation of the active thin structure is realized with a set of two double layer piezoceramic transducers. The MIDÉ QuickPack QP45n double layer actuators are mounted on both fixed sides of the beam and receive the same electric signal with opposing polarity. The high voltage driving signal is provided through a MIDÉ 1225 capacitive precision amplifier with a 20× V/V gain on the signal input. The active control signal is emulated by a pseudorandom binary voltage signal, having its low level at 0 V and high level at 4 V, resulting in 80 V impulses to the transducers. The PRBS signal has a 10 bit/s base sample (in order to excite the first resonance of the beam) and is generated by an Agilent 3351A signal generator unit. In addition to the measured position inputz(k), the extended Kalman filter requires the input signal u(k) as well; therefore the output from the signal generator is directly connected to the AI1 analog input port of the prototyping board and then rescaled to the amplified levels within the EKF algorithm.The low-cost electronics prototyping board employed in this work is an Arduino Mega 2560 Revision 3 device, featuring a simple 8-bit Atmel ATmega 2560 microcontroller unit (Figure2). The MCU used in this board has a clock speed of 16 MHz and a 256 kB large flash memory for storing user algorithms [49], enabling its online external use with the Simulink product range. The detected fault signal e(k) is passed on to the DO9 digital output of the board with an LED connected to it for signal indication. Although this is only a binary information relating to the occurrence of a mechanical fault, the parameter estimates may also be output using a pulse-width modulated signal or encoded into several digital outputs. The board is powered through the USB port, which also serves as the means for programming and online data transfer for logging and storage. Programming of the board and data logging is realized using a laptop computer with an Intel Core i3 CPU running at 2.4 GHz and 4 GB operating memory. The computing power of the design computer has no effect on the speed of the proposed algorithm, as it is running in a stand alone mode on the Arduino Mega 2560 prototyping board.


## 2.2.2. Experimental Software

The proposed estimation and fault detection algorithm was developed within the Mathworks Matlab and Simulink framework. The software was implemented in Matlab and Simulink 8.0 (R2012b), running under the Microsoft Windows 7 × 64 operating system. In addition to the programming role of this software implementation setup, the prototyping board was also providing the estimates and other signals through the USB connection in external mode. There was an approximately 4 s communication delay between the prototyping board and the Simulink data logging and visualizing interface. Even though this portion of the experiments is not shown on figures here, the proposed algorithm was still fully functional on the MCU during this initialization period.The simplified Simulink scheme of the software implementation is featured in Figure3. Communication between the analog input and digital output ports of the Arduino Mega 2560 board and the real-time logging capability was ensured using the Simulink Target for use with Arduino Hardware (v 2.0) suite, providing access to the peripherals of the MCU and a possibility to deploy Simulink schemes onto the hardware. The 10-bit in-processor A/D converter [49] provides a decimal integer reading of the voltages, which is converted using type conversion blocks. The fault signal is transferred into the digital port using the corresponding output block. Internal state estimates are displayed and saved using a scope block, with its logging capability enabled.Figure 3
Block scheme of the algorithm running on the microcontroller unit.In addition to the readily available standardized Simulink blocks, the core of the algorithm featured previously was developed in the Matlab M-script language and then automatically transcribed into a C/C++ version for the embedded MCU. The custom algorithm block noted as EKF in Figure3 contains means to manipulate voltage inputs to physical units, the extended Kalman filter providing state and parameter estimates, and the simple algorithm generating the fault signal based on the calculated damped first resonant frequency. This M-script code is automatically recompiled into embedded C/C++ code and downloaded to the target hardware using the Embedded Coder (v 6.3) and Simulink Coder (v 8.3). The compilation itself is performed by the Microsoft Visual C++ 2010 compiler.


## 2.3. Experiment Design and Settings

### 2.3.1. Experiment Design and Procedures

The experiments featured in this paper were designed to demonstrate the functionality of the proposed estimation and fault detection algorithm and the viability of its real-time implementation on low-cost microcontroller hardware. Experiments were lasting 80 s with the pseudorandom binary driving signal supplied to the actuators during the entire duration of the test. After initializing the experiments with completely identical tuning and settings, the estimated dynamic states and the parameters converged to the values in agreement with the mechanical properties of the nominal beam. The majority of this starting convergence is not indicated in this paper, due to the previously mentioned communication issue between the MCU and the laptop computer, which is performing data logging and visualization.In addition to a nominal beam configuration with the beam itself and the heavy bolt-washer-nut combination, 6 types of different weights ranging from 16 g to 345 g were added to the end of the beam approximately at the ~40 s time mark (see Figure4). The added weights emulated the effects of mechanical failures and changes with different severity. As these weights were placed manually, there may be a ±1-2 s variance in this time between the individual tests. The weights and individual experiment types are summarized in Table 1. The added weights changed the overall dynamic behavior of the beam and yielded parameters which diverged from the nominal values. If the approximate of the first damped resonance f^(k) left the ft=±0.1 Hz tolerance band of the nominal first resonance at fn=2.27 Hz, the alarm signal e(k) was switched on and this fact was logged along with the rest of the parameters. The first resonant frequency of the thin structure in its nominal configuration was determined by supplying the piezoceramic transducers with a PRBS signal and performing a position measurement. This time domain measurement was then analysed using the fast Fourier transform (FFT) to extract the first fundamental resonant frequency.Table 1
Summary of experimental scenarios.
Number
Weight
Change
Nominal
Color
1
0 g
0%
✓
Blue
2
16 g
3%
✗
Red
3
50 g
10%
✗
Green
4
99 g
20%
✗
Black
5
195 g
40%
✗
Magenta
6
246 g
50%
✗
Orange
7
345 g
70%
✗
GreyFigure 4
Estimated dynamic states for the experiment with the 246 g added weight.


### 2.3.2. Experimental Parameters and Settings

All experiments were initialized with the initial state estimate of(28)x^a0+=0010011E-410.5T,denoted in the experiments with an upside down black triangle located at zero time. Likewise, the initial covariance estimates remained the same throughout the experiments with (29)P0+=1E-50000001E-40000001E+50000001E-20000001E-50000001E+0.The values in the P0+ matrix provide a compromise between a realistic guess of the variance in the initial state estimate error and an incremental tuning of the EKF algorithm.The covariance of the single measurement signal was estimated based on the square of the maximal amplitude resolution of the A/D input of the prototyping board in meters, yielding ~2.5 V at 10 bits [49] or in the terms of the measurement covariance (30)R=R=1E-5.The greatest freedom in tuning the behavior of the algorithm can be achieved by manipulating the continuous-time process noise covariance matrix. The process noise matrix was tuned for best performance based on incremental testing and the expected range of changes in the physical variables and then fixed for all experiments featured in this paper as (31)Qt=5E-3000001E+1000001E-1000001E-12000003E-1,with its discrete-time analogy recomputed by Q=Qt/Ts.Global sampling times were set asTs = 0.1 s, while the continuous model was numerically integrated using the Euler method with a dt=0.02 s time step.


## 2.3.1. Experiment Design and Procedures

The experiments featured in this paper were designed to demonstrate the functionality of the proposed estimation and fault detection algorithm and the viability of its real-time implementation on low-cost microcontroller hardware. Experiments were lasting 80 s with the pseudorandom binary driving signal supplied to the actuators during the entire duration of the test. After initializing the experiments with completely identical tuning and settings, the estimated dynamic states and the parameters converged to the values in agreement with the mechanical properties of the nominal beam. The majority of this starting convergence is not indicated in this paper, due to the previously mentioned communication issue between the MCU and the laptop computer, which is performing data logging and visualization.In addition to a nominal beam configuration with the beam itself and the heavy bolt-washer-nut combination, 6 types of different weights ranging from 16 g to 345 g were added to the end of the beam approximately at the ~40 s time mark (see Figure4). The added weights emulated the effects of mechanical failures and changes with different severity. As these weights were placed manually, there may be a ±1-2 s variance in this time between the individual tests. The weights and individual experiment types are summarized in Table 1. The added weights changed the overall dynamic behavior of the beam and yielded parameters which diverged from the nominal values. If the approximate of the first damped resonance f^(k) left the ft=±0.1 Hz tolerance band of the nominal first resonance at fn=2.27 Hz, the alarm signal e(k) was switched on and this fact was logged along with the rest of the parameters. The first resonant frequency of the thin structure in its nominal configuration was determined by supplying the piezoceramic transducers with a PRBS signal and performing a position measurement. This time domain measurement was then analysed using the fast Fourier transform (FFT) to extract the first fundamental resonant frequency.Table 1
Summary of experimental scenarios.
Number
Weight
Change
Nominal
Color
1
0 g
0%
✓
Blue
2
16 g
3%
✗
Red
3
50 g
10%
✗
Green
4
99 g
20%
✗
Black
5
195 g
40%
✗
Magenta
6
246 g
50%
✗
Orange
7
345 g
70%
✗
GreyFigure 4
Estimated dynamic states for the experiment with the 246 g added weight.


## 2.3.2. Experimental Parameters and Settings

All experiments were initialized with the initial state estimate of(28)x^a0+=0010011E-410.5T,denoted in the experiments with an upside down black triangle located at zero time. Likewise, the initial covariance estimates remained the same throughout the experiments with (29)P0+=1E-50000001E-40000001E+50000001E-20000001E-50000001E+0.The values in the P0+ matrix provide a compromise between a realistic guess of the variance in the initial state estimate error and an incremental tuning of the EKF algorithm.The covariance of the single measurement signal was estimated based on the square of the maximal amplitude resolution of the A/D input of the prototyping board in meters, yielding ~2.5 V at 10 bits [49] or in the terms of the measurement covariance (30)R=R=1E-5.The greatest freedom in tuning the behavior of the algorithm can be achieved by manipulating the continuous-time process noise covariance matrix. The process noise matrix was tuned for best performance based on incremental testing and the expected range of changes in the physical variables and then fixed for all experiments featured in this paper as (31)Qt=5E-3000001E+1000001E-1000001E-12000003E-1,with its discrete-time analogy recomputed by Q=Qt/Ts.Global sampling times were set asTs = 0.1 s, while the continuous model was numerically integrated using the Euler method with a dt=0.02 s time step.


## 3. Results and Discussion

### 3.1. Dynamic State Estimates

A typical experiment is shown in Figure4, where the nominal beam configuration is altered by adding the 246 g weight to the free end at time t=41.4 s, thereby simulating a mechanical failure of the structure. The results logged during the experiments (from top to bottom) are the estimated and measured position x^1,z≈q, the estimated and calculated velocity x^2,Δz≈q˙, and the measured voltage input to the transducers u. The PRBS input signal is supplied constantly, regardless of the change.Both dynamic state estimates start from the common initial state ofx^0+=00T, shown as black triangles in the figure. As it has been previously mentioned, the initial portion of the measurement cannot be logged, due to the properties of the hardware implementation. The difference in the position and velocity estimates before and after the emulated failure is visually observable. The amplitude of the estimated position and speed signals becomes smaller after the failure, while due to the decrease in the first resonant frequency the oscillations are less frequent, which is also visually noticeable in the figure.The estimated positionx^1(k) is compared to the measured signal z(k). As it is evident from the figure, the estimates agree with the measurement very well. The estimate error x~1=|x^1(k)-z(k)| was calculated to be x~1=3.8E-5±3.5E-5 m for this demonstration case. As direct speed measurements are not available, the estimated speed x^2(k) is compared to a speed signal calculated from the position measurement using [38] (32)Δzk≈zk+1-zk-12Ts.The difference between the computed and estimated velocity is larger than that in the case of the position. The EKF algorithm tends to overestimate the velocity, in comparison with the calculated values. However, one has to note that the calculated velocity signal is hardly exact and cannot be considered as the conventionally true velocity. It is not less likely that the estimates are actually closer to the true velocity than the calculated values. Keeping this in perspective, the velocity estimate error x~2=|x^2(k)-Δz(k)| was calculated to be x~2=4.9E-3±4.1E-3 ms−1 for this demonstration case.


### 3.2. Parameter Estimates and Failure Detection

The estimated first resonance frequency with the instances of fault detection (vertical dashed lines) for the different failure emulation cases is illustrated in Figure5 along with the estimated model parameters. The figure shows (from top to bottom) the estimates of the first damped resonance frequency f^, the dynamic mass of the structure m^=x^6, the dynamic equivalent stiffness coefficient k^=x^3, the dynamic equivalent damping b^=x^4, and the voltage-force conversion coefficient c^=x^5. Again, the hardware communication with the Arduino board prevents the first 3.8–4.5 s of the measurement to be logged; however, the real-time algorithm running on the MCU is fully functional at this time and the parameters estimates are converging from the initial estimates (triangles at t=0 s) to the nominal values.Figure 5
Estimated first resonance and model parameters.After a brief convergence period, the frequency estimates tend to the neighborhood of the nominal frequency within ±0.015 Hz (maximal mean value, betweent = 20–30 s), well within the tolerance band of the failure detection limits. This is also true for the estimated parameters; the nominal mass, stiffness, damping, and conversion coefficients are located in the same relative region for all experimental scenarios. The most variation in the nominal case is visible in the conversion coefficient, which is harder to estimate in general [33].After adding the different weights to emulate the failure, the parameters start to diverge from the nominal levels and so do the estimated resonances. As soon as the estimated resonance frequency leaves the preset tolerance band, the fault signal is activated. The time of detection depends entirely on the severity of the fault and varies between 0 s for the largest weight to ~36.5 s for the smallest (see Table2). The estimated frequency naturally corresponds to the amount of weight added, and it varies between the nominal frequency f^=2.26 Hz to an average of f^=1.81 Hz for the largest weight (see Table 2).Table 2
Summary of estimated results (mean estimates att = 50–60 s).
Number
Weight
Detected
f
^ (Hz)
m
^ (kg)
k
^ (Nm−1)
b
^ (Nsm−1)
c
^ (NkV−1)
1
0 g
✗
2.26
0.42
83.8
1.7
0.5
2
16 g
76.5 s
2.23
0.42
83.2
1.8
0.5
3
50 g
49.2 s
2.20
0.44
83.8
1.8
0.4
4
99 g
46.1 s
2.11
0.47
82.9
1.7
0.4
5
195 g
44.4 s
1.98
0.56
86.9
1.9
0.4
6
246 g
41.7 s
1.91
0.60
86.5
2.2
0.3
7
345 g
39.9 s
1.81
0.61
79.0
1.9
0.2The estimated parameters react to the change in mechanical properties on various levels (see Table2). Computed as mean values between the times t = 50–80 s, the mass parameter varies the most as one would expect, ranging from m^=0.42 kg to m^=0.61 kg for the different fault emulation cases. The stiffness and damping coefficients do not vary as much for the individual cases; average stiffness estimates k^ stay within 10% of the nominal stiffness, while the shift in the damping coefficient is also moderate. The effect of the actuators tends to diminish as greater and greater weights are placed at the end of the beam. This is demonstrated in Figure 5, where the estimated voltage conversion parameter c^ starts to decrease with the increasing weight as the actuators are capable to provide less force with more weight.While the first resonant frequency estimates are stable in all experimental scenarios, the parameter estimates reveal possible stability issues with the heaviest weight (Experiment number 7, grey). Although the initial convergence after the change is meaningful, the parameters begin to decline to levels that are not representative for the change. This may cause numerical instability in the estimation algorithm, which—in case only a binary fault signal is required—does not pose a problem, as it is likely that the change will be detected much earlier than the instability occurs. This divergence of parameters is largely due to the fact that the conversion coefficientc^ will be understandably compromised, if extensive mechanical faults are present.It is important to note that the identified physical parameters are meaningful for the dynamic model of the structure and are not necessarily identical to the parameters one would obtain using exact first-principle measurements. Although the estimated parameters have a physical meaning and sense, they are evaluated in a way to match the measured position to the behavior of the model, in accordance with the general EKF formulation. For example, a change in the estimated mass is rather a dynamic change and one cannot expect to obtain the exact weight change. The nominal beam-weight configuration with no additional change ism0=0.487 kg, which is estimated to a lower value of m^=420 kg. This is to be expected and normal as the algorithm deals with the parameters only on the level of the dynamic behavior of the structure. Also, if only one physical parameter is changed—as it is the case in this work—its effect is translated to a change in all parameters.The global sampling times used in this experiment were too fast for the external mode Simulink simulation and caused a computation overflow, meaning that the time necessary to evaluate the algorithm along with the communication and logging procedures was longer than the allocated sampling time. The overflow detection was enabled in the simulation and indicated by lighting an LED connected to the output port DO13 of the Arduino Mega 2560 prototyping platform. This overflow posed only a problem in the case where the algorithm was running with the data logging facilities enabled. In case the fault detection scheme ran on the MCU in a stand-alone mode, the sampling periods were strictly enforced. Data logging is only necessary to present the experimental results in this work; therefore it is not considered to be a serious limitation of the application. Moreover, the estimated parameters can be output in computationally more efficient ways than the use of USB communication. The proposed fault detection algorithm was able to run down to a sampling rate ofTs=0.08 s on the Atmel ATmega 2560 MCU. Increasing the number of integration steps has the greatest effect on computation load, as a dt=Ts/20 s was unable to run without overflows with a Ts=0.1 sampling.


### 3.3. Verification of the Frequency Estimates

The verification of the mean EKF frequency estimates was carried out by comparing the mean observed first resonancef^ in each experimental scenario to the first resonant frequency obtained from the power spectral density (PSD) Sxx(f) of the measured displacement signal. The mean frequency estimate was computed from the EKF results after the additional weights have been placed on the beam, that is, between 40–80 s of the measurements. The power spectral density of the displacement measurement signal z(k) was computed on the same interval. The spectral estimate was acquired using the Yule-Walker parametric autoregressive analysis of the displacement signal, assuming an 8th order autoregressive model of the PSD. The results of this comparison are demonstrated in Figure 6, where the PSD from the displacement is marked by the thicker solid lines, while the vertical thin dashed lines indicate the average EKF estimates. From this, it is clear that the difference between the EKF estimates and the first resonances judged by the raw signal itself match very well.Figure 6
Comparison of power spectral density frequency estimates with average EKF estimates.To quantify this verification procedure in more detail, Table3 lists the periodogram estimates, EKF estimates, and their absolute and percentual difference. The absolute difference of the frequencies is well below |Δf|=0.05 Hz in all experimental cases, meaning a difference of <1.5% in all the listed cases.Table 3
Verification experimental results (mean EKF estimates att = 40–80 s).
Experimental setting
0 g
16 g
50 g
99 g
195 g
246 g
345 g
Periodogram estimate (Hz)
2.289
2.239
2.214
2.139
2.015
1.940
1.841
EKF estimate (Hz)
2.261
2.238
2.209
2.131
2.014
1.945
1.835
Abs. difference (Hz)
0.028
0.001
0.005
0.008
0.001
0.005
0.006
Percentual difference (%)
1.22
0.04
0.023
0.37
0.05
0.26
0.33Note that comparing the direct coefficient estimates (m^,b^,k^,c^) to the true or assumed values would be inconsequential in this case. While the frequency estimates of the EKF and the real resonant frequency of the structure can be directly related and verified, the simple model assumed in (17) does not allow a direct assessment of the observer accuracy. As it has been previously noted, the coefficients of the model represent the dynamic equivalent of the system, not a physically accurate mechanical system. The estimation of the dynamically equivalent m^(k),b^(k),k^(k), and c^(k) is useful for structural detection and monitoring or as online identification in adaptive control, but naturally, it is absurd to expect them to fully represent the beam from the mechanical standpoint.


## 3.1. Dynamic State Estimates

A typical experiment is shown in Figure4, where the nominal beam configuration is altered by adding the 246 g weight to the free end at time t=41.4 s, thereby simulating a mechanical failure of the structure. The results logged during the experiments (from top to bottom) are the estimated and measured position x^1,z≈q, the estimated and calculated velocity x^2,Δz≈q˙, and the measured voltage input to the transducers u. The PRBS input signal is supplied constantly, regardless of the change.Both dynamic state estimates start from the common initial state ofx^0+=00T, shown as black triangles in the figure. As it has been previously mentioned, the initial portion of the measurement cannot be logged, due to the properties of the hardware implementation. The difference in the position and velocity estimates before and after the emulated failure is visually observable. The amplitude of the estimated position and speed signals becomes smaller after the failure, while due to the decrease in the first resonant frequency the oscillations are less frequent, which is also visually noticeable in the figure.The estimated positionx^1(k) is compared to the measured signal z(k). As it is evident from the figure, the estimates agree with the measurement very well. The estimate error x~1=|x^1(k)-z(k)| was calculated to be x~1=3.8E-5±3.5E-5 m for this demonstration case. As direct speed measurements are not available, the estimated speed x^2(k) is compared to a speed signal calculated from the position measurement using [38] (32)Δzk≈zk+1-zk-12Ts.The difference between the computed and estimated velocity is larger than that in the case of the position. The EKF algorithm tends to overestimate the velocity, in comparison with the calculated values. However, one has to note that the calculated velocity signal is hardly exact and cannot be considered as the conventionally true velocity. It is not less likely that the estimates are actually closer to the true velocity than the calculated values. Keeping this in perspective, the velocity estimate error x~2=|x^2(k)-Δz(k)| was calculated to be x~2=4.9E-3±4.1E-3 ms−1 for this demonstration case.


## 3.2. Parameter Estimates and Failure Detection

The estimated first resonance frequency with the instances of fault detection (vertical dashed lines) for the different failure emulation cases is illustrated in Figure5 along with the estimated model parameters. The figure shows (from top to bottom) the estimates of the first damped resonance frequency f^, the dynamic mass of the structure m^=x^6, the dynamic equivalent stiffness coefficient k^=x^3, the dynamic equivalent damping b^=x^4, and the voltage-force conversion coefficient c^=x^5. Again, the hardware communication with the Arduino board prevents the first 3.8–4.5 s of the measurement to be logged; however, the real-time algorithm running on the MCU is fully functional at this time and the parameters estimates are converging from the initial estimates (triangles at t=0 s) to the nominal values.Figure 5
Estimated first resonance and model parameters.After a brief convergence period, the frequency estimates tend to the neighborhood of the nominal frequency within ±0.015 Hz (maximal mean value, betweent = 20–30 s), well within the tolerance band of the failure detection limits. This is also true for the estimated parameters; the nominal mass, stiffness, damping, and conversion coefficients are located in the same relative region for all experimental scenarios. The most variation in the nominal case is visible in the conversion coefficient, which is harder to estimate in general [33].After adding the different weights to emulate the failure, the parameters start to diverge from the nominal levels and so do the estimated resonances. As soon as the estimated resonance frequency leaves the preset tolerance band, the fault signal is activated. The time of detection depends entirely on the severity of the fault and varies between 0 s for the largest weight to ~36.5 s for the smallest (see Table2). The estimated frequency naturally corresponds to the amount of weight added, and it varies between the nominal frequency f^=2.26 Hz to an average of f^=1.81 Hz for the largest weight (see Table 2).Table 2
Summary of estimated results (mean estimates att = 50–60 s).
Number
Weight
Detected
f
^ (Hz)
m
^ (kg)
k
^ (Nm−1)
b
^ (Nsm−1)
c
^ (NkV−1)
1
0 g
✗
2.26
0.42
83.8
1.7
0.5
2
16 g
76.5 s
2.23
0.42
83.2
1.8
0.5
3
50 g
49.2 s
2.20
0.44
83.8
1.8
0.4
4
99 g
46.1 s
2.11
0.47
82.9
1.7
0.4
5
195 g
44.4 s
1.98
0.56
86.9
1.9
0.4
6
246 g
41.7 s
1.91
0.60
86.5
2.2
0.3
7
345 g
39.9 s
1.81
0.61
79.0
1.9
0.2The estimated parameters react to the change in mechanical properties on various levels (see Table2). Computed as mean values between the times t = 50–80 s, the mass parameter varies the most as one would expect, ranging from m^=0.42 kg to m^=0.61 kg for the different fault emulation cases. The stiffness and damping coefficients do not vary as much for the individual cases; average stiffness estimates k^ stay within 10% of the nominal stiffness, while the shift in the damping coefficient is also moderate. The effect of the actuators tends to diminish as greater and greater weights are placed at the end of the beam. This is demonstrated in Figure 5, where the estimated voltage conversion parameter c^ starts to decrease with the increasing weight as the actuators are capable to provide less force with more weight.While the first resonant frequency estimates are stable in all experimental scenarios, the parameter estimates reveal possible stability issues with the heaviest weight (Experiment number 7, grey). Although the initial convergence after the change is meaningful, the parameters begin to decline to levels that are not representative for the change. This may cause numerical instability in the estimation algorithm, which—in case only a binary fault signal is required—does not pose a problem, as it is likely that the change will be detected much earlier than the instability occurs. This divergence of parameters is largely due to the fact that the conversion coefficientc^ will be understandably compromised, if extensive mechanical faults are present.It is important to note that the identified physical parameters are meaningful for the dynamic model of the structure and are not necessarily identical to the parameters one would obtain using exact first-principle measurements. Although the estimated parameters have a physical meaning and sense, they are evaluated in a way to match the measured position to the behavior of the model, in accordance with the general EKF formulation. For example, a change in the estimated mass is rather a dynamic change and one cannot expect to obtain the exact weight change. The nominal beam-weight configuration with no additional change ism0=0.487 kg, which is estimated to a lower value of m^=420 kg. This is to be expected and normal as the algorithm deals with the parameters only on the level of the dynamic behavior of the structure. Also, if only one physical parameter is changed—as it is the case in this work—its effect is translated to a change in all parameters.The global sampling times used in this experiment were too fast for the external mode Simulink simulation and caused a computation overflow, meaning that the time necessary to evaluate the algorithm along with the communication and logging procedures was longer than the allocated sampling time. The overflow detection was enabled in the simulation and indicated by lighting an LED connected to the output port DO13 of the Arduino Mega 2560 prototyping platform. This overflow posed only a problem in the case where the algorithm was running with the data logging facilities enabled. In case the fault detection scheme ran on the MCU in a stand-alone mode, the sampling periods were strictly enforced. Data logging is only necessary to present the experimental results in this work; therefore it is not considered to be a serious limitation of the application. Moreover, the estimated parameters can be output in computationally more efficient ways than the use of USB communication. The proposed fault detection algorithm was able to run down to a sampling rate ofTs=0.08 s on the Atmel ATmega 2560 MCU. Increasing the number of integration steps has the greatest effect on computation load, as a dt=Ts/20 s was unable to run without overflows with a Ts=0.1 sampling.


## 3.3. Verification of the Frequency Estimates

The verification of the mean EKF frequency estimates was carried out by comparing the mean observed first resonancef^ in each experimental scenario to the first resonant frequency obtained from the power spectral density (PSD) Sxx(f) of the measured displacement signal. The mean frequency estimate was computed from the EKF results after the additional weights have been placed on the beam, that is, between 40–80 s of the measurements. The power spectral density of the displacement measurement signal z(k) was computed on the same interval. The spectral estimate was acquired using the Yule-Walker parametric autoregressive analysis of the displacement signal, assuming an 8th order autoregressive model of the PSD. The results of this comparison are demonstrated in Figure 6, where the PSD from the displacement is marked by the thicker solid lines, while the vertical thin dashed lines indicate the average EKF estimates. From this, it is clear that the difference between the EKF estimates and the first resonances judged by the raw signal itself match very well.Figure 6
Comparison of power spectral density frequency estimates with average EKF estimates.To quantify this verification procedure in more detail, Table3 lists the periodogram estimates, EKF estimates, and their absolute and percentual difference. The absolute difference of the frequencies is well below |Δf|=0.05 Hz in all experimental cases, meaning a difference of <1.5% in all the listed cases.Table 3
Verification experimental results (mean EKF estimates att = 40–80 s).
Experimental setting
0 g
16 g
50 g
99 g
195 g
246 g
345 g
Periodogram estimate (Hz)
2.289
2.239
2.214
2.139
2.015
1.940
1.841
EKF estimate (Hz)
2.261
2.238
2.209
2.131
2.014
1.945
1.835
Abs. difference (Hz)
0.028
0.001
0.005
0.008
0.001
0.005
0.006
Percentual difference (%)
1.22
0.04
0.023
0.37
0.05
0.26
0.33Note that comparing the direct coefficient estimates (m^,b^,k^,c^) to the true or assumed values would be inconsequential in this case. While the frequency estimates of the EKF and the real resonant frequency of the structure can be directly related and verified, the simple model assumed in (17) does not allow a direct assessment of the observer accuracy. As it has been previously noted, the coefficients of the model represent the dynamic equivalent of the system, not a physically accurate mechanical system. The estimation of the dynamically equivalent m^(k),b^(k),k^(k), and c^(k) is useful for structural detection and monitoring or as online identification in adaptive control, but naturally, it is absurd to expect them to fully represent the beam from the mechanical standpoint.


## 4. Conclusions

This paper proposed an algorithm framework for detecting mechanical faults in thin active cantilever beams and other smart structures with dominant first resonance frequencies. This structural health monitoring algorithm was then implemented on a low-cost embedded microcontroller unit.The presented experimental results verify the feasibility of the hardware implementation and demonstrate the functionality of the fault detection algorithm. The emulated mechanical faults were detected in a matter of seconds. In addition to the digital fault signal, the algorithm is capable of identifying the equivalent parameters of the single degree of freedom driven mass-spring-damper representing the dynamics of the active structure.Upon contrasting the estimated first resonance frequencies, it has been found that they match very well to the true frequencies assessed by the power spectral density estimate of the displacement measurement signal.The embedded 8-bit microcontroller unit considered in this work enables sampling speeds down toTs = 0.08 s. This suggests that the use of the proposed fault detection method with the Atmel ATmega 2560 MCU is limited to thin active structures with a largely dominant first resonance frequency of up to ~3-4 Hz. Of course, this work used a machine-transcribed real-time code for the experiments; therefore this is the worst case scenario for the usability of the method on low-priced hardware.Commonly available 8-bit microcontrollers like the Atmel ATmega328 [50] or Texas Instruments C2000 [51] series with an even lower price range are likely to run the proposed algorithm without issues. In case the EKF algorithm is directly transcribed to a lower-level language (C++ or assembler), it is safe to assume that either the complexity of the algorithm can be increased in order to estimate more than one resonant frequency, or the sampling speed may be increased for stiffer structures.Excluding the cost of other components—such as a MEMS semiconductor accelerometer instead of a precise laser system as used here—the unit price of the Atmel ATmega328 embedded microcontroller is currently under 2 EUR for quantities over 100 pieces [52], rendering the proposed fault detection system viable for mass production and deployment.



---

*Source: 506430-2015-05-11.xml*