Nonlinear dynamic analysis of asymmetric bistable energy harvesters

Nonlinear vibration energy harvesting systems can potentially increase the power collected from the kinetic energy available in their operating environment since they usually can recover energy in broadband frequencies compared to their linear counterpart. However, these systems have a high degree of complexity, sensitivity to slight variations of the parameters and the initial conditions, and may present multiple solutions. For these reasons, it is rare for the designer to have a deep understanding of the dynamic behavior of this type of nonlinear oscillator. This situation is even more peculiar when geometric imperfections from the system’s manufacturing process are present, as they can signiﬁcantly inﬂuence the energy recovery process. In-tending to ﬁll this lack of understanding about general aspects of the nonlinear dynamics of these kind of systems, the present paper presents a broad numerical investigation of local and global characteristics of the underlying dynamical systems using bifurcation di-J

agrams and basins of attraction.Bifurcation analysis is performed by exploring the broad spectrum of a harmonic signal, going from low to high amplitude and frequency of excitation.Basins of attraction analysis based on 0-1 test for chaos is proposed as an efficient statistical technique to identify chaotic and periodic solutions.Different levels of asymmetry are investigated, and a particular situation is defined and analyzed when a value of the sloping angle where the system is attached compensates for the asymmetry of the quadratic term.The result shows the different solutions defined by excitation forces and initial conditions, indicating the best scenario for increasing the power output.The adverse effects of the asymmetries are presented.However, we also demonstrated that it is possible to around this behavior using the sloping angle to compensate for the asymmetric influence, obtaining good performance as a symmetric system.

