Time-domain model and optimization for single-axis kinetic energy harvesters driven by arbitrary non-harmonic excitation

Energy harvesting is employed to extend the life of battery-powered devices, however, demanding applications such as wildlife tracking collars, the operating conditions impose size and weight constraints. They also only provide non-harmonic mechanical motion, which renders much of the existing literature inapplicable, which focuses on harvesting energy from harmonic mechanical sources. As a solution, we propose an energy harvesting architecture that consists of variable number of evenly-spaced magnets, forming a fixed assembly that is free to move through a series of evenly-spaced coils, and is supported by a magnetic spring. We present an electromechanical model for this architecture, and evolutionary optimization process that finds the model parameters which describe the time-domain behaviour observed in ground truth measurements. The resulting model can predict the time-domain behaviour of the energy harvester for any configuration of the proposed architecture and for any mechanical excitation. We also propose an optimization process that, using the electromechanical model, optimizes the energy harvester configuration to maximize the power delivered to a resistive load. The resulting optimized harvester design is specific to the particular kind of non-harmonic mechanical excitation to which it will be exposed. To demonstrate the effectiveness of our proposed model and optimization procedure, we constructed four energy harvesters, each with different configurations, and compared their measured behaviour with that predicted by the model, given an excitation that approximates footstep-like motion. We show that the model predictions were consistently within 25% of the RMS load voltage. We then synthesize an optimal energy harvester using the proposed optimization process. The resulting optimal design was constructed and tested using the same footstep-like excitation, and delivered an average power of 1.526 mW to a 30Ωload. This is a 2.8-fold improvement over an unoptimized reference design. We conclude that our proposed behavioural model and optimization process allows the determination of energy harvester designs that are optimized for a non-harmonic and specific input excitation.


Introduction
Energy harvesting is a well-established approach to supplement and extend the life of battery-powered devices in a variety of operating environments and applications [1][2][3]. This is particularly relevant for devices that have restrictive size and weight constraints, such as wearable electronics [4][5][6] and animal tracking collars [7][8][9][10][11][12]. In the case of wildlife in particular, tracking collars have the additional complication that the replacement of batteries is time-consuming, expensive, and carries particular health risks for the animal, which must be captured before the battery can be replaced [13][14][15]. The size and weight constraints also present an optimization problem: given a particular volume and weight budget, how can harvested energy be maximized? This combination of highly restrictive size and weight constraints and the impracticality of replacing the batteries of wildlife tracking collars was the primary motivation for the energy harvester modelling and optimization we present in this work.
A typical means of harvesting energy is by converting mechanical energy into electrical energy using an electromagnetic transducer [16]. Mechanical energy sources can themselves be divided into two distinct categories: harmonic motion (vibration) [17] and non-harmonic motion [18]. Energy harvesters that rely on harmonic energy sources are by far the most commonly reported in the literature. Since the excitation is harmonic, analytical models with closed-form solutions can be used [19][20][21][22][23]. Optimizing and testing harmonic energy harvesters can also be performed more easily, using simulated and empirical frequency analysis, albeit with the latter requiring carefully calibrated testing equipment [23][24][25][26].
However, the motion exhibited by people or wildlife rarely provides a mechanical excitation that can be considered harmonic [27][28][29]. Since the inspiration for this work is the harvesting of energy from the motion of wildlife, we deliberately focus on energy harvesters that are excited by non-harmonic sources of mechanical energy. In contrast to harmonic energy harvesters, non-harmonic sources of mechanical energy pose a unique set of challenges. Primarily, dynamic non-linear systems that are excited by non-harmonic input do not have closed-form solutions. Instead, they require the numerical solution of differential equations in order to model their behaviour [30][31][32]. Without a closed-form solution, the analytical optimization of design parameters becomes impractical, and hence an alternative approach is required. Although some forms of non-linear optimization have been reported, these suffer from a very limited set of tunable parameters, are entirely specific to a particular, pre-selected input excitation, or are impractical due to their complexity [16,33,34].
We therefore set out with two goals in this work. Firstly, to develop a system model that can predict the behaviour of an energy harvester for an arbitrary input excitation that is not constrained to be harmonic. Secondly, to develop an accompanying optimization process, that improves on the limitations of the current state-of-the-art by being relatively simple to execute, having a larger number of tunable design parameters and allowing the input excitation to be chosen freely.
To achieve this we propose significant extension to a model and optimization process reported in the literature that utilizes steady-state motion as a proxy for the optimization of the actual dynamic non-linear behaviour of the microgenerator device [35]. To do this, we reconsider the same microgenerator architecture, but instead define its behaviour with an extended system model so that we can predict the full dynamic electrical and mechanical behaviour of an electromagnetic energy harvester in the time-domain, given an arbitrary specific choice of input-excitation. We begin by selecting some design parameters, by direct inspection, to conform to real-world constraints. Then, we describe an evolutionary process that allows a separate set of model parameters to be optimized in order to minimize the difference between the model's time-domain predictions and empirical ground truth measurements. Finally, we propose an optimization method that, using the timedomain system model, can optimize the remaining design parameters to maximize the power delivered to a load by the energy harvester. Notably, this optimization can be performed for an arbitrarily specified input excitation, without requiring additional data collection. It also encompasses a larger set of design parameters than previously reported in a single work, encompassing the number of magnets, number of coils, number of vertical and horizontal coil windings, as well as the magnet and coil spacing.
We demonstrate the effectiveness accuracy of our proposed approach by using it to predict the time-domain electrical and mechanical behaviour of a number of reference devices that are driven by a particular example of free-form non-harmonic excitation, and subsequently performing a comparison with the measured behaviour of physically-constructed devices. Finally, we use our proposed optimization process to design an optimal microgenerator, for the same specific non-harmonic input excitation, and show how this optimized design vastly outperforms all other non-optimized reference device under consideration.

Microgenerator design
The kinetic microgenerator consists of one or more copper coils wound around a circular tube, as shown in figure 1. The tube has a wall thickness l th and a total length l L . All coils are identical, with the windings assumed to have a square packing configuration. Thus, we consider each coil as n w windings wide and n z windings tall, leading to a coil width w c and coil height l c . When there are multiple coils, a distance l ccd separates their centres. The windings of adjacent coils are always in opposite directions. Inside the microgenerator body, a fixed magnet attached to the floor of the tube suspends a cylindrical levitating magnet assembly. We use axially-magnetized neodymium permanent magnets with a diameter w m and a height l m throughout. The levitating magnet assembly consists of one or more identical magnets of this type. When this assembly consists of multiple magnets, iron spacers separate the magnets and adjacent magnets have an opposite orientation so that their facing poles match. We place coils in such a way so that, at rest, the magnet assembly is positioned entirely below the lowermost coil. We assume the microgenerator body is sufficiently long for the magnet assembly to be free to move beyond the uppermost coil.
We attach the microgenerator to a mobile fixture, such as an animal's leg. As the fixture moves, the levitating magnet assembly inside is free to move, thus inducing a changing flux in the coil(s) which, in turn, induces an emf. We can harvest this emf and use it to power electronics.

Analytical electromechanical system model
In this section, we first derive an expression for the analytical electromechanical model in terms of its components. We then derive expressions for the behaviour of each component before finally combining the component expressions with the analytical electromechanical model to produce a single set of differential equations that describe the complete analytical electromechanical system.

