Analytical Modelling and Optimization of a Piezoelectric Cantilever Energy Harvester with In-Span Attachment

In this paper, the location of masses and of a piezoelectric patch for energy harvesting reported onto a vibrating cantilever beam is studied and optimized. To this aim, a genetic algorithm is adapted and utilized to optimize the voltage amplitude generated by the piezoelectric patches by choosing attachment mass, attachment mass moment of inertia, attachment location, piezoelectric patch location and force location on the beam as parameters. While an analytical approach is proposed to evaluate the voltage amplitude, a multi-layer perceptron neural network is trained by the derived characteristic matrix to obtain an approximate function for natural frequencies based on the attachment parameters. The trained network is then used in the core of genetic algorithm to find the best optimization variables for any excitation frequency. Numerical simulation by COMSOL Multiphysics finite element software validates the calculated voltage by analytical approach. The optimization method successfully matches the natural frequency of the beam with the excitation frequency which therefore maximizes the output energy. On the other hand, the superiority of the optimized design over the conventional configuration in harvesting the energy at high frequency excitation is also approved.


Introduction
Energy harvesting by piezoelectric materials has attracted lots of interests during the last years due to its high power density and architectural simplicity [1]. Piezoelectric materials are generally used in the vibration based energy harvesting specially at small scales [2]. Among different types of structures which have been developed to scavenge the energy from the ambient vibrations sources [3][4][5][6][7], the conventional and still widely used configuration in this case is cantilever based piezoelectric energy harvester which is deeply investigated in several researches by Erturk et al. [8][9][10][11]. In order to maximize the harvested energy from a piezoelectric cantilever beam, previous works considered geometry and parameters optimization on the piezoelectric cantilever beam [12][13][14][15][16] by using topology optimization [17] or by using interval techniques [18], or added a tip attachment as a proof mass to match the fundamental natural frequency of the beam to the excitation frequency [19,20]. However, in some applications with high excitation frequencies like automotive engines, industrial machinery or micro dimension applications the excitation frequency can exceed the fundamental natural frequency of the beam. Then, the best location for the attachment to have maximum piezoelectric voltage may not be the tip of the beam due to vibrational mode shapes. For this case, Erturk et al. [21] investigated the higher modal energy harvesting with a unimorph beam without attachment while the possibility of having the attachment in-span of the beam can be considered. On the other hand, studies on energy harvesting from high frequency excitation with having the mass in-span of the cantilever beam received a restricted attention. Researches in this area did not consider modelling [22] or the piezoelectric electromechanical coupling effect [23]. Furthermore, in all of the researches mentioned above, attached mass on the beam is modelled as a lumped mass and the effects of mass moment of inertia on the harvested energy are not investigated.
In piezoelectric cantilever configuration, it is usual to have piezoelectric layers which cover the whole beam from clamped side to free end [9,10,[24][25][26]. This configuration is mostly convenient for the microscale dimensions. However, in the mesoscale dimensions this configuration suffers from low power density since the most strain occurs close to the clamped side [4]. Besides, in excitation frequency close to higher vibrational modes, charge cancellation may happens and continuous electrodes should be avoided [21]. Because of these reasons some researches prefer to have a piezoelectric patch mounted to a passive (non-piezoelectric) beam with the length of the piezoelectric patch being smaller than that of the beam [4,19,27,28]. Generally, in these cases, the piezoelectric patch are mounted near to clamped side of the beam to take the advantage of the maximum strain. However, for excitation frequencies higher than the first natural frequency, again the clamped side may not be the best place for the piezoelectric patch.
In this paper, for a general case of a given passive cantilever beam with mounted piezoelectric patch and several in-span attachments as shown in Figure 1, the voltage amplitude of the piezoelectric patch is found by analytical approach by extending the existing vibrational analysis in the literature [29][30][31][32]. In this vibrational analysis, the method of sectioning the beam between the attachments are utilized to find the response of the beam due to excitation. As such, the model proposed for the voltage amplitude is general for any number of attachment. On the other hand, in order to optimize the harvested energy from an external source of vibration, optimization is applied on a particular case of a beam with one in-span attachment. The derived equation for piezoelectric voltage is used as objective function of the optimization algorithm. In this case, the employed optimization variables are: attachment mass, attachment mass moment of inertia, attachment location, piezoelectric patch location and external force location. In order to deal with this multi parameter optimization, genetic algorithm (GA) is used. However the attachment parameters do not appear in the voltage equation while affecting the natural frequencies. To handle this, one can use the Rayleigh's quotient method to calculate the fundamental natural frequency of a cantilever beam with one in-span attachment analytically [33]. But, for higher natural frequencies the Rayleigh's quotient method is almost impossible to be used since there is no information about the shape functions of the beam for higher modes and using polynomial function as shape function [34] can produce significant error especially with in-span attachment on the beam. Therefore, here a two steps optimization is suggested.
First, a multi-layer perceptron (MLP) neural network is trained in order to obtain approximate functions for the natural frequencies based on the attachment parameters. Second, the trained network is used in the evaluation process of the genetic algorithm to form a neural network based genetic algorithm which can find the best combination of the optimization variables to match one of the natural frequencies of the beam to the excitation frequency and maximize the piezoelectric voltage. Both GA and MLP are already explored in the literature. However, utilization of MLP in the evaluation process of the GA has two advantages: first, it can increase the speed of GA algorithm significantly which is a known bottleneck of the GA algorithm in the literature. The second advantage is that MLP can have very high accuracy in calculation of the natural frequencies if the training data is sufficient. These advantages will be explained in details in the text.
To validate the analytical modelling and proposed optimization method of this paper, commercial finite element software COMSOL Multiphysics (COMSOL, Inc., Burlington, NJ, USA) is used.
The performance of optimization method in matching the natural frequency of the structure with the excitation frequency is also assessed.
The rest of the paper is organized as follows-in Section 2, vibrational analysis is performed and the piezoelectric voltage due to external excitation is found. In Section 3, first, neural network fitting approach and genetic algorithm are described separately, then their combination in maximizing the voltage amplitude based on the excitation frequency is presented. Section 4 is devoted to simulation results in which the piezoelectric voltage from the proposed approach is compared with conventional method in the literature. Finally, conclusions are discussed in Section 5.

Modelling
For a beam with several in-span attachments and a piezoelectric patch which is mounted on the beam as shown in Figure 1, the regular method for deriving the equation of motion in transverse vibration is to separate the beam between the attachments. In this case, the equations of motion for each section by considering electrical coupling is in the following form [9]: In (1), w n (x, t) is the transverse displacement, n is the number of each segment, EI is the flexural rigidity of the beam, µ is the mass per unit length, C is the damping constant, F(t) is the external force, x F is the application location of the external force and δ is the Dirac delta function. The coupling term ϑ which comes from the electrical circuit can be written as [9], where, EI p is the flexural rigidity of the piezoelectric patch, d 31 is the piezoelectric coupling coefficient, h p and b p are piezoelectric thickness and width respectively while h b is the beam thickness. h c is the distance from the top of piezoelectric layer to the neutral axis. It should be noted that the term specifies the end points of the piezoelectric patch on the beam; and for the sections of the beam where there is no piezoelectric patch, this term is equal to zero. This term is due to the fact that induced voltage in the piezoelectric patch generates point moments at its boundaries which affects the mechanical response of the beam. Now, to find the piezoelectric voltage due to external mechanical excitation, the response of partial differential Equation (1) should be found. It is worthwhile to mention that the thickness and length of the piezoelectric patch is much less than the beam. As such the rigidity of the piezoelectric patch is negligible in comparison to the beam. In addition, the effects of the attachment's geometry on the beam's structural dynamic are highly superior to the effects of the piezoelectric patch geometry. Therefore, in the following modal analysis of the beam, the geometry of the piezoelectric patch is neglected in comparison to the beam and attachment.

Free Vibration Analysis
By following the Ansatz separation, solution of (1) is expressed in the following form: ϕ n (x) is the mode shape of the beam and q(t) is the generalized time-dependent coordinate which satisfies In (4), ω is the natural frequency of the system. In order to find the natural frequencies and mode shapes of the beam with several in-span attachments, the external force and damping constant in (1) are considered to be zero. As such, the equations of motion for each segment of the beam in short circuit condition (v(t) = 0) is converted to the following eigenvalue problem [35,36]: Now based on (5), the exact solution for the eigenfunctions which are the mode shapes of the beam, is in the following form [35,36]: ϕ n (x) = A n sin βx + B n cos βx + C n sinh βx + D n cosh βx, (6) in which, In (6), A n , B n , C n and D n are unknown constants that can be found by applying proper boundary and continuity equations. The clamped end boundary condition is written in the following form: For the attachment point of each segment, the continuity of deformations and equilibrium of forces and moments can be written as M n is the mass of the attachment and J n is the related moment of inertia of the attachment. Finally, for the free end with attachment, boundary conditions are By applying (8)-(10) to (6), the following characteristic equation is formed: whereJ 4N×4N is the characteristic matrix which is the function of β. On the other hand, β is a function of ω based on the explicit expression mentioned in (7). In (11), P 4N×1 is the matrix of the mode shapes constants, To find the nontrivial solution of (11), the determinant of the characteristic matrixJ 4N×4N should be zero which forms the frequency equation and numerical approach should be used to find the natural frequencies. On the other hand, by setting the determinant of characteristic matrix to zero, characteristic Equation (11) will become undetermined. In this case, by integrating the normalization condition, the set of equations become solvable and the coefficients in vector P 4N×1 are found. In this method, orthogonality condition can be used to find the perfect normalization. The generalized orthogonality condition which results in a decoupled ODE of motion for a beam with in-span attachments is expressed in the following form [31,32]: This orthogonality condition has a companion form as, In (13) and (14), (r) and (s) are the numbers of mode. Now, based on (13), the normalization condition can be written for the case when (r = s) and in this case the right hand side of (13) is equal to 1. After integrating the normalization condition to the set of (8)-(10), the coefficients in vector P 4N×1 in (12) are found and the mode shapes can be found for each natural frequencies. Then, the displacement for each point of the beam with considering the beam segments can be found by using the expansion theorem in the following form:

Forced Vibration Analysis
By substituting (15) in (1) and multiplying each side of the equation by ϕ (s) n and integrating over the length of the section, the following equation is derived: Equation (16) is true for each section of the beam. By summing these equations for all sections in the entire length of the beam, the following familiar form of ordinary differential equation (ODE) of motion is derived:M∆ +C∆ +K∆ =F −V, (17) in which,M In (17),M rs is the element of the mass matrix,C rs is the element of the damping matrix,K rs is the element of the stiffness matrix,F r is the modal vector's element of external forces,V is the vector of the electromechanical coupling voltage, ∆ is the vector of the time dependent coordinates and H is the Heaviside function. Now, by applying the orthogonality conditions and normalization in (13) and (14), the mass and stiffness matrices can be rewritten in the following form: where δ rs is the Kronecker delta symbol. Based on (25), mass and stiffness matrices are diagonal while the damping matrix is not. For solving the equation of motion in (17) analytically, it is better to have a set of decoupled ODE which is impossible with non-diagonal damping matrix. However, the damping of a system is a model parameter and it is not a physical parameter [31]. Therefore, it is possible to assume a modal damping matrix for the system in the following form: where ζ (r) is the modal damping constant. It is possible to define the numerical value for this modal damping constant by using the Rayleigh damping theory which is common in the Finite Element Method (FEM) software. The Rayleigh damping can be defined as where α, β are the Rayleigh's damping coefficients. By substituting (25) and (26) in to (27), the modal damping can be calculated as Now, by considering diagonal damping matrix, the decoupled set of ODE of motion can be written in the following form: (29) is found by the electrical circuit equation with mechanical coupling as follows [9]: and In (30), R l is the resistive load, L p is the piezoelectric length, ε S 33 is the piezoelectric permittivity constant at constant strain by assuming plan-stress assumptions and h pc is the distance of the piezoelectric patch center in thickness direction to the neutral axis. In order to solve (30) and (29) first, external harmonic force is modelled by F(t) = F 0 e jΩt , where Ω is the excitation frequency. By considering a linear electromechanical system, the response of v(t) is also harmonic in terms of v(t) = V 0 e jΩt while V 0 and F 0 are the amplitudes of the voltage and external force respectively. Now, the particular solution of non-homogeneous differential equation in (29) is: By substituting (32) in (30) the following equation for voltage amplitude is derived: in which, τ c is the time constant of the electrical circuit. The voltage equation in (33) is similar to one which is reported by Erturk et al. [9]. However, the difference lies on calculating the natural frequencies (ω (r) ) and mode shapes (ϕ (r) ) which are based on modal analysis of a beam with in-span attachments. After finding the amplitude of voltage, the displacement of each point of the beam with the help of (15) and (32) can be found by the following equation: As can be seen in (33), voltage amplitude is a complex number which has an absolute value and a phase angle. In the framework of this paper, just the absolute value is important for optimization.

Optimization
It is desired to have the maximum possible of piezoelectric voltage amplitude in (33) for any excitation frequency. The geometrical optimization on the beam and piezoelectric has been studied before [12,37]. Therefore, by considering constant geometrical parameters for beam and piezoelectric patch and considering just one attachment on the beam, the optimization variables that can influence the voltage amplitude in (33) are piezoelectric patch location (l p1 ), force location (x F ), attachment location (l 1 ), attachment mass (M) and attachment mass moment of inertia (J). Among these optimization variables, force location and piezoelectric patch location directly appear in the voltage amplitude (33). However, the other three optimization variables related to attachment do not directly appear in the voltage amplitude. Instead, they affect the continuity conditions in (9), and only these latter affect the characteristic matrix in (11) which itself affects the natural frequencies and mode shapes that appear in the voltage amplitude (33). On the other hand, there is no closed form expression between the natural frequencies and the attachment parameters. In fact, there is just closed form expression for first (fundamental) natural frequency of the beam with in-span attachment based on the Rayleigh's quotient method [33]. But, for higher natural frequencies the Rayleigh's quotient method cannot be used since there is no information about the shape functions without knowing the natural frequencies. As such, performing optimization based on the attachment parameters is a challenge that will be addressed in the next section.

Multi Layer Perceptron Neural Network
In this section, neural network fitting algorithm in MATLAB software is utilized for the purpose of finding approximate functions for the natural frequencies based on the optimization variables. The algorithm is based on the multi-layer perceptron (MLP) which is a class of artificial neural network. MLP consists three layers of artificial neurons-input layer, output layer and hidden layers. There is just one input layer and one output layer in the MLP network and the number of neurons in input and output layers are equal to the number of input and output variables. However, there can be different numbers of hidden layers in MLP and the main method to increase the performance of the MLP is to define the adequate number of hidden layers which can be done with the trial and error approach.
Having this background, it is now desired to have a MLP network that can get attachment parameters and gives the natural frequencies. To do so, the MLP network should be trained with the training data which consist of input data and target data. Input data are optimization variables regarding the attachment including attachment location, attachment mass and attachment mass moment of inertia. Target data consist of natural frequencies related to the input data. In this paper maximum number of natural frequencies and mode shapes for modal analysis is considered to be 4 which is sufficient enough to model a cantilever even under high excitation frequencies. Therefore, each set of input data consists of 3 optimization variables, which has a set of target data consist of 4 natural frequencies.
The MLP network should have the ability to get any combination of the attachment parameters and give the related natural frequencies. To form the training data for this MLP network, 20 different attachment masses in the domain of (0 < M < µL), 20 different attachment mass moments of inertia in the domain of (0 < J < 0.0042 µL 3 ) and 20 different attachment locations in the domain of (0 < l 1 < L) have been chosen to form 8000 combinations of attachment parameters. For each of these combinations, there are 4 natural frequencies which should be found by the numerical approach on the characteristic matrix in (11). Therefore, input training data is a matrix with 3 rows and 8000 columns while the target data is a matrix with 4 rows and 8000 columns. Number of hidden layers is 50 and the Bayesian Regularization method has been used to train the network.
After successful training of the network by obtaining natural frequencies for finite set of variables, a continuous function is approximated as a black box which can get any variable in the predefined domain and gives the related natural frequencies immediately as shown in Figure 2 and it is not necessary to calculate the natural frequencies by numerical approach on the characteristic matrix in Equation (11). This will boost the speed in the evaluation procedure of the GA optimization which will be discussed later.

Finding Analytical Expression for Mode Shape Constants
By using MLP neural network, the natural frequencies in (33) are found as functions of the attachment parameters. But, in voltage amplitude (33), mode shapes are also functions of the attachment parameters since attachment parameters affect the continuity (9) and mode shapes constants in (12). Finding analytical expressions for the mode shape constants based on (8)-(10) are almost impossible due to huge size of the expression for each of the constants. Alternative approach proposed by Naguleswaran [36,38] is used here to decrease the number of the unknown constants in (12).
For the case when the beam has just one in-span attachment, some algebraic simplification can be done on the mode shapes by defining two coordinate systems for two sections of the beam as shown in Figure 3. The reference of the second coordinate system has been placed at the tip of the beam and the points in the second coordinate system are shown by "˜".
By applying the clamped end boundary conditions mentioned in (8) to Section 1 of the beam, the following result related to the mode shape constant of the first section is obtained: Now, by considering the second coordinate system for the second section of the beam and considering no tip mass for the beam, the free boundary condition can be written for the second section of the beam in the following form: By applying (38) to the mode shape of the second section of the beam, (39) is obtained: Using (37) and (39), mode shapes for Segment 1 and 2 can be written as By defining two coordinate systems, there are 4 unknown constants instead of 8 unknown constant. These 4 unknown constants can be found by rewriting (9) for the new coordinate system, Now, to find each constant, first B 2 is considered to be 1. Then, the first three conditions of (41) are solved analytically with the help of MAPLE software to find the remaining three constants as functions of the natural frequencies and attachment parameters. The code which is written in MAPLE software is mentioned in the Appendix A. On the other hand, in this case there is no normalization in the mode shapes and the mass matrix in (17) is not an identity matrix while it is still diagonal. Therefore, (29) without normalization will be in this form: M rrq (r) (t) +M rr 2ζ (r) ω (r)q(r) +M rr (ω (r) ) 2 q (r) (t) = By considering (42) instead of (29), the mass matrixM rr enters the piezoelectric voltage amplitude (33) and this latter is rewritten in the following form: The value ofM rr is found analytically as function of the mode shape constants and attachment specifications by Maple software and it is reported in Appendix A.
In the approach of the last two sections, natural frequencies and mode shapes have been found as functions of the attachment specifications.

Genetic Algorithm Optimization
Now, it is desired to find the optimal voltage amplitude from (43) by finding the optimal position of the attachment, the attachment mass, mass moment of inertia, external force location, and piezoelectric patch location. As such, genetic algorithm is utilized in this paper, which has a great performance in dealing with multi variable optimization problems. Genetic algorithm is an optimization method that works based on the evolution theory and it is completely a simulation of real life in nature. By choosing the voltage amplitude (43) as fitness function to be optimized, the mathematical framework of GA optimization method works in the following steps: Pairs of individuals are selected as parents with Roulette wheel method. The ratio of individuals selected as parents to the overall individuals is 0.8.

•
Two selected parents give birth to two offspring (two new possible solution) in the crossover procedure.

•
Then chromosomes of some offspring will be changed in mutation procedures. The chromosomes for the optimization problem are the optimization variables. The percentage of offspring who experienced the mutation to the whole population is called mutation percentage which is considered here to be 10 percent.

•
The last step is the evaluation procedure in which individuals with lowest fitness value will be replaced by offspring. These newly born offspring with remaining individuals form the next generation. To calculate the fitness value of each individual, MLP network and analytical expression for the mode shape constants should be used which will be explained in the next section.

•
Newly formed generation will undergo the same procedure of the previous generation.
MATLAB GA toolbox performs the aforementioned steps and stops after a restricted number of times.

Neural Network Based Genetic Algorithm (NN-GA)
In the evaluation step of genetic algorithm, the fitness value of each individual should be determined. As has been shown in Figure 4, each individual has 5 chromosomes in which 3 of them are related to the attachment. MLP neural network calculates the natural frequencies of the system based on the attachment parameters. Then, mode shapes are calculated by the analytical expressions for the mode shape constants. Finally, with the help of mode shapes, natural frequencies, force location and piezoelectric patch location, the fitness value of each individual are calculated by using the voltage amplitude (43).
Using MLP neural network in the evaluation process of GA will increase optimization speed significantly by eliminating the cumbersome numerical procedures to find the natural frequencies related to each individual. Therefore, it is possible to choose the GA parameters (population number, mutation percentage, ...) different to those mentioned above to compare the results and to avoid trapping in the local optimums.

Simulation and Results
In this section, first the effect of each optimization variables on the voltage amplitude is investigated solely. It should be noted that the effects of the electrical circuit parameters on the voltage is not investigated here. Therefore, for constant electrical circuit parameters particularly the resistance, the optimization of voltage is equivalent to optimization of power. Investigations are performed on a cantilever beam with an in-span attachment which has a cubic shape with length (L M ), thickness (H M ) and width (W M ). All the features and dimensions of the system are reported in Table 1. In next part of results, the neural network based genetic algorithm is used to find the best optimization variables based on the excitation frequency of the external force. The natural frequencies of the optimized structures are calculated to asses the performance of the optimization method in matching one of the natural frequencies of the structure with the excitation frequency. Then, by simulating the optimized structure in the COMSOL FEM platform, the natural frequencies and voltages which are derived by analytical method are compared with ones obtained by FEM simulation in COMSOL. To investigate the performance of the optimization method, the voltages of the optimized structures are also compared with the classical configuration in which a lumped mass is attached at the tip of the beam and a piezoelectric patch mounted on clamped side of the beam and the results are discussed at the end.

Discussion
In Figure 5, the effects of patch location vs. the attachment location have been illustrated on the voltage amplitude. As can be seen in this figure, it is obvious that when the excitation frequency is near to first natural frequency, the best piezoelectric patch location and attachment location are near to the clamped end. The reason for piezoelectric patch location is that in first vibrational mode shape it has the most possible strain close to clamp side of the beam. The reason for optimum attachment location is that the force is also applying at the end of the beam and in the next figure it will become clear that having the attachment and force at the same point of the beam should be avoided. When the excitation frequency approaches to the second natural frequency, then the best placement for the piezoelectric patch and attachment is near to the midpoint of the beam. On the other hand, for the excitation frequency near to higher natural frequencies, there are many local optimums for the mass and piezoelectric patch location. . In Figure 6, the effects of attachment mass vs. force location on the voltage amplitudes are shown. Based on this figure, attachment mass does not have optimums for voltages. It means that tuning the mass for matching the natural frequency with the excitation frequency is enough. On the other hand, for the excitation frequencies close to first natural frequencies the best force location is at the tip of the beam while for excitation frequencies close to the higher natural frequencies the optimum force location is close to the mode shape apexes. Another important point is that for the excitation frequencies close to higher natural frequencies, when the force and attachment are at the tip of the beam, increasing the mass decreases the voltage amplitude. As such, when the attachment is at the tip of the beam the optimum location for force is not at the tip of the beam. In Figure 7, the effects of attachment mass moment of inertia and force location on the voltage amplitude are shown. When the excitation frequency is equal to the first natural frequency, there is no optimum for the mass moment of inertia. When the excitation frequency reaches to second natural frequency, the behavior of voltage amplitude is highly dependent on the force location. On the contrary, for third and fourth natural frequency it can be deduced that increasing the mass moment of inertia will increase the voltage generally. It will be shown in the next section that mass moment of inertia will increase the bandwidth of accessible natural frequencies to match with the excitation frequency.

MLP Neural Network Fitting Function
In Figure 8, all possible accessible range of natural frequencies for the classical setup of cantilever with tip attachment and for the proposed setup of this paper is shown. Figure 8a illustrates the natural frequencies when the mass moment of inertia is zero and the attachment is at the tip of the beam. Figure 8b illustrates the natural frequencies for all 8000 combination of mass, mass moment of inertia and location. As can be seen, by considering in-span attachment and mass moment of inertia, a higher frequency bandwidth becomes accessible for the system of cantilever beam and attachment. For example, when excitation frequency is equal to 100 (Hz), 300 (Hz) or 600 (Hz), without considering in-span attachment and its mass moment of inertia, it is impossible to match the natural frequency of the system with the excitation frequency. On the other hand by increasing the mass moment of inertia these frequencies become accessible as it is shown in Figure 8b. Based on Figure 8, for some excitation frequencies such as 20 (Hz), 100 (Hz), 140 (Hz), 300 (Hz), 370 (Hz), 600 (Hz) and 750 (Hz), there are infinite combinations that match the natural frequencies of the system to the excitation frequencies. Now an interesting question is: which of them produce more voltage amplitude? Furthermore, some other excitation frequencies such as 50 (Hz) and 220 (Hz) are impossible to match with the natural frequencies of the system even by considering the mass moments of inertia. So, in these cases is it better to eliminate the attachment? These are the questions which will be addressed in the optimization section.

Optimization
The trained MLP neural network is a continuous function which gets 3 inputs and gives 4 outputs. Therefore, plotting this network based on the inputs and outputs in a 3D plot is impossible. But, the mean squared error after 1000 training epochs reached to 10 0 which is completely satisfying since the scale of the natural frequencies as the target data are close to 10 2 (rad/s) for the first natural frequency and 5 × 10 3 (rad/s) for the fourth natural frequency. To check the fitting performance of the network, it is possible to consider one of the inputs constant and plot the results for variations of the two other inputs. This is done in Figure 9, in which the four natural frequencies of the beam are illustrated based on the variation of the attachment location and attachment mass. It can be seen that the data obtained from the trained network fits the discrete training data with high accuracy. The discrete data plots for the first and second natural frequency are similar to plots obtained by Low [29]. With a combination of an MLP neural network and genetic algorithm, now we have an optimization algorithm which gets the excitation frequency and gives the optimum possible optimization variables to maximize the voltage of the piezoelectric patch due to external excitation. Now, by setting the desired excitation frequency, NN based GA will choose the optimization variables based on the properties defined in Table 1. For each frequency specified in Figure 9, the best and mean fitness values in the generations are illustrated in Figure 10. It should be noted that the fitness values are obtained in the function evaluation step of GA for the output voltage multiplied to a negative sign to convert the maximization to a minimization problem. Based on Figure 10, it is obvious that less than 200 generations was enough for convergence. In Table 2, for different excitation frequency, the optimization variables which have been chosen by NN based GA are reported. Each optimized set of variables named as config (excitation frequency). In Table 3, the natural frequencies for optimized structures are calculated by analytical method and FEM modal analysis by COMSOL. To model the beam and in-span attachment, 3D environment of COMSOL is used. It should be noted that the cubic attachment and beam should not be in complete contact. Otherwise, the continuity of the beam will be affected by the rigidity and thickness of the attachment. Whereas, in continuity Equation (9), only the attachment mass and mass moment of inertia is taking in to consideration. To remedy a very small rigid cylindrical connector is modelled in COMSOL to mount the attachment on the beam. Furthermore, the same modal damping of analytical modelling is chosen in the COMSOL 3D environment.
As can be seen in Table 3, the optimization method successfully matched one of the natural frequencies of the structure with the excitation frequency, when the excitation frequency is in the accessible frequency bandwidth of the structure. It is worthwhile to mention that when the natural frequencies of the structure are calculated by analytical Euler-Bernoulli beam theory, one of the natural frequencies of the structure is exactly matches the excitation frequency with 10 −1 percent error. However, when the natural frequencies of the structures are calculated by the FEM analysis, the error will be increased to 10 0 percent. This increase in error comes from the fact that in 3D FEM analysis of the beam specially with attachment, there are other degrees of freedom like chord wise bending or torsional degree of freedom that are coupled with the transverse bending of the structure. This will deviate the natural frequencies of the beam in 3D FEM analysis from those of Euler-Bernoulli beam theory. To compare the piezoelectric voltages of the optimized design and classical configuration, the frequency response plot for each configuration around its excitation frequency are illustrated in Figure 11.  Table 2 Several points can be deduced by the frequency response plot in Figure 11. First, it is obvious that the piezoelectric voltages derived by analytical method and FEM analysis are in good agreement. In fact, in COMSOL the Rayleigh damping is defined with same coefficient which is very important in terms of verification of the results. As can be seen in the figure, the amplitude of the voltage even at the peaks are in good agreement with the analytical calculations.
The next important point that can be seen in the frequency response plot is that, the voltage of the optimized design is highly superior to the classical configuration even for the cases when the excitation frequency is in the accessible bandwidth of the classical configuration(like 140 HZ, 370 HZ and 750 HZ). For other frequencies like 100 Hz, 300 Hz and 600 Hz which is outside the accessible bandwidth of the classical configuration, obtained voltages for optimized configurations are extremely superior to the classical configuration thanks to consideration of in-span attachment and its mass moment of inertia, which can convert the excitation frequency to resonance frequency.

Discussion
By using the evolutionary optimization algorithms like GA, it is always possible to be in the local optimums. To reduce the possibility of local optimums, GA can be applied several times with different parameters (population of individuals, mutation rate, etc.). In general, this procedure is very time consuming. However, by introducing the MPL in the core of evaluation process, the GA is performing fast enough to be applied several times with different parameters. On the other hand, even without applying the GA for several times, still the improvement of obtained results in comparison to the classical approach of tip attached cantilever beam is significant.
The optimization methodology is applied on a particular case of a beam with one in-span attachment. However, with the proposed analytical approach for the calculation of the piezoelectric voltage, the optimization methodology can be applied on a beam with several in-span attachment as well. On the other hand in this case more training data will be required for the MLP neural network.
Finally, while the optimization methodology proposed in this paper is for piezoelectric energy harvesting scope, it could be worth to explore its application to piezoelectric sensing scope. Most of the existing miniaturized piezoelectric sensors studies are focused on on fabrication technology and on their integration, but also on their combination with piezoelectric actuation in order to form one single and same structure for both actuation-sensing (named as self-sensing) [39][40][41]. Perspective works could therefore explore the use of optimization methodologies, among which the one in this paper, to design sensing or self-sensing structures with piezoelectric elements.

Conclusions
In this paper, a neural network based genetic algorithm has been proposed to maximize the voltage amplitude of a piezoelectric patch mounted on a cantilever beam with in-span attachment. The effects of attachment location, attachment mass, attachment mass moment of inertia, piezoelectric patch location and force location on the beam have been illustrated on the piezoelectric voltage amplitude. It has been shown that when the excitation frequency of a beam exceed the fundamental natural frequency, conventional configuration of cantilever beam with tip attachment and clamp side mounted piezoelectric patch does not provide the optimum voltage amplitude. The optimization results demonstrate that for excitation frequencies higher than fundamental natural frequency the piezoelectric voltage amplitude with optimization variables suggested by the proposed NN-GA has significant superiority over the conventional configuration.
Future works would extend the proposed approach to 2D structures like plates or circular diaphragms to optimize the harvested energy from the harmonic pressure. The work will also focus on the fabrication and realization of the optimized structures. Future works also include the fabrication of one optimized design from the proposed technique and perform extensive experimental tests on it. Funding: This work has been supported by the national CODE-Track project (ANR-17-CE05-0014-01, Control theory tools for optimal design of piezoelectric energy harvesters devoted to birds tracking devices). This work has also been partially supported by the Bourgogne Franche-Comté region project COMPACT.