Introduction
Overusing electronic mobile devices that depend on batteries requires techniques to provide greater autonomy.Energy harvesting techniques can be an essential step in this technological advance.They consist in producing electrical energy from available secondary energy sources, such as vibration, heat, wind, and radiofrequency, which are often neglected.The applications of this technology range from large-scale energy generation to small device batteries, such as embedded electronic sensors in the Internet of Things (IoT) appli-cations, microelectromechanical systems (MEMS), and nanoelectromechanical systems (NEMS), etc. [18,31].
This technology proved to be promising, especially in low energy demand applications.Because of the immense possibilities for development, these systems have been studied in several recent works.In Abdelkareem et al. [1], a regenerative energy harvest from automotive suspension is presented, while Phillips [29] investigates the topic of nanoscale sensors.Karami and Inman [14] proposed energy harvesters from heart vibration to power the pacemaker, and Catacuzzeno et al. [4] proposed electrical energy harvested from cell membrane potential to power wireless communication.Novel physical layouts of energy harvesting devices are proposed and discussed in [15,41,40].Other recent high-impact works in the scientific community can also be seen in [23,33].
The present work deals with a nonlinear bistable piezo-magnetic-elastic energy harvester presented by Erturk et al. [9] to convert vibration energy into electricity.Linear energy harvesting devices, despite their dynamic simplicity, typically provide reduced output power in regions away from the resonance, limiting their application possibilities.On the other hand, nonlinear systems, as first demonstrated by Cottone et al. [6] and Erturk et al. [9] can potentially increase the energy harvested over a wide frequency band.Several works present an in-depth investigation of the dynamics of this system, mainly to find improvements in energy efficiency.Daqaq et al. [8] presented a review of the role of nonlinearities in the transduction of energy harvesters and their benefits for different types of excitations, conditions, and potential functions.Mann and Owens [22] investigated the presence of co-existing solutions and a broadening in the frequency response for some levels of base excitation.While Tan et al. [32] indicated an optimal performance analyzing the coefficients of the linear and cubic elastic external forces.The response to different forms of excitation is studied in [16], in particular, harmonic noise and stochastic white noise.Erturk and Inman [10] investigated the high-energy orbits over a range of excitation frequencies.Costa et al. [5] investigated the influence of cubic polynomial nonlinearities subjected to harmonic excitation.The nonlinear dynamics of this system are studied in [21,28] to determine the influence of the external force parameters on the system response.While [19] studied the stationary Gaussian white noise influence.Other works investigated this nonlinear vibration energy harvesting from different points of view, using nonlinear piezoelectric coupling [34], rectifies circuit [2], time-frequency analysis [35] and global sensitivity analysis [24,25].Some studies proposed an enhancement of this system [7,30] and others explored systems with trior multi-stable configurations [17,20,42].
Although numerous efforts have been dedicated to bistable energy harvesters, most research focuses on symmetric systems.However, constructing a perfectly symmetric system is arduous or almost impossible.Several environmental conditions and manufacturing processes yield asymmetries.The uncertainty of the geometry and material properties are examples of asymmetry source that, despite the accuracy of the current manufacturing machines, it is unavoidable.The operation can also lead to asymmetries, for instance, interactions with mechanical loads, magnetic fields, etc.In the bistable devices, the not identical magnets and their positions are examples of this feature.The authors in [12,13] characterized the asymmetric potentials as the nonlinear restoring force using a quadratic term.Wang et al. [36] investigated the multiple solutions of nonlinear asymmetric potential bistable energy harvesters under harmonic excitations, demonstrating a negative influence on output performance.Wang et al. [37] studied asymmetric potentials on the response probability and electrical outputs of bistable energy harvester driven by Gaussian white noise.The self-gravity force caused by the sloping angle of the plane, where the system is attached, is another source of asymmetry.This condition is often found in many engineering applications.The authors in [3,39] studied this asymmetry effect when the energy harvesting device is applied to human motion that the leg motion caused a swing over a certain angle.While Wang et al. [38] employed a sloping to enhance the harvesting performance by compensating the asymmetric potentials.
In order to implement the described applications, the underlying dynamical system must be studied in detail.Analyzing the system's responses according to its forcing configurations and initial conditions is an essential step in its characterization and the design of more efficient models.However, to the best of the authors' knowledge, no comprehensive analysis is available in the literature.It is possible to find only local assessments or configurations of interest delimited from the system dynamics.Thus, the present work differs from past investigations by offering comprehensive analyses of the broad spectrum and dynamic behaviors, also taking into account a more sophisticated model with asymmetries.
The objective of this work is to explore the system's dynamics under low values of amplitude and excitation frequency to high ones and characterize the dynamic behavior of each one.The analysis is performed from a time series of displacement, velocity, and voltage.Bifurcation diagrams are calculated to examine the dynamic behavior at wide observation windows.In addition, the basins of attraction are presented with the corresponding attractors to investigate the incidence of chaos and multiple solutions concerning the initial conditions.In this paper, we proposed the basins of attraction calculations by statistical analysis of the electrical voltage using 0-1 test for chaos.This tool has a lower computational cost than the Lyapunov exponent method because it is unnecessary to solve computationally expensive eigenvalue problems.On the other side, it solves computationally simple correlation functions.Therefore, we conducted a comprehensive study for symmetric and asymmetric bistable energy harvesting, exploring different amplitude and frequency of excitation.Besides, a particular condition is established and analyzed when the asymmetries cancel themselves.
The rest of this manuscript is organized as follows.Section 2 presents the symmetric and asymmetric dynamical system of interest.Bifurcation analysis by amplitude and frequency of excitation of these systems are presented in Section 3. In Section 4, the basins of attraction, based on 0-1 test for chaos, are conducted to analyze the multiple solutions for the initial conditions.Finally, Section 5 presents the main conclusions.

