Differential Evolution in Waveform Design for Wireless Power Transfer

: The technique of transmitting multi-tone signals in a radiative Wireless Power Transfer (WPT) system can signiﬁcantly increase its end-to-end power efﬁciency. The optimization problem in this system is to tune the transmission according to the receiver rectenna’s nonlinear behavior and the Channel State Information (CSI). This is a non-convex problem that has been previously addressed by Sequential Convex Programming (SCP) algorithms. Nonetheless, SCP algorithms do not always attain globally optimal solutions. To this end, in this paper, we evaluate a set of Evolutionary Algorithms (EAs) with several characteristics. The performance of the optimized multi-tone transmission signals in a WPT system is assessed by means of numerical simulations, utilizing a simpliﬁed Single Input Single Output (SISO) model. From the model evaluation, we can deduce that EAs can be successfully applied to the waveform design optimization problem. Moreover, from the presented results, we can derive that EAs can obtain the optimal solutions in the tested cases.


Introduction
A lot of research has been recently carried out on the Internet of Things (IoT) concept. IoT is an evolutionary approach that is promising to make a great impact on everyday life and business activities. It is believed that its wide deployment will also bring economic growth, but this requires the advancement of various enabling technologies [1]. One of these technologies is Radio Frequency (RF) Wireless Power Transfer (WPT), which is used for powering up low-power devices without the need of batteries or cables [2][3][4].
In a traditional radiative WPT system, a transmitter device generates a continuous electromagnetic field, which transmits power across the medium to a receiver device (rectenna). However, the transmission of pulsed electromagnetic fields that demonstrate higher peak-to-average power ratio (PAPR) has proved to be advantageous. Pulsed wave signals are multi-tone signals that exploit the nonlinear behavior of the receiver's rectifying circuit and increase its RF-DC conversion efficiency [5,6].
When the propagation between the transmitter and the receiver is characterized by a multipath channel function, Channel State Information (CSI) should be taken into account for the design of multi-tone signals. Towards this direction, multi-tone waveforms constructed by Time Reversal (TR) processing are proposed for WPT in [7,8]. TR waves are designed in a predetermined way, according to the CSI, and are capable of focusing on the receiver with the highest peak voltage. Consequently, they exploit both the multipath propagation channel and the nonlinear behavior of the receiver's The paper is organized as follows. In Section 2, we describe the WPT system model, including the structure of the multi-tone waveforms, the propagation channels, the rectenna circuit model, and the selected objective function. Section 3 depicts the tested EAs and the modifications needed for their application in this optimization problem. Finally, In Sections 4 and 5, we outline the derived results and the conclusions of the paper, accordingly.

Wireless Power Transfer System Model
We consider a Single Input Single Output (SISO) WPT system model. This model corresponds to the simplest utilized WPT system, consisting of a transmitter, a multipath propagation channel, and a receiver rectenna (Figure 1). The wireless transfer of the power is obtained by generating and transmitting an RF signal which, after propagating across the medium, is received and rectified by the receiver rectenna (rectifier + antenna) to supply an electrical load. Hence, the objective is to design a waveform s(t), such that the DC power at the output of the rectenna (P out ) is maximized.