Equations of motion
Consider figure 2(a), which shows the free-body diagram and electrical circuit of the microgenerator. The force experienced by the levitating magnet (due to the fixed magnet below it) is F mag , the frictional force opposing motion is F fric , the electrical damping force due to the opposing magnetic field (induced by the current flowing through the coil) is F damp , and the force the magnet assembly experiences when it collides with the top of the microgenerator body is F top . The external force acting on the tube and causing it to accelerate is F acc . We can derive an expression for the mechanical component of the system model in terms of these forces. Unless otherwise indicated, all derivatives will be with respect to time.
With reference to figure 1, recall that y mag is the position of the center of the moving magnet (with mass m mag ) and y tube is the position of the microgenerator tube (with mass m tube ). Now let: The microgenerator consists of a tube with total length l L , wall thickness l th , bottom thickness l bth and inner radius r t . We wind one or more coils, each of height l c and width w c , around the outside of the tube. We use the same copper wire to wind all the coils, but the direction of the windings alternates for adjacent coils. Each coil consists of n z windings in the vertical direction and n w windings in the horizontal direction. The distance between the centers of adjacent coils is l ccd . We fix a magnet to the bottom of the inside of the tube. Above this, we place a levitating magnet assembly, with its bottom-most pole matching the upper pole of the fixed magnet. The magnet assembly consists of one or more identical permanent magnets, each with width w m and height l m . When the assembly consists of multiple magnets, we match adjacent poles and separate adjacent magnets with an iron spacer of width w m such that the spacer separates the centers of the constituent magnets by length l mcd . We use an identical magnet for the fixed magnet at the bottom of the tube. When the levitating assembly consists of a single magnet, spacers are not needed. We use y tube as a reference point to the absolute position of the microgenerator body relative to some fixed stationary point at which analysis begins. Similarly, y mag references the absolute position of the center of the lowermost magnet in the magnet assembly, relative to the same fixed point as y tube . then, Note that y diff represents the displacement of the center of the lowermost magnet in the assembly, measured relative to the bottom of the tube, as shown in figure 1. We can further expand equations (3) and (4) by adding an expression for the emf  induced in the coil, which is itself a function of y diff . Let f be the flux linkage between the magnet assembly and the coil, as a function of the relative magnet position y diff , then from Faraday's Law: To describe the analytical electromechanical model, we thus need to find expressions for F mag , F fric , F damp , F top , and f.

3.2.
The microgenerator body length, l L Consider the microgenerator design shown in figure 1, which shows the dimensions of the individual microgenerator components, and figure 3, which shows the position of the magnet assembly at rest and at its maximum height. If the microgenerator body is to allow the magnet assembly to move entirely through all the coils, we must define the height of the microgenerator body l L as: where l bth is the thickness of the bottom of the tube, l th is the thickness of the tube everywhere else, l hover is the distance between the fixed magnet and the magnet assembly at rest,  l is the distance between the magnet assembly and the lowermost coil at rest. It is also the distance between the magnet assembly and the uppermost coil at the magnet assembly's maximum height. Parameters c and m are the number of coils and magnets, Figure 2. Free-body diagram of the mechanical system (a) and circuit diagram of the corresponding electrical system (b). (a) The mechanical system model consists of a levitating magnet with mass m m whose center is at position y m , suspended inside a tube with position y t . F mag is the force experienced by the levitating magnet due to the fixed magnet, F fric is the opposing force experienced due to friction and F damp is the force experienced due the opposing magnetic field induced by the current flowing through the coil within the electrical system. F acc is the external driving force that acts on the tube. F top is the force the magnet experiences when it collides with the top of the microgenerator body. (b) A variable voltage source and two resistors model the electrical system. The open-circuit voltage produced by the levitating magnet is emf. The closed-circuit current is i cc . The internal resistance of the coil is R c and the load resistance is R l . V l is the voltage over the load. respectively, while l ccd and l mcd are the distances between the centers of adjacent coils and centers of adjacent magnets, respectively, and parameters l c and l m are the lengths of each individual coil and magnet, respectively.
From equation (7), note the disproportionate effect of m, l mcd , and l m (which are all parameters related to the magnet assembly) on l L : it is spatially expensive to have a large magnet assembly.
3.3. The magnetic spring force, F mag Due to the pole-matching, the levitating magnet experiences a repulsive force from the fixed magnet at the bottom of the microgenerator body, forming a magnetic spring. We approximate the force experienced by the levitating magnet due to the magnetic spring by a sum of three decaying exponentials, which can approximate the exponential decay inherent between two magnetic poles: where s i is the size of each exponential and is an element of S = [s 1 , s 2 , s 3 ] and a i is the rate of decay of each exponential and is an element of A = [a i , a 2 , a 3 ]. We find the elements of S and A using a curve-fitting procedure described in section 5.
3.4. The mechanical friction force, F fric We approximate the mechanical loss due to friction by modelling it as a constant damper that is proportional to the mass of the magnet assembly: where b m is the damping coefficient, which must be experimentally determined, and m m is the mass of the magnet assembly.
3.5. The rebounding impact force, F top We model the rebounding impact force experienced by the levitating magnet as it encounters the top of the microgenerator body as a joint spring-damper system that is not physically fixed to the levitating magnet. This can be adequately represented by the following relationship: Figure 3. We arrange the coils so that the bottom edge of the lowest coil is a distance  l from the upper edge of the magnet assembly at rest. In this situation the magnet assembly hovers at a distance of l hover above the fixed magnet. We achieve this by deliberately positioning the center of the bottom coil c c to satisfy this condition. Furthermore, the magnet assembly body must be sufficiently long so that the magnet assembly can move through all the coils and, at its maximum position, provide a distance  l between the upper edge of the top coil and the lower edge of the magnet assembly. The maximum possible displacement of the magnet assembly is y diff,max . F H y y Q y y b y 10 where Q is the strength of the spring, y diff,max is the maximum magnet assembly displacement, relative to the tube, before it impacts the top cap of the tube. Note that we relate y diff,max to the length of the microgenerator body with the following expression: Also, in equation (10), b f is the damping constant, which models the energy lost during the impact. We determine the value of b f using a parameter search process, described in section 5. H is the Heaviside step function of the form: which ensures that F top acts on the levitating magnet only when it impacts the top of the microgenerator body.
To model a rebounding impact with the top of the microgenerator body that has little 'compression', we select the spring strength Q to have a large value, Q = 1e7. The larger the spring strength, the less compression occurs, approximating a solid material. Note that a large value of Q does not influence the speed with which the magnet assembly rebounds after impact (this is, instead, determined by b f ), as the law of conservation of energy is not violated.
3.6. The flux linkage, f In order to calculate the electrical damping force F damp and the emf, we need to model the magnetic flux linkage f passing through the coil. The magnetic flux linkage is the amount of magnetic flux, produced by the magnetic field of the magnet assembly, that passes through the windings of the coil. We model the flux linkage as a function of the relative displacement between the magnet assembly and tube y diff , the number of windings in the axial direction n z , the number of windings in the radial direction n w and our particular choice of magnet. We assume a single choice of magnet throughout our described process.
To develop our expression for the flux linkage, we begin by modelling it for a single-coil, single-magnet configuration. We then systematically extend this to a generalized model that we can apply to any number of coils and number of magnets.
3.6.1. Single-coil, single-magnet configuration The simplest configuration of our microgenerator has one coil and one magnet. We have not been able to determine a closed form expression for the true value of f as a function of coil parameters n z and n w , and believe this is not possible because the relationship is dependent on a number of hidden parameters, such as the magnetic field strength vector and magnetization of the magnet assembly, as well as the field's interaction with each winding in the coil. Hence, we will model the flux linkage as a weighted sum of k smooth kernels, with the kernel centers evenly distributed about y y 0, , and sampled over evenly distributed points p ä y diff . Modelling the flux linkage in this way provides two advantages. Firstly, such superpositions of smooth kernels have the ability to approximate a wide variety of smooth functions. Secondly, and more importantly, if we can find the relationship between n z and n w and the kernel approximation (for a particular choice of magnet), we can then predict the flux linkage as a function of n z and n w . We explore this particular concept further in section 6, where we construct a model that can predict the flux linkage from n z and n w for our particular choice of magnet.
Let us define each kernel as a function similar to that of a Gaussian radial basis function. If we define the center of a kernel as k c and recall that we sample each kernel over evenly spaced points p, then we can express an individual kernel as: In equation (13) k c controls the kernel location (which is set differently for each kernel so as to evenly distribute the kernels over y y 0, ), D controls the kernel magnitude and s controls the kernel width. We treat D and s as constants that we select, and thus omit D and s as parameters.
Next we combine the kernels into matrix K with dimensions (k, n p ), where k is the number of kernels and n p is the number of linearly spaced points in p. If we take the dot product of K with weight vector w k , we can approximate the flux curve with the resulting vector Φ: where we assume that our choice of w k produces a weighted sum that approximates a smooth curve with a single peak at p c , which corresponds to the center of the flux curve. We can easily achieve the desired p c by linearly shifting the discrete sampling points p and kernel centers k c of the kernels present in matrix K.
Thus, by carefully manipulating w k , we can approximate different flux curve shapes. We show an example of this with p c = 5 in figure 4.
Additionally, when the levitating magnet's center aligns with the center of the coil, the total flux linkage reaches a peak, since at this position the flux entering and exiting the coil is equal, resulting in a net flux of zero. We can place the peak of f by setting p c equal to the location of the coil's center f(p c = c c ) for the single-coil, single-magnet (c = 1, m = 1) case. This allows us to model the position of the coil, relative to the position of the moving magnet.
We deliberately model the flux linkage as a function of the relative displacement y diff between the magnet assembly and the tube, as this allows us to calculate the induced emf, it's derivative (velocity) and the flux linkage between the magnet and the coil given equation (6). Later, this is further utilized to couple the mechanical and electrical systems.
Note that equation (14) produces a vector of points, while we require a continuous function to later solve the set of equations that describe the electromechanical system. We therefore use a linear interpolation function, Γ, that accepts the resulting weight kernel sum f and a single point p ä y diff as arguments, and returns the interpolated flux value at this point p to define a new function f 0 : Now, given a kernel weight vector w k , and a flux curve center p c , we can use equation (15) to estimate the instantaneous flux linkage given the relative position y diff between the magnet assembly and the coil for the case of c = 1 and m = 1.