Nonlinear dynamical systems
The dynamical system of interest in this work, illustrated in Fig. 1, is the bistable piezo-magneto-elastic energy harvester adapted from [24].This electromechanical system comprises a clamped-free ferromagnetic elastic beam in a vertical configuration.The beam is attached to a rigid base where a pair of magnets is placed on the lower part.Two piezoelectric laminae are placed on the beam's highest part and connected to a resistive circuit.The rigid base is periodically excited by an external source, which, together with the magnetic field generated by identical magnets, induces large amplitude vibrations.Once this motion is perceived by the piezoelectric laminae, the mechanical energy of vibration is converted into electrical power that is dissipated in the resistor.This description describes the symmetric bistable energy harvester at Fig. 1a.The asymmetric bistable energy harvester is presented in Fig. 1b.We introduced two features that distinguish the symmetric model: (i) the harvester is attached to a plane with a sloping angle ϕ; (ii) the magnets are not identical, inducing asymmetric potential energy.These features allow obtaining a sophisticated model, which can assume asymmetries arising from the manufacturing and operation process.
The symmetric bistable energy harvester studied is presented by Erturk, Hoffmann, and Inman [9].The  behavior of this dynamical system evolves according to the following dimensionless initial value problem where t denotes time; x is the modal amplitude of oscillation; v is the voltage in the resistor; ξ is the damping ratio; f is the rigid base oscillation amplitude; Ω is the external excitation frequency; λ is a reciprocal time constant; the piezoelectric coupling terms are represented by χ, in the mechanical equation, and by κ in the electrical one; x 0 represents the modal initial position; ẋ0 is the modal initial velocity; v 0 denotes the initial voltage over the resistor.The upper dot is an abbreviation for time derivative.These parameters are all dimensionless.Cao et al. [3] firstly proposed the asymmetric bistable energy harvester.They introduced a quadratic term on the nonlinear restoring force to model the asymmetry from the magnets' eccentricity as described by works [12,13].The sloping plane's effect is modeled as the harvester's self-gravity force.So, the following initial value problem describes the governing lumpedparameter equation of motion for the dimensionless asymmetric bistable energy harvesting model where δ is a coefficient of the quadratic nonlinearity, p is the equivalent dimensionless gravity of ferromagnetic beam, and ϕ is the sloping angle.
The power output at time t is given by

Optimal sloping angle
The asymmetry can be a harmful condition for the energy harvesting process, as pointed out by the authors in [36,37].In order to avoid this undesired configuration, an optimal sloping angle value is calculated.This optimal value is used to compensate for the asymmetry in the potential functions [38].Assuming that the derivative of nonlinear restoring force of the model, , is zero, such that its stable solutions are given by To compensate the asymmetry condition and make the system symmetric, the sum of the nonlinear restoring force for the stable solutions of Eq.( 8) need to be zero, F r (x 1 ) + F r (x 2 ) = 0.In this way, the optimal sloping angle value is obtained as 3 Bifurcation analysis Numerical simulations are carried out considering the symmetric and asymmetric bistable energy harvester, assuming both increasing and decreasing amplitude and frequency values sampling ordering, here referred like forward and backward.Sampling analysis intervals for forcing parameters are set as 0.1 ≤ Ω ≤ 1.4 for the frequency and 0.01 ≤ f ≤ 0.3 for the amplitude.The remaining model parameters are assumed to be ξ = 0.01, χ = 0.05, λ = 0.05, κ = 0.5, p = 0.59 and δ = 0.15 for initial conditions as (x 0 , ẋ0 , υ 0 ) = (1, 0, 0).The parameter ϕ is set to 35 • to analyze a strong asymmetry effect.The optimal value obtained by the Eq.( 9) are also analyzed, which for the δ mentioned previously, Forward and backward bifurcation diagrams are, respectively, the sweeping up and down the parameter of interest, and they are typically depicted on a cool and warm colors scale.Figure 2 presents a brief schematic of the methodology of the bifurcation diagram.The diagrams are built considering the system steady-state condition, here defined as the last 10% of the response time series.For each excitation amplitude or frequency condition, sampled from a set of 1200 regularly spaced values, taken from the analysis intervals, the initial conditions problem are updated with the last system state of the previous solution, considering 1000 forcing cycles.
To explore the excitation amplitude as the parameter of interest, the frequency observation window is set from Ω = 0.1 to Ω = 0.9, regularly increased of 0.1, while the excitation amplitude, on the frequency as the parameter of interest, is set from f = 0.019 to f = 0.275, regularly increased of 0.032.Nonlinear dynamics effects are investigated.The symmetric model is compared with the asymmetric models.The power forward and backward bifurcation diagrams are presented, from which the voltage-time series are sampled from dynamics regions of interest, arbitrarily set by diagrams visual inspection.The bifurcation diagrams are plotted in 3-dimension to a broad visualization, but they are shown individually in Supplementary Material.The voltage-time series are also shown in Supplementary Material.

