Monte Carlo-Based Procedure for Determining the Maximum Energy at the Output of Accelerometers

: The solutions presented in this paper can be the basis for mutual comparison of di ﬀ erent types of accelerometers produced by competing companies. An application of a procedure based on the Monte Carlo method to determine the maximum energy at the output of accelerometers is discussed here. The ﬁxed-point algorithm controlled by the Monte Carlo method is used to determine this energy. This algorithm can only be used for the time-invariant and linear measurement systems. Hence, the accelerometer nonlinearities are not considered here. The mathematical models of the accelerometer and the special ﬁlter, represented by the relevant transfer functions, are the basis for the above procedure. Testing results of the voltage-mode accelerometer of type DJB A / 1800 / V are presented here as an example of an implementation of the solutions proposed. Calculation of the energy was executed in Mathcad 14 program with the built-in Programming Toolbar. The value of the maximum output energy determined for a speciﬁed time interval corresponds to the maximum integral-square error of the accelerometer. Such maximum energy can be a comparative ratio just like the accuracy class in the case of instruments used for the static measurements. Hence, the main analytical and technical contributions of this paper concern the development of theoretical procedures and the presentation of their application on the example of a real type of accelerometer.


Introduction
The basic condition for mutual comparison of different types of accelerometers [1][2][3] is that the accelerometers under test have the same frequency bandwidth and it is necessary to develop the model of special filter [4] which is the reference for such an analysis. The frequency bandwidth of the filter must be the same as the operation frequency range of the compared accelerometers.
The reference of energies at the outputs of the accelerometer and the filter is necessary because, from the point of view of the practical implementation of measurements, only the value of energy for the operation range of the accelerometer, not for the full frequency bandwidth, is interesting [5]. Such comparison is carried out by subtracting the energy value at the outputs of both accelerometers. The difference between these energies is defined below as the reference energy and, for the purpose of determining its maximum value, it is convenient to apply the algorithm discussed in Section 3.
The reference energy of an accelerometer corresponds to the integral-square error [6] of the sensor related to this type of error at the output of the special filter. Developing the procedure determines the maximum value of the integral-square error (maximum reference energy) and seems to be the justified proposal for the comparison of different accelerometers.
The basis for the implementation of the task described above is the mathematical model of the accelerometer determined through its parametric identification [7][8][9][10][11][12], as well as the mathematical model of the special filter. The model of the accelerometer is defined by the second-order transfer function [1][2][3], whereas the model of the filter should be a high-order (a range of 10 to 15 order is

Accelerometer Modeling
The model of the accelerometer is defined by using the transfer function, where S a and β are the accelerometer voltage sensitivity and damping factor, respectively, while ω 0 = 2π f 0 , where f 0 denotes the non-damped natural frequency. Parametric identification of the accelerometer model is directed towards determining the parameters of the transfer function given by Equation (1) together with associated uncertainties [8,9]. The first step in this subject is the determination of the measurement points for both frequency responses. This is achieved by applying the back-to-back method commonly used in the measurement practice. The block diagram for this type of identification method is shown in Figure 1. The reference accelerometer is used here as the comparison signal source, while the data acquisition card together with the LabVIEW software are employed for the control of the shaker and recording of the measurement data.

Accelerometer Modeling
The model of the accelerometer is defined by using the transfer function, where and are the accelerometer voltage sensitivity and damping factor, respectively, while = 2 , where denotes the non-damped natural frequency. Parametric identification of the accelerometer model is directed towards determining the parameters of the transfer function given by Equation (1) together with associated uncertainties [8,9]. The first step in this subject is the determination of the measurement points for both frequency responses. This is achieved by applying the back-to-back method commonly used in the measurement practice. The block diagram for this type of identification method is shown in Figure 1. The reference accelerometer is used here as the comparison signal source, while the data acquisition card together with the LabVIEW software are employed for the control of the shaker and recording of the measurement data.
The sensitivity of the tested accelerometer is determined on the basis of the relation, where , , and are the sensitivity of the reference accelerometer, the voltage output of the system for signal conditioning at the output of the accelerometer, and the voltage output of the reference, respectively. Based on the measurement points of the frequency responses of the accelerometer, the estimates , , and are associated with the parameters , , and which occur in the transfer function given by Equation (1), as well as the associated uncertainties , , and , which are determined by applying the WLSP [3,13]. This procedure employs the MCP which utilizes the Box-Muller and the Wichmann-Hill pseudo-random number generator with normal and uniform distributions, respectively. The minimum number of Monte Carlo trials (MCts) should be determined based on the relation, where denotes the assumed coverage probability [13].

Algorithm for Determining the Maximum Energy (ME)
The optimization algorithm for determining the ME is based on the integral equation, where Figure 1. Block diagram of system for back-to-back identification of accelerometer.
The sensitivity S a of the tested accelerometer is determined on the basis of the relation, where S r , V a , and V r are the sensitivity of the reference accelerometer, the voltage output of the system for signal conditioning at the output of the accelerometer, and the voltage output of the reference, respectively. Based on the measurement points of the frequency responses of the accelerometer, the estimates S a , f 0 , and β are associated with the parameters S a , f 0 , and β which occur in the transfer function given by Equation (1), as well as the associated uncertainties u S a , u f 0 , and u β , which are determined by applying the WLSP [3,13]. This procedure employs the MCP which utilizes the Box-Muller and the Wichmann-Hill pseudo-random number generator with normal and uniform distributions, respectively. The minimum number M of Monte Carlo trials (MCts) should be determined based on the relation, where p denotes the assumed coverage probability [13].

Algorithm for Determining the Maximum Energy (ME)
The optimization algorithm for determining the ME is based on the integral equation, Energies 2020, 13, 1552 where while x(t) is any signal at the input of the accelerometer, ξ is the Lagrange multiplier [5], k a (t) and k f (t) are the impulse responses of the accelerometer and the filter, respectively, and T denotes the time of accelerometer testing. The impulse responses k a (t) and k f (t) are determined as the inverse Laplace transforms of the corresponding transfer functions K a (s) and K f (s). In turn, the transfer function K a (s) is obtained by way of the parametric identification of the accelerometer, while the transfer function K f (s) is convenient to be defined by the low-pass Butterworth filter given by, In Equation (6), L is the order of the filter, f c is the filter's cut-off frequency, which is equal to the frequency bandwidth of the accelerometer, and j denotes the imaginary number.
The solution to Equation (4) determines the signal x m (t) that yields the maximum value of the coefficient ξ for which the ratio is satisfied, where E(x) and E(y) are the energies of any input and output signals x(t) and y(t), respectively. The error y(t) is defined by, where and are the output signals of the accelerometer and filter, respectively. Both the signal x(t) as well as the signal x m (t), maximizing the coefficient ξ, are limited in time T and magnitude a. Thus, we can write, The denominator and nominator of Equation (7) are represented by, and Substituting Equations (9) and (10) into Equation (13) and taking into account the relationship between the impulse responses, we finally have When the signal x(t) is limited in magnitude M (i.e., it has a rectangular shape), then we can consider the condition given by, where sgn denotes the signum operation. Hence, the signal maximizing the energy E(y) is represented by, and ξ max is the maximum value of ξ(x, y) defined by Equation (7). The procedure for determining the input signal that maximizes the energy E(y) is based on the FPA [5] and includes the following stages:
Determine the initial input signal x 0 (τ), where ξ 0 is the initial Lagrange multiplier, assumed in advance by the formula,

3.
Determine the maximizing signal on the basis of the iterative algorithm, given by, where J denotes the number of assumed iterations. It should be emphasized that for each iteration step i, the Lagrange multiplier ξ i+1 must be determined based on Equation (7). 4.
Calculate the energy E(y) by using Equation (15).
The signal x m (t) that yields the maximum value of the quotient ξ(x, y) is the one that maximizes the energy E(y). This ME denoted by E max (y) is obtained by substituting the signal x m (t) into Equation (15). Then, we have The energy E max (y) is the accelerometer output energy, referred to as the output energy of filter, and is valid for the accelerometer frequency bandwidth. Figure 2 shows the block diagram of the algorithm described above applied for determining the ME at the output of accelerometers. The switch S has the position 1 for i = 0, and it has the position 2 when i > 0. During the particular iterations i = 1, 2, . . . , J−1 of the algorithm above, the energies E y i and E y i−1 are compared with each other, and the value of the current ME and the corresponding input signal limited in magnitude are stored in the memory. This signal can be then processed by using specialized IT tools, e.g., neural network [30]. The signal producing the ME has the property that any other signal contained in its limitations can only yield less energy than the energy value caused by the maximizing signal [4].  The limited signal ( ) in Figure 3 has seven time switches denoted by , , …, . The number of these switches and the times at which they occur are determined by the FPA, and they constitute the shape of the signal that maximizes the energy at the output of the accelerometer.

Monte Carlo-Based Procedure for Determining the ME
On the basis of obtaining both estimates , , and of the accelerometer model and the associated uncertainties , , and , the values of the parameters denoted by , , and that produce the energy ( ) should be determined. This can be carried out for the rectangular ranges shown in Figure 4 (the lines slanting to the right), by applying the MCP with the number of trials determined using Equation (3), where: and corresponds to the estimates , , and , respectively, in the case of three ranges shown in Figure 4. The final solutions of the ME and corresponding signal are delivered for i = J − 1. It should be noted that this algorithm converges quickly because the maximum number J of iterations does not exceed 30.
The signal producing the ME has the property that any other signal contained in its limitations can only yield less energy than the energy value caused by the maximizing signal [4]. Figure 3 shows examples of the maximizing signal x m (t) limited simultaneously in time and magnitude and three other signals x 1 (t), x 2 (t), and x 3 (t) contained in its limitations denoted by ±a (magnitude) and T (time). The signal producing the ME has the property that any other signal contained in its limitations can only yield less energy than the energy value caused by the maximizing signal [4].  The limited signal ( ) in Figure 3 has seven time switches denoted by , , …, . The number of these switches and the times at which they occur are determined by the FPA, and they constitute the shape of the signal that maximizes the energy at the output of the accelerometer.

Monte Carlo-Based Procedure for Determining the ME
On the basis of obtaining both estimates , , and of the accelerometer model and the associated uncertainties , , and , the values of the parameters denoted by , , and that produce the energy ( ) should be determined. This can be carried out for the rectangular ranges shown in Figure 4 (the lines slanting to the right), by applying the MCP with the number of trials determined using Equation (3), where: and corresponds to the estimates , , and , respectively, in the case of three ranges shown in Figure 4. The limited signal x m (t) in Figure 3 has seven time switches denoted by t 1 , t 2 , . . . , t 7 . The number of these switches and the times at which they occur are determined by the FPA, and they constitute the shape of the signal that maximizes the energy at the output of the accelerometer.

Monte Carlo-Based Procedure for Determining the ME
On the basis of obtaining both estimates S a , f 0 , and β of the accelerometer model and the associated uncertainties u S a , u f 0 , and u β , the values of the parameters denoted by S max a , f max 0 , and β max that produce the energy E max (y) should be determined. This can be carried out for the rectangular ranges shown in Figure 4 (the lines slanting to the right), by applying the MCP with the number of trials determined using Equation (3), where: (21) and α corresponds to the estimates S a , f 0 , and β, respectively, in the case of three ranges shown in Figure 4.  The trials of the MCP executed in three ranges were carried out by employing a uniform distribution. Hence, it is recommended to use the Wichmann-Hill pseudo-random number generator [28], according to the guidelines included in the JCGM supplement 1 [13]. The horizontal ranges of a uniform distribution in Figure 4 are constrained by the values of uncertainty associated with the parameter estimates of the accelerometer, while the vertical ranges of this distribution are contained in the range 〈0, 〉, where = 1/ ( ) is inversely proportional to the range of the parameter estimate variability. Figure 5 shows the block diagram of the MCP for determining the ( ). This procedure includes three main sub-blocks: 1. Back-to-back procedure for identification and determination of the accelerometer parameters and associated uncertainties. 2. Selection of the special filter which is the reference for calculation of the energy at the output of accelerometer. 3. Execution of the FPA.  The trials of the MCP executed in three ranges were carried out by employing a uniform distribution. Hence, it is recommended to use the Wichmann-Hill pseudo-random number generator [28], according to the guidelines included in the JCGM supplement 1 [13]. The horizontal ranges of a uniform distribution in Figure 4 are constrained by the values of uncertainty associated with the parameter estimates of the accelerometer, while the vertical ranges of this distribution are contained in the range 0, γ , where γ = 1/r(α) is inversely proportional to the range of the parameter estimate variability. Figure 5 shows the block diagram of the MCP for determining the E max (y). This procedure includes three main sub-blocks: 1.
Back-to-back procedure for identification and determination of the accelerometer parameters and associated uncertainties.

2.
Selection of the special filter which is the reference for calculation of the energy at the output of accelerometer.

3.
Execution of the FPA. The trials of the MCP executed in three ranges were carried out by employing a uniform distribution. Hence, it is recommended to use the Wichmann-Hill pseudo-random number generator [28], according to the guidelines included in the JCGM supplement 1 [13]. The horizontal ranges of a uniform distribution in Figure 4 are constrained by the values of uncertainty associated with the parameter estimates of the accelerometer, while the vertical ranges of this distribution are contained in the range 〈0, 〉, where = 1/ ( ) is inversely proportional to the range of the parameter estimate variability. Figure 5 shows the block diagram of the MCP for determining the ( ). This procedure includes three main sub-blocks: 1. Back-to-back procedure for identification and determination of the accelerometer parameters and associated uncertainties. 2. Selection of the special filter which is the reference for calculation of the energy at the output of accelerometer. 3. Execution of the FPA.   The input parameter for the MCP is the number of MCts assumed in advance according to Equation (3), while the input parameters for the FPA are: number of iterations J, value of time T, quantization step size ∆ s of the calculation, values of the signal limit a, and initial Lagrange multiplier ξ 0 .

Results
The results of testing the accelerometer of the type DJB A/1800/V are presented below. Typical parameters for the accelerometer, taken from the corresponding datasheet [29], are as follows: The B&K8305 charge mode accelerometer was applied as the reference. The parameters of this accelerometer, contained in the datasheet [31], are: charge sensitivity (±10%): 1.08 pC/g, -amplitude (±10%) and phase (±1 • ) range: 0.2 Hz-1 kHz.
The results of the back-to-back identification are obtained based on the measurement of the amplitude and phase responses carried out for 42 frequencies within the range 30 Hz-6 kHz. This range was selected in such a way as to exceed the value of the resonant frequency since it has a significant influence on the value of the estimates f 0 and β.
The measurement points of both frequency responses, obtained for the above frequencies (the first column denoted by F), are contained in the second and third columns (denoted by A and Φ) of Table 1. Table 1. Measurement points of the frequency responses together with expanded uncertainties (k = 2) for the accelerometer DJB A/1800/V.   The measurement points above are registered by the NI-6221 data acquisition card with a sampling rate equal to 100 kS/s and the computer programs developed in the LabVIEW software.
The expanded uncertainties u(A) and u(Φ) which are contained in the fourth and fifth columns of Table 1 were determined based on Table 1 included in the ISO 16063-21:2003(E) standard (7). These uncertainties are valid for the coverage factor k = 2. The accelerometer model responses (solid lines) obtained by utilizing the WLSP are shown in Figure 6. The covariance matrices involved in the basis kernel of the WLSP used for calculating the estimates of the accelerometer model parameters and applied for determining the uncertainties associated with these estimates are determined by using the MCP. The number M of MCts, which is equal to 2 × 10 5 , was taken as the minimum acceptable value obtained based on Equation (3). The coverage probability p equal to 0.95 was assumed in advance. The measurement points above are registered by the NI-6221 data acquisition card with a sampling rate equal to 100 kS/s and the computer programs developed in the LabVIEW software.
The expanded uncertainties u(A) and u( ) which are contained in the fourth and fifth columns of Table 1 were determined based on Table 1 included in the ISO 16063-21:2003(E) standard (7). These uncertainties are valid for the coverage factor = 2.
The accelerometer model responses (solid lines) obtained by utilizing the WLSP are shown in Figure 6. The covariance matrices involved in the basis kernel of the WLSP used for calculating the estimates of the accelerometer model parameters and applied for determining the uncertainties associated with these estimates are determined by using the MCP. The number of MCts, which is equal to 2 × 10 , was taken as the minimum acceptable value obtained based on Equation (3). The coverage probability equal to 0.95 was assumed in advance. The expanded uncertainties, u(A) and u(Φ), are marked in Figure 6 by the solid vertical lines within the measurement points of both responses. The graphs of both responses are presented by using a logarithmic scale due to the wide range of considered frequencies.
The estimates , , and of the parameters of the accelerometer model and the associated uncertainties , , and determined by using the WLSP are tabulated in Table 2. Table 2. Estimated parameters and associated uncertainties obtained by using the weighted least squares procedure (WLSP).

[V/(m/s 2 )] [Hz]
[-] The expanded uncertainties, u(A) and u(Φ), are marked in Figure 6 by the solid vertical lines within the measurement points of both responses. The graphs of both responses are presented by using a logarithmic scale due to the wide range of considered frequencies.
The estimates S a , f 0 , and β of the parameters of the accelerometer model and the associated uncertainties u S a , u f 0 , and u β determined by using the WLSP are tabulated in Table 2. Table 2. Estimated parameters and associated uncertainties obtained by using the weighted least squares procedure (WLSP). The model responses shown in Figure 6 are obtained by relative approximation by using the functions A( f ) and ϕ( f ) for the parameters tabulated in Table 2, where Energies 2020, 13, 1552 10 of 13 and The estimates and associated uncertainties from Table 2 were utilized to determine the parameters S max a , f max 0 , and β max giving the maximum value of the reference energy determined by the FPA controlled by the MCP. The input parameters for the FPA were selected as: J = 30, T = 2 ms, ∆ s = 2 µs, and ξ 0 = 1, and the value of signal limit a = 1.037 V was assumed as equal to the estimate of the voltage sensitivity S a of the considered accelerometer. The time T was obtained from the steady state of the impulse response k(t) calculated based on Equation (14). This steady state was assumed to be equal to 5% of the maximum deviation (peak response) from zero of this impulse response. The tenth-order (L = 10) analog Butterworth filter was adopted and is described by Equation (6). The filter cut-off frequency f c is equal to the frequency response of the considered accelerometer and is 1 kHz. The values of both the gain factor a and the signal limit a are equal to 1.03 V (i.e., the value of the estimate S a of the accelerometer voltage sensitivity). The upper limit of the filter order results from the computational limitation of the MathCad program which was applied for calculation of the impulse response in the way of determining the inverse Laplace transform on the basis of the relevant transfer function K f (s). Figure 7 shows the impulse responses k(t), where the maximum deviations of these responses from the steady state value are equal to 1764. The number of MCts for the procedure controlling the FPA was 2 × 10 5 , similar to the WLSP. Based on parameters and uncertainties given in Table 2 The estimates and associated uncertainties from Table 2 were utilized to determine the  parameters  , , and giving the maximum value of the reference energy determined by the FPA controlled by the MCP. The input parameters for the FPA were selected as: = 30, = 2 ms, ∆ = 2 μs, and = 1, and the value of signal limit = 1.037 V was assumed as equal to the estimate of the voltage sensitivity of the considered accelerometer. The time was obtained from the steady state of the impulse response ( ) calculated based on Equation (14). This steady state was assumed to be equal to 5% of the maximum deviation (peak response) from zero of this impulse response. The tenth-order ( = 10) analog Butterworth filter was adopted and is described by Equation (6). The filter cut-off frequency is equal to the frequency response of the considered accelerometer and is 1 kHz. The values of both the gain factor and the signal limit are equal to 1.03 V (i.e., the value of the estimate of the accelerometer voltage sensitivity). The upper limit of the filter order results from the computational limitation of the MathCad program which was applied for calculation of the impulse response in the way of determining the inverse Laplace transform on the basis of the relevant transfer function ( ). Figure 7 shows the impulse responses ( ), where the maximum deviations of these responses from the steady state value are equal to 1764. The number of MCts for the procedure controlling the FPA was 2 × 10 , similar to the WLSP. Based on parameters and uncertainties given in Table 2  The relationship between the maximum reference energy ( ) and the time of the considered accelerometer is linear [5] for the steady state of the impulse response ( ) (i.e., for > 2 ms). Thus, it is possible to determine the functional relationship between this energy and any time higher than 2 ms. Hence, the value of energy ( ) should be determined for any time higher The relationship between the maximum reference energy E max (y) and the time T of the considered accelerometer is linear [5] for the steady state of the impulse response k(t) (i.e., for T > 2 ms). Thus, it is possible to determine the functional relationship between this energy and any time T higher than 2 ms. Hence, the value of energy E max (y) should be determined for any time higher than 2 ms, and then the relevant formula can be easily determined based on the linear regression. For example, if T = 10 ms, then we have E max (y) = 34.115 mV·m 2 and u(E max ) = 0.012 mV·m 2 , which are obtained for m = 93,515 and J = 9. Based on the two values of the energy E max (y) and the values of time T, we have: Energies 2020, 13, 1552 11 of 13 Figure 8 shows the signal x m (t) with two reversals producing the E max (y) at the output of the considered accelerometer obtained for T = 2 ms.   The reversals of the signal ( ) are 0.68 ms and 1.49 ms.

Conclusions
The paper presents a proposal for the application of the FPA controlled by the MCP to determine the maximum reference energy at the output of accelerometers. Block diagrams together with a detailed description of the proposed solutions are presented.
The reference energy of the accelerometer should be understood as its output energy decreased by the output energy of the special filter (modeled by the low-pass type) with the cut-off frequency corresponding to the upper limit of the accelerometer frequency response. The lower limit of this response is neglected in this paper due to its low value (0.2 Hz).
Additionally, the MCP was used for parametric identification of both the parameters associated with the transfer function of the accelerometer and the associated uncertainties. Then, the value of the maximum reference energy was determined along with the associated uncertainties at the output of the accelerometer of the type DJB A/1800/V, considered as an example of the accelerometer. This was performed by checking the full range of the parameter variability of the accelerometer model by value-associated uncertainties. This is possible to perform only by using the MCP based on the pseudo-random number generator with a uniform distribution. This necessity results from the fact that the energy above can be obtained for the values of parameters contained within the variability ranges and not, for example, at their ends or in their middles. Only in this way it is possible to accurately determine the maximum reference energy, which can be the basis for mutual comparison of different types of accelerometers. The presentation of the relevant procedures together with a practical example giving the possibility of such comparison is the main contribution of this paper.  The reversals of the signal x m (t) are 0.68 ms and 1.49 ms.

Conclusions
The paper presents a proposal for the application of the FPA controlled by the MCP to determine the maximum reference energy at the output of accelerometers. Block diagrams together with a detailed description of the proposed solutions are presented.
The reference energy of the accelerometer should be understood as its output energy decreased by the output energy of the special filter (modeled by the low-pass type) with the cut-off frequency corresponding to the upper limit of the accelerometer frequency response. The lower limit of this response is neglected in this paper due to its low value (0.2 Hz).
Additionally, the MCP was used for parametric identification of both the parameters associated with the transfer function of the accelerometer and the associated uncertainties. Then, the value of the maximum reference energy was determined along with the associated uncertainties at the output of the accelerometer of the type DJB A/1800/V, considered as an example of the accelerometer. This was performed by checking the full range of the parameter variability of the accelerometer model by value-associated uncertainties. This is possible to perform only by using the MCP based on the pseudo-random number generator with a uniform distribution. This necessity results from the fact that the energy above can be obtained for the values of parameters contained within the variability ranges and not, for example, at their ends or in their middles. Only in this way it is possible to accurately determine the maximum reference energy, which can be the basis for mutual comparison of different types of accelerometers. The presentation of the relevant procedures together with a practical example giving the possibility of such comparison is the main contribution of this paper.
Funding: The presented research was conducted at the Faculty of Electrical and Computer Engineering, Cracow University of Technology, and was financially supported by the Ministry of Science and Higher Education, Republic of Poland (grant no. E-31/2020).

Conflicts of Interest:
The authors declare no conflicts of interest.