This work deals with the solution of a non- convex optimization problem to enhance the perfor- mance of an energy harvesting device, which involves a nonlinear objective function and a discontinuous con- straint. This optimization problem, which seeks to nd a suitable con guration of parameters that maximize the electrical power recovered by a bistable energy har- vesting system, is formulated in terms of the dynamical system response and a binary classi er obtained from 0-1 test for chaos. A stochastic solution strategy that combines penalization and the cross-entropy method is proposed and numerically tested. Computational exper- iments are conducted to address the performance of the proposed optimization approach by comparison with a reference solution, obtained via an exhaustive search in a re ned numerical mesh. The obtained results il- lustrate the e ectiveness and robustness of the cross- entropy optimization strategy (even in the presence of noise or in moderately higher dimensions), showing that the proposed framework may be a very useful and pow- erful tool to solve optimization problems involving non- linear energy harvesting dynamical systems. Keywords energy harvesting nonlinear dynamics  cross-entropy method evolutionary algorithm  non-convex optimization 0-1 test for chaos 1 Introduction Energy harvesting is a process of energy conversion in which a certain amount of energy is pumped from A. Cunha Jr Rio de Janeiro State University { UERJ, Institute of Mathe- matics and Statistics, Rio de Janeiro, Brazil ORCID: 0000-0002-8342-0363 E-mail: americo.cunha@uerj.bran abundant source (e.g. the sun, wind, ocean, envi- ronmental vibrations, electromagnetic radiation, etc.) into another system that stores and/or use this en- ergy for its self-operation [24, 48, 59]. This type of technology has a broad eld of applicability, that in- cludes: (i) large-scale devices such as ocean thermal en- ergy converters [38], magnetic levitators [52], etc.; (ii) middle size applications like electromechanical systems [26, 37], sensors/actuators [5, 23, 61], living trees [66] etc.; (iii) small or very small systems, namely medical implants [47], micro/nano electro-mechanical systems (MEMS/NENS) [44, 58], graphene structures [42], bio cells [8], etc. Enhance the amount of energy collected by an en- ergy harvester is a key problem in the operation of these kind of devices, being the object of interest of sev- eral research works among the last decade, with e orts divided essentially in two fronts: (i) approaches with great physical appeal, which seek to explore complex geometric arrangements and nonlinearities in a smart way to improve the performance of the system[1, 3, 18, 21, 32, 33, 50, 62, 68]; (ii) works that look at this task from a more mathematical perspective, conducting theoretical analysis in mathematical models [51], us- ing control theory to optimize the underlying dynamics [2, 20, 60, 65] or formulating and solving non-convex op- timization problems to nd the optimal system design [27, 39, 40, 46, 67]. The approaches that seek to optimize energy har- vesters exploring their physics, in general, try to take advantage of a key (and very particular) characteristic of the dynamics, often neglected in a naive look, in favor of better performance. Although they can produce very ecient solutions for certain types of systems, they lack generality, which makes them, in a sense, an \art", as 2 Americo Cunha Jr it requires a high level of knowledge about the system behavior by the designer's part. On the other hand, approaches based on numerical optimization have a good level of generality, however, they come up with theoretical and computational di- culties inherent to the solution of non-convex optimiza- tion problems de ned by objective functions that are constructed from observations of the nonlinear dynam- ical systems. These diculties are faced not only when dealing with energy harvesters, but they can also be seen in the optimization of several types of nonlinear dynami- cal systems, generally being circumvented with the use of computational intelligence-based algorithms [43, 49], since traditional gradient-based methods have no guar- antee of nding global extremes in the absence of con- vexity [6, 7, 45]. Among the computational intelligence algorithms available in the literature, the most frequently used for optimization in nonlinear dynamics include genetic al- gorithms (and their variants), particle swarm optimiza- tion, di erential evolution, arti cial neural networks (and other machine learning methods), etc [49]. As they are global search methods, they are often successful in over- coming the (local) diculties faced by gradient-based methods, at the price of losing computational eciency. But the loss of eciency is not the only weakness of these global methods, most of then need to have the underlying control parameters tuned for proper func- tioning. This task is usually done manually, in a trial and error fashion, which is not at all practical, as the meaning of many of these parameters is often not in- tuitive. This peculiarity practically \condemns" these tools to be used in the black-box format, without much control by the user, which may induce signi cant losses in performance and accuracy. However, the computational intelligence literature has at least one global search algorithm that is robust and simple, for which the control parameters are very intuitive, it is known as the cross-entropy (CE) method , proposed by R. Rubinstein in 1997 [53, 54, 56] for rare events simulation. Soon after it was realized that it could be a very appealing global search technique for challenging combinatorial and continuous optimization problems [36, 57]. It is a sampling technique, from the family of Monte Carlo methods, which iteratively at- tacks the problem, re ning the candidates to be the solution according to a certain optimality criterion. Surprisingly, the nonlinear dynamics literature has been neglecting the CE method, although it is a rela- tively known technique in the combinatorial optimiza- tion community. To the best of the author's knowledge, there are few studies in the open literature applyingthe CE method for numerical optimization problems involving nonlinear mechanical systems [15, 16, 17, 25, 34, 63], which leaves space for new contributions in this line. In particular, the development of a simplis- tic and robust optimization framework, easily customiz- able, which can be used by naive users without major performance losses. In this context, this work deals with the numerical optimization of nonlinear energy harvesting systems, trying to nd a strategy to increase the performance of a bistable energy harvesting device subjected to a periodic excitation. For this purpose, a non-convex op- timization problem, with a nonlinear objective func- tion and discontinuous constraint, is formulated and a stochastic strategy of solution that combines penaliza- tion and the CE method is proposed and numerically tested. As this method has a theoretical guarantee of convergence and evidence in the literature proving its e ectiveness [36, 57], the proposed cross-entropy frame- work for optimization of dynamical systems has the potential to be a very general and robust numerical methodology. The rest of this manuscript is organized as follows. Section 2 presents the energy harvesting device of inter- est and the underlying dynamical system. An optimiza- tion problem associated with this nonlinear system is de ned in section 3, and the CE method strategy of so- lution is presented in section 4. Numerical experiments, conducted to test the e ectiveness and robustness of the stochastic solution strategy, are shown in section 5. Fi- nally, in section 6, nal remarks are set out. 2 Nonlinear dynamical system 2.1 Physical system The energy harvesting system of interest in this work is the (sinusoidal excited) piezo-magneto-elastic beam presented by Erturk et al. [22], which is based on the (stochastically excited) inverted pendulum energy har- vester proposed by Cottone et al. [9]. An illustration of this bistable energy harvesting system can be seen in Figure 1, where it is possible to see that it consists of a vertical xed-free beam made of ferromagnetic mate- rial, a rigid base and a pair of magnets. In the beam upper part, there is a pair of piezoelectric laminae cou- pled to a resistive circuit. The rigid base is periodi- cally excited by a harmonic force, which, together with the magnetic force generated by magnets, induces large amplitude vibrations. The piezoelectric laminae convert the energy of movement into electrical power, which is dissipated in the resistor. Enhancing the performance of a bistable energy harvesting device via the cross-entropy method 3 Ferromagnetic beamPiezoelectric materialRigid baseResistor Magnets N Nx Figure 1 Illustration of the bistable piezo-magneto-elastic energy harvesting system. Although this system dissipates energy, instead of using it to supply some secondary system, it is a typical prototype of a piezoelectric energy harvesting system. In an application of interest, the resistor is replaced by a more complex electrical circuit, which stores (and possibly manipulates) the voltage delivered by the me- chanical system. 2.2 Initial value problem The dynamic behavior of the system of interest is de- scribed by the following initial value problem x+ 2_x1 2x 1x2  v=fcos ( t); (1) _v+v+_x= 0; (2) x(0) =x0;_x(0) = _x0; v(0) =v0; (3) whereis the damping ratio; is the piezoelectric cou- pling term in mechanical equation; is a reciprocal time constant;is the piezoelectric coupling term in electri- cal equation; fis the external excitation amplitude; is the external excitation frequency. The initial condi- tions arex0, _x0andv0, which respectively represent, the beam edge initial position, initial velocity and the initial voltage over the resistor. Also, tdenotes the time, so that the beam edge displacement at time tis given by x(t), and the resistor voltage at tis represented by v(t). The upper dot is an abbreviation for time-derivative, and all of these parameters are dimensionless.2.3 Mean output power Since the objective in this work is to maximize the amount of energy recovered by an energy harvesting process, the main quantity of interest (QoI) associated with the nonlinear dynamical system under analysis is the mean output power P=1 TfT0ZTf =T0v2( )d ; (4) which is de ned as the temporal average of the instan- taneous power v2over a given time interval [ T0;Tf]. This QoI plays the role of objective function in the op- timization problem de ned in section 3. 2.4 Dynamic classi er Due to dynamical system nonlinearity, the steady-state dynamical response (over a given time interval) of the energy harvesting device may be chaotic or regular (non- chaotic), such as illustrated in Figure 2, which shows typical voltage time-series for this kind of bistable os- cillator. 01000 2000 3000 4000 5000−2−1.5−1−0.500.511.52 timevoltage (a) regular dynamics 01000 2000 3000 4000 5000−2−1.5−1−0.500.511.52 timevoltage (b) chaotic dynamics Figure 2 Two typical voltage time-series for the bistable energy harvesting system. The time-series in (a) has a regu- lar steady-state dynamics, while in (b) a chaotic steady-state dynamics is observed. To distinguish between the chaotic and regular dy- namic regimes, the 0-1 test for chaos by Gottwald and 4 Americo Cunha Jr Melbourne [28, 29, 30, 31] is employed. This test, which is based on an extension of the dynamical system to a two-dimensional Euclidean group [4], uses a binary classi erKto identify the dynamic regime of oper- ation underlying the system of interest. This classi- er is constructed from a dynamical system observa- tion (time-series) (t). In fact, it is sucient to use = (t1); (t2);; (tN) , a discrete version of the observable (t) that is obtained through a numerical integration (sampling) process. Receiving the discrete observation (time-series) as an input, the analytical procedure of the 0-1 test con- sists of the following steps: 1. A real parameter c2[0;2) is chosen and, for n= 1;2;;N, the discrete version of (t) is used to de ne the translation variables pn(c) =nX j=1 (tj) cos (jc); (5) qn(c) =nX j=1 (tj) sin (jc): (6) 2. Then, for n= 1;2;;N, the time-averaged mean square displacement of the dynamics trajectory in (pc;qc) space is computed by Mn(c) = lim N!11 NNX j=1~Mn(c); (7) where ~Mn(c) = pj+n(c)pj(c)2+ qj+n(c)qj(c)2: (8) 3. Finally, de ning the mean-square vector Mn= (M1;M 2;;Mn); (9) and the temporal mesh vector tn= (t1;t2;;tn); (10) the dynamic classi er is constructed through the correlation Kc= lim n!1cov(tn;Mn)p var(tn)var(Mn); (11) where cov(;) and var() respectively denote the covariance and variance statistical operators.It can be proved that Kc2f0;1g, beingKc= 0 for regular dynamics and Kc= 1 for chaotic dynam- ics [30, 4]. For numerical implementation purposes, the above procedure is performed several times, for several randomly chosen values of c2[0;2), and considering the voltage time-series as the system observation, i.e., (t) =v(t). The discrete version of (t) is constructed through temporal integration via the standard 4th or- der Runge-Kuta method. The limit processes in Eqs.(7) and (11) are replaced by the condition n N, and the classi erKis calculated as the median of Kcrealiza- tions, i.e.,K=median (Kc). Indeed, for a careful done numerical simulation one has K0 orK1. See [4, 31] for further details. 3 Non-convex optimization problem 3.1 Problem de nition The objective of this work is to nd a set of parame- ters that maximize the mean output power dissipated in the resistor. For this purpose, the electromechani- cal system properties ( , ,and) sound as natural choices for the design variables, once they are intrin- sic to the energy harvester embedded physical char- acteristics, something a designer might want to mod- ify for better performance. However, the present work uses, rst, the dynamical system excitation parameters, namely the amplitude fand the frequency , as de- sign variables. Although this choice sounds quite arti- cial at rst since, in general, the designer does not have control of the level (or frequency) of vibration in which the energy harvester is operating, it is also of interest in the optimization of these systems to know in which vibration settings the device can recover more energy. As the numerical methodology developed here is completely general, it can be easily applied in the cases where the other parameters are used as the de- sign variables, but by using fand to test the pro- posed methodology, an advantage is gained in terms of physical intuition about underlying physics, as the lit- erature presents better knowledge about the e ect of these two parameters on the behavior of this bistable energy harvesting system. Many combinations of ( f; ) lead the dynamical system to operate in the chaotic regime, an undesirable condition for electric power usage in principle, since the process of rectifying the chaotic electrical signal, which is necessary to enable its use in applications, can con- sume a signi cant part of the available energy. Although it is possible to explore the chaotic dynamics in favor of greater eciency of the system [2, 18], as done by the author and collaborators in [20], through the use Enhancing the performance of a bistable energy harvesting device via the cross-entropy method 5 of chaos control techniques, this is not the approach followed in the present work, which seeks to develop a very general numerical framework. Thus, as not every pair ( f; ) is an acceptable choice for the optimal design, it is necessary to impose a con- straint that ensures the regularity (non-chaoticity) of the system dynamics. Taking advantage of the 0-1 test for chaos, described in section 2.4, the constraint to ensure a regular (non-chaotic) dynamic regime can be formulated as K= 0. Note that the optimization prob- lem is extremely nontrivial, once the constraint presents jump-type discontinuities since K2f0;1g. 3.2 Constrained problem formulation In an abstract way, one can formulate the constrained nonlinear optimization problem described above as nd a feasible vector of design variables x?that maximize a certain objective function x2D7!S (x)2R, i.e., x?= arg maxx2DS(x); (12) where the set of admissible (feasible) parameters is de- ned by the bounded region D= xjxminxxmax andG(x) = 0 : (13) In this context, it is straightforward to see that the design variables vector is x= (f; ), the objective function isS(x) =P(f; ), the binary constraint is G(x) =K(f; ), and the limits of the design variables arexmin= (fmin; min) and xmax= (fmax; max). In the case of a problem with more variables, or with di erent objective function and/or constraints, the adaptations to be made are straightforward. 3.3 Penalized problem formulation A penalized version of this constrained nonlinear opti- mization problem is introduced here in order to facili- tate the computational implementation of the solution algorithm. In this formulation the constraint G(x) = 0 is replaced by the weaker condition G(x)" 1, once in practice the best one has is K0. In this way, the problem de ned by (12) is replaced by the penalized problem which seeks a pair x?such that x?= arg maxx2D~S(x); (14) where the set of feasible parameters is now de ned by D= xjxminxxmax ; (15)and the penalized objective function is given by ~S(x) =S(x) max 0;G(x)" : (16) The penalty parameter is heuristically chosen, being = 10 the value used in all simulations reported in this work. In the same way, "= 1=10 is adopted. 4 The cross-entropy method 4.1 The general idea of CE method The key idea of the CE method is to transform the given non-convex optimization problem into an \equivalent" rare event estimation problem, that can be eciently treated by a Monte Carlo like algorithm. The only re- quirement is that the problem has a single solution, i.e., the global extreme is unique. In this framework, the feasible region is sampled ac- cording to a given probability distribution chosen by the user, and the low-order statistics (mean and standard deviation) of these samples are used to update the opti- mum point estimation, and the stopping criteria metric (respectively). This is a two-step iterative process: 1.Sampling : First, the feasible region is sampled ac- cording to a given probability distribution, and then the objective function is evaluated in each one of these samples. The information embedded in these samples are used in the algorithm adaptive process; 2.Learning : A special subset of these samples, dubbed theelite sample set , is de ned by the samples that produced the highest values for the objective func- tion. The parameters of the probability distribution are updated using statistics obtained from this elite sample set, modifying the given distribution in a sense that tries to make it as close as possible to a Dirac delta centered on the global optimum. Throughout this process, the distribution mean pro- vides the approximation for the global optimum. Its up- date is done in the sense of moving the distribution cen- ter towards the optimization problem optimum, while the standard deviation is reduced, thus \shrinking" the distribution around its central value. The continuous composition of these e ects of translation and \shrink- ing" characterizes the process in which the distribution is\shaped"towards a point mass Dirac distribution cen- tered on the optimal. 4.2 Theoretical framework Without loss of generality, suppose that the problem of interest is to maximize an objective function x2 6 Americo Cunha Jr D7!S (x)2R, as stated in Eqs.(12). If the penalized formulation from Eq.(14) is adopted, just consider x2 D7! ~S(x)2Rinstead ofS(x). Denote the global optimal by x?and the correspond- ing maximum value by ?=S(x?), i.e., ?= maxS(x) for all x2D. Using this global maximum as a reference value, and considering a randomized version of the vector x, de- noted by X, it is possible construct the random event S(X) ?, which represents the scenarios where the random variable S(X) is bigger or equal to the (de- terministic) scalar value ?. In other words, this ran- dom event considers the possibility of choosing, ran- domly, points in the feasible region Dthat produce val- ues greater than or equal to the global maximum. Note that this random event is related to the optimization problem de ned in Eq.(12), it can be thought of as its randomized version. Since ?is the maximum value of S(x), no realiza- tion of Xcan produceS(X)> ?and only X=x? can makeS(X) = ?. Therefore, the probability of this random event is zero, i.e., P S(X) ? = 0. Relaxing the reference value for a scalar < ?, one has the random event S(X) , which de nes a rare event if  ?, for whichP S(X) 0. In his 1997 seminal work [53], R. Rubinstein con- ceived the CE method as a tool to eciently estimate such type of probabilities, associated with events lo- cated at the distribution tail. Sometime later [54], due to the connection between the random event and the optimization problem described above, he noticed that this rare event estimation process could be used as a global search method to optimize the given objective function. Notice that, with the aid of the expected value op- erator Efg, and the indicator function 1A(x) =( 1;ifx2A 0;ifx=2A;(17) the probability of the random event S(X) can be written as P S(X) =En 1S(X) o ; (18) which provides a practical way of estimating its value. Once the right side of Eq.(18) is the mean value of the random event 1S(X) , it can be easily computed through the sample mean of Nssamples of X, i.e., En 1S(X) o 1 NsNsX k=01S(Xk) ; (19) where the realizations Xk(k= 1;;Ns) are drawn according to g(x;v), the probability density functionofX, parametrized by the hyper-parameters vector v. Very often, v= (;), where andare the the mean the standard deviation vectors of X, respectively. The idea of the CE method is to approximate the so- lution of the underlying optimization problem by solv- ing the rare event probability estimation from Eq.(18). For this purpose, it employs a multilevel approach, that generates an optimal sequence of statistical estimators for the pair ( ;v), denoted by ( b `;bv`), such that b `a:s:! ?andg(x;bv`)a:s:!(xx?); (20) i.e., the reference level tends with probability 1 to the maximum value, and the family of distributions g(;v) goes (almost sure) towards a point mass distribution, centered on the optimization problem global optimum. This sequence of estimators is optimal in the sense that it minimizes the Kullback-Leibner divergence be- tween 1S(X) andg(;v) [57]. More concretely, the feasible region Dis sampled withNsindependent and identically distributed (iid) realizations of the random vector X, drawn from its densityg(x;v). For each of these samples, the objective function is evaluated, generating the sequence of values S(X1);S(X2);;S(XNs). Then, an elite sample Et= Xk:S(Xk)b t is de ned by lumping the Ne< Nspoints that better performed, i.e., those which produced the highest values forS(x). Note that this elite set is de ned in terms of the maximum value statistical estimator, given by b t=S(NsNe+1): (21) The hyper-parameters vector vis updated using the maximum likelihood estimator so that, with the aid of the elite setEt, it is written as bvt= arg maxvX Xk2Etln g(Xk;v) : (22) Depending on the distribution chosen for X, the stochastic program from Eq.(22) needs to be solved numerically. However, for the distributions of the ex- ponential family, which includes the Gaussian and its truncated version, this estimator can be calculated in a analytic way, with each component of mean and devia- tion given by the formulas bt=X Xk2EtXk Ne; (23) and bt=vuuutX Xk2Et(Xkbt)2 Ne; (24) Enhancing the performance of a bistable energy harvesting device via the cross-entropy method 7 respectively. Although this process has a theoretical guarantee of converging to a point mass distribution centered on the global optimum, in computational practice it is com- mon to see the updated distribution numerically de- generating before it\reaches the target". Sometimes the standard deviation decreases very quickly, causing the distribution to \shrinks" in a region far from the global optimum. In this scenario, only samples far from the op- timum point are drawn, producing poor estimates for the optimization problem solution. At the theoretical level, where it is possible to sam- ple an in nite number of times, this pathological sit- uation is bypassed, because at some moment a sam- ple on the tail is drawn, moving the distribution cen- ter out of the \frozen region", as this outlier has a lot of weight in the mean estimate. However, in computa- tional terms, as any robust implementation requires a maximum number of iterations (levels), the algorithm may stop with the distribution center within a \frozen region". This problem may be solved through the use of a smooth updating scheme for the hyper-parameters bt:= bt+ (1 )bt1; (25) bt:= tbt+ (1 t)bt1; (26) t=   11 tq ; (27) where the smooth parameters are such that 0 < 1, 0:8 0:99 and 5q10 [35, 56], and the estimations at tandt1 are obtained by solving the Eq.(22), which de nes the nonlinear program that gives rise to the vector v. 4.3 The computational algorithm The geometric idea of the CE method is described at the beginning of this section, being complemented by the theoretical formalism presented in the previous sub- section. The compilation of these ideas in the form of an easy to implement computational algorithm is pre- sented below: 1. De ne the number of samples Ns, the number of elite samples Ne<Ns, a convergence tolerance tol, the maximum of iteration levels `max, a family of probability distributions g(;v), an initial vector of hyper-parameters bv0forg, and set the level counter `= 0; 2. Update the level counter `=`+ 1; 3. Generate a total of Nsindependent and identically distributed (iid) samples from g(;bv`1), denoted byX1;;XNs;4. Evaluate the objective function S(x) at the samples X1;;XNs, sort the results S(1)S (Ns), and de ne the elite sample set Etwith theNepoints which better performed; 5. Update the estimators b `andbv`with aid of the elite sample set, using order and maximum likeli- hood statistic estimators, given by the Eqs.(21) and (22), respectively. If necessary, apply the scheme of smooth updating; 6. Repeat the steps (2) up to (5) of this algorithm while a (standard deviation dependent) stopping criterion is not met. For instance, max fg<tol. A schematic representation of this algorithm, illus- trating all the stages of the sampling and learning phases, can be seen in Figure 3. 4.4 Remarks about CE method Among the several characteristics that make this rela- tively novel technique interesting, the following can be highlighted: {Simplicity { very intuitive algorithm with few con- trol parameters ( Ns,Ne,`maxandtol), each of then with a very clear interpretation; {Robustness { theoretical results ensure that, under typical conditions, the method is guaranteed to con- verge if the problem has a single global extreme; {Eciency { the method typically presents fast con- vergence in comparison with classic global search meta-heuristics, such as genetic algorithms; {Generality { the method does not require any reg- ularity of the objective function and can be applied to almost any type of non-convex optimization prob- lem (even non-di erentiable or discontinuous); {Extensibility { the theory is general, in principle it can be applied to problems of any nite dimension, the computational cost and the \curse of dimension- ality" being the limiting factors in practice; {Easy implementation { the algorithm can be im- plemented with a few lines of code in a high level programming language. The mathematical development associated with the formalism described above is relatively non-trivial, but it has been well established throughout the rst decade of the method, including theorems that strictly estab- lish the conditions where the method has guaranteed convergence. These details are suppressed from this pa- per as they are outside the scope of the journal. But for the interested reader, the following