Amplitude analysis
Forcing amplitude on the forward and backward bifurcation diagrams under different frequencies for the symmetric and asymmetric nonlinear models are displayed in Fig. 3.
The symmetrical nonlinear power forward and backward bifurcation diagrams, Fig. 3a, reveal a wealth of chaotic dynamics responses, providing regular multiple period regions and discontinuities.We can quickly note an improvement in recovering power for high frequencies.Regular dynamics is prevalent in the early portion of the amplitudes interval for Ω ≤ 0.7, both on forward and backward diagrams, showing the lowest output power values.Higher power outputs are noted for frequencies from 0.5 onwards on the backward case, with a smooth output power increase starting at about f = 0.015.
Multiple period regions emerge on backward diagrams for the highest amplitudes at Ω = 0.5.For higher frequencies, the last portion of the amplitudes interval, i.e., for the highest amplitude values, the diagrams remain regular and unchanged both on forward and backward sweeping, which can be an interesting operational condition for practical applications.The output power reaches a significant steady-state value about f ≥ 0.2 and Ω = 0.6 onwards, disregarding of sweeping direction.Despite this, the lower frequencies lead to chaotic behavior for the same amplitudes in both cases.
The forward diagrams differ from the backward ones for high frequencies, revealing massive dynamic regions of chaotic behavior for Ω = 0.8 and Ω = 0.9 on the forward case, depicted by blurred areas.The analysis of such diagrams for these frequencies also reveals multiple periodic behaviors between the chaotic dynamic regions for both cases and even before them, for Ω = 0.9, as seen from the output voltage waveforms depicted in Supplementary Material.
The asymmetrical nonlinear power forward and backward bifurcation diagrams for ϕ = 35 • are shown in Fig. 3b.The system presents a unique periodic solution for all amplitude and frequency ranges.In addition, the power output is low for all excitation conditions, indicating poor behavior for energy harvesting.The recovered power of the asymmetric case is less than the symmetric case in the whole range of amplitude and frequency.The strong asymmetry is a physical reason for this behavior, requiring more external energy for the system to perform the snap-through motion and harvest more energy.
The bifurcation diagram for the asymmetric system using the optimal sloping angle, ϕ opt = −4.95• , are shown in Fig. 3c.Multiple period regions for the highest amplitudes and Ω ≥ 0.5 are noticed.There are chaotic regions for Ω ≤ 0.4 and Ω ≥ 0.8.The behavior is similar to the symmetric case, and the chaotic behaviors occur almost in the same conditions of f and Ω.Thus, for this sloping angle value, the recovered power is higher than ϕ = 35 • .