Multiple-coil, single-magnet configuration
We now extend the flux model to allow for multiple-coils, but still assume a single magnet in the assembly. As shown in figure 1, we wind adjacent coils in opposite directions and have their centers offset by the distance l ccd . This ensures that, as the magnet assembly moves through the coils, the total flux linkage between the magnet and all coils adds constructively, resulting in an increased induced emf.
Using equation (14) and the knowledge that we space the coils uniformly, we can define an expression for the total flux linkage f c1,m=1 for the multi-coil, single-magnet configuration by adding the flux linkage contributions due to the individual coils for the spacing between the coil centers l ccd : where c c references the position of the center of the lowermost coil, and c is the number of coils in the configuration. We cannot make the choice of coil center distance l ccd arbitrarily. Consider figure 5, which shows the individual flux linkage due to each coil in a microgenerator with two coils and one magnet (c = 2, m = 1). Note how the flux linkage adds constructively when the gradient of the flux linkage for both coils has the same sign or zero for one of them (regions shown in blue). In contrast, destructive superposition occurs when the gradient of the flux linkage in the coils have opposite signs (regions shown in red). As our goal is to maximize the emf  d dt = f , we must maximize the constructive superposition and minimize the destructive superposition by selecting an appropriate value of l ccd . Figure 6 shows the RMS of the induced emf as a function of l ccd for the same microgenerator configuration shown in figure 5, which we calculate by computing the RMS of the derivative of equation (16) (i.e. the emf) with respect to y diff . We see that, when l ccd is too small, RMS is reduced due to destructive superposition. Of particular interest is that an optimal l ccd value exists that results in a maximum RMS. We therefore select this optimal value * l ccd to maximize the RMS of the emf: = is the emf induced across all coils for the multiple-coil, single magnet microgenerator configuration: 3.6.3. Single-coil, multi-magnet configuration Next, we extend our flux model to allow the microgenerator to have multi-magnets, but a single coil. This is similar to the multi-coil, single-magnet configuration described in the previous section. As shown in figure 1, adjacent magnets in the magnet assembly have opposite polarity, and we separate their centers by a distance l mcd . As for the single-magnet multi-coil microgenerator configuration, this ensures that the flux linkage of each magnet with the coil superimposes constructively, resulting in a larger induced emf.  = that maximizes the RMS of the emf  of the two-coil single magnet microgenerator. Any microgenerator where c > 1 should thus use this optimal value of l ccd . When the value of l ccd is too small, it diminishes the resulting RMS due to destructive superposition of the flux linkage. Increasing l ccd , we find the optimal point * l ccd where we maximize the constructive superposition. If we further increase l ccd , no superposition of the flux linkage occurs, causing the RMS to approach a stable value that is lower than that of the maximum RMS.
Recall from figure 1 and equation (1) that y diff is the displacement of the center of the lowermost magnet in the assembly above the bottom of the tube.
We define an expression for the flux linkage between all m magnets and the single coil f c=1,m1 by adding the terms representing the contribution due to each magnet: where c c is the location of the center of the coil, and the spacing between the magnet centers is l mcd . Note here that we subtract the flux curve offset for additional magnets from the location of the coil, since y diff describes the position of the lowermost magnet. This is in contrast to equation (16).
As with the multi-coil, single-magnet configuration, it is again important to select an optimal distance l mcd between the magnet centers, so as to maximize the constructive superposition and minimize the destructive superposition of the flux linkages due to the individual contributions by the magnets. As for the case of multiple coils and a single magnet, selecting a value for l mcd that is too small results in a reduced power output, due to the destructive superposition of the flux linkage. In this particular case, there is also a lower-bound on l mcd due to the practical difficulty of assembling magnets in a pole-matched configuration. We denote this lower-bound as l mcd,min .
By inspecting equations (16) and (19) we see that there is no difference, from the perspective of the resulting flux curve, between the configurations (c 1, m = 1) and (c = 1, m 1). Therefore, as in equation (17), there must again exist an optimal value of * l l mcd mcd = that maximizes the RMS of the emf induced for the single-coil, multi-magnet case: is the emf induced across the coil for the single coil, multi-magnet case:

Multi-coil, multi-magnet configuration
Finally, we derive the flux model for the general case in which the microgenerator has multiple coils as well as multiple magnets. As discussed in section 3.6.3, when a magnet assembly with multiple magnets moves through a single coil, it produces a flux curve described by the superposition of m curves as modelled by equation (19). When there are both multiple coils and multiple magnets, the produced flux curve will be a superposition of c curves, as described by equation (16), where each of these c curves is itself a superposition of m curves, as described by equation (19). Therefore, in this general case, we derive the microgenerator's flux curve from following double summation: From equation (22) we see that both the coil-center distance l ccd and the magnet-center distance l mcd affect the flux curves in the same way. Therefore, their optimal values must be equal: As an example of a multi-coil, multi-magnet flux curve, we show f c=2,m=2 (y diff ), the case of two coils (c = 2) and two magnets (m = 2), in figure 7, together with the four individual flux curves superimposed by equation (22). Note the alternating polarity of each individual flux curve, due to the alternating polarity of adjacent magnets and the alternating direction in which we wind adjacent coils. In this particular example, the optimal coil center * l ccd , and magnet center * l mcd distance are used. Note that, the values of l ccd and l mcd are independent of the number of coils c and number of magnets m. As a result, we only need to find either * l ccd or * l mcd for one of two cases: (c = 2, m = 1) or (c = 1, m = 2), respectively. We can then infer the other distance parameter using equation (23).