Waveforms Structure
The waveform signals that have been utilized in this paper incorporate multiple sub-carriers evenly-spaced in the frequency domain. These signals are pulsed waves of a periodic form with a pulse period of T o = 1 ∆ f , considering that the frequency space between two adjacent tones is ∆ f . Each transmitted signal can be expressed in the time domain by: s(t) = N ∑ n=1 a n e jθ n e j2π f n t = N ∑ n=1 a n cos(2π f n t + θ n ) (1) where n is a positive integer, N ≥ 1 is the total number of sub-carriers, and a n , θ n , f n are the amplitude in Volts, phase in radians, and frequency in Hz, of the nth sub-carrier respectively. We organize the a n , θ n , f n values in vectors of the following form: According to Parseval's theorem, if S( f ) is the Fourier transform of s(t), then the average power of s(t) in Watts can be calculated as: The definition of the waveform's central frequency ( f c ), bandwidth (B), and pulse period (T o ) determines the total number of sub-carriers and the associated f n frequencies of the signal. The minimum and maximum frequency tones in the spectrum are limited by f min = f c − B 2 and f max = f c + B 2 , accordingly. The resulted frequencies f n of the sub-carriers correspond to the multiples of 1/T o in the interval [ f min , f max ].
To compute the received signal at the input of the rectenna, the frequency response function of the propagation across the medium is required. For a selected frequency f n , this can be expressed by: where b n is the amplitude and ψ n is the phase of the propagation channel's frequency response, respectively. These values can be also organized in vectors b and ψ accordingly. If we presume that the propagation channel is linear time-invariant, the incident wave at the input of the rectenna r(t) is given by: a n b n e jθ n e jψ n e j2π f n t = N ∑ n=1 a n b n cos(2π f n t + θ n + ψ n )

Propagation Channels
A realistic testbed of the propagation channel was fabricated in a laboratory. Since the use of CSI is more practical in a multipath propagation environment, we developed complex channels by placing a rectangular box of highly reflective surfaces and various metal obstacles between the transmitter and the receiver. A set of VNA measurements were conducted to characterize four different propagation channels, based on the four different sites of the receiver antenna, respectively ( Figure 2). The obtained data by the VNA measurements were utilized to compute b and ψ depending on the waveform characteristics. Top view (schematic representation) of the fabricated realistic testbed for conducting VNA measurements in a complex multipath propagation channel. Details of the schematic view include the three-dimensional box interposed between the transmitter and the receiver and the highly reflective metal objects placed randomly inside the box (tinfoil has adhered to its inner walls). The sites of the receiver antenna that correspond to the propagation channels are also marked (in red color) and numbered.

Rectenna Model
The utilized rectenna model [9,10] is depicted in Figure 3. It comprises an antenna, a single diode rectifier, and a load. The receiver antenna model also includes a voltage source v s (t) and a series impedance (R s = 50 Ω). The rectifier's input impedance R in , as well as its input voltage v in (t) is also depicted in Figure 3. In order to simplify the rectenna model, a lossless antenna and a perfect matching between the antenna and the rectifier (R s = R in = 50 Ω) is considered. Moreover, the antenna noise and the diode's series resistance are omitted. The value of the load is set to R L = 1.6 kΩ and the value of the capacitance C is equal to C = 50 T o R L , resulting in a mitigation of the output fluctuation v out .  [9,10]. The same circuit depicting the equivalent impedance R in due to the rectifier and the load is also shown (right).

Objective Function
The optimization problem in hand is to design a waveform s(t) (find the corresponding a and θ) that maximizes the DC power delivered to the load R L , for a specific constraint in the average power of the designed waveform. For this problem, we apply the objective function that is utilized in [10], which is based on Kirchhoff's circuit laws of the equivalent circuit depicted in Figure 3. The assumption of a lossless antenna, perfect matching between the antenna and the rectifier, ideal diode and negligible fluctuation of the output voltage, led to a tractable nonlinear rectenna model and made the derivation of the objective function possible [9,10]. According to [10], the formulation of the objective function can be expressed as: which is subject to: a n ≥ 0 and Based on [10], the maximization of (6) with respect to (7) maximizes the DC power delivered to the load. In (6), we set η = 1.05 as the ideality factor of the diode and V t = 25.86 mV its thermal voltage. Also, we denote as P t in (7), the power constraint imposed on the designed waveform in Watts.
For the numerical computation of the integral in (6), we apply the 2-point Newton-Cotes formula [10]. To this end, we set the sampling frequency at f s = 20 f c . Considering that the sampling period is given by ∆t = 1/ f s , the number of sub-intervals by Q = T o /∆ t , and Q is rounded to the nearest integer that is ≤ Q, (6) derives to: an bn cos(2π fn q∆t+θn +ψn ) ηV t The optimal selection of the phases θ n in the transmitted signal are defined as θ n = −ψ n [9]. Consequently, if we choose a specific propagation channel and a set of f c , B, and T o values, the optimization problem confines to obtain the optimal value of the vector of amplitudes a that maximizes (8) with respect to the constraint in (7).

