Surrogacy-Based Maximization of Output Power of a Low-Voltage Vibration Energy Harvesting Device

Featured Application: The analyzed system is a nonlinear electromagnetic resonant vibration energy harvester designed with the intention of application in a low-cost battery management system of a wireless sensor network. The developed idea of maximization of the output power generated by the system makes it more practical due to noticeable improvement of the rating and simple manufacturing. Abstract: The coreless microgenerators implemented in electromagnetic vibration energy harvesting devices usually su ﬀ er from power deﬁciency. This can be noticeably improved by optimizing the distribution of separate turns within the armature winding. The purposeful optimization routine developed in this work is based on numerical identiﬁcation of the turns that contribute most to the electromotive force and the elimination of those with the least contribution in order to reduce the internal impedance of the winding. The associated mixed integer nonlinear programming problem is solved comparatively using three approaches employing surrogate models based on kriging. The results show very good performance of the strategy based on a sequentially reﬁned kriging in terms of the ability to accurately localize extremum and reduction of the algorithm execution time. As a result of optimization, the output power of the system increased by some 300 percent with respect to the initial conﬁguration.


Introduction
There has been considerable interest in small converters harvesting energy from various sources available in the environment, such as industrial machinery, civil structures, and even humans or animals. These systems are applicable in sensor networks or monitoring systems [1][2][3][4][5][6][7][8][9] that require autonomy of power supply. The use of small generators of electric energy in such systems makes it possible to significantly increase an overall reliability, to eliminate wiring and reduce maintenance service activity over the whole system [3,[10][11][12][13][14].
Most of the attention in the literature is focused on systems harvesting energy from mechanical vibrations. The minigenerators applied in these systems use various physical phenomena as the means of energy conversion. The most commonly implemented are piezoelectric, magnetostrictive, magnetoelectric, electrostatic and electromagnetic converters [2,15,16]. This work concerns the latter system due to its reasonable power per unity of mass ratio and low manufacturing cost compared to the other types [11,13].
The coreless (iron-free) magnetic circuit is a common choice for application in electromagnetic vibration energy harvesting converters due to its simple design guidelines [17]. It is well known that these systems suffer from a low value of induced electromotive force (emf), and thus generated power. Their very low volt-per-turn ratio often makes the impedance of the armature coils too high for the electronic boost of the induced voltage to any practical value. Here, we propose a routine that improves the performance of such a system by means of providing an optimal distribution of individual turns within the coil. This allows for a reduction of the impedance to an extent where the microgenerator connected to the electronic converter develops considerable output power.
In the design of electromagnetic vibration energy harvesters, the electromagnetic part is modelled using either numerical or analytical models [8,[17][18][19]. In this work, the model segment of the routine involves determination of useful electromagnetic quantities by combining a 3D integral method, analytical formulas with a steady-state circuit simulation considering the electronic converter. For the purposes of the optimization process, it is not necessary to consider the whole system and its complicated kinematics resulting from the connection of the generator parts with elements of mechanical support [20]. As a consequence, the motional effects are modelled as simply as possible, assuming that the displacement of the moving element of the converter is sinusoidal in time, which means that the system operates at a prescribed frequency and magnitude.
The optimization segment aims at the maximization of the power generated by the whole system. The considered optimization task is a mixed integer nonlinear programming (MINLP) problem [21]. Due to the considerable evaluation time of the objective function, the kriging surrogate model [22] is involved in the process. Three solution strategies are comparatively examined with respect to their cost and accuracy characteristics. The optimized system is successfully built and tested. The new contributions of this work are:

•
Elaboration of an efficient and uncomplicated routine for the maximization of output power for a class of microgenerators for vibration energy harvesting devices; • Comparative analysis of features for three solution strategies for the associated MINLP problem using a surrogate model with particular attention focused on the strategy using sequential refinement of the computer experiment grid [23].

System and Its Mathematical Model
The CAD diagram of the considered microgenerator, variations of potential energy, and electronic converter are depicted in Figure 1a-c.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 13 power. Their very low volt-per-turn ratio often makes the impedance of the armature coils too high for the electronic boost of the induced voltage to any practical value. Here, we propose a routine that improves the performance of such a system by means of providing an optimal distribution of individual turns within the coil. This allows for a reduction of the impedance to an extent where the microgenerator connected to the electronic converter develops considerable output power.
In the design of electromagnetic vibration energy harvesters, the electromagnetic part is modelled using either numerical or analytical models [8,[17][18][19]. In this work, the model segment of the routine involves determination of useful electromagnetic quantities by combining a 3D integral method, analytical formulas with a steady-state circuit simulation considering the electronic converter. For the purposes of the optimization process, it is not necessary to consider the whole system and its complicated kinematics resulting from the connection of the generator parts with elements of mechanical support [20]. As a consequence, the motional effects are modelled as simply as possible, assuming that the displacement of the moving element of the converter is sinusoidal in time, which means that the system operates at a prescribed frequency and magnitude.
The optimization segment aims at the maximization of the power generated by the whole system. The considered optimization task is a mixed integer nonlinear programming (MINLP) problem [21]. Due to the considerable evaluation time of the objective function, the kriging surrogate model [22] is involved in the process. Three solution strategies are comparatively examined with respect to their cost and accuracy characteristics. The optimized system is successfully built and tested. The new contributions of this work are: • Elaboration of an efficient and uncomplicated routine for the maximization of output power for a class of microgenerators for vibration energy harvesting devices; • Comparative analysis of features for three solution strategies for the associated MINLP problem using a surrogate model with particular attention focused on the strategy using sequential refinement of the computer experiment grid [23].

System and Its Mathematical Model
The CAD diagram of the considered microgenerator, variations of potential energy, and electronic converter are depicted in Figure 1a The emf in armature coils 1 is induced as the result of the relative motion of magnets 3 along the y axis direction. Magnets 2 are stationary with respect to the armature and are used to develop an auxiliary magnetic force that brings nonlinear mechanical resonance into the system. The converter belongs to a class of nonlinear electromagnetic vibration energy harvesters whose insightful analysis including kinematic and dynamic effects can be found in [11] and [20,25]. Figure 1b shows variations of potential energy vs. relative displacement of permanent magnets (No. 3 in Figure 1a). The senses of magnetization of the magnets in Figure 1a determine the softening action of the magnetic force within the feasible range of displacement (between −5 mm and +5 mm), and two stable points around zero displacement (see Figure1b). With the variations of potential energy due to magnets (magnetic energy) and that due to auxiliary cantilever spring (spring energy) there is one static stable point at zero displacement. The frequency range where the system produces electric power is approximately equal to 10 Hz.
Due to lack of core element, the distribution of magnetic flux over the coil of the microgenerator shown in Figure1a is too nonuniform to assume that the induced emf can be calculated considering the coil as a block. For this reason, the optimization routine assumes that the emf is estimated for each turn individually. This apparently simple computation causes, however, a burden that needs further organization. Consider a single racetrack-shaped coil in Figure 2. The emf in armature coils 1 is induced as the result of the relative motion of magnets 3 along the y axis direction. Magnets 2 are stationary with respect to the armature and are used to develop an auxiliary magnetic force that brings nonlinear mechanical resonance into the system. The converter belongs to a class of nonlinear electromagnetic vibration energy harvesters whose insightful analysis including kinematic and dynamic effects can be found in [11] and [20,25]. Figure 1b shows variations of potential energy vs. relative displacement of permanent magnets (No. 3 in Figure 1a). The senses of magnetization of the magnets in Figure 1a determine the softening action of the magnetic force within the feasible range of displacement (between −5 mm and +5 mm), and two stable points around zero displacement (see Figure 1b). With the variations of potential energy due to magnets (magnetic energy) and that due to auxiliary cantilever spring (spring energy) there is one static stable point at zero displacement. The frequency range where the system produces electric power is approximately equal to 10 Hz.
Due to lack of core element, the distribution of magnetic flux over the coil of the microgenerator shown in Figure 1a is too nonuniform to assume that the induced emf can be calculated considering the coil as a block. For this reason, the optimization routine assumes that the emf is estimated for each turn individually. This apparently simple computation causes, however, a burden that needs further organization. Consider a single racetrack-shaped coil in Figure 2. δ-air-gap between wires; w, h-assumed width and height of coil walls, respectively.
As shown in the figure, the turns over a cross section of a uniformly wound coil are distributed one after another and then layer after layer in order to fill the assumed area w × h and are numbered repeatedly from left to right. Because the number of turns Nt has to be considered as a design variable, it should be noticed that it is not uniquely attributed to an electromagnetic quantity, such as the emf. It is only unique if it is attributed to turns ordered with respect to the value of the quantity. This attribution will be carried out as follows.
For the system in Figure 1, the emf is induced only by the z-th component of the magnetic flux density, whose distribution due to the i-th magnet can be determined from: where 0 μ is the vacuum permeability, is the magnetic dipole density, The useful quantity is, however, a derivative of turn-linked magnetic flux not the flux itself, as the former is in proportion to the turn-induced emf. For the j-th turn, it can be conveniently calculated using Green's theorem: where t S is the area enclosed by a single turn and t L is the path along turn. It should be noted that the use of the final result of Equation (2) considerably reduces the time required for numerical integration and eliminates unwanted "noisy" variations of the quantity due to variations of turn dimensions that appear using surface integration over t S and numerical differentiation with respect to y. After evaluation of (2) for each turn of the cluster in Figure 2, the results are ordered with respect to value by sorting and numbering from 1 to Nt. This makes the turn ordinal number a unique design variable. At this instance, it should also be noted that: • If, for example, the call from the search algorithm includes a request for Nt = 100, the above described model returns the values of quantity (2) determined for the first 100 turns with the highest contribution to the emf.

•
Quantity (2) is determined without displacement of magnets 3 from the origin of the system of coordinates. It can be deduced that this position of the magnets corresponds with the maximum value of the turn-induced emf.

•
An important observation that makes the results of the optimization practical regards the fact that the ordered turns are not spread within the coil cross section, but they form an irregular cluster. This will be demonstrated further. As shown in the figure, the turns over a cross section of a uniformly wound coil are distributed one after another and then layer after layer in order to fill the assumed area w × h and are numbered repeatedly from left to right. Because the number of turns N t has to be considered as a design variable, it should be noticed that it is not uniquely attributed to an electromagnetic quantity, such as the emf. It is only unique if it is attributed to turns ordered with respect to the value of the quantity. This attribution will be carried out as follows.
For the system in Figure 1, the emf is induced only by the z-th component of the magnetic flux density, whose distribution due to the i-th magnet can be determined from: where µ 0 is the vacuum permeability, Q mi = M i · n is the magnetic dipole density, M i and n are the magnetization and surface-normal vector, respectively, r and r are the source and observation point, respectively, and S i+ and S i− are the magnet surfaces parallel to the xy plane. The useful quantity is, however, a derivative of turn-linked magnetic flux not the flux itself, as the former is in proportion to the turn-induced emf. For the j-th turn, it can be conveniently calculated using Green's theorem: where S t is the area enclosed by a single turn and L t is the path along turn. It should be noted that the use of the final result of Equation (2) considerably reduces the time required for numerical integration and eliminates unwanted "noisy" variations of the quantity due to variations of turn dimensions that appear using surface integration over S t and numerical differentiation with respect to y. After evaluation of (2) for each turn of the cluster in Figure 2, the results are ordered with respect to value by sorting and numbering from 1 to N t . This makes the turn ordinal number a unique design variable. At this instance, it should also be noted that: • If, for example, the call from the search algorithm includes a request for N t = 100, the above described model returns the values of quantity (2) determined for the first 100 turns with the highest contribution to the emf.
• Quantity (2) is determined without displacement of magnets 3 from the origin of the system of coordinates. It can be deduced that this position of the magnets corresponds with the maximum value of the turn-induced emf.

•
An important observation that makes the results of the optimization practical regards the fact that the ordered turns are not spread within the coil cross section, but they form an irregular cluster. This will be demonstrated further.
Determination of power P generated by the system requires calculation of the coil impedance and the emf (sum of quantity (2) through j) vs. y for the given N t . The resistance and inductance are calculated using analytical formulas [26,27], taking into account the active and end-turn segments of the coil (Figure 2) where r c,i , σ Cu are inner radius of i-th turn and copper conductivity, respectively. The end-turn inductance of i-th turn L r,i , inductance of active segment of turn L z , end-turn mutual inductance of i-th and j-th turns M r,ij and mutual inductance of active segments of i-th and j-th turns M z,ij are calculated from formulas where K(κ) and E(κ) are complete elliptic integrals of the first and second kind and Please note that the above formulas are valid even when L = 0, which means that the racetrack coil becomes a solenoid. The purpose of optimization is thus not only to optimally distribute turns within the coil cross-section but also decide whether a racetrack-like shape is better or worse than a solenoidal one.
The power P is determined at a single frequency assuming a sinusoidal variation, a constant magnitude of displacement of magnets 3 and considering the power converter in Figure 1c. The steady-state time-periodic circuit simulation using interpolated variation of ∂λ/∂y vs. y is carried out for this purpose using program Plecs [28]. The diagram of the computer model which implements all of the above described equations is shown in Figure 3. only a slight impact on power. The resolution of the second variable is low as only six different wire diameters are taken into account. Moreover, the sensitivity of power to variation in D is considerable. It should also be noticed that this variable should be converted from enumeration to an integer.
Due to the considerable time required for the evaluation of power using the computer model presented in Section 2, the function P(X) is replaced by the surrogate model. The response R of the computer model is evaluated at nodes of the Latin hypercube sampling-based design (LHS) grid ) ( ) ( : Finally, the surrogate model ) (X P is created using the kriging approximation )) ( , ( ) ( : The three following strategies are employed for the solution of problem (10): Strategy I: Problem (10) is solved as its stands using nonlinear programming with an "integer" constraint equation in the form 0 ) ( . This can be considered as the most ad hoc approach that can easily be implemented using many existing nonlinear search algorithms. Here we use the combination of global → local searches using genetic [29] and the interior point algorithms [30] one after another. The choice of a genetic algorithm as a global search method results from its ability to reliably localize the region with global extremum. The interior point algorithm is applied

Optimization Strategies and Results
The considered optimization regards a 4D problem that is carried out in the following mixed integer design space X = (N t , E(D), L, R l ) (see Figures 1c and 2), with N t being an integer and E(D) the enumeration of wire diameters. The variable N t is considered continuous because its resolution is high (practical value is of order of a few hundreds) and rounding to the closest integer value has only a slight impact on power. The resolution of the second variable is low as only six different wire diameters are taken into account. Moreover, the sensitivity of power to variation in D is considerable. It should also be noticed that this variable should be converted from enumeration to an integer.
The problem definition is maximize P(X) subject to Due to the considerable time required for the evaluation of power using the computer model presented in Section 2, the function P(X) is replaced by the surrogate model. The response R of the computer model is evaluated at nodes of the Latin hypercube sampling-based design (LHS) grid X LHS : P(X LHS ) ⇐ R(X LHS ) . Finally, the surrogate model P(X) is created using the kriging approximation K : P(X) = K(X LHS , R(X LHS )).
The three following strategies are employed for the solution of problem (10): Strategy I: Problem (10) is solved as its stands using nonlinear programming with an "integer" constraint equation in the form round(X 2 ) − X 2 = 0. This can be considered as the most ad hoc approach that can easily be implemented using many existing nonlinear search algorithms. Here we use the combination of global → local searches using genetic [29] and the interior point algorithms [30] one after another. The choice of a genetic algorithm as a global search method results from its ability to reliably localize the region with global extremum. The interior point algorithm is applied as a local search method due to its high precision and the ability to prevent violation of inequality constraints.
The strategy can be explained by the following pseudocode.
Generate LHS grid and round X 2 Perform computer experiment Create 4-D surrogate model Optimize 4-D surrogate model with "integer" constraint The results obtained for different numbers of the LHS grid points are summarized in Table 1.

Strategy II:
Relaxation of integer variable. With only one integer variable considered, the solution strategy becomes a simplification of the branch-and-bound algorithm [21] with far fewer sequences required to find the global extremum. For the solution of each 3D subproblem, we use once again the genetic and interior point algorithms. The following pseudocode explains the strategy. The assumed resolution of the LHS grid is limited to 200 points, as shown in Table 1, to provide enough accurate localization of the extremum even for the case with four variables accounted for. The results obtained from the current solution strategy are given in Table 2.  (10) is solved using a sequentially refined surrogate model. The strategy assumes starting from a low-resolution surrogate model and employs an adaptation of a design grid around the extremum to provide a reduction of execution time involving the creation of the surrogate model. The adaptation is performed repeatedly by adding extra LHS samples around the extremum unless the relative variation of the extremum is lower than the assumed error tolerance.

Strategy III: Problem
At each adaptation step the algorithm adds ten extra LHS points. The term "around extremum" means that these points are generated in the same design space, although with the boundary limited to 20 percent (in each direction) of the value of solution in previous step.
The search itself is carried out once again using the genetic and interior point algorithms one after another. The integer variable is constrained using the "integer" constraint equation, as in strategy 1. The pseudocode below explains the strategy: The results for different numbers of coarse LHS grid points are summarized in Table 3. The synthesis of all results leads to the following conclusions: • The most direct strategy (1) did not provide as high a value of power P as the other strategies even for a relatively fine LHS grid. The reason resides in the flatness of the objective function around the extremum. • Strategy 2 provides nearly the same values of power P as strategy 1, though exploration of the whole design space becomes extremely expensive. • Strategy 3 demonstrates the best performance with respect to the mentioned criteria. The highest value of power P found by the strategy confirms the flatness of the objective function, which can be correctly explored using only very fine design grids. It can also be seen from Table 3 that starting from a finer LHS grid can reduce the number of adaptation sequences. The differences between the values in the second and third row are related with the assumed error tolerance equal to 0.01 (see the pseudocode of the strategy III).
The maximum power is given in the second row of Table 3. The corresponding optimum solution is X opt = (N t = 746, D = 0.4 mm, L = 10.5 mm, R l = 1824 Ω). Figure 4 displays the distribution of wires over the cross section of a single coil of the twin-coil armature corresponding with the above solution. It should be noted that this result can easily be accomplished in the physical system, as shown in Figure 4c.   (Figure 4a) and for that with the optimized racetrack coils (Figure 4b). As shown in Figure 5, the output power increased as much as by 300%. The slight (14 percent) differences between the results in Table 3 and the maximum power measured for the optimized system are most likely due to complex (non-straight line) motional effects of magnets 3 on the induced emf in the physical harvester, which is a specific feature of the considered system that we identified and analyzed in [20], and is also due to inaccuracy of manufacturing of physical coils where the turn distribution does not match exactly one presented in Figure 4b. Though, the values of total coil resistance and inductance Rc and Lc measured using the laboratory multimeter agree to those calculated from formulas (3)-(9) within 2 and 6 percent, respectively (see Table 4).  Figure 2) before (a) and after optimization (b). Dashed line denotes assumed area w× h with r c = 2 mm, and w, h and δ equal to 22, 33 and 0.01 mm, respectively. Colormap displays distribution of quantity (2) in V·s/m. c) Assembled armature using a 3d-printed coil carcass including formed bobbin and flat lid so as to obtain the distribution of turns shown in Figure 4b. The coils in the physical system are connected in parallel. Figure 5 displays the laboratory test-stand for validation of the obtained results. The comparison of powers measured for the systems with initially designed cylindrical coils with area w × h (dashed line in Figure 4) filled uniformly with 2000 turns (Figure 4a) and for that with the optimized racetrack coils (Figure 4b). As shown in Figure 5, the output power increased as much as by 300%. The slight (14 percent) differences between the results in Table 3 and the maximum power measured for the optimized system are most likely due to complex (non-straight line) motional effects of magnets 3 on the induced emf in the physical harvester, which is a specific feature of the considered system that we identified and analyzed in [20], and is also due to inaccuracy of manufacturing of physical coils where the turn distribution does not match exactly one presented in Figure 4b. Though, the values of total coil resistance and inductance R c and L c measured using the laboratory multimeter agree to those calculated from formulas (3)-(9) within 2 and 6 percent, respectively (see Table 4).

Boost converter
Optimized system mounted on shaker   Figure 5c presents the measured characteristics of output power vs. frequency for the initial and optimal designs. The maximum power is reached at 35.5 Hz and is equal to 6.46 mW. Due to nonlinear magnetic force caused by interaction between moving magnets (magnets 3 in Figure 1a) and stationary magnets (magnets 2 in Figure 1a), the proposed energy harvester operates in the wide range of vibration frequency (around 10 Hz) providing output power sufficient to charge small Li-Ion battery in a power supply of a wireless sensor network. In comparison to reference systems of similar size operating at the same value of acceleration our converter has very similar range of operation frequency (7.8 Hz in [31] and 10 Hz in [32]), although it provides significantly higher output power (1.18 mW in [31] and 2.5 mW in [32]). The latter confirms the applicability and novelty of the proposed method for maximization of the output power of such systems.

Discussion
This work has shown a new possibility to significantly increase the output power for a class of electromagnetic vibration energy harvesting devices with its mathematical framework. The term "class" refers here to systems employing the coreless electromagnetic minigenerators that to this moment had to be designed with high number of turns detracting from their ability to deliver practical value of the output power. The overall idea of optimal placement of coil turn would also be applicable to other classes of converters, although the model for determination of circuit parameters outlined is Section 2 may not be applicable directly.
It should be noted that the method was successful due to intentional choice of structure of the electronic converter which has an ability to boost the input voltage from values lower than 0.05 V, and which we selected testing earlier a few other structures. However, for such an excellent performance it requires a very small impedance of the coil connected to the input terminals. Without the boost converter the voltage generated by the optimized system would be even smaller than that of the reference systems mentioned in the preceding section.
It must be admitted that, generally, an optimization strategy based on the refined surrogate model, which we identified in this work as the best among the three considered strategies, is a known idea [23]. However, it was used here in a different way and validated with respect to its potential to solve the mixed integer nonlinear optimization problem.
The successful optimization prompts the authors to modify the whole designing procedure of the converter so as to take account for more design criteria such as frequency bandwidth and overall mass. An application of the optimization strategy based on the refined surrogate model in solution of such a multi-objective optimization problem is a very challenging task.