3.7.
Coil center distance, c c From figure 3, we can see that it is possible to calculate the coil center distance, c c , from the magnet length l m , the magnet center distance l mcd (in the case of multiple magnets), the magnet assembly hover height at rest l hover , and the chosen distance l ò between the top edge of the magnet assembly and the lower edge of the first coil. This ensures that the magnet assembly has little coupling with the coil(s) at rest, thus increasing the possibility of realizing a large rate of change of flux, and thus emf, when the microgenerator experiences motion. From figure 1 and 3, we can express the coil center distance c c as follows: We can calculate the value of l hover using the magnetic spring force from equation (8) and the weight of the magnet assembly: We must know the coil resistance R c in order to calculate the electrical damping force, F damp , and therefore we must define a model for it. We first model the length of the copper wire used in the coil, l coil . Given the assumed square-packing configuration of the windings and the parameters shown in figure 8, we can express the length of the wire as: Figure 7. Here we show a particular case with c = 2 coils and m = 2 magnets as an example. We sum the four individual flux linkages (first four plots), each given by equation (14), as described by equation (22) to produce the final flux curve profile of the microgenerator (bottom plot).
where n w is the number of windings in the radial direction, n z is the number of windings in the axial direction, l th is the tube wall thickness, r t is the inner radius of the microgenerator tube, and r c is the radius of the copper wire used for the winding. Once we have found l coil , we can calculate the total resistance of the coil: where β is the wire gauge's resistance per unit length, which has a unit of Ω/mm.
3.9. The electrical damping force, F damp As a magnet moves through a coil, it induces an emf. If there is a closed circuit, as shown in figure 2(b), Lenz's Law states that a current will flow in the direction that induces a magnetic field which opposes the magnetic field inducing the current (in this case the magnet assembly). This opposing magnetic field exerts a force on the magnet assembly that is proportional in magnitude to the current flowing through the coil. We therefore model the electrical damping force as directly proportional, with constant b e , to the current i cc flowing through the coil: Given the electrical circuit in figure 2(b), we can calculate the closed-circuit current i cc as: Substituting equations (6) and (31) into (30) gives the final expression for the electrical damping force: 10. The complete system of equations By substituting equations (8) to (32), (22) and (10) into (3), (4) and (6), we obtain a set of equations that models the joint mechanical and electrical system: Figure 8. A detai/t view of a cross section of the microgenerator body and a wound coil. The microgenerator tube has an inner radius of r t and a wall thickness of l th . The wire used for the winding has a radius of r c . The coil consists of n z windings in the vertical direction and n w windings in the radial direction.
Using equations (33) to (35), we can numerically solve for y t , y m , and f as functions of time, given initial conditions and input acceleration force F acc . We can then use the solutions for y t , y m and f to calculate the complete predicted electrical and mechanical behaviour of the microgenerator.
There are three categories of parameters that we require in order to solve equations (33) to (35): a set of parameters that are chosen by inspection, a set of parameters that are optimized to minimize the difference between model predictions and empirically-measured (or simulated) data, and a set of variable parameters that we can optimize to maximize the power that is generated.
The first category consists of parameters that we choose, by inspection, at the outset to conform to practical constraints. These include the specific magnet type and shape (controlled by l m and w m ) used in the magnet assembly and magnetic spring, the inner radius r t of the tube, the tube wall thickness l th , coil wire radius r c (and thus also the coil resistance per unit length β), the load resistance R l , the magnet assembly to coil distance l ò and the microgenerator tube length l L . We discuss the parameters belonging to this category in greater detail in section 4.
The second category of parameters consists of those that describe inherent physical properties of the microgenerator behaviour. These parameters are optimized to minimize the difference between the system model predictions and empirical measurements using an evolutionary process in conjunction with a number of physical reference devices. These include the magnet spring parameters (S, A), the mechanical friction coefficients b m , the electrical damping constant b e , and the impact force damping constant b f . The process of discovering these parameters is discussed in section 5.
The third category of parameters consists of those that can be varied to alter the behaviour of the microgenerator, but without affecting the predictive accuracy of the model. We optimize these parameters to discover microgenerator designs that maximize the power delivered to the load. These include the number of coils c, number of magnets m, coil-center distance l ccd , and magnet-center distance l mcd . The base flux waveform f 0 is a special case that also falls into this category, because it is dependent on the number of coil windings in the axial direction n z and radial direction n w (in addition to the choice of magnet), both of which we optimize to maximize load power. As a result, we give special attention to how we obtain f 0 from n z and n w in section 6.
Although not strictly an optimization parameter, we include the center height of the lowermost coil c c in this parameter category also, since the coil center height c c depends on parameters m and l mcd (equation (24)) and we therefore need to update its value each time we change m and l mcd during the optimization process. We discuss the optimization process for this third category of parameters in section 7.
The benefit of equations (33) to (35) is that, once we determine the first two categories of parameters (chosen parameters, and discovered parameters), we can predict the behaviour of the microgenerator for any input excitation F acc and for any set of parameters that fall into the optimization category. This allows us to optimize the design of the microgenerator for any specific excitation.
It is critical to note that the parameters of the second and third categories (found parameters and optimized parameters) will be specific to the first category of parameters (those that we chose). As a result, if we modify any of the hand-chosen design parameters, we must repeat both the parameter discovery (including acquiring new measurement data) and parameter optimization processes for valid results. This is a direct result of the wide range of inter-dependencies between the parameters that are inherent to this microgenerator architecture, and contributes considerably to the difficulty of optimization.

Parameter selection
The microgenerator design includes a number of parameters that do not describe inherent behaviour in the electromechanical system and that are therefore not part of the parameter discovery or optimization processes, and are chosen arbitrarily. We provide these chosen parameters and their values in table 1, and we consider them as constant values for the remainder of this paper.
We choose identical N38-grade neodymium magnets for the magnet assembly and the fixed magnet at the bottom of the tube. The magnets are cylindrical with a length l m = 10 mm and diameter w m = 10 mm. We therefore select the inner radius of the tube to be r t = 5.5 mm, to provide sufficient clearance for the magnet assembly to move freely. We selected the tube wall thickness to be l th = 1 mm, as this was the minimum that we could achieve using an additive 3D printing process with PLA as the printing material. We select the bottom of the tube to be l bth = 3 mm thick. We selected American Wire Gauage (AWG) 35 copper wire to use for the coils, which has a radius of r c = 0.0715 mm and a resistance per unit length of β = 1079 × 10 −6 Ω/mm. We select the distance between the uppermost magnet edge and the lower edge of the bottom coil to be l ò = 5 mm. Finally, we select our target load to be a constant resistive load of 30Ω.
Note again that all subsequent calculations, such as the parameter discovery in section 5 and parameter optimization in section 7 are entirely specific to this particular choice of parameter values. If the value change, we must repeat all subsequent calculations to achieve correct results.

Parameter discovery
Having chosen and fixed some parameters in section 4, we now turn to describing a process to find those parameter values in equations (33) to (35) that must be determined by measurement or simulation. These are the magnetic levitation force F mag parameters S = [s 1 , s 2 , s 3 ], A = [a 1 , a 2 , a 3 ], the magnet assembly levitation height l hover , which we require to determine the coil center height c c , the mechanical damping coefficient b m that models the mechanical losses due to friction, the constant b e that models the coupling between the electrical damping force F damp and the coil current, and the damping coefficient b f required to model the impact force F top .
We show an overview of the parameter discovery process in figure 9. First, as an independent step, we determine the magnetic levitation force parameters S and A using least squares regression and FEA data. Then, we physically construct a number of reference microgenerator devices, with the same values selected in section 4, and measure each reference devices' behaviour and output experimentally, including the acceleration of the microgenerator tube y tube ¢ . We then use an evolutionary algorithm to find the remaining unknown parameters b m , b e and b f . The algorithm iteratively samples new possible values for b m , b e and b f based on previous sampled values, which we substitute into our model equations (33) to (35) together with the magnetic levitation force parameters S and A. We then solve these parameterized model equations using the measured microgenerator tube accelerations y tube  as the input excitation. Next we update the best parameter values if they produce a lower loss than the current best. We repeat this process for a large number of iterations, resulting in set of parameters that best predict the measured behaviour of the reference microgenerator devices.
Note that, while solving the system equations during the parameter search process, we replace the base flux linkage curve f 0 with flux linkage curves obtained using FEA. This prevents potential inaccuracies in the base flux linkage curve model, as given in equation (15), from propagating through the system, resulting in an inaccurate selection of parameter values.