Evolutionary Algorithms
In a typical EA like DE, an initial random population of candidate solutions to the underlying optimization problem is evolved by applying various mathematical operators to the population's individuals. Then, the generated solutions that demonstrate better fitness to their predecessors substitute them. The same procedure of generating new solutions and performing selection is repeated until a stopping criterion is fulfilled. The fittest individual that has survived through this process is the proposed solution to the optimization problem.

Differential Evolution
The iterations in DE algorithm are called generations (g) and the termination criterion can be set by selecting a maximum number of generations (g max ). Alternatively, the termination criterion can be a maximum number of objective function evaluations (FE max ). The parameter vectors are called population vectors or individuals and are denoted by x i,g , where i = 1, 2, . . . , NP and g = 1, 2, . . . , g max . NP is called population size and it is one of the three control parameters of DE. Each vector x i,g consists of D elements that correspond to the decision variables of the optimization problem and are denoted by x j,i,g , where j = 1, 2, . . . , D. Consequently, the dimensionality of the optimization problem is D. The population vectors in generation g are usually organized in a population matrix denoted by P x,g , which is given in (9).
The first step in DE algorithm is the initialization of NP individuals. The population vectors of the first generation are initialized randomly inside the feasible search space. Consequently, when the decision variables are subject to some constraints, these constraints should be taken into account. After the initialization phase, the algorithm enters the main loop, which consists of three different operations called mutation, crossover, and selection.
Mutation is the operation that generates mutant vectors. A mutant vector is created by perturbing one of the population vectors with the weighted difference of two others. In this way, DE exploits information from different individuals to direct the search. During mutation, each x i,g serves once as a target vector and a corresponding mutant vector v i,g is generated by: where y i,g is called the base vector. In the classic DE version it is y i,g = x r1,g . Indices r1, r2 and r3 are randomly chosen every time a mutant vector is generated, under the condition that r1, r2, r3 ∈ {1, 2, . . . , NP} and r1 = r2 = r3 = i. The scale factor (SF), denoted by F, is one of the control parameters of DE. It is always > 0 and while there is no upper limit, it is rarely > 1 [23]. We should note that (10) describes the mutation strategy in the classic version of DE, but many other different mutation strategies have been introduced in the literature. Crossover is the operation responsible for the diversity enhancement of the population. During crossover, a new trial vector u i,g is generated for each target vector by combining some parameters of v i,g with some parameters of x i,g . This operation is controlled by a control parameter called Crossover Rate (CR), which may take values in the range [0, 1]. For the generation of a trial vector, a random number in the range (0, 1) is generated for each j. If this number is smaller than or equal to CR, the trial vector inherits its jth parameter from the corresponding mutant vector, otherwise, it retains the jth parameter of the target vector. Moreover, a random index j rand ∈ {1, 2, . . . , D} is generated for each i and the new trial vector inherits its jth = j rand parameter from the corresponding mutant vector. This is to ensure that u i,g takes at least one element from v i,g . The crossover operation is described by: At this point, we should note that various other crossover strategies have been also introduced in the literature.
Selection is the operation of choosing between a trial vector and the corresponding target vector. The chosen vector "survives" and continues in the next generation. The selection process is very simple in DE as the vector with the lowest fitness value will survive in a minimization problem or the vector with the largest fitness value will survive in a maximization problem. For a maximization problem, the selection operation is described by:

CoDE, L-SHADE, Jaya Algorithms
In this subsection, we are going to highlight the most important characteristics of CoDE, L-SHADE, and Jaya. We are not going to give a full description of the algorithms since their basic structure is similar to DE.
CoDE is an algorithm that is based on DE and utilizes the distinct attributes of different trial vector generation strategies and control parameter values. The term trial vector generation strategy refers to the selected mutation and crossover strategies. The authors in [20] utilized other researchers' experiences of applying the DE algorithm in various problems to choose three trial vector generation strategies and three pairs of control parameter values (F, CR) that demonstrate distinct advantages. Each of them is suitable for different kinds of problems and different stages of the optimization process. The novelty in CoDE is that each one of the trial vector generation strategies is combined randomly with one of the (F, CR) pairs to produce three different trial vectors for each target vector. The three trial vectors compete with each other before competing with the target vector. This gives CoDE its adaptive capabilities and enables its effectiveness for a variety of optimization problems. One may say that CoDE uses adaptive control both for the control parameters and the trial vector generation strategies. The only control parameter value that needs to be adjusted is the population size.
L-SHADE is also based on the DE algorithm and it automatically adapts its control parameter values during the search process. These values, which differ between individuals, are indicated as F i , CR i , and are adapted according to the control parameter values of the trial vectors that have successfully survived selection in past generations. L-SHADE features a modified mutation scheme that directs the population towards the top p × NP individuals, where p ∈ (0, 1). Also, an optional external archive contains former target vectors that have failed to survive selection and they are used in a way that enhances the diversity of the current population. Finally, L-SHADE incorporates a feature called Linear Population Size Reduction (LPSR), meaning that the population size NP is continuously reduced from generation to generation according to a linear function. Hence, L-SHADE automatically adjusts not only the F i , CR i control parameters, but the population size as well. The control parameters that need to be provided in L-SHADE are the mutation scheme parameter (p), the size of the external archive (A), the size of the historical memory (H), the initial population size (NP init ), and the final population size (NP min ).
Jaya is not classified to a DE variant; nonetheless it is also a population based heuristic algorithm. It has only one mechanism for perturbing the population vectors and generating new ones. The basic concept is that the generated solutions move towards the individual with the best fitness and away from the one with the worst. Jaya is a very simple algorithm and its only control parameter is the population size.

Modifications for Application in Waveform Design
Here we describe the necessary modifications made in the original versions of the algorithms to apply them in the waveform design for the WPT problem. These include the initialization of the population and the enforcement of boundary constraints.
In the optimization problem of waveform design for WPT, the population vectors of any EA coincide with the a vectors, while the dimensionality of the problem is equal to the waveform's total number of sub-carriers (D = N). The population vectors of the first generation are initialized randomly inside the feasible solution space. According to (7), the minimum value of a n is 0. On the other hand, the maximum value is √ 2P t , considering a scenario where all the available power would be allocated in only one of the sub-carriers. Hence, each element of the population vectors initially takes a random value between 0 and √ 2P t . Then, the vectors that violate the transmit power constraint are normalized so that the average transmitted power for each one of them is P t before entering the main loop of the algorithm. Considering that each x i,g represents a vector of amplitudes a, the normalization is accomplished by multiplying the invalid vectors with a term Z i that is calculated individually for each one of them: The invalid vectors generated during mutation and crossover operations in DE are also substituted by vectors that lie inside the feasible solution space. This process is applied in two steps right after the crossover operation in a generation basis. In the first step, the negative valued elements of each trial vector are substituted by the respective elements of the corresponding base vector multiplied by a uniformly distributed random number in the range (0, 1): The random number rand j (0, 1) is generated individually for each element substitution. In the second step, the trial vectors that violate the transmit power constraint are normalized in the same way as described earlier: Algorithm 1 summarizes the pseudocode of DE in the WPT optimization problem. The rest of the algorithms are modified and applied in the same manner.

