Multi-Objective Optimization of Hybrid Renewable Energy System Using an Enhanced Multi-Objective Evolutionary Algorithm

Abstract: Due to the scarcity of conventional energy resources and the greenhouse effect, renewable energies have gained more attention. This paper proposes methods for multi-objective optimal design of hybrid renewable energy system (HRES) in both isolated-island and grid-connected modes. In each mode, the optimal design aims to find suitable configurations of photovoltaic (PV) panels, wind turbines, batteries and diesel generators in HRES such that the system cost and the fuel emission are minimized, and the system reliability/renewable ability (corresponding to different modes) is maximized. To effectively solve this multi-objective problem (MOP), the multi-objective evolutionary algorithm based on decomposition (MOEA/D) using localized penalty-based boundary intersection (LPBI) method is proposed. The algorithm denoted as MOEA/D-LPBI is demonstrated to outperform its competitors on the HRES model as well as a set of benchmarks. Moreover, it effectively obtains a good approximation of Pareto optimal HRES configurations. By further considering a decision maker’s preference, the most satisfied configuration of the HRES can be identified.


Introduction
The increasing consumption of fossil fuels coupled with environmental degradation has resulted in the growth of eco-friendly renewable energy resources [1].The hybrid renewable energy system (HRES) that integrates various renewable energy resources together has become increasingly popular [2].The optimization of HRES is a multi-objective problem (MOP) in nature.To be specific, multiple objectives have to be considered simultaneously [3] such as minimizing the lifetime system costs, the carbon emissions, maximizing the system reliability, and the utilization of renewable energies.
Given multifarious variables and non-linear models, the single-objective algorithms are hard to find the global optimum.In recent decades, multi-objective optimization methods, due to their good performance on obtaining the trade-offs for MOPs, have been extensively applied to the optimal design of HRESs [4].The most popular approach proposed to solve these kinds of problems is a multi-objective evolutionary algorithm (MOEA), which has been used to design or schedule many hybrid power systems [1,[5][6][7][8][9].For example, Dufo-López and Bernal-Agustín [7] applied the strength Pareto evolutionary algorithm (SPEA) to optimize the cost, pollution and unmet load of the HRES.In [10], the non-dominated sorting genetic algorithm (NSGA) [11] was introduced to the optimization of HRES with two objectives, i.e., minimizing the total cost and pollutant emissions.Following those Pareto-based algorithms, the MOEA based on decomposition (MOEA/D) [12] was also used in this research area for computational efficiency and scalability.It has been applied in [13] to optimize an HRES in which the cost, CO 2 and SO 2 emissions are minimized and the output power is maximized.
However, most of the existing studies tend to merely consider renewable resources as an isolated-island mode (the HRES responds to the load demand on itself).In fact, the HRES can also connect to the power grid (that is, the HRES exchanges power with the public grid [14,15]), which can further enhance the reliability and the economic benefit of HRES.In this study, the multi-objective design of HRES that operates in both isolated-island and grid-connected modes is modelled and simulated.The contributions are that (i) a comprehensive set of components are considered in the design of HRES.Taking the batteries as an example, the shortened lifespan of batteries resulted from the frequent charge-discharge cycles is calculated; (ii) this paper describes the HRES that can operate in both isolated-island and grid-connected modes.In addition, the optimal operational strategy for each mode is provided; and (iii) a simple yet effective algorithm, namely, MOEA/D with a localized penalty-based boundary intersection (LPBI) method (denoted as MOEA/D-LPBI) is proposed for the multi-objective design of HRES in both modes.The MOEA/D-LPBI can eliminate the need of penalty value setting in the original MOEA/D-PBI, and significantly outperforms MOEA/D-PBI on both benchmarks and the multi-objective design of HRES.
The remainder of this paper is presented as follows.The optimization model and the control strategy of HRES are presented in Section 2. Section 3 introduces the LPBI method and the derived algorithm, MOEA/D-LPBI.Section 4 examines the performance of MOEA/D-LPBI on benchmarks and demonstrates its performance on the design of HRES.Section 5 concludes the paper and identifies future studies.

Mathematical Model of the HRES
In general, an HRES consists of photovoltaic (PV) panels, wind turbines, the battery storage, diesel generators and some accessory devices.The considered HRES in which the power grid can be also included is illustrated in Figure 1.The details about the optimization model and the operational mechanism are elaborated in the following subsections.

Optimization Model
As a piece of important generating equipment, the PV panel has been widely used because of its safety and no pollution.Assuming that the PV array combines N s of PV panels in series and N p strings in parallel, the maximum output power can be calculated as Equation (1) [16]: where the factor F loss takes into account the losses of power due to shadows, dirt, etc.The total number of PV units (N pv ) is available by N s multiplying N p .The calculation method to obtain the electric current I pv and the voltage U pv is employed from the previous work [17].
The wind turbine consists of the wind tower and the generator, whose output power can be calculated as follows [18]: with: where C p is a coefficient obtained by the maximum wind power dividing the actual power output, ρ is the air density, A wg is the cross section of the rotor, and P wgr denotes the rated power of the wind turbine.V c , V r and V f are the cut-in, the rated and the cut-off wind velocity, respectively.In addition, the wind velocity v at the (actual) height of H wg is converted by the anemometer data (v r ) at the reference height of H r as Equation (3).Note that the height of wind tower is retained within The batteries are frequently applied in an HRES as the storage device.The quantity of electric charge and discharge is dependent on not only the load demand but also the state of charge (SOC).The SOC of the battery storage at the simulation time t is obtained by Equation (4): where P bat (t) is the input/output power (positive during charging and negative during discharging) of batteries.∆t is each simulation time step that is assumed to be an hour.The round-trip efficiency η bat is defined as 80% for charging and 100% for discharging models.Additionally, N bat , C bat and V bat denote the number of batteries, the nominal capacity and the nominal voltage of each battery, respectively.The battery model retains the SOC between the lower limit (SOC min ) and the upper limit (SOC max ) to ensure safety.Correspondingly, the minimum and the maximum battery power are also defined respectively, as shown in Equations ( 5) and (6): Since the number of the charge-discharge circles and the depth of discharge in each circle both affect the lifespan of the batteries, the "rainflow"-based method described in [7] is adopted to estimate the actual lifetime of the battery storage Li f e bat .When the working hours reach the longevity (Li f e bat ), new batteries will replace those old ones.
The diesel generator usually acts as the back-up power and its fuel consumption F cons is identified [7] as Equation (7): where P dgr and P dg are the rated and actual output power, respectively.Note that the latter is not allowed to be larger than the former.γ 1 and γ 2 are coefficients of fuel consumption.
In addition to the above components, some accessory devices, such as the inverter and the rectifier, are also included.Their efficiencies are 90% and 95%, respectively.
As for the objectives in this HRES, it is always expected that the system cost and the output pollution can be minimized.In addition, the load demand can be met.Due to the fact that the HRES can be in either grid-connected mode or isolated-island mode, the optimal design of HRES could be different.Specifically, for the isolated mode, the cost of system (C s ), fuel emissions (F e ) and the probability of unsatisfying the load demand (P u ) are to be minimized; for the grid-connected mode, the cost of system (C s ), fuel emissions (F e ) and the proportion of non-renewable energy (P nre ) are to be minimized.Afterwards, we present the definitions of these objectives.
First, we assume that the lifespan of the HRES is determined by the lifetime of PV panels.Since the batteries have shorter service life and the frequent charge-discharge pattern further shortens the time, its replacement is also considered in this work.Overall, the total cost is comprised of the following: • Payment of purchasing and installing the PV panels, wind turbines, batteries, diesel generators, inverters and rectifiers, denoted as C initial .Among them, the cost of the wind turbines includes that of the wind generator and the tower, and the expense in tower is positively related to its height.• Cost of repair and maintenance of the devices, denoted as C repair .
• Expenditure on the fuel consumed through the lifespan, denoted as C f uel , which is defined by unit price (C uni f uel ) multiplying the quantity of fuel consumption (F cons ).• Cost of replacing the old batteries with the new ones, denoted as C replace .It is calculated by , where C unitBat and N replace are the unit price of a battery and the number of replacement, respectively.• Cost/profit of exchanging power with the public power grid (positive for buying power and negative for selling power), denoted as C grid .
The costs considered above are indicated as Equation ( 8): where C grid is zero when the system operates in the isolated-island mode.Second, we adopt the total amount of CO 2 through the lifespan (T) to measure the pollutant emissions as CO 2 occupies the largest percentage of all emissions in fuel combustion: where E f is the emission factor depending upon the property of the diesel generator and T(1 ≤ t ≤ T) is the maximum simulation time.The amount of emission accumulates every time the device is powered [19].Third, the probability of failing to satisfy all of the load (P u ) is employed to measure the system reliability in the isolated-island mode, which is calculated by Equation ( 10): with: where the available power supply E avail (t) and the load demand E load (t) are compared for each simulation time, and 1 is returned to Bl (boolean variable) if the former is smaller than the latter and vice versa.Thus, P u refers to the proportion of hours of power shortage during the lifetime, within [0, 1].
In the grid-connected mode, P u is always zero since the deficit power can be bought from the power.Thus, the proportion of using non-renewable energy (P nre ) is considered instead, and is defined in Equation ( 12): Overall, the design of an HRES is an MOP where three objectives are to be minimized as shown in Equation ( 13): min F = min (C s , F e , P u ) , the isolated-island mode, min (C s , F e , P nre ) , the grid-connected mode, ( subject to:

Systematic Planning of Operation Mechanism
The energy flow in the HRES is illustrated in Figure 2 with two modes considered.Mode I is the isolated-island mode, and mode II is the grid-connected mode.• Mode I: in this mode, the system is not connected to the power grid, thus the A module in Figure 2 is not considered.The energy produced by the PV panels and the wind turbines is directly provided for the direct-current (DC) load and flows to the alternating-current (AC) load through the inverter.If the provided energy exceeds the total demand (DC load plus AC load), the surplus energy will be saved in the battery storage after meeting the load.On the contrary, if it cannot satisfy the load demand, the batteries will discharge based on the SOC.If there still exists unmet load, the diesel generators start as the emergency power-supply.Note that the power from the diesel generators is alternating current and the part used to supply DC load is converted by the AC-to-DC rectifier.• Mode II: in this mode, when the produced renewable energy is more than the demanded quantity and the batteries have been charged to the maximum, the surplus energy will be sold to the power grid for the profit.On the contrary, if the produced energy cannot meet the load demand, the method of obtaining the supplement is determined by the electricity price (C ep ).Providing low C ep , all of the deficit power will be bought from the public power grid and supplied to the load.On this condition, the B module does not operate.However, if the C ep is high, first the battery storage and the diesel generators are used as supplement power.Once these devices still cannot cover the gap, the required power will be bought from the power grid.To describe it more clearly, the choosing strategy when the renewable energy is insufficient is shown in Algorithm 1.
Algorithm 1: Decision strategy for the source of supplementing deficit energy in mode II.To simulate the operational process, the flowchart is shown in Figure 3.In this flowchart, P(t) and P L (t) are the total power produced by renewable resources (PVs plus wind turbines) and the load demand at the simulation time t, respectively.At the beginning, the hourly solar irradiation, the hourly wind velocity, the temperature and the load demand during T time are input into this model.

Multi-Objective Optimization Using MOEA/D-LPBI
As the multi-objective design of the HRES features non-linearity, mixed variables, etc., it is necessary to adopt an effective algorithm to solve this model.Amongst the state-of-the-art MOEAs, the MOEA based on decomposition (MOEA/D) [12] has gained great popularity, and thus is adopted.It decomposes an MOP into a set of subproblems by a scalarizing method with evenly distributed direction vectors, and then optimizes these subproblems in a collaborative manner [20,21].

The Localized PBI Method
The penalty-based boundary intersection (PBI) method is a frequently-used scalarizing method [22,23], which is written as follows: where θ(θ ≥ 0) is a user-definable penalty value.In addition, λ and z * are the direction vector and the ideal point, respectively.Different penalty values lead to different performances of the PBI method.Specifically, the PBI with a smaller θ emphasizes more on reducing the d 1 distance while the PBI with a large θ is inclined to curtail the d 2 distance.In other words, the PBI is sensitive to the value setting of θ.In addition, the PBI with a θ value offers a fast convergence yet cannot solve the problems with non-convex Pareto fronts (PFs) [24].
To enhance the convergence speed of the PBI method, while handling non-convex PF shapes in the meantime, a localized PBI (LPBI) is proposed [25].In the LPBI method, the newly generated solution is merely compared with solutions that are inside the same hypercone, as shown in Figure 4, instead of the entire population.Given N direction vectors, the objective space can be divided into N sub-regions (i.e., hypercone) by the LPBI method.In any hypercone, e.g., the i-th hypercone, the center line is assumed to be the associated direction vector λ i .Correspondingly, the apex Φ i (i ∈ {1, 2, . . . ,N}) is obtained by averaging those angles between the λ i and its m closest direction vectors (where m is the number of objectives, see Equation ( 16)): where φ k i is the angle between the k-th closest direction vector and λ i .For two vectors, e.g., λ 1 and λ 2 , the cosine value of the intersection angle φ is available by Equation (17): For each λ, only solutions whose objective vectors are inside the associated hypercone are chosen.Taking Figure 4 as an example, the solution for λ 1 is merely selected from x 1 and x 2 .Even though the θ value is small, e.g., θ = 0, the objective vectors of the obtained solutions are still not quite far away from their reference direction lines (λ i ) (see the solid dots (•)).To be exact, the angles formed by those objective vectors and their corresponding direction lines are no more than Φ 2 even in the worst case.
Overall, the LPBI method is responsible for obtaining an optimal solution within each hypercone.Ideally, the LPBI will find different optimal solutions along different direction vectors, providing a set of diversified solutions regardless of the PF shapes.Specifically, the LPBI is implemented as follows: the scalar values of solutions inside the hypercone are calculated by Equation ( 15), while those outside the hypercone are set to g PBI = ∞.

MOEA/D-LPBI
In this subsection, we incorporate the LPBI method into the MOEA/D framework.The pseudo code of the derived algorithm, denoted as MOEA/D-LPBI, is shown in Algorithm 2. Initially, a population of N individuals (solutions) is randomly generated, S = {x 1 , x 2 , . . . ,x N } where x i is the i-th individual.Also a set of evenly distributed direction vectors is generated, denoted as {λ 1 , λ 2 , . . . ,λ N }.Each individual, e.g., x i , is randomly paired-up with a direction vector, e.g., λ i .For each direction vector, its T b neighbours (measured by the Euclidean distance) are found denoted as B(λ i ), where T b is a user-defined positive integer.Their associated neighbouring solutions are denoted as B(x i ).Then, the apex of the hypercone for each direction vector is computed by Equation ( 16).

end
Following that, the population evolves for maxGen generations (lines 5 to 22).In lines 7-11, the parent set Q of x i is probabilistically selected from the neighbours B(x i ) or the whole population S, and is used to generate an offspring x new i by the differential evolution (DE) and polynomial mutation (PM) operators [26].Lines 12-14 merge the current population and the newly generated individuals, denoted as JointS, and obtain the angle matrix φ sw .Note that the proposed algorithm is implemented within a (µ + µ) elitism framework [27].For each vector, lines 15-19 select the solution whose PBI scalar values is the smallest amongst the JointS according to the framework.At the end of each generation, the offline archive archiveS is updated with the newly obtained solutions based on the Pareto-dominance relation (line 20).
In terms of the time complexity of the derived algorithm (MOEA/D-LPBI), calculation of the scalar values of all individuals on all direction vectors runs at O (N × N) , where N is the population size.Additionally, identifying the minimal PBI values from 2N values runs at O(N).Therefore, the overall time complexity of MOEA/D-LPBI is O N 2 .In comparison, if the conventional method is adopted instead, all of the possible values of the variables need to be tried.Assuming that the precision accuracy is 0.01, the total complexity is O(100 2 × N pv × N wg × N dg × N bat × H high × α), which is far larger than that of MOEA/D-LPBI (α is a variable in the PV model).

Experimental Results
MOEA/D-LPBI is applied to solve the multi-objective HRES model after it has been examined on a set of MOP benchmarks (see Appendix A).
Prior to the optimization, the required parameter settings in the above equations are provided in Tables 1 and 2. Among them, Table 1 shows the initial investment cost and the annual repair expenses of all the devices.In addition, the expenditure in replacing old batteries is similar to the initial cost, that is, C unitBat equals C initial (Battery).All of the other parameters of HRES are summarized in Table 2. Readers can refer to [17] for more details.
In addition, the data of the solar radiation intensity, the wind velocity and the environmental temperature over a year is considered as input data.They are assumed to be constants for every hour time step.Their distributions during a year as well as the annual load are shown in Figure 5.In this experiment, T is the total number of hours in a year (also the maximum simulation time) and is calculated as T = 24 × 365 = 8760.To obtain the best configuration of the HRES in the two modes, MOEA/D-LPBI is applied to handle the model.In addition, not only does the original MOEA/D-PBI serve as a reference, but also two classic algorithms, i.e., preference-inspired coevolutionary algorithms using goal vectors (PICEA-g) [17,28] and MOEA/D with the weighted Tchebycheff (MOEA/D-TE) [24], are compared with the proposed algorithm afterwards.
Related parameter settings of these algorithms are described below: • Algorithm runs: the maximum generation is set to 50.Assuming that each PV string connects four units in series, the number of N pv is a multiple of four.• MOEA/D parameters: the neighbourhood size T b is set to be 10 and the probability of selecting in the neighbourhood is set as δ = 0.8.In the reference algorithm, the replacement size nr is 2. • Penalty values: the considered penalty values are θ = {0, 0.05, 0.1, 0.5, 1, 5, 20, 100}.
Simulations of the HRES in two modes are respectively operated.Note that the electricity price during peak period (8:00 to 22:00) and trough period (22:00 to 8:00 on the next day) in mode II is 0.2$/kWh and 0.1$/kWh, respectively.The MOEA/D-LPBI under various penalty values and its competitors are applied to obtain a set of Pareto optimal HRES configurations.To quantitatively compare these algorithms, the hypervolume (HV) is adopted to evaluate their performance.This indicator is one of the most frequently used performance metrics in literature [29], the larger value indicating the better performance.It is calculated as follows: first, all the objective values of solutions are normalized within [0,1], and then the HV metric measures the region enclosed by the normalized non-dominated solutions and a reference point z re f , which is set to [1.1, 1.1,. . ., 1.1].The mean HV results of the obtained solutions from 31 runs are shown in Table 3.The results from Table 3 reveal that MOEA/D-LPBI obtains more competitive performances than the traditional framework for almost all penalty values.Moreover, MOEA/D-LPBI with θ = 20 performs best in both modes.The corresponding Pareto optimal solutions are plotted in Figure 6.Provided the suitable penalty setting (i.e., θ = 20), the MOEA/D-LPBI is further compared with another two competitors.Similarly, their mean HV values are illustrated in Table 4 with the Wilcoxon-ranksum two-sided comparison [30] procedure at 95% confidence level employed.Through quantitative comparison, the proposed algorithm (with θ = 20) is also shown to be superior than the classic versions PICEA-g and MOEA/D-TE.
Amongst the obtained Pareto optimal solutions by the MOEA/D-LPBI, a decision maker can choose the preferred one.For example, • in mode I, if the decision maker prefers to meet the energy demand 100% and limit the emission within 250 kg, then one can set P u = 0% and F e < 250.Then, the three filtered are shown in Table 5. Amongst these solutions, the one that produces the minimum cost is selected, i.e., the third one (in bold).• in mode II, if the decision maker prefers to minimize the fuel emission and the utility of non-renewable energy, the solution having minimum F e and P nre is selected, i.e., the first one (in bold).Additionally, it can be observed that the diesel generators are not used in this case.According to the selected results, the final optimal configurations of the HRESs in the two modes are confirmed.Their processes of supplying power during a year are simulated as shown in Figures 7 and 8.It can be observed that the configuration is different in the two modes and the power grid takes the place of the diesel generators in mode II.

Conclusions
This paper studies the multi-objective optimal design of hybrid PV-wind-diesel-battery system for the power-supply, considering that three objectives are simultaneously to be minimized in the specialized operational mode.In mode I (isolated), the system relies on the battery storage and the diesel generator to supplement deficit power; in mode II (grid-connected), it prefers buying power from the power grid to using the diesel generators (under the provided parameters).By a MOEA/D using localized PBI, denoted as MOEA/D-LPBI, a set of Pareto optimal configurations of the HRES in each mode is obtained, from which a decision maker can choose the most adequate one.In addition, Figure A1 demonstrates that MOEA/D-LPBI exhibits excellent performance and is robust to the PFs, that is, the solutions nearly cover the whole PF region and are close to the PF regardless of the choice of θ.However, the performance of MOEA/D-PBI varies significantly and is inferior to the improved version for all problems.

Figure 1 .
Figure 1.Illustration of the hybrid renewable energy system (HRES) with power grid.PV: photovoltaic.

Figure 2 .
Figure 2. Illustration of the energy flow in the HRES.

Figure 3 .
Figure 3. Flowchart of system simulation.(a) the isolated-island mode; and (b) the grid-connected mode.SOC: state of charge.

17 x j ← arg min x j ∈JointS C ij ; 18 Replace 20 Update
x i in S with x j in JointS; 19 end archiveS with S based on Pareto-dominance relation; 21 gen ← gen + 1;

•
Individuals: the population size is N = 100.Each candidate solution consists of six variables in the form: [N pv |N wg |N bat |N dg |H wg |α], which, respectively, represent the number of four main components, the height of the wind tower and the inclination angle of the tilted PV panel.

Figure 6 .
Figure 6.Optimal solutions for HRES in different modes.(a) mode I; and (b) mode II.
Buy all the deficit energy from the public power grid 1 if C ep is low then 2

Table 1 .
The initial cost (per unit) and the annual repair expenses of the system components in HRES.

Table 2 .
Parameter specifications for the environment and the system components.

Table 3 .
The mean hypervolume (HV) results of MOEA/D-LPBI and the original version with different θ for HRES.The symbol "+", "−" or "=" means that the counterpart is statistically better than, worse than or comparable to MOEA/D-LPBI.The best HV value for each problem is marked in boldface.

Table 4 .
The mean HV results and the standard deviations (in the parentheses) of MOEA/D-LPBI, PICEA-g and MOEA/D-TE for HRES.Just as above, their comparable results are marked by "+", "−" or "=".

Table 5 .
Selected results using MOEA/D-LPBI with specified reference.