Magnetic levitation force parameters S and A
We find the magnetic levitation force parameters of S = [s 1 , s 2 , s 3 ] and A = [a 1 , a 2 , a 3 ] by fitting equation (8) to data acquired by FEA simulation for the specific magnet shape and type that we selected in section 4.
We use the software package FEMM [36] to carry out the FEA force simulation. The axisymmetric test setup shown in figure 10. The simulation consists of an axially-magnetized magnet with a diameter of 10 mm and a height of 10 mm placed at a reference position. We specify the magnet material as N38 grade neodymium. The magnet dimensions and material match those of the magnets chosen for the physical microgenerator device. We position an identical magnet directly above the reference magnet, but rotated so that its bottom pole matches the top pole of the reference magnet. We translate top magnet upwards in 1 mm increments (starting at a distance of 0 mm and ending at 60mm) from the reference magnet, and calculate force experienced by the top magnet as a function of the distance between the two magnets. Using this FEA force data, we find the parameters S and A using a least squares optimization [37,38]. For the particular magnet under consideration, we found that the parameter values that produce the best fit are S * = [5.134, 26.291, 0.1256 × 10 6 ], and A * = [0.122, 0.352, 11.609]. A comparison between the FEA simulation and the curve-fit model of the magnetic levitation force in figure 11 shows strong agreement between the proposed model and the results of the FEA simulation.

5.2.
Parameters b m , b e and b f Next, we describe the parameter search process used to find values for the mechanical damper coefficient b m , the electromechanical damper coupling factor b e and the damping coefficient b f of the impact force F top . This process requires empirical, measured ground truth data of the relative magnet assembly position y diff , the induced emf  and the corresponding driving acceleration y tube  , as well as our results for the magnetic spring parameters S and A, discussed in section 5.1.

Reference devices and ground truth data collection
We require ground truth data in order to find parameters b m , b e and b f . To achieve this, we physically construct and assemble n devices = 4 possible reference microgenerator devices (A, B, C and D), each with different sets of chosen parameters as shown in table 2. It is important to have sufficiently many reference devices with differing designs to sufficiently explore and generalize over a wide design space. As a result, we configured the reference devices with substantial differences in their design parameters. For example, each device's coil has a different number of windings N, coil length l c , coil width w c , coil center c c and coil resistance R c . All other parameters match those chosen in section 4 and given in table 1.
We 3D-printed the reference device bodies using PLA through an additive process, and wound the coils using a speed-controlled lathe.
We indicate the experimental setup in figure 12. We mounted a reference microgenerator to one end of a base plate and a camera, capable of high-framerate video, to the opposite end. We also mounted an accelerometer (ADXL345) to the base plate, directly behind the microgenerator. We connected the microgenerator coil to a Schottky diode full-wave bridge rectifier, whose rectified output powers a static resistive electrical load. This load voltage was attenuated using a voltage divider, and connected to the ADC input of a microcontroller (Teensy 3.5). We measured the voltage attenuation factor to be approximately 2.92: 1. Furthermore, each reference microgenerator body was manufactured with an axial groove down the sides, thereby exposing the magnet assembly to the camera. We programmed the microcontroller to record both the accelerometer output and the scaled load voltage, while the camera recorded video footage of the magnet Figure 9. An overview of the process that we use to find parameters that describe inherent physical properties of the microgenerator. Note again that this process is specific to the the second category of parameters, i.e. the parameters that are optimized to minimize the difference between the system model predictions and experimentally-measured or simulated data. First, we use FEA data to find the magnetic spring parameters S and A. In our particular case, we construct four references devices, and take ground truth measurements for the relative magnet assembly position y diff , the induced emf  and the acceleration applied to the microgenerator tube y tube  . We then solve the set of parameterized equations (33) to (35) for new, sampled values of the mechanical damping coefficient b m , the electromechanical damping coupling factor b e and the impact force damping coefficient b f , using an evolutionary algorithm. We use the measured tube acceleration from experimental measurements as the input excitation. After solving the parameterized model equations, we compute the loss between the behaviour predicted by the simulation and the measure behaviour, and update our set of best parameters with the newly-suggested parameters if they result in a lower loss. After a large number of such iterations, we are left with the an optimized parameter set * * assembly moving through the microgenerator body. This experimental setup allowed us to record both the mechanical system behaviour and the electrical system behaviour of the reference microgenerators simultaneously.
To generate a ground truth sample, the camera and microcontroller recorded simultaneously at 120 FPS and 1600 Hz respectively. The base plate was moved in an approximately footstep-like motion for a single cycle. As illustrated in figure 13, the motion consists of an upwards acceleration across a distance of approximately 15 cm, a brief pause, and a subsequent downward acceleration followed by an impact. Figure 13(b) shows a typical such accelerometer measurement in the y-direction. We repeated this process n input = 7 times for each device.
After data collection, we labeled the position of the magnet assembly for each frame in the video recording. We also descale the ADC voltage measurements, using the voltage divider ratio, to calculate the true total load voltage. We use an Savitzky-Golay filter, which is a type of filter that is both simple to implement and easy to tune [39], which was empirically-tuned to reduce noise that was present on the output of the accelerometer.

Search process
The parameter search procedure is shown in figure 14. This process combines the ground truth measurements and our system of model equations to iteratively find the best set of parameter values for b m , b e and b f . The search process begins by sampling new values for parameters b m , b e and b f provided by the (1+1) evolutionary algorithm [40], as implemented in the Nevergrad computational library [41]. We deliberately chose to use an evolutionary algorithm to avoid the combinatorial explosion that occurs when using a exhaustive gridsearch with fine-grained parameter values. The sampled parameter values, along with magnetic spring parameters S and A and fixed parameters from table 1, are then substituted into the system of model equations (33) to (35), for each device. Next, we solve the model equations for each device using each of the Figure 10. The axisymmetric FEMM simulation setup which consists of two identical N38 grade neodymium magnets. We translate the top magnet's position upwards in 1 mm increments (starting at a distance of 0 mm and ending at 60mm) in order to calculate the force experienced by the top magnet as a function of the distance between the two magnets. n acc = 7 corresponding measured acceleration as the input excitation using an ODE numerical solver of type Explicit Runge-Kutta method of order 5(4) [42], implemented in the SciPy library [38]. The solution provides estimates for the relative magnet position y diff and the load voltage V l , the latter of which is then used to calculate a loss relative to the corresponding ground truth measurement of the load voltage V l . We define our loss as the mean normalized load power difference across each reference device and every corresponding input excitation: where  i is the mean normalized load power difference for the ith reference device: Here, P ĵ is the mean power dissipated in the load for the jth input excitation for the ith reference device, and P j is the corresponding ground truth measurement. The mean load power is calculated as follows: After calculation, we return the loss to the evolutionary algorithm, which updates the current best parameter set * * b b , m e and * b f if we find the loss is better than the current best. We repeated this process for n evo = 2000 iterations, in an effort to converge to a parameter set that minimizes the loss function. Figure 11. Comparison of the force experienced by the levitating magnet as predicted by FEA simulation and the proposed magnet spring model, given by equation (8). We observe strong agreement between the FEA simulation and the proposed magnetic spring model.