Application of Algorithms
In this section, several algorithms (DE, CoDE, L-SHADE, Jaya, and SCP-QCLP) are applied in the waveform design for the WPT problem. The objective is to design multi-tone signals with f c = 910 MHz and B = 100 MHz. We set T o = 20, 40, 80, 160, 320 ns, which lead to waveforms with N = 2, 4, 8, 16, 32 sub-carriers respectively. The propagation channels examined are the four channels described in Section 2.2. Consequently, 20 different cases are assessed. The power constraint in the designed waveform is set to P t = −30 dBm. This choice leads to a received signal r(t) with average power that varies from −14.4 dBm to −20.6 dBm, depending on the channel and the waveform. We consider this as a reasonable choice taking into account that in [9,10] the examined transmitted signals delivered with a received signal power of −20 dBm and ≤ −11.67 dBm respectively. Moreover, according to [24], a received signal of 10 µW to 100 µW is enough to power a modern wireless low power device. We note that the power of the finally transmitted signal is higher than −30 dBm since the utilized propagation channels include an amplifier and antenna before its transmission. Each algorithm is applied 50 times per case. For each obtained solution a, we calculate the rectified DC power using the appropriate expressions derived by the rectenna model of Figure 3 and a bisection method as described in [10]. For this calculation we set the sampling frequency at f s = 100 f c in order to increase the accuracy in the numerical computation of the integral. From the obtained data, we compute the average rectified DC power and the corresponding standard deviation over 50 runs for each algorithm.
Before the execution of each of the applied EA algorithms, the control parameter settings, as well as the termination criteria used in each case, are determined. Table 1 lists the termination criteria and the optimal control parameter values (F, CR) for DE. These values are obtained by conducting a parametric study for the 5 cases of the second propagation channel, whereas the population size is set to NP = 10D. For the CoDE algorithm, the population size was set equal to D, except for the cases of D < 6, where it is set to NP = 6. This limitation derives from the minimum required number of population vectors in CoDE, which is 6. In L-SHADE, we set the mutation scheme parameter to p = 0.11, the size of the external archive to A = 1.4NP, the size of the historical memory to H = 5, the initial population size to NP init = 18D, and the final population size to NP min = 4. Finally, in Jaya, we set the population size to NP = 10D.

Performance Comparison
The average rectified DC power for the waveforms designed by the algorithms under test is listed in Table 2. We can easily derive that all algorithms yield very similar results. Considering the performance of SCP-QCLP as the reference, we will compare the other algorithms with it. In DE, we observe a maximum percentage decrease of 2.96% in channel 1 for T o = 320 ns. This is a slight performance deterioration that could be explained by the fact that the tuning of F and CR parameters was based only on the characteristics of the second propagation channel. However, taking into consideration the overall outcome, the selected values of F and CR give satisfactory results for all the channels under test. Jaya's performance is just slightly deteriorated for N = 32 if we take into account the maximum percentage decrease of 0.66% in channel 4. CoDE and L-SHADE generally score values of a very slight improvement in the rectified DC power compared to SCP-QCLP, with a maximum of around 0.06% in channel 4 for T o = 40 ns. However, this improvement is trivial and it can be attributed to the value of the termination criterion applied in SCP-QCLP mode (it was set to ∈= 10 −3 as in [10]).   Table 3 lists the standard deviation values of the rectified DC power over 50 runs for each algorithm as they presented in Table 2. We can easily conclude that these values are quite low compared to the corresponding average ones. Note that SCP-QCLP should normally demonstrate zero standard deviation in any of the cases since it is not a stochastic algorithm and it is developed to converge always to the same solution. However, the observed values for SCP-QCLP are not exactly zero, due to the roundoff error of floating-point arithmetic. Tables 4 and 5 present the best and worst values of rectified DC power detected by each algorithm. We can derive that the best values are quite similar to the worst, as expected from the corresponding low standard deviation. Also, we should remark that EAs did not detect any solution in any of the cases to demonstrate more than 0.06% increase in the rectified DC power compared to SCP-QCLP.
We should point out that the SCP methods used in [9,10] make use of some techniques to transform the original non-convex objective function into a convex one that approximates the problem and can be solved iteratively until convergence. Also, in these methods, the search starts from only one point in the solution space, which could easily lead to obtaining a locally optimal solution. On the other hand, in this paper, we employ evolutionary algorithms using the original non-convex objective function. EAs initially sample randomly the solution space and utilize a population of solutions to guide the search. Furthermore, they are equipped with mechanisms, like the crossover, to enhance the population's diversity and prevent convergence to a local optimum. However, convergence to a globally optimal solution is still not guaranteed. We should also note that the four EAs tested in this paper have different characteristics and are suitable for different types of problems. Considering all the above, the results indicate that the obtained solutions may be the globally optimal ones. Table 3. Summarized results (standard deviation of the presented average rectified DC power in Watts of Table 2) for DE, CoDE, L-SHADE, Jaya, and SCP-QCLP WPT modes. The rest of the parameter values (runs, propagation channels, multi-tone waveforms) are equivalent to the values listed in Table 2.

