{'fraction_non_alphanumeric': 0.026615340038908895, 'fraction_numerical': 0.010099755784593734, 'mean_word_length': 4.879654417133122, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimised graded metamaterials for energy harvesting via reinforcement learning

A reinforcement learning approach to design optimised graded metamaterials for energy harvesting is described. Through the proximal policy optimisation algorithm, the reinforcement agent is trained to optimally set the lengths and the spacing of an array of resonators. The design optimisation problem is formalised in a Markov decision problem by splitting the optimisation procedure into a discrete number of decisions. The a priori physical knowledge of the problem has been exploited to constrain the space of possible configurations by adopting few continuous variables to describe the resonator arrangements. A preliminary analytical investigation has been performed to characterise the dispersive properties of the analysed system by treating it as a locally resonant system. The outcomes of the optimisation procedure confirms the results of previous investigations, highlighting both the validity of the proposed approach and the robustness of the systems of graded resonators when employed for energy harvesting. The role of the resonator spacing is shown to be secondary with respect to the resonator lengths or, in other words, with respect to the oscillation frequencies of the resonators. However, it is also demonstrated that reducing the number of resonators can be advantageous for energy harvesting. The outcomes related to the joint optimisation of the resonator lengths and spacing, thanks also to the adaptive control of the analysis duration, overcome significantly the performance of previously known systems by working almost uniquely on enlarging the time in which the harvester oscillations take place without amplifying these oscillations. The proposed procedure is suitable to be applied to a wide range of design optimisation problems in which the effect of the design choices can be assessed through numerical simulations.


Introduction
In the last two decades metamaterial concepts have witnessed an increasing popularity to control the propagation of waves across much of physics and engineering, with multiple realizations in electromagnetism [1,2,3], acoustics [4,5] and elasticity [6]. In the context of elasticity, considerable effort has been devoted to the investigation of novel mechanisms to manipulate elastic waves for numerous applications of technological relevance such as nondestructive evaluation [7,8], vibration isolation [9,10], seismic protection [11,12,13] and cloaking [14,15,16] to name a few. Concurrently, optimal design solutions have been developed [17], together with advanced modelling methods for complex metamaterials [18,19,20]. Metamaterials are often combined with multi-physics materials, leveraging energy conversion phenomena between mechanical deformations and, for instance, electrical stimuli via piezoelectric coupling. In this context, multifunctional metamaterials have been recently proposed for energy harvesting purposes, thanks to their ability to simultaneously provide vibration isolation and mechanical energy enhancement [21,22,23,24,25]. Whilst local resonators in elastic metamaterials allows for strong energy enhancement, the harvested power is noticeable only close to the bandgap frequency and energy confinement is accompanied by strong scattering effects [26].
To overcome such limitation, a versatile way to obtain inherent broadband and low scattering designs is based on graded metamaterials, which are structures incorporating the gentle variation of resonating elements. The term graded refers to a smooth variation of a particular parameter of the local resonators along space (conventionally the resonance frequency), which enables spatially varying effective properties of the medium. These devices take advantage of local band gaps to control wave propagation; array guided waves slow down as they transverse the array with different frequency components localising at specific spatial positions, resulting in the so-called rainbow effect. Originated in electromagnetism using axially non-uniform, linearly tapered, planar waveguides with cores of negative index material [27], the rainbow effect has been extended to acoustics [28,29,30] and elasticity [31], with multiple realizations for trapping [32] or mode conversion [33]. Graded metamaterials have thus attracted increasing attention due to their ability to manipulate waves by confinement over some spatial region along the structure, enabling wideband vibration attenuation. Within this framework, a number of graded metamaterials have been proposed for energy harvesting, using elastic beams with graded resonators [34,35,32,36].
Even if several works have been done on graded arrays, the definition of an optimal spatial modulation of the medium properties is still an unsolved problem. In acoustics, a stronger sound enhancement has been demonstrated in exponentially chirped crystals rather than in linearly chirped crystals [30]. In elasticity, several grading laws have been compared, and different performances both in terms of vibration isolation and energy harvesting have been reported [37,38,39,40,36]. In this article, a general procedure for the grading optimisation (both in terms of frequency and spacing) is proposed. In doing so, we opt for an elastic beam with lumped resonators, as commonly done in the literature on graded metamaterials. This allows us to define a relatively simple 1D proof of concept framework, which could be further generalised to 2D and 3D. Specifically, the structure under study is a rainbow based metamaterial consisting of a graded array of resonant rods connected to an elastic waveguide [34,35]. As the resonance frequencies are determined by the resonator lengths, we have treated the resonator lengths and spacing as the main design parameters.
A Markov Decision Process (MDP) has been used to formalise the optimal design search as done in [41]. Reinforcement Learning (RL) [42] has been then employed to solve the MDP. The use of RL and of the MDP formalisation concepts have been preferred: (i) to gradient based methods, as these approaches would have been negatively affected by the non linear dependence between the metamaterial harvesting capacities and the design parameters; (ii) to genetic algorithms [43], as they suffer from an high computational cost; (iii) to particle swarm optimisation [44], due to the difficulty of imposing constrains on the design parameters [45]. Moreover, RL enjoys the theoretical advantage of handling possible source of stochasticity affecting the optimal design search. It is our aim to investigate this aspect in future work.
In order to provide information to the RL agent, the Finite Element Method (FEM) has been employed on an Euler-Bernoulli beam with lumped resonators to evaluate the response of different arrangements of resonators. In [46], a RL agent exploited experimental data to discover active control strategies for drag reduction and to gain an insight on the problem, showing that best design solutions often triggers particular physical phenomena. Similarly, in [47] rareevent designs of 3D printed multi-material metamaterials were found out by learning a map from the space of design parameters to the space of mechanical properties through deep learning.
Another aspect of interest is the representation strategy handling the space of the design parameters. Similarly to what was done by [48] in structural shape optimisation, the possible resonator arrangements have been described through the coordinates of few points by exploiting B-spline interpolation. In this way, a large number of configurations and possible modifications of the resonator arrangement could been obtained by playing on a limited number of variables. This approach has been suggested by the physics of the problem. Paying attention to the representation strategy is one of the points of contact between the proposed approach and a family of strategies for automatic system design termed computational design synthesis [49]. A second even more evident similarity is the strong automation of the design process.
In the light of the above, the main contribution of this work is the combination of: the MDP formalization of the design optimisation problem; the constraining of the design space based on physical consideration; the use of an actorcritic RL algorithm for solving the MDP. The potentiality of the method is demonstrated by presenting optimised configurations for a 1D rainbow based metamaterial in which variable resonator lengths or variable resonator spacing, or variable resonator lengths and spacing have been respectively considered. Interestingly, a different number of resonator bars has resulted from the performed optimisations, showing that employing the maximum allowed number of resonators can be useless or even detrimental. Obtained configurations can be interpreted in a physical sense in accordance with previous works [34,35], but revealing aspects of novelty. As it is in the following be shown, the best grading rule privileges the extension of the harvester oscillation time with respect to the amplification of the wave excitation along the guide. The methodology can be straightforwardly applied to other design optimisation problems in which design choices can be evaluated through numerical simulations. A certain robustness to the algorithm hyperparameters, here almost identically employed in the three optimisation cases, further increases the attractiveness of the method.
The reminder of the paper is arranged as follows. The proposed approach is illustrated in Sec. 2, together with an introduction to MDP and RL. The discussion and motivation of the representation describing the possible design arrangements, exploiting few points and B-spline interpolation, is also included. In Sec. 3, the optimisation outcomes concerning the cases featuring variable resonator lengths (Sec. 3.2), and variable resonator lengths and spacing (Sec. 3.3) are presented, together with physics based interpretations justifying their improved performance. These results are preceded in Sec. 3.1 by an insight on the mechanical problem obtained by working out the dispersion relation for the rainbow based metamaterial when represented as a waveguide with resonators. Final considerations on the proposed methodology and on the obtained resonator configurations are contained in Sec. 4, together with a discussion on the future developments of this work.

Methodology
The operation workflow for design optimisation is reported in Fig. 1. While the need of constraining the design space will be addressed in the following, it is considered that if no a priori knowledge of the final design is available, the optimisation of a complex mechanical system is usually split into a sequence of decisions, producing a set of configurations defining specific design descriptions of the systems. In the following, decisions are termed actions, and configurations are named states, while a sequence of actions and configurations leading to a final design is defined trajectory. In particular, the th configuration , together with the corresponding reward , is obtained by modifying  entities interacting with each other, namely the agent and the environment. Actions are said to be taken by an agent, while states and rewards are said to be related to an environment.
For the case at hand, a sketch of this interaction is reported in Fig. 2. A reward measures how well the optimisation task is pursued in a certain state. Here, the state of the mechanical system is defined by the resonator arrangements, while actions modify the resonator lengths and/or spacing. The reward is defined by the sum over the analysis duration of the harvested elastic energy , a quantity strictly related to the electricity obtained by coupling the harvester with a piezoelectric material.
It is now useful to introduce some quantities that will be exploited by the employed actor-critic RL algorithm. First, we observe that the interest is not to maximise the immediate reward at a generic < , but to set a strategy to achieve large rewards at the end of the optimisation process for = . For this reason, the notion of expected return is employed. It consists in evaluating how good is to end up in after actions and then adopting a certain strategy to take the subsequent decisions. For a finite MDP featuring < ∞, the expected return is usually defined as: Second, the concept of policy is introduced. The policy defines the strategy guiding the agent decisions. If a deterministic policy is employed, a unique action is associated to each state; if is stochastic, a Probability Density Function (PDF) over the set of possible actions is returned. An example of PDF is the one related to a Gaussian distribution. A suitable policy allows to solve the MDP by setting a certain sequence of actions. Here, the wanted policy is deterministic, because no uncertainties affect the design at the considered stage. However, during the optimisation it is useful to employ stochastic policies to guarantee the exploration of the action and state spaces.
Last, the concept of value function is discussed. The value function ( ) of a state (here intended as a random variable, while indicates a possible realisation of at ) under a policy is defined by exploiting the notion of expected return as where denotes the expected value under or, in other words, the expected value computed starting from and following thereafter.
By considering as input space the one defined by the combination of the state and action spaces, two functions, namely the action-value function ( , ) and the advantage function ( , ) connected to ( ) are introduced where is an action, here intended as random variable, while is its possible realisation at . Specifically, ( , ) quantifies the expected return of when is taken and then is adopted; ( , ) is the advantage function. Within a MDP, treating ( , ) is often preferred than handling ( ) because the estimation of the advantage function is less affected by variance than the estimation of the value function [50].
The concept of value function is important because by evaluating if ( ) ≥ ̄ ( ) in every state, it is possible to conclude that a policy is better or equal than a policȳ . The policy gradient method aims to learn the best policy * ( | ) [51]. In RL, * is searched by exploiting indications coming from a large number of agent-environment interactions. In contrast, action-value methods do not explicitly look for the best policy focusing, first, on estimating ( ) and, secondly, in reconstructing the best policy by picking up the states with the largest value functions.
In this work, policy gradient methods have been preferred [42]. First, they allow for asymptotically approaching a deterministic policy starting from a stochastic policy. Secondly, they have the capacity of automatically learning appropriate levels of exploration. Moreover, they enjoy the theoretical advantage of expressing the effect of a policy change on the value functions without computing derivatives with respect to the state distribution. This is guaranteed by the policy gradient theorem [52], which is at the basis of the employed RL algorithm. Last, policy gradient methods easily handle continuous action spaces.
The convenience of constraining the design space exploiting physical knowledge is now discussed. If a numeric vector collecting the resonator lengths is adopted to represent the state, the Markovianity of the sequential decision process will be fulfilled as all the information necessary to define the reward and to plan the next action will be contained in the state representation. If we set as action the modification of a single resonator height, and if we  fix the number of resonators and the number of possible levels for the resonator lengths, a discrete number of states will be treated. Discretising the resonator length will not limit the number of obtainable configuration because, on one hand, the discretisation level can be arbitrary increased, while on the other hand the tolerances always present in a manufacturing process agrees with this schematisation. This state definition promises to easily solve the optimisation process, e.g. through a brute force approach consisting in investigating all the states of the system. However, proceeding in this way is generally impossible for the exploding number of combinations. Things do not improve by adopting a RL approach, given that modifying the resonator height one by one does not produce large changes in the rewards, providing weak indications to the agent on how to improve the starting configuration. In other words, if we look at the tackled optimisation task as a gradient based function optimisation, agent explorations corresponds to function evaluations in flat regions.
A more convenient definition of the system state has come from adopting a limited number of continuous variables to describe the possible resonator arrangements. This has been made possible by constraining the design space including the most physically meaningful configurations. From a ML perspective, the variance of the optimisation process has been reduced by introducing a bias based on a previously formed understanding of the problem. The adopted constrain consists in enforcing smooth graded patterns for the resonator lengths in agreement both with theoretical [34] and experimental [35] results. The state continuous variables are used to fix the position of a few interpolation points defining the envelope = ( ) of the resonator lengths, as shown in Fig. 3. Cubic B-splines are exploited for the interpolation. The agent modifies the coordinates of the interpolation points. The sequence of actions is fixed: in the th episode of the agent-environment interaction, the agent always sets the position of a certain interpolation point. Thanks to the concept of expected return, defining a priori an order of operation does not affect the optimisation outcome. The effectiveness of the proposed state and action definitions will be further assessed by showing that increasing the number of state variables does not necessary benefit the optimisation procedure.
To treat continuous policy and action spaces, the state value functions have been approximated by parameterised functions as where ∈ ℝ is a vector of tunable weights. Similarly, associating a PDF featuring a Gaussian distribution to the policy, deterministic function approximators have been used to parametrise the dependence on the state of the policy mean and standard deviation as where here is just the number ≈ 3.1415 …; and are two function approximators relying on the parameter vector ∈ ℝ . [m] [m] The way in which the MDP is solved by RL is now explicitly addressed. In particular, the Proximal Policy Optimisation (PPO) algorithm [53] has been employed to approximate the optimal * finally solving to the MDP by tuning ∈ ℝ . The PPO algorithm belongs to a subset of policy gradient methods named actor-critic approaches. In these methods, state value functions ( for PPO) are used to assign credit to the agent actions. This requires the algorithm both to approximate the advantage function ( , ) ≈ , , and to calibrate the policy parameters. Two fully connected Neural Networks (NNs), ruled respectively by and , have been used for this goal. The NN differentiability has been exploited within the backpropagation algorithm [54] to maximise the PPO objective by tuning via gradient ascend with Adam [55], where: is an hyperparameter usually set to 0.2;̂ is the empirical mean of the advantage function over trajectories; , , is the advantage function based on the current critic; old | , old is the policy run to collect the trajectories. "min" requires to take the minimum (pessimistic bound) between a clipped and an unclipped objective. The unclipped objective is the one that, under a certain constrain, is minimised by the trust region policy optimisation algorithm [56] that represents the basis of PPO.
Importance sampling is employed due to the difference between the updated policy and the one used to generate trajectories. The clipping operation "clip" removes the incentive for changing rapidly the policy making the ratio between | , and old | , old moving outside the [1 − , 1 + ] interval. In the following, this ratio will be indicated by . The clipping operation allows the definition of the following piecewise probability distribution for , , > 0 and < 1+ , or , , < 0 and > 1− , whose empirical mean is required to compute  . The need of cyclically update , , points out the actor-critic scheme of the PPO. The parameters modelling the critic are tuned via gradient descend by minimizing a loss function defined as Algorithm 1 reports the complete PPO algorithm employed to solve the MDP with which the design process is formalised. Rewards are computed through FE simulations. Additional details on PPO can be found in [56]. As it can be noted, a major element of interest of RL with respect to supervised learning is the possibility to collect data during the training, potentially online.

Results
Aim of the work is to show how RL can be applied to metamaterial-based energy harvesters, highlighting the potentialities of the method as a design optimisation tool in mechanics. Differently from the supervised and the unsupervised paradigms, RL does not require to assemble a data set a priori. Agent-environment interactions are collected on the fly by simulating wave propagation in the system through the Finite Element Method (FEM).
Both the resonators and the waveguide, depicted in where is the resonator length. The metastructure is excited through the wave packet depicted in Fig. 4. The frequency content of the wave packet coincides with the angular frequency of the harvester first bending mode ℎ = 17.67 MHz, obtained by substituting the length ℎ = 5.028 ⋅ 10 −4 in Eq. 9 to . The same type of excitation was provided through an electrodynamic shaker in an experimental setting by [35]. The wave number = 3470 rad/m is determined through the following dispersion relation by setting = and ℎ = , see [26]. Consequently, the wave length is equal to 1.8 ⋅ 10 −3 m. To avoid wave reflections, two absorbing layers with length abs = 5 = 9.1 ⋅ 10 −3 m have been implemented at the extremities of the waveguide as in [57]. Comparable absorbing properties were experimentally obtained through acoustic black holes [58,59,35]. Due to the presence of the absorbing layers, the exciting force has been applied at = 9.1 ⋅ 10 −3 m. Different FE discretisations have been employed for the considered optimisation cases. However, every discretisation has employed five FEs in between  Table 1 Metasurface mechanical data.
neighbouring resonators to catch localised effects due to the resonator interaction. Resonators have been modelled using a lumped spring-mass representation. Wave propagation has been simulated for = 1.25 ⋅ 10 −5 s for the first two analyzed optimisation setting. The second case, involving the grading of both the resonator lengths and spacing, will require further discussions to set . A Rayleigh damping featuring a ratio roughly of 0.2% for the first 40 structural modes has been considered. Mechanical data are collected in Tab. 1.
The input force generates two wave fronts, one propagating left (towards points with smaller ), and the other propagating right. The left propagating wave is immediately damped out by the absorbing layer, while the other wave propagates until it reaches the second absorbing layer. The wave travels with group velocity equal to 1.018 ⋅ 10 4 m∕s determined by substituting Eq. 10 in = ∕ .
In Sec. 3.1, the dispersion relation for a local resonant waveguide is analytically derived to get a mechanical insight on the optimisation problem. Outcomes concerning the optimisation procedure are reported: in Sec. 3.2 with respect to the grading of the resonator lengths; in Sec. 3.3, with respect to the resonator lengths and spacing. For sake of clarity, interpolation points will be differently labelled: with -IP ("IP" refers to interpolation point), when they refer to the resonator lengths; with -IP, when they set the resonator spacing.

Dispersion relation for 1D rainbow-based metamaterials
A dispersion relation expresses the dependence between and for a wave propagating in a certain periodic medium. Provided the grading is gentle enough, the global behaviour of a rainbow metamaterial is deduced from the local dispersion curves of the constituent elements [31]. This allows to compute the dispersion relation for a transverse wave in the rainbow metamaterial by formulating the infinitesimal equilibrium equation of the periodic problem, according to the Euler-Bernoulli beam theory where ( , ) is the waveguide displacement field, and by employing a lumped mass-stiffness model to include the presence of resonators.
To avoid boundary effects, the length of the waveguide has been assumed to be infinite. By assuming that the motion of the resonators is dependent only on the motion of the underlying waveguide, it is possible to write: where: ( ) is the displacement of the th mass of the local resonance system; is the lumped stiffness of the spring that connects the lumped mass with the underlying waveguide; is the distance between two resonators; is the Dirac delta.
To solve the system, the solution of the equation of motion of the th resonator is needed. This is formulated as follows: granted that the stiffness = , the mass = , and the distance are assumed the same, thanks to the locally periodic behaviour of the graded metamaterial.
In this work, = , while = 2 , where has been computed as in Eq. (9). In case of identical resonators, the solution obtained by [60] is: where = 2 is the reciprocal lattice vector, and is just the number.
The summation in Eq. (14) is needed to model the Bragg effect that develops near to the boundary of the First Brillouin Zone (FBZ). If the local resonance band gap is sufficiently separated from the Bragg one, i.e. they arise at different frequencies, this effect can be neglected as it marginally contributes in defining the shape of the dispersion relation, leading to simplify Eq. (14) as In turn, this equation can be reshaped to show the dependence of over the lattice spacing, obtaining: for a fixed frequency , stiffness and mass of the resonators. It is interesting to see that  contains all the contributions that define the dispersion relation for the infinite beam, while  contains the contribution brought by the resonators. The distance is kept out of  to show how it affects .
Eq. (16) highlights which are the design parameters to be optimised for a proper grading: the resonator stiffness , appearing in the formulation of ; the distance in between the resonators. From a geometric point of view, setting is allowed by setting the resonator length . In light of that: Sec. 3.2 will investigate how to optimally grade the resonator lengths or spacing ; Sec. 3.3 how to optimally combined both and

Optimisation of the resonator lengths or spacing
The optimisation procedure detailed in Sec. 2 has been applied to set the lengths, and therefore the lumped mass and stiffness properties, of 25 resonators. Resonators have been spaced by roughly ∕11, within the range ∕15, ∕2 for which resonator interaction is expected at the subwavelength scale [61]. The precise value depends on the employed discretisation. Specifically, the waveguide has been discretised through 376 Euler Bernoulli FEs [62] with length equal to 0.0344 ⋅ 10 −3 m where the resonators are placed, and equal to 0.344⋅10 −3 m outside that zone. The integration time step Δ has been set to 3 ⋅ 10 −9 s. The harvester coincides with the 18th resonator from the left. Its height has been fixed to ℎ .
Four points have been used to define the envelope = ( ) of the resonator lengths. In Fig. 3, the four points have been depicted by star markers. Specifically, the 3rd -IP has always been placed at the tip of the harvester, while the 1st and 4th points have been forced to match the tips of the first and of the last resonators, respectively. Therefore, = 4 continuous variables have been exploited, coinciding in the -plane with the ordinates of the 1st and 4th -IPs, and with the ( , ) coordinates of the 2nd -IP. The agent has been asked for setting these continuous variables by taking the following sequence of = 4 actions: (i) set the coordinate of the 1st -IP; (ii) set the coordinate of the 4th -IP; (iii) set the abscissa and then the ordinate (iv) of  Table 2 Optimisation of the resonator lengths: description of the agent actions. The IPs are the ones reported with cyan markers in Fig. 3.
the 2nd -IP. The order of operations is exemplified in Fig.  5 and schematically reported in Tab. 2.
The configuration space has been constrained on the basis of the physical understanding of the system acquired in [34,35]. No other a priori knowledge has been exploited. For example, no notion has been employed to set the system configuration at the start of each trajectory. Indeed, all resonators have been set equal to ℎ at each trajectory start, as depicted in the top-left plot in Fig. 5.
The -IP ordinates have been allowed for varying in the 0, max interval, where max ℎ = 9.156 ⋅ 10 −4 m. The rationale behind this value is that if a sinusoidal force with frequency equal to max is used to excite a single degree of freedom system centered at ℎ , a 10% attenuation will be produced, with max computed as in Eq. (9) by setting = max . The coordinate of the second -IP have been allowed to move between the first resonator position and the harvester position. A description of the agent actions is reported in Tab. 2. Cubic B-splines have been exploited to define ( ). If a resonator length results smaller than ℎ ∕20, this resonator is removed from the design. This threshold value has been set to avoid the ill conditioning of the stiffness matrix caused by very high stiffness of the lumped resonators. In this way, the procedure can possibly reduce the number of resonators.
A NN featuring two fully connected layers with just 32 neurons and = 1348 tunable weights has been employed to approximate the advantage function where  − ∶ ℝ → ℝ is the transformation induced by a fully connected layer receiving a vector of length in input and returning a vector with length in output. A similar NN architecture, exploiting = 1282 parameters, has been used to calibrate the policy parameters Performance has not been improved either by adding hidden layers or by increasing the number of neurons. A number = 32 of trajectoris have been collected for each estimate of the PPO objective  expressed in Eq. (6). It has been verified that a smaller deteriorates the performance of the  Table 3 Optimisation of the resonator lengths: PPO algorithm approximating functions and hyperparameters.
procedure, and that a larger increases the computational effort without improving the procedure outcome. The total number of agent-environment interactions has been fixed to = 100, 000 FE simulations. As it will be assessed by Fig.  6, this number has guaranteed the convergence to a quasideterministic policy obtained by greatly reducing , with respect to , for some state . The learning rate has been set to 5 ⋅ 10 −4 after a trial and error procedure. The employed function approximators and hyperparameters are collected in Tab. 3. The stable baselines library [63] has been exploited for the implementation of the PPO algorithm.
In Fig. 5, the sequence of actions leading to the configuration with the highest reward is reported. As it can be observed, the 2nd -IP has been moved very close to the 1st -IP, so much that the resulting curve could have been obtained by interpolating 3 points only. This further confirms that the adopted state description has reduced the space of possible designs but including the most effective configurations for energy harvesting. The convergence to a deterministic policy is assessed in Fig. 6, where the evolution of the episode final reward during the agent training is depicted. The graph ordinate is normalised with respect to the reward collected for the waveguide equipped with the harvester only, and it is plotted against the number = 1, … , of agent-environment interactions. The ratio ∕ quantifies the advantage of adopting a certain metastructure design with respect to the configuration featuring a single resonator. The deterministic policy obtained after agent-environment interactions generates a configuration, depicted in Fig. 7d, featuring a ∕ ratio equal to 3.537. This ratio is slightly smaller (by ≈ 3.5%) than the ∕ = 3.669 ratio associated to the configuration obtained through the action sequence depicted in Fig. 5. In other words, the best generated configuration has been produced before the total number of policy updates have been performed. The performance of the optimisation procedure is evaluated not just by comparing with , but also with the reward featuring the linear grading rule shown in Fig. 7b, as proposed in [34].
As it can be seen, the best RL discovered configuration largely enhances the harvesting capability with respect to . This enhancement has been obtained constraining the design space as unique a priori knowledge guiding the agent actions. The RL agent has been able not only to recover what previously known, but also to overcome the performance of  the linear grading resonator arrangement by 4.71%. Despite the limited increase, no further major improvements are expected by playing only on tuning the resonator lengths, as also shown by the work of [38] when comparing different linear and parabolic grading rules. The lower performance of the finally obtained deterministic policy with respect to the best discovered configuration does not mean that the MDP solution has been unsuccessful or meaningless. On the contrary, the utility of the proposed approach has been confirmed: by the convergence of the PPO algorithm to a suboptimal policy generating a resonator arrangement close in terms of ∕ to the best configuration found during the algorithm explorations; by the fact that the best discovered configuration has been produced after the reduction of the space of possible configurations due to the first policy updates. Moreover, arriving to a suboptimal policy is the usual outcome of a RL procedure unless of a dramatic increasing in computing time [42].
Some theoretical arguments discussing how to theoretically find an optimal grading are reported to further assess the outcome of the optimisation procedure. A possible way to approach the problem is to think that the propagating energy of the wave is progressively reflected by the resonators. The matter becomes to define the intensity of such reflection and to find a grading law that minimizes them. With this goal, we consider the waveguide as made of a series of layers of different homogeneous material whose properties are given by the dispersion relation stated in Eq. 14. Therefore, it is possible to define a reflection coefficient between each layer, so between each resonator cell and the preceding one, in the following way: and a transmission coefficient that is: Where is the wavevector associated to the th layer, and a resonator cell is defined by a resonator and by the attached waveguide portion. These equations are the same used in optics, and they stem from imposing both displacements and force continuity at the interfaces between the layers. Physically, we are interested in the square of the transmission coefficients from one layer to the next, namely transmittance and reflectance. In these terms, the optimisation problem can be thought as finding which sequence of layers maximises the wave intensity on the desired final layer of the target resonator, i.e the one that maximises the product of all the layer transmittance. The best RL configuration (Fig. 7c) partially fulfills such condition, with an average variation of Δ ≃ 174 rad∕m for the initial portion of the array. On the other hand, this approach fails close to the target resonator, since it is not able to predict the real wavefield generated by the multiple scattering and emissions of the resonators. In particular, as explained in the last section, the effectiveness of rainbow metastructures is strongly related to the interaction time between the wave and the harvester. This means that, the maximum peak power provided by a proper tuning of the transmittance, may not be the optimal solution to harvest energy in a given period of time. The deterministic context in which the optimisation problem has been set allows for considering the best configuration discovered during the exploration phase as the outcome of the procedure. As the number of agent-environment interactions run up that point has been quite limited (≈ 5000), the method seems suitable to investigate mechanical problems involving more time demanding FE simulations. Relying on the configuration related to the policy finally determined by the PPO algorithm would have been necessary if any source of stochasticity, e.g. related to the loading conditions, had been included in the FE simulations. In that case, the best configuration could have not be generated by the sequence of actions leading to the highest during the policy exploration phase, as that value of may have been connected to special loading conditions.
The interest is now directed to optimise the spacing in between the resonators. Eq. (16) directly shows the dependence  Table 4 Optimisation of the resonator spacing: description of the agent actions. The -IPs are reported with cyan markers in Fig. 8a.
of the wavevector with respect to just a change in distance in the lateral resonators. This is also true considering that an enlargement of the cell is reducing the influence of the resonators over the dispersion relation, while smaller spacing between resonators straightens the contribution of .
A RL agent has been asked for optimally ruling the resonator spacing for energy harvesting. The minimum spacing has been set to ∕16, the maximum spacing to ∕2. A discrete number of spacing levels ( ∕16, ∕8, 3 ∕16 …) has been allowed by increasing the minimum distance by ∕16. A number = 6 of continuous variables setting the coordinates of four -IPs in the ( , ) plane and ruling the envelope = ( ) has been employed to constrain the design space. The coordinates have been allowed to vary continuously, operating the discretisation of the spacing levels in a second time. Tab. 4 resumes the continuous variables possible range of variations and starting values clarifying the ordering of the agent actions. Concerning the resonator lengths, a linear grading rule has been assumed. The same tuning parameters employed for the PPO algorithm in Sec. 3.2 have been exploited, still modelling , , and calibrating the policy parameters , , , through two different NNs featuring two hidden layers with 32 neurons.
The best configuration discovered by the RL agent is reported on the top of Fig. 8. The number of resonators has been reduced from 35 to 23 by enlarging the spacing of the resonators placed before the harvester. On the bottom, a spectrogram representation [64] has been employed to depict the outcome of the wave propagation obtained for this resonator arrangement. The predictions of Eq. 16 agree well with the outcomes of the FE simulation even by simultaneously modifying the distance and the ratio ∕ .
The ratio ∕ quantifying the performance of the optimisation procedure is equal to 3.600. It implies just a 1.5% performance improvement with respect to the case featuring a linear grading rule with 35 resonators. This can can be addressed considering from Eq.14 that a change in height gives a more significant change in the wavevector with respect to a change in distance. Therefore, both analytical arguments and the outcomes of the optimisation procedure support the conclusion that playing on the resonator spacing  is not as effective as playing on the resonator lengths for sake of energy harvesting. However, this investigation has also revealed that employing the maximum allowed number of resonators is unnecessary for getting the best performance, therefore promising savings in the material and in the manufacturing costs.

Optimisation by combining resonator lengths and spacing
The joint optimisation of the resonator lengths and spacing has been tackled considering for the waveguide the same geometry of Sec. 3.2. The agent actions reported in Tab. 5 have been used to modify the value of = 10 continuous variables by operating on: the coordinates of four -IPs belonging to the ( , ) plane and ruling the envelope = ( ) of the resonator lengths; the coordinates of four -IPs belonging to the ( , ) plane and ruling the envelope  Table 5 Optimisation of the resonator lengths and spacing: description of the agent actions. The -IPs are reported with cyan markers in Fig.; the -IPs with orange markers in  , , and to calibrate the policy parameters , , , , respectively, have been doubled to account for the larger dimension of the configuration space. The learning rate has been set to 1 ⋅ 10 −3 to speed up the training in view of the greater number of actions performed by the agent.
A first optimisation attempt in which the analysis time has been set to 1.5 ⋅ 10 −5 has led to contradictory outcomes highlighting an unusual decreasing profile of the lengths of the resonators preceding the harvester, as shown in Fig. 9a. This result disagrees both with what reported in Sec. 3.2 and with the literature, see [34,35]. It has been necessary to inspect the evolution of the harvester  to explain the reason behind this unexpected outcome. As shown in Fig. 9d, the analysis has been too short to damp out the oscillations of the harvester when the linear grading rule of Fig. 9c has been employed. Consequently, the agent discovered a resonator configuration capable of exciting stronger but for a shorter time the harvester. Computing the ratio * ∕ * for = 2.0 ⋅ 10 −5 for the obtained configuration and for a resonator arrangement ruled by a linear grading confirms what has been suggested, highlighting the better performance of the linear grading rule.
It must be said that previously reported outcomes have not been affected by the adopted , as in Sec. 3.2 was set to a sufficiently large value, while the effect of the resonator spacing on the harvesting capabilities of the system has been shown to be limited in any case. However, the need of handling the possibly different time interval necessary to damp out the oscillations of the harvester, and the possibly different interaction time between the propagating wave and the resonators justifies the need of employing an adaptive criterion to set . On light of that, the next analyses have   been stopped only when the  has been lower than 5 ⋅ 10 −20 J for at least 500 times in row. The count has been prescribed to start when the threshold value has been overcome for the first time.
The outcome of the optimisation procedure for an adaptive is reported in Fig. 10a. As in Sec. 3.2, the lengths of the resonators preceding the harvester monotonically increases and features a negative concavity ′′ ( ). However, the slope of the grading is much lower than the one depicted in Fig. 7c. Moreover, working on the spacing has enabled to take out some resonators from the guide. Despite the limited discrepancies, the mechanism exploited by the obtained configuration to increase  during analysis is largely different from the one featuring the linear grading of Fig.  9c. This reference resonator arrangement increases the time interval Δ lin  in which the harvester oscillates with respect to a configuration featuring the harvester without other resonators, but primarily focuses on amplifying the harvester oscillations. On the other hand, the obtained resonator array works almost uniquely enlarging the oscillation time interval Δ RL  , so that Δ RL  > Δ lin  as clearly shown in Fig.  10b. finally leading to enhance by 12.8% the elastic energy exploitable for energy harvesting. This increase has been possible despite harvester oscillations even smaller than the ones related to the configuration featuring just the harvester without other resonators. The time-spatial representation of 10c further elucidates the matter by depicting how the wave packet propagates in the waveguide for the two cases. For both configurations, the wave travelling to the left has been immediately damped out by the absorbing conditions, while the wave propagating to the right has been converted in resonator oscillations. After a while, resonator oscillations have turned back to a wave packet left propagating in the waveguide. The oscillation conversion has taken place after Δ lin  for the blue curve depicting the wave propagation for the linear grading configuration, and after Δ RL  for the red curve referring to the optimised grading. The inequality Δ RL  > Δ lin  has been allowed by the repeated wave reflections between the harvester and the neighbouring resonators, as highlighted by the several  picks reported in red in Fig.  10b.

Conclusions
A reinforcement learning approach to optimise graded metamaterials for energy harvesting has been proposed. The procedure leverages on constraining the design space on the basis of the knowledge of the physical behaviour of the system, and on the use of a policy gradient actor-critic algorithm to train the reinforcement learning algorithm. Results related to the optimisation of the grading rule of the resonator lengths and spacing highlight the potentiality of the method. The procedure can be applied to to other design optimisation problems in which the effect of a design choice can be numerically simulated. A preliminary analytical investigation has been carried out by obtaining the dispersion relation for a local resonant metamaterial considering both variations in the resonator oscillating frequencies and spacing. The predictions of this dispersion relation have been found in agreement with numerical outcomes.
Specifically, results on the resonator lengths optimisation have confirmed previous works reported in the literature. The procedure has also highlighted the robustness for energy harvesting of this type of mechanical system with respect to limited changes in the resonator configuration.
Results on the resonator spacing have demonstrated the secondary effect of this design parameter with respect to the resonator lengths, however showing the possibility of reducing the number of resonator in the system still obtaining comparable, if not better, performances. Finally, the outcomes of the optimisation procedure when applied to the joint calibration of the resonator lengths and spacing with an adaptive rule to set the analysis duration has further highlighted the potentiality of the method by producing a significant increase in the harvesting capability with respect to a reference optimised configuration featuring a linear grading. The obtained configuration has exploited a harvesting mechanism different from this reference configuration by focusing almost uniquely on enlarging the time in which the harvester oscillations take place without amplifying them. These outcomes have confirmed the fundamental importance that the analysis duration plays on the optimisation outcomes.
In future works, the proposed reinforcement learning procedure empowered by physically based constrains of the design space will be applied to 2D and 3D resonator systems. Attention will be also devoted to design optimisation problems involving sources of uncertainty, for example the one related to the frequency content of the wave excitation, as well as the possibility of harvest energy on the entire set of resonators. In blue: evolution of  for the reference configuration ∕ = 3.870. In red: evolution of  for the best RL discovered configuration. ∕ = 4.365. (c) Wave propagation time-spatial outcomes. Blue scale: outcomes for the reference configuration (see Fig. 9c. Red scale: outcomes for the best RL discovered configuration.