Search results
The search process found a parameter set of * 3.108 f = that minimized the loss function.
To evaluate these parameter values, we again executed the procedure shown in figure 14, this time using * * b b , m e and * b f , and record the predicted relative magnet displacement y diff and predicted load voltage V l . We then evaluate our model according to a number of criteria.  figure 13(a), after smoothing. In this case, the impact occurred at 4.0 s. Figure 12. The experimental setup used to take ground truth measurements. We attached the microgenerator to one end of the base plate and a camera that can record high framerate video to the other end. An accelerometer was attached to the base plate directly behind the microgenerator. We rectified the output of the microgenerator coil using a Schottky full wave bridge rectifier, which we attached to a static resistive load. We measured the load voltage using the ADC input of the microcontroller, via a resistive voltage divier.
To evaluate temporal accuracy, we calculate the dynamic time warping (DTW) distance between the measured and the computed relative magnet distances and load voltages. This gives an indication of how well we are able to model the mechanical component and electrical component, respectively, of our reference microgenerators.
The DTW algorithm computes the degree to which two temporal signals can be made to align in time [43,44], if we allow a certain degree of stretching and compression of the time axes. As a result, the smaller the DTW distance, the more similar two signals are in terms of their time behaviour. The literature reports DTW as a more robust measure of similarity than Euclidean (and related) distance metrics in multiple applications [45].
In addition to the DTW, we also calculate the normalized RMS difference between the measured and computed load voltages. This measure indicates how accurately the system model predicts the eventual overall power delivered to the load.
To ensure comparable results, we constrain the computed and the measured signals to the same length, and use interpolation to achieve the same fixed sampling rate for both. For DTW, we normalize the calculated distance by dividing it by the number of samples in the ground truth signals. This allows for ground truth measurements of different lengths to be compared. The normalized DTW distance between the computed and measured relative magnet distances is and * b f ) that minimizes  is indicated in a different color and with arrows. We observe that every measure generally decreases as the loss decreases. This is significant, as it indicates that our parameter search process minimizes the loss (which is the difference between the calculated and measured load power as per equation (36)  b f . We substitute these candidate values, in conjunction with the magnetic spring parameters S and A and fixed device parameters from table 1, into the system of model equations equations (33) to (35) for each reference device, of which there are n dev = 4. Using the measured ground truth accelerometer measurements (n acc = 7 per device) as input excitation, we solve each instantiated system of model equations using a numerical solver. The results are then compared to the corresponding ground truth measurements and we calculate the mean loss using equation (36). We compare this loss to the current best, and, if a set of candidate parameters results in a lower loss, we update the current best set of parameters * * b b , m e and * b f . power predictions, but inaccurate time-domain behaviour predictions. Figure 15 indicates that this is not the case: our proposed search process (in conjunction with our proposed model) finds parameter values that lead to accurate time-domain behaviour prediction, which in-turn results in accurate load power estimations. Of interest is that the minimum value of y y DTW , diff diff (ˆ) does not perfectly coincide with the minimum loss . This indicates that, of all the possible parameter sets that were sampled, some lead to more accurate estimation of the mechanical time-domain behaviour, but worse estimation of the load power. This indicates the possible existence of some inaccuracies of the proposed model (possibly due to simplifying assumptions) or in the ground truth data collection process (possibly due to noise), or both. Notably, however, the difference between the parameter sets that lead to the best value of y y DTW , diff diff (ˆ) and the best value of  is small, such that it does not significantly affect our ability to estimate y diff .
We inspect the results for the best parameter set in * (ˆ) for each reference device. We see generally good agreement between the predicted and measured results, indicating that the parameters of * * b b , m e and * b f found in the previous section are both accurate and able to generalize to varied microgenerator designs. In particular, note that the predicted load voltage RMS is almost always within 25% of the measured result. This is noteworthy given the free-form nature of the experimental setup, and suggests that we can be reasonably confident in the model's ability to accurately predict the behaviour and load power of the energy harvester with parameters * * b b , m e and * b f .

Mapping from coil features to flux curves
We now turn to the last component we require before our final power optimization process can take place: we must model the mapping from the coil features to the flux curve defined in equation (14).

Flux curve model
The relationship between coil parameters (n z and n w ), magnet properties and the resulting base flux curve f 0 is not easily determined. The calculation is typically performed using finite element analysis (FEA), which is computationally expensive and therefore not suitable to an iterative application. We propose an alternative method, which only uses FEA initially to generate a training dataset. This dataset is subsequently used to train a linear model for the computation of f 0 that is inexpensive to apply for arbitrary values of n z and n w , including those not in the training set. Recall the model given by equation (14), which approximates the flux curve f by a weighted sum of k kernels across n p points, defined as matrix K, by multiplying with a weight vector w k . Since we fix K, we can consider the elements of weight vector w k as trainable parameters.
If we attempt to make a prediction of the weight vector w k as follows: and * b f ) that minimizes the loss  is indicated. We observe that every measure generally improves as the computed evolutionary loss decreases.
we can define a linear model that does so given some coil features x: This allows us to modify the model in equation (14) to predict the flux curve directly from the coil features: x is an f-dimensional row vector, where f is the number of features, and W w is a matrix of feature weights with dimensions ( f, k).
We can extend this to the general case for n sets of coil features X = [x 1 , x 2 ,...,x n ] to produce a corresponding set of n flux curve predictionsF: where K still has dimensions (n p , k), X has dimensions (n, f ), and W w still has dimensions ( f, k). The resulting outputF has dimensions (n, n p ). We treat the number of kernels k, the shape parameter s and the kernel magnitude D as model hyperparameters. Note that the only trainable component remains the matrix of feature weights W w .
To find W w , we first generate a flux curve dataset with FEA. Using ANSYS Maxwell, we simulate the n flux linkage curves Φ at n p discrete points for a range of y diff values for a single magnet passing through a single coil for all possible combinations of a set of n z values {n z1 , n z2 ,...} and set of n w values {n w1 , n w2 ,...}. We simulate the magnet moving at a constant velocity of 0.35 m/s, starting 80 mm above the center of the coil and stopping 80 mm below the center of the coil. All design parameters that are chosen, remain fixed at the values discussed and chosen in section 4, and shown in table 1. We select the range of windings in the axial direction to be n z = {1, 6,K,200} and the range of windings in the radial direction to be n w = {1, 6,...,200} for a total of n p = 1681 combinations. It is important to note that both the FEA dataset and resulting trained weights W w will be entirely specific to these same chosen parameters, and we must repeat the process if any of these parameters are chosen differently. We use least squares regression to find an optimal set of kernel weights * W k , that maps K to Φ with minimum error: where W k is the trainable feature weights matrix, and where SSE is the element-wise sum of the squared errors for two matrices: We then perform a subsequent least squares regression to determine the feature weights * W W that best map the coil features X to the optimal kernel weights * W k : (43) then produces a linear model that predicts a set of flux curvesF from a set of coil features X: This approach to modelling the flux curves from the coil features eliminates the need for computationallycostly FEA simulation during the optimization process, and replaces it with the FEA synthesis of a training set that is only carried out once, before performing the microgenerator optimization. In this way it decouples the optimization process from the FEA simulations. The result is an end-to-end optimization process that requires measured input excitations and a trained flux curve model as input, and outputs an energy harvester design that has been directly optimized for the chosen design parameters in section 4 and the input excitation in question.

Flux curve model evaluation
We trained the model on n = 840 randomly-selected samples (out of n p = 1861 possible combinations), where each sample is a n z and n w value pair. We used the remainder of the dataset as a testing set. We randomly choose an additional development set of n = 100 samples from the training set, which we use to optimize the hyperparameters with the Nevergrad optimization library [41]. We transform every n z and n w value-pair into the coil features matrix X, where X is a fifth-degree polynomial feature matrix [46], that consists of all polynomial combinations of every n z and n w value-pair, resulting in f = 21 features. The hyperparameter optimization process found that the lowest average mean-absolute-error (MAE) across the development set was obtained for k = 53 kernels, a kernel magnitude of D = 1.883, and a shape parameter of s = 0.00485.
We show a plot of the MAE for every flux curve for small values of n z and n w in figure 17 and for larger values in figure 18. We show the plots separately for clarity. Every flux curve and predicted flux curve sets are normalized together prior to calculating the MAE to account for the large differences in magnitude in the resulting of the flux curves.
The figures show that we observe small errors across almost the entire space, with the test set reporting an average MAE of 0.1821 × 10 −6 . To put this result into perspective, the average MAE is approximately six orders of magnitude smaller than the mean flux in the training set, which is 0.052 87 Wb. Poorer performance is only observed for small values of n z and n w , particularly when both n z and n w 6. This behaviour is also observed in figure 19, which shows a selection of ground truth FEA and predicted flux curves for different values of n z and n w . Here we show that for n z and n w 6 the flux curve model does not make useful predictions, apparently struggling when the magnitude of the flux curve is particularly small. However, for all other cases, the linear model provides an accurate approximation to the FEA curves. Importantly, the model's inaccuracy at particularly small values of n z and n w does not diminish its utility. Coils with n z and n w 6 correspond to small flux curves with small gradients that do not produce a large emf (and therefore power) as per Faraday's law. As a result, coils with n z ä {1, 6} ∧ n w ä {1, 6} can simply be excluded from the optimization without affecting the outcome of the optimization process.