Convergence Plots
The assessment of the derived results so far has exhibited that among the tested EAs, no one stands out for yielding superior solutions. The next step for the evaluation of their performance is to check which one of them converges within a lower number of function evaluations. This can be achieved by examining the convergence plots of the algorithms in each one of the cases. These plots are generated by processing the already obtained data and depict the mean fitness value over 50 runs as a function of objective function evaluations (Figures 6-9).
By observing the convergence plots for all channels and waveforms, we could say that CoDE algorithm converges slightly faster or at least equally fast in most of the cases compared to the other competitors. It seems that the adaptive nature of CODE makes it suitable for different problems. An exception is detected in Channel 1 for B = 100 MHz and T o = 40 ns. The superiority of CoDE's convergence speed is more apparent for a larger number of sub-carriers (N = 16, N = 32). On the other hand, Jaya performs satisfactorily for N = 4 and N = 8, yet poorly for N = 16 and N = 32, probably due to its greedy vector generation strategy. Finally, DE and L-SHADE demonstrate quite similar behavior, as they require a larger number of function evaluations for convergence than CODE and Jaya for N = 4 and N = 8, while they converge faster than Jaya and slower than CoDE for N = 16 or N = 32.

Conclusions
In this paper, we evaluated the application of EAs in the problem of optimal waveform design for WPT systems. In detail, we applied various algorithms (DE, CoDE, L-SHADE, and Jaya) for the optimization of the transmitted multi-tone waveforms. A SISO WPT system was utilized and a multipath propagation channel was implemented. From the derived results, EAs seem to successfully obtain the optimal solutions to this non-convex problem, since all four tested algorithms converged to very similar results. However, the differences between the obtained optimal solutions from EAs and the SCP-QCLP algorithm are indiscernible. As a result, we can conclude that the objective function's landscape for the given propagation channels is not very complicated. Consequentially, the complexity of EAs was not beneficial compared to the SCP-QCLP algorithm for the given optimization problem. Also, the results increase the confidence that the SCP approach can provide optimal waveforms despite its simplicity. Finally, CODE exhibited a slightly better convergence speed compared to the other EA algorithms.
The advantage of using EAs is that they are very flexible and can be used in conjunction with any objective function. To this end, a more realistic model of the rectifier, designed in some simulation software, could be integrated into the objective function. In that case, a charge pump circuit could be used instead of the single diode rectifier. For the given optimization problem, the disadvantage of EAs compared to SCP-QCLP is the longer execution time due to their stochastic nature.
It would be interesting to examine the application of EAs in waveform design for more complex WPT systems (more challenging objective functions). One scenario would be to examine waveform design in systems that comprise multiple transmitter antennas and/or multiple users like in [9,17]. However, in that case, it is expected that the computational time for the calculation of the objective function, and consequently the total execution time, would increase. Another scenario could be to apply EAs in systems that require the joint optimization of WPT and some other characteristic, like SNR [15] or information transfer rate [16]. That would require the modification of EAs by integrating some mechanisms for the enforcement of the additional constraints.
The optimized waveforms demonstrate high PAPR that is beneficial for the receiver rectenna's RF to DC conversion efficiency, but not for the amplifier used in the transmission. This is because the behavior of common amplifiers is nonlinear for these types of signals. Consequently, in a real-world application, an additional constraint on PAPR should be imposed. In that case, since there is a strong limitation on the possible s n values, it would be interesting to examine if optimizing both s n and φ n using EAs could provide some advantages.