Frequency analysis
As in the amplitude case, Fig. 4 presents the bifurcation diagrams sweeping the frequency for the symmetric and the asymmetric models.The bifurcation diagrams are shown individually in Supplementary Material with the voltage-time series.
Once again, the symmetric model, Fig. 4a, provides a rich behavior regards the chaos and period multi- plicity occurrence, both in the forward and backward sweeping analysis.Due to the model parameters set, a negligible amount of power can be observed for f = 0.019 and f = 0.051 even for the higher frequencies.Despite the multiple periodic responses being pointed by the backward bifurcation (about Ω = 0.62 and Ω = 0.75), the power is still negligible.
The analysis of backward sweeping of the diagrams from f = 0.083 onwards suggests that such bifurcations evolve to chaotic patterns, as seen in about 0.6 ≤ Ω ≤ 1.2.In the same regions, the forward sweeping of the bifurcation diagrams shows higher regular power outputs, especially for f ≥ 0.115, with a smooth increasing response from Ω = 0.4 onwards.A discontinuities region spread on the early frequencies, both on forward and backward diagrams emerge for amplitudes about f = 0.115 and 0.147.
There are significant differences between the profiles of forwarding and backward sweeping of the bifurcation diagrams for f = 0.179.In the early portion of the frequencies analysis interval for f ≥ 0.179, there are regions with chaos incidence, both in sweeping forward and backward cases.Although, one can also note some narrow regions with multiple periods, both in sweeping forward and backward cases.The analysis of the respective time series in the Supplementary Material reveals details of such dynamics and corroborates the general observations.
The bifurcation diagrams sweeping the frequency for the asymmetric model when ϕ = 35 • are shown in Fig. 4b.The chaotic regions disappeared, and the dynamic behavior is periodic for all excitation values.Some specific values conducted period response to multiple frequencies, but the motion with one frequency is dominant.There are conditions where the power output is higher, mainly for higher frequencies.As seen in Fig. 3b, the strong asymmetry damaged the energy harvesting process.
The bifurcation diagrams sweeping the frequency for the asymmetric system using the optimal sloping angle, ϕ opt = −4.95• , are shown in Fig. 4c.The chaotic regions appeared for the frequencies near 0.8 and the low frequencies with high amplitude values.The high power output is seen when the amplitude value is higher than 0.115, and the frequency is higher than 0.4.For this value of the sloping angle, the behavior is similar to the symmetric case.

Basins of attraction
In order to give a global qualitative overview of the harvester's dynamics, this section presents a numerical study where the different attractors that can coexist with the underlying dynamical system are identified.This numerical study integrates the dynamics for different sets of initial conditions, randomly chosen, and maps the possible attractors in which the phase space trajectory can accumulate.
In this paper, we propose a novel approach to obtaining the basins of attraction, where the 0-1 test for chaos is used to characterize the dynamics, determining if the system evolves to the chaotic or regular response [11].This tool is based on time series and has a better computational performance compared to the Lyapunov exponent.Despite presenting more complex mathematics and requiring a steady-state regime, the 0-1 test is computationally more efficient because it handles an approach to statistical divergences.On the other side, the Lyapunov exponent requires phase space reconstruction and eigenvalue calculation, which is usually a computationally intensive task.Although the Lyapunov exponent can define the strength of chaos by the maximum exponents even in a transient regime and determine the bifurcation points, being more informative, the 0-1 test is cheaper and a more general test.
In this technique, the time series x(t) is mapped into a pair of transformed coordinates defined by the Euclidean extension x(t j ) cos (j c), (10) q n (c) = n j=1 x(t j ) sin (j c), (11) which c is a random value uniformly drawn in the support [0, 2π) and n = 1, 2, . . ., N .For regular behavior, the p n and q n coordinates show bounded motion; otherwise, for chaotic behavior, the p n and q n coordinates asymptotically behave as Brownian motion.Figure 5 show the projections of time series in the p×q space for two different values of the excitation amplitude.For periodic behavior in Fig. 5a, the projection in p×q space is bounded.However, for chaotic behavior in Fig. 5b, the projection is not limited, and the system has a diffuse motion.
In order to analyze the diffusive (or non-diffusive), the mean square deviation is calculated by so that, the classifier is obtained by the correlation  where Cov and V ar are covariance and variance operators.K c is calculated for several c values, and the median is obtained.Then, one classifies the dynamic behavior K > 0.8 as chaotic behavior, K < 0.2 as regular behavior and K < 0.8 & K > 0.2 as an inconclusive behavior.Figure 6 presents a schematic of the methodology of the 0-1 test to classify the dynamic behavior.
This process is repeated for a grid of initial conditions to obtain the basins of attraction as presented in Fig. 7.When the dynamic behavior is classified as a chaotic response, we scheduled to color the node with a grey color.When the dynamic behavior is classified as a regular response, we introduce one more step to identify the attractors and distinguish the different regular orbits.The green color is a pattern color to indicate the regular high-energy orbit that is the best scenario for the harvesting process.The red and blue colors indicate regular low-energy orbits.Other colors, e.g., magenta and cyan, are also used to indicate regular behaviors different from red and blue.