Optimization
Up to this point, we have developed an electromechanical model in section 3, selected our first category of parameters (parameter that we choose by inspection) in section 4, optimized the our second category of parameters (in order to minimize the difference between the system model predictions and the empiricallymeasured data) in section 5, and modeled the relationship between the coil parameters and resulting flux curves in section 6. We now turn to the third category of parameters: those that we optimize in order to maximize the power delivered to the load.
The remaining parameters in this third optimization category are given in table 3. They are: the number of coils c, the number of coil windings in the axial direction n z , the number of coil windings in the radial direction n w , the number of magnets m and, where applicable, the coil center distance l ccd and magnet center distance l mcd . Note again that all other parameters have either been selected in section 4 or found using empirical data in section 5.
As our objective function, we will maximize the mean power that the microgenerator delivers to a load across a number of input accelerations n acc : where n t is the number of timesteps. A visual overview of the optimization process is presented in figure 20. We perform an exhaustive search for each possible combination of parameters c, m, n w , n z and, where applicable, l ccd and l mcd , in order to find a combination of optimization parameters that result in a microgenerator design that delivers the most power to the load.
For each parameter combination, we estimate a base flux curve 0 f from coil features n z and n w using the curve model in equation (47). If multiple coils or multiple magnets (i.e. c > 1 or m > 1) are under consideration, then the optimal coil center distance * l ccd and/or magnet center distance * l mcd is also calculated using equations (17) and (20). Next, we calculate the coil center height (c c ) of the coil (or the lowermost coil when c > 1) using equation (24). This allows us to compile a parameterized equation set by substituting * * l l c m c , , , , ccd mcd c and 0 f as well as the best parameter set * * b b , m e and * b f , found in section 5, into the system model equations equations (33) to (35). The parameterized equation set is then solved for a range of accelerometer measurements (n acc = 15) as input. The accelerometer measurements are taken in an identical way as described in section 5.2.1. The parameterized equation set is solved using the same ODE numerical solver as in section 5.
Each solution produces an estimate for the instantaneous power that the microgenerator will deliver to the load, P load for a particular acceleration input and parameter set. The result is n = 15 solutions per single parameter set, one for each of the accelerometer measurements. We calculate the mean load power P load across this set of solutions, which gives an aggregate estimation of the performance of this particular parameter set and    Table 3. Parameters that are treated as optimizable. Each of these parameters are optimized using the search process presented in figure 20, in order to maximize the mean power P load delivered to the load.
Parameter Description c

Number of coils m
Number of magnets in the magnet assembly n z Number of coil windings in the axial direction per coil n w Number of coil windings in the radial direction per coil l mcd Distance between the centers of each subsequent magnet l ccd Distance between the centers of each subsequent coil the resulting hypothetical device design. We repeat this process for each possible parameter set of c, m, n w and n z . We select the best design by choosing the parameters that deliver the greatest mean power to the load, denoted as * * * * * c m n n l , , , , w z ccd and * l mcd .

Optimization results
We executed the procedure shown in figure  We solve each possible parameter combination for n acc = 15 different accelerometer measurements as input for the acceleration of the microgenerator tube y tube  , resulting in 576 240 solutions total.
We show the mean power delivered to the load P load for each optimization parameter combination in figure 21, separated into four 'bins', each representing a different {c, m} grouping. It is evident that values of n z and n w that are too small result in low power output, since the coil captures little flux. Similarly, for n z and n w values that are too large we also observe low power output, since in this case both the rate of change of flux is small (since a large portion of the magnet assembly's magnetic flux is captured by the coil(s) no matter the assembly's relative position) and the total coil resistance matches poorly the load poorly. Between these two Figure 20. The optimization procedure follows an exhaustive search across all combinations of the values for the parameters n z , n w , c and m that we consider. Using the coil parameters n z and n w , we synthesize a base flux curve 0 f using the linear model in equation (47). If the design involves a multi-coil and/or multi-magnet design, we calculate the optimal spacing between the coils and/or magnets (l ccd and l mcd , respectively). We substitute 0 f , and optionally l ccd and l mcd , into the parameterized system of equations represented by equations (33) to (35) that describe the system behaviour, as well as the optimal values found for the mechanical, electrical and impact damping coefficients * * b b , m e and * b f , respectively in figure 9. The set of parameterized equations is then solved for multiple accelerometer measurements, each solution producing an estimate of the power delivered to the load for the respective input excitation. We select the parameter set that produces the greatest mean power, across all accelerometer inputs, as the best design. extremes, however, there exists an optimum that provides a maximum power output. Furthermore, figure 21 indicates that this optimal point is also the global optimum for a particular {c, m} combination, and that the contour of the resulting load power 'surface' appears smooth. We also see that microgenerators with multiple coils are optimal when coils are shorter (smaller n z ) and narrower (smaller n w ) than microgenerator architectures with fewer coils. This is likely due to the increased resistance resulting of multiple coil configurations. The optimization also appears to prefer coils that are taller than they are wide (i.e. n z > n w ). This is likely due to the particular choice of magnet and its corresponding magnetic field used in our microgenerator architecture.
We show the optimal parameter combination, and corresponding power P load delivered to the load for each {c, m} combination in table 4. We find that the best configuration delivers a mean power P 1.526 10 W load 3 =´to the 30Ωload across all input excitations (n = 15) with parameter combination c = 1, m = 2, n z = 68, n w = 20 and l mcd = 17 mm. Furthermore, we find that the best design delivers slightly more power than for the best c = 2 and m = 2 configuration, suggesting that adding both additional coils and magnets does not always guarantee more power. Figure 21 shows that optimal parameters not only exist, but achieve vastly greater power output than most other designs that we considered in the exhaustive optimization set. As an example, let us consider the c = 1 , m = 1 microgenerator as a 'naive' reference design. Even with the best possible coil parameters for n z and n w , our optimum c = 1, m = 2 devices achieves ≈ 2.8x the load power output of this reference design. Here we note that the 'naive reference' is taken as the optimum c = 1, m = 1 device after optimization. Non-optimal choices of n z , n w leader to even lower output. This demonstrates the benefit of our proposed optimization process. The complete set of design parameters for the optimal microgenerator design is shown in table 5.

Optimized physical device evaluation
Having found our optimal parameters in the previous section, we now physically construct and test the optimal microgenerator design using the experimental setup already shown shown in figure 12 and described in section 5.2.1. By comparing the output of the physically constructed devices with the results computed as part of our optimization we can assess the accuracy of our model equations, as well as the actual power output produced by our optimal microgenerator device. Our physical optimal microgenerator is constructed in the same was as the reference devices, described in section 5.2.1. The microgenerator body is 3D-printed using PLA using an additive process, and the coil's wire wound using a slowly-turning speed-controlled lathe. As before, we include an axial groove into the microgenerator body to expose the magnet assembly to the camera. A photograph of the manufactured optimal microgenerator device is shown in figure 22. =´for a parameter combination of c = 1, m = 2, n z = 68, n w = 20 and l mcd = 17 mm.
Next, we followed the same test procedure we used for the references devices in section 5. Measurements were repeated the experiment for n input = 7 input excitations while measuring the load voltage V l , the acceleration experienced by the microgenerator body y tube  , and the relative position of the magnet assembly y diff . The accuracy of our system of model equations equations (33) to (35) was assessed by calculating the same measures prescribed in section 5, only this time for the optimal microgenerator device. These are the normalized Dynamic Time Warping (DTW) distance between the predicted and measured relative magnet assembly displacement y y DTW , diff diff (ˆ), the normalized DTW distance between the predicted and measured load voltage values V V DTW , load load (ˆ) and the normalized RMS difference between the predicted and measured load voltage The results for the physical optimal microgenerator device are shown in figure 23, alongside the results calculated for the references devices earlier in section 5. The optimal microgenerator labeled as device O and shown in a separate colour. We can see from figure 23 that our system model is able to predict the behaviour of our optimal microgenerator device with similar accuracy as the reference devices. This indicates that our system model is sufficiently robust to generalize to unseen designs.
Additionally, in figure 24, we show the predicted and measured instantaneous relative magnet assembly displacement (y diff and y diff , respectively) in the time domain for the optimal microgenerator device for each of the input excitations. Also here, we observe a high level of agreement, with the model predicting the same dynamics present in the measurements. This corresponds with our earlier aggregate results shown in figure 23. Figure 25 shows the predicted and measured instantaneous load voltage (V load and V load , respectively) for the optimal microgenerator device for the same input excitations. We again observe a high level of agreement. We believe that this is, at least in part, due to the explicit coupling between the mechanical and electrical aspects of the system of model equations. Notably, we accurately predict both the peak values and waveform shapes, when compared to the ground truth measurements.
It is noteworthy that this high level of predictive accuracy is possible given the free-form nature of the experimental process (as discussed in section 5.2.1), which results in a high degree of potential variation in the input excitations. This suggests that the system of model equations are general and robust to some significant Table 4. Optimal combination of parameters n z , n w , l ccd and l mcd for each combination of {c, m} and the corresponding mean power P load , over n acc = 15 input excitations, that is delivered to the 30Ωload. c m n z n w l ccd (mm) l mcd (mm) P load (W)  extent (i.e. not specific to any particular input excitation). The implies is that we could hypothetically apply the process described in this work in the general sense, and use it to design optimal microgenerator devices for arbitrary sources of input excitation. However, despite the positive outlook, further investigation and testing for different applications and microgenerator designs are required to substantiate this claim. Finally, we assess the benefit of the optimization process with regards to an energy harvester's ultimate goal: powering a load. Using our ground truth measurements and predictions, we compute and compare the predicted and measured power dissipated in the 30Ωresistive load (P and P, respectively) for both the reference and optimized microgenerator devices. We calculate the mean load power by taking the RMS of the load voltage over an 8 second period, which encompasses the entire duration of a single input excitation, and then use the RMS to calculate the mean load power using equation (50): We plot the results in figure 26, which shows the box-plot of the predicted and measured mean load power (P and P, respectively) for all reference devices and the optimized microgenerator device. We observe that both the direct measurements and the system of model equations indicate that the optimized microgenerator, O, delivers the most power to the load. What is notable, however, is that the ranking in terms of the power delivered to the load by the system model is identical to that measured in practice.
We observe that the system of model equations slightly underestimates the load power for devices A, B, C and slightly overestimates the load power for devices D and O. The definitive cause of this discrepancy is difficult to pinpoint, and could be due to the simplifying assumptions used in the components in the system model equations (33) to (35), inaccuracies in the empirically-derived parameters A, S, b m , b e and b f from the parameter discovery process described in section 5, or from inaccuracies in the coil curve prediction model equation (43). These inaccuracies can propagate and accumulate throughout the modelling process. Despite this potential for error Figure 22. A photograph of the manufactured physical optimal energy harvester, with design parameters given by table 5. The groove reveals the magnet assembly consisting of m = 2 magnets, separated by a spacer, as well as the fixed magnet at the bottom of the device. accumulation and compounding, the predictions made by the system of model equations remain highly representative of the practical measurements. Figure 26 also demonstrates the practical benefit we derive from the proposed optimization process. The ground truth measurements indicate that, with a measured mean load power of P 1.059 10 W 3 =´-, the optimized microgenerator O delivers ≈ 8.825x more power to the load than device A, ≈ 4.9x more than device B, ≈ 10x more than device C and ≈ 1.94x more than device D. This is noteworthy when we consider that all the microgenerator devices under test are of similar size and volume. We thus observe that the proposed optimization process is able to produce a microgenerator design that delivers 2-10 times more power to the load than a number of 'naive' reference designs, despite having the same volumetric-budget. Given that the optimization itself is specific to the particular input excitation (and thus can target arbitrary input excitation), we can consider the resulting benefit of the proposed optimization as substantial in application.

Conclusion
In this work, we have presented a single-axis microgenerator architecture that consists of a fixed assembly of one or more magnets levitated by a magnetic spring within a tube around which one or more coils are wound. We then proposed a system of model equations that describes the joint mechanical and electrical behaviour of this microgenerator architecture when driven by a particular but arbitrary non-harmonic input excitation. This system model includes a number of different sets of parameters, each determined by means of a different (ˆ) for the previously-seen reference devices (A, B, C, D) and our optimal microgenerator device (O) found by the optimization process. Notably, the system models predicts the behaviour of the optimal device to a similar degree of accuracy as the reference devices. Figure 24. The predicted and measured instantaneous relative magnet assembly displacement for the optimal energy harvesting device across a number of input excitations. We observe a high level of agreement, particularly with regards to the dynamic behaviour. process: the first set is chosen by direct inspection, the second set is optimized to minimize the difference between time-domain predictions and empirical ground truth data, and the third set is optimized in order to maximize the generated power.
The first parameter set, chosen by direct inspection, includes the magnet length l m , magnet diameter w m , the inner tube radius r t , tube wall thickness l th , tube bottom thickness l bth , coil radius thickness r c (and therefore coil resistance per unit length β), coil-magnet gap l ò and the target load R l .
The second parameter set, optimized to minimize the difference between the time-domain prediction of the model and empirical data, includes the mechanical damping coefficient b m , the electromechanical coupling constant b e and the mechanical spring damping constant b f , since these are characteristics of the particular physical construction. To determine these parameters for a particular microgenerator architecture, we proposed  a method that makes use of an evolutionary algorithm to minimize the difference between the system model's predictions and corresponding ground truth measurements for four distinct real-world reference devices. The measurements in question were the acceleration of the input excitation and the subsequent induced mechanical and electrical behaviour of each reference device. We demonstrated this procedure by means of physical freeform excitation that approximates the motion of a human footstep. The evolutionary algorithm samples different candidate values of b m , b e and b f , eventually converging to a set of values that achieves the smallest difference between measured and simulated behaviour. Using the values of b m , b e and b f determined in this way, we analyzed the accuracy of the system model and could report a high degree of agreement between the predicted and experimental behaviour of all four reference devices. The predicted RMS load voltage V load was typically within 25% of the measured value.
Finally, we proposed an optimization process, that can be performed for any input excitation, and that optimizes the final set of parameters c, m, n z and n w to maximize the power delivered to the chosen load. The process consists of a grid search across a range of possible c, m, n z and n w values, where we instantiate the system model using the previously-determined parameters b m , b e and b f , and use this to estimate the power delivered to the load. This procedure includes a method that allows the estimation of a flux curve directly from coil properties, since multiple flux curves are required during the optimization. This estimate was achieved using a linear-weighted sum of discrete kernels, and eliminates the need for computationally-expensive FEA simulations. We solved the fully-instantiated system model for 15 different input excitations and selected the parameter set that produced the greatest average load power across all input excitations as the optimal design.
To practically evaluate our computational method, we constructed this optimal device, and tested it using the same footstep-like excitation. We found that this physical realization of the optimal device provided a mean load power of 1.526 × 10 −3 W over the duration of the excitation to a 30Ω target load. This represents a ≈ 2.8fold increase over the 'naive' single-coil, single-magnet reference device, which served as our baseline. We also found that the system model predicted the behaviour of the optimal device to a similar degree of accuracy as the reference devices, indicating that our system model is a good approximation across a wide range of different energy harvester designs. Moreover, we noted that our system model is able to to correctly rank each device under test according to the power delivered to the load.
We conclude that our proposed system model and associated optimization procedures provide a good approximation to the true mechanical and electrical behaviour of a physically-constructed energy harvester, and that the framework allows the design of an energy harvester providing maximized power output given physical size and weight constraints for specific and non-harmonic input excitation.