Symmetric bistable energy harvester
When Ω = 0.1, two regular low energy attractors are observed, with low amplitude orbits concerning displacement and velocity.As the excitation frequency increases, a narrowing of the basins is clear, as well as chaotic behavior.There is a predominance of a specific attractor, red (Ω ∈ {0.3, 0.5}), blue (Ω = 0.4) and green (Ω ∈ {0.6, 0.7, 0.8}) and, in boundary regions, erosions are identified through fractals.In particular, at Ω = 0.4, there are cyan and magenta attractors, but they are imperceptible in the basins of attractions.Furthermore, near the Ω = 0.8, the basins narrow as they drift away from the equilibrium points.At Ω = 0.8, there are only high energy, periodic and chaotic orbits.Three types of attractors are associated with the dynamical system: a chaotic, and two regular attractors, with low or high energy orbits.
Figures 10 and 11 shows, respectively, the basins of attractions and the attractors for the amplitude of excitation f = 0.115 and the frequency of excitation Ω going from 0.1 to 0.9.Under these conditions, it is notable that chaotic behavior is present at almost all excitation frequencies, except for the Ω = 0.3, Ω = 0.6, and 0.7.The low energy orbits, blue and red attractors, can be found, mainly at low frequencies.At Ω = 0.3, a yellow basin appears with high energy orbits with more than one frequency.For Ω = 0.5, besides the green attractor emerges, the blue and red orbits present two frequencies.The green attractor, which has a high energy orbit, is also frequent and dominant as the frequency increases.At Ω = 0.8, the basins are entirely covered by high energy orbit, demonstrating the desired situation for the energy harvesting process.At Ω = 0.9, there is a vast presence of chaos, along with the green basin and, scattered through chaos, the blue and red basins with two frequencies.
The following analysis concern the response of the system when the amplitude of excitation varies by 0.019 ≤ f ≤ 0.275 with 0.032 intervals and the frequency is fixed at Ω = 0.8.This frequency value is defined, in particular, due to the excellent performance in the occurrence of green basins.Figures 12 and 13 displays the basins of attractions and the respective attractors.For f = 0.019, four regular attractors with pulverized regions express the sensitivity to the initial conditions.Attractors display the orbits, all of which have a low energy measure.As the amplitude of excitation increases to f = 0.051, the green basin appears, but the purple           attractor is predominant.The results for f = 0.083 and f = 0.115 have already been discussed previously.The primary feature to note is that, as the amplitude growsand more energy enters the system -the more energetic the attractors become.For f = 0.115, the responses are qualitatively the same, and the amplitude of excitation is so high that the only possible basin is the most energetic or the chaotic attractor.This behavior is noted for higher amplitudes of excitation.

Asymmetric bistable energy harvester
The analysis is done for the asymmetric bistable model when the ϕ = 35 • .Figures 14 and 15 displays the basins of attractions and the attractors.The frequency of excitation is fixed at Ω = 0.8, which is the frequency of the best performance for the symmetric model, and its amplitude of excitation varies by 0.019 ≤ f ≤ 0.275 with 0.032 intervals.The initial displacement × initial velocity plane is defined in the region (x 0 , ẋ0 ) at the restriction plane to v 0 = 0.For f = 0.019, the red basins, low energy attractor, is prevalent throughout the plane.The orbit is so small that the attractor looks like a point, which is the worst condition for energy harvested visualized in the results of this work.For f ∈ {0.051, 0.083}, there are three solutions, green, red and magenta basins.The energy of the red attractor increases, and it is possible to visualize its orbit.The green attractor is visually asymmetric, and the positive equilibrium point has a deeper well than the negative one.The magenta basins arise near the negative equilibrium point, and its attractor is a homoclinic orbit.This basin vanishes when the f = 0.115 and part of this solution becomes a yellow basin with two periods solution.For f = 0.147, there are two solutions, and the green basin is dominant, but the other solution presents low energy.As the amplitude of excitation increases, the red basins become more frequent, and the green basins almost disappear.In particular for f = 0.243, a purple basin arise with several frequencies.The chaotic behavior is not visualized for the asymmetric model in any condition investigated in this work.
The blue basins, low energy orbit around the negative equilibrium point, did not appear in any case.This indicates that the solution is attracted to the other equilibrium point, which is physically explained by the deep of the energy well of this equilibrium point, which is shallow.This way, when the system exceeds the poten- tial barrier, the system performs a snap-through tion.On the other hand, when the initial conditions are near to negative equilibrium point, even with low external energy, the system goes to a positive one, and the solution is attracted to a red orbit.Note that the blue basins often appear near the negative equilibrium point for the basins of attractions of the symmetric model (Fig. 12).However, this region is filled by red basins (or basins with a snap-through motion) for the asymmetric model.
The basins of attractions of the asymmetric model for ϕ = 35 • demonstrate unusual behavior.The low orbit response is dominant as more external energy is inserted into the system.The best energy harvesting performance in this angle value is not for the high excitation amplitude.This system can harvest more energy at f = 0.147, where the green basin is visually dominant.However, in all scenarios, the symmetric model performs better than the asymmetric model.Therefore, the asymmetry is an undesired configuration for energy harvesting.

Optimal asymmetric bistable energy harvester
The basins of attractions were also obtained to the asymmetric bistable model for the optimal slopping plane, ϕ opt = −4.95.Figures 16 and 17 shows the basins of attractions and the attractors, respectively.The frequency of excitation is fixed at f = 0.8, while its amplitude of excitation varies by 0.019 ≤ f ≤ 0.275 with 0.032 intervals.The initial displacement × initial velocity plane is defined in the region (x 0 , ẋ0 ) ∈ [−3 ≤ x 0 ≤ 3] × [−3 ≤ ẋ0 ≤ 3] at the restriction plane to v 0 = 0.
For f = 0.019, there are four solutions colored with blue, cyan, red, and magenta basins with low energy response.The green basin arises for f = 0.051, and the system has five solutions.Besides the red and blue basins, a chaotic solution appears when f = 0.083.In this value of the amplitude of excitation, the green attractor disappears.For f = 0.115, two solutions have snap-through motions; one is chaotic, while another is regular.However, the chaotic solution is imperceptible, and the green basins fill the entire region.As the amplitude of excitation increases, the chaotic solution disappears, and the system has only the green attractor as a solution.
It is essential to observe that the attractors are symmetric for this value of the slopping angle.This in-     , where it is near the optimal angle, the system has two basins, green and purple.The purple basin is a multiple periods response.For ϕ = 0 • , the asymmetry of the system is driven only by the δ coeffi-cient, and the system has two solutions, green and red.The red basin is a low-energy orbit around the positive equilibrium point.As the ϕ increases, the asymmetric also increases, and the system has green and red basins.In addition to these basins, for ϕ = 35 • , there is a small yellow basin with a homoclinic orbit.
The asymmetry is inverted after ϕ = −5.For negative values of the sloping angle, the blue basin concentrates at negative values of the displacement.When the value of ϕ approximates the optimal value, a part of this region goes toward positive displacement values.After the optimal value, the blue basin is replaced with red, and the blue basin attracts no more responses.As the value of ϕ increases, the red solution is concentrated around the positive equilibrium point.
Figure 20 shows the ratio of the occupied zone of each basin of attraction on the domain of analysis ( (x 0 , ẋ0 ) | − 3 ≤ x 0 ≤ 3 and − 3 ≤ ẋ0 ≤ 3 ) for the symmetric, asymmetric (with ϕ = 35 • ), and optimal asymmetric models for the different amplitude of excitation when Ω = 0.8.We plotted the ratios of the occupied zone by the stacked bar, and each color corresponds to the basins of attractions as defined previously.We can compare the behavior of the basins for the different models and excitation amplitude.The symmetric (SM)  and optimal asymmetric (OAM) models present almost the same ratios of the occupied zone for each basin.As the amplitude increase, the green basin is dominant, the best scenario for harvesting energy.For the asymmetric model with ϕ = 35 • (AM35), the red basin is dominant, and as increase the amplitude of excitation, the green basin achieves 50% of the occupied zone and drastically drops.For all values of the excitation amplitude, the asymmetric with ϕ = 35 • model presents a smaller region of green basins than the others.

Conclusions
This paper numerically investigates the dynamic behavior of the bistable energy harvesters.Symmetric and asymmetric models were investigated.A quadratic term for the nonlinear restoring force and a sloping angle in the plane where the system is attached was introduced to consider the asymmetries.Two asymmetric conditions were studied, when the system has a strong influence of asymmetry and at a specific value of the sloping angle that the asymmetry in the potential function is canceled.
A comprehensive analysis was performed for a broad range of values of the amplitude and frequency of ex-citation and sloping angle.Forward and backward bifurcations diagrams of the amplitude and frequency indicated multiple solutions.Basins of attraction calculations, based on 0-1 test for chaos, were proposed to investigate the attractors for different initial conditions.This analysis demonstrated multiple solutions for initial conditions and exhibited asymmetries of the basins.The dynamic behaviors of the asymmetric model were analyzed, and its negative effect on harvesting was proved.The performance of the optimal angle demonstrated that the system has similar dynamic behavior to the symmetric bistable, indicating that it is possible to around the negative effect of asymmetry using the slope angle.Therefore, this paper extensively analyzes the nonlinear energy harvesting system, providing crucial information on their dynamic behavior and performance.The best excitation conditions to harvest more energy were defined.The comprehension of how initial conditions attract the solutions and the asymmetries affect the harvesting performance were exposed.

Figure 1 :
Figure 1: Illustration of the two piezo-magneto-elastic energy harvesting systems analyzed in this work.

Figure 2 :
Figure 2: Schematic representation of the methodology used to construct the bifurcation diagrams.

Figure 3 :
Figure 3: Power bifurcation diagram sweeping forward (blue scale) and backward (red scale) amplitude of excitation of the symmetric and asymmetric bistable energy harvesters.

Figure 4 :
Figure 4: Power bifurcation diagram sweeping forward (blue scale) and backward (red scale) frequency of excitation of the symmetric and asymmetric bistable energy harvesters.

Figure 5 :
Figure 5: Dynamic extension in p × q space.The bounded projection is shown for f = 0.115 and the diffuse projection for f = 0.083.

Figure 6 :
Figure 6: Schematic of the 0-1 test for chaos, which classifies the dynamic behavior as regular or chaotic based on the statistical behavior of an Euclidian extension of the original dynamics.

Figure 7 :
Figure 7: Schematic representation of the methodology used to compute the basins of attractions.The 0-1 testfor chaos is executed for each initial condition.When the motion is chaotic, we schedule to catalog the grid with a grey color.If the motion is regular, there is one step to attractors identification, and the grid is cataloged in a color scheme for different attractors.

Figure 20 :
Figure 20: Relative area of the basin with respect to the domain of initial conditions for the symmetric model (SM), asymmetric model with ϕ = 35 • (AM35), and optimal asymmetric model with ϕ = −4.95• (OAM) under different amplitude of excitation for Ω = 0.8.