Multi-Objective Optimization for Coordinated Day-Ahead Scheduling Problem of Integrated Electricity-Natural Gas System With Microgrid

This paper presents a multi-objective optimization algorithm for coordinated day-ahead scheduling problem of integrated electricity-natural gas system with microgrid (IENGS-M). Mathematically, the day-ahead scheduling of IENGS-M is formulated as a multi-objective optimization problem considering multitudinous constraints. In order to solve the problem efficiently, we introduce an acceleration of differential evolution, Lévy search strategy and a treatment mechanism to multitudinous and complex constraints into the original Non-dominated Sorting Genetic Algorithm-III (NSGA-III). Furthermore, a decision making method based on a fuzzy function approach is used to determine a final optimal solution from the Pareto-optimal solutions. Simulation studies are carried out on a modified IEEE 39-bus system and 15-node gas system to verify the effectiveness of the modified NSGA-III (MNSGA-III), in comparisons with the NSGA-II and NSGA-III. The simulation results show that the Pareto-optimal solutions obtained by MNSGA-III has better convergence performance and diversity than the NSGA-II and NSGA-III.


I. INTRODUCTION
WITH the growing concerns over the climate change and energy security around the world, more and more attention is being paid to the development and utilization of the integrated energy system (IES). As the increasing penetration of renewable energy sources, the interconnections of distributed renewable energy will bring a huge impact on the operation of electricity network and gas network, and the coordination of IES will face severe challenges [1], [2].
Over the last decades, a number of methods have been developed to solve the optimization problems. In general, these methods can be divide into mathematicsbased optimization methods and heuristic methods [3], [4]. Mathematics-based optimization methods including linear programming [5], non-linear programming [6], [7], mixed The associate editor coordinating the review of this manuscript and approving it for publication was Huai-Zhi Wang . integer linear programming [8], mixed integer non-linear programming [9], [10]. Some of the mathematics-based optimization methods can solve this problem very fast. However, these methods have some disadvantages when dealing with non-linear, non-smooth and non-convex day-ahead scheduling problem. Moreover, these methods are sensitive to the choice of the initial values and may be fail into local optimum due to the improper initial points. Heuristic methods have aroused intense interest due to their versatility, flexibility, and robustness in seeking the global optimum solution, such as GSO [11], PSO [12], DE [13], GSA [14], and etc. Attempts have been made to apply these algorithms, including machine-learning algorithms, to tackle large-scale and economically important optimisation problems, such as optimal power flow and economic dispatch in power systems [15].
In recent years, NSGA-III, which proposed by Deb et al. [16], [17], has been show perfect performance for a number of different types of benchmark functions and multi-objective problems. What's more, Mkaouer et al. proposed for the first time a scalable search-based software engineering approach based on this newly proposed evolutionary optimization method NSGA-III, which addresses software engineering problems when a large number of objectives are to be optimized [18], [19]. The authors of reference [20] presented an extended NSGA-III algorithm to select new generation by introducing the dominance relationship criterion based on constraint violation. The superior quality and good distribution of the Pareto-optimal solutions show the NSGA-III can find an efficient alternative for optimizing multi-objective economic emission hydro-thermal-wind scheduling problem. Reference [21] proposed an orthogonal design-based non-dominated sorting genetic algorithm-III (ONSGA-III) to optimize the lockage co-scheduling of Three Gorges Dam and Gezhou Dam. An improved non-dominated sorting genetic algorithm is proposed by Chen to solve reservoir flood control operation problem [22]. Reference [23] proposed improved NSGA-III based on objective space decomposition for enhance the convergence of NSGA-III and tested it on number of many-objective optimization problems with state-of-the-art algorithms. In 2016, Bhesdadiya et al. employed the NSGA-III algorithm to solve multi-objective combined economic emission dispatch problem [24]. The experimental and practical results all show perfect performance of the NSGA-III. However, there is still an insufficiency in NSGA-III algorithm, which emphasized the diversity but neglected the speed of convergence and the ability of exploitation.
Inspired by the NSGA-III algorithm, this paper proposes a MNSGA-III algorithm by introducing an acceleration of differential evolution, Lévy search strategy to mitigate the imperfections of original NSGA-III algorithm. The acceleration of differential evolution is applied to accelerate the searching process and enhance the diversity at the same time. Compared with random walk, Lévy search strategy is more efficient, due to it can escape away from the local optimum to find well-distributed Pareto-optimal solutions in search process. In the proposed MNSGA-III algorithm, we use the Lévy search strategy to improve the ability of exploitation. Moreover, a treatment mechanism is applied to solve multiobjective optimization problem which contains multitudinous and complex constraints. To verify the efficiency of the proposed algorithm, the MNSGA-III is tested on the day-ahead scheduling problem of IENGS-M system, and the results obtained by the MNSGA-III are compared with those by other reported methods.
The rest of this paper is organized as follows: Section II presents the multi-objective operation optimization for the IENGS-M. Section III introduces the framework of the MNSGA-III in details. A fuzzy function approach is presented in Section IV. Simulation results and comparisons are given in Section V. Finally, the last section draws the conclusion of this paper.

II. PROBLEM FORMULATION A. FRAMEWORK OF IENGS-M
The framework of IENGS-M system considered in this paper is shown in Fig. 1. As shown in the figure, the IENGS-M consists of electricity network, natural gas network, distributed heating and cooling stations (DHCs). Specifically, the electricity network and natural gas network are connected by the gas-fired generators. Furthermore, microgrids are connected to the electricity network to embedded distributed renewable energy resources. Consequently, the aim of the multiobjective operation problem of IENGS-M is to determine an optimal energy scheduling over a day-ahead period of time that minimizes the operation cost and emissions subject to several constraints simultaneously.
The electricity part includes electricity network, electricity load, and the electricity supply system-microgrids. Each microgrid consists of wind turbine (WT) generators, photovoltaic generators (PV), diesel generators(DG), which is shown in Fig. 2. The gas part includes gas supply system, gas network, and gas load.

B. OBJECTIVE FUNCTIONS
In this paper, two objective functions in the multi-objective operation of the IENGS-M are defined and minimized at the same time which can be formulated as follows: (1) VOLUME 8, 2020 J 1 and J 2 represent the total operation cost and the pollution of the IENGS-M system, respectively. h(X , Y ) and g(X , Y ) denote the constraints.

1) TOTAL OPERATION COST
The total operation cost of the IENGS-M includes the operation cost of electricity network and the operation cost of gas network. At the same time, the voltage offset is taken into consideration to make electricity network operate in a safe state. The objective function is shown in (2) which represents the total operation cost of the IENGS-M, while (3), (4) and (5) represent the cost of electricity network, the cost of gas network and the penalty term of voltage set, respectively. The cost of electricity network considers the start-up and showdown cost of the wind/PV/battery/generation units. The ϕ(x) in the penalty term of voltage set is described by (6).
T is the total scheduling period. N 1 is the total number of wind/PV/battery/generation units. N 2 is the number of DG. N 3 is the number of gas sources. u j,t is the status of wind/PV/battery/generation unit j at time t. x en,j is the power output of jth wind/PV/battery/generation units at time t. B j is the operation cost of jth generation units. S j,start and S j,shut are the start-up and shut-down cost for jth generation units. Q j W ,t is the amount of jth gas sources at time t. Cost gas,j,t is the gas purchase cost of supplier i at time t. n is number of electricity network nodes. V i,t is the voltage magnitude of node i at time t. V ref i is the anticipant voltage magnitude of node i which is taken as 1.0 p.u. in this paper, δV i is the permitted maximum voltage deviation of node i which is set 0.05 p.u.

2) POLLUTION
The second objective function is the pollution of the IENGS-M. Four of the most pollutants are considered in this objective function: CO, CO 2 , SO 2 , NO x . The mathematical formulation of the pollutant function is described in (7)- (9).
where P i,t denotes the output power of ith generation units at time t. G is the number of the types of pollutant emissions. ψ g represents the amercement of the gth type of atmospheric pollutant emission. E i,g elec and E i,g gas represent the emission rate of gth type pollutant emitted by i unit in electricity and natural gas network, respectively.

III. MODIFIED NSGA-III ALGORITHM A. NSGA-III ALGORITHM
The NSGA-III, which proposed by Deb et al. [16], [17] utilizes the fast non-dominated sorting concept, elitism, a set of well-distributed reference points to update the obtained solutions of diversity preservation. The detail description of the NSGA-III steps are given in [16] and the main procedure of NSGA-III is briefly described below.

1) POPULATION INITIALIZATION
At the first step, SN numbers of D-dimensional individuals are generated randomly. Every single individual in the initial population of the solutions is also randomly generated using uniform random numbers ranging over the feasible limits of each control variable by using the following expression: where D is number of optimization parameters. x min j and x max j are the lower and upper limits of the jth dimension, respectively. rand(0, 1) is a uniformly distributed random number between 0 and 1.

2) GENERATE REFERENCE POINTS ON A HYPER-PLANE
The predefined set of reference points are used to ensure diversity of the obtained solutions. In this paper, Das and Dennis's systematic approach that produces points on a normalized hyper-plane is used. If p divisions are considered along each objective, the total number of reference points (H ) in an M -objective problem is given by usually, the number of H is directly related to the desired number of trade-off points.

3) TOURNAMENT SELECTION OF PARENT POPULATION
At the ith generation, the parent population is P t and its size is SN. Then the parent population P t is divided into different non-domination levels (F 1 , F 2 ,. . . , F l ) based on nondominated sorting. In addition to front rank, a new parameter called crowding distance is also calculated for each individual which is a measure of how close the individual is to its neighbors in the objective function. The larger value of the crowding distance is, the better diversity of individual is. Then pre-offsprings are selected from the parent population by applying Binary Tournament Selection based on nondominated operator (≺ n ). Non-dominated operator is based on front level and crowding distance as follows: 1) The individual with lower front rank value is greater than the other regardless of crowding distance, and it's selected.
2) The individual is selected with the larger crowding distance when two individuals located in the same front.

4) CROSSOVER AND MUTATION
The offspring population (Q t ) is created after the operation of Simulated Binary Crossover and mutation to the preoffspring. Thereafter the combined populations of parents and offspring (S t = P t Q t ) are sorted again based on non-domination.

5) NORMALIZATION OF OBJECTIVE FUNCTION
Before associating the combined populations with the reference points, the objective functions of the combined populations need to be normalized so they have an identical range. Firstly, calculate the ideal point obj min by finding the minimum corresponding objective value. Secondly, convert to the new objective value by subtracting the objective value obj with the ideal point obj min . Thirdly, calculate the extreme points by identifying the solutions that make the achievement scalarization function minimum. Note that the achievement scalarization function here is formed with the objective value obj and the weight vector close to the corresponding objective value. Then, calculate the intercept points a obj . These points are the interception between the corresponding objective values with its corresponding extreme points obtained from the last step. Lastly, normalize the objective values with the following equation: when a hyper-plane is constructed, it will yield M i=1 norm(obj.i) = 1, where M is the number of objective functions.

6) SELECTION BASED ON REFERENCE POINTS
After normalizing all the objective values in S t , it need to associate them with the reference points. Firstly, the perpendicular distance between individual in S t and each of reference lines, which is constructed by joining the ideal point with the reference point, is calculated. Each individual in S t is then associated with a reference point having the minimum perpendicular distance. Finally, the niche-preservation operation which plays an important role in maintaining the diversity of the solutions is used to select members from F l . For the jth reference point, the niche count ρ j is the number of the individuals form S t /F l that are associated with the jth reference point. The reference point having the minimum niche count is identified and the member from the last front F l which is associated with the minimum niche count is included in the final population. The niche count of the identified reference point is increased by one and the procedure is repeated to fill up population S t+1 .

B. MODIFIED NSGA-III ALGORITHM 1) ACCELERATION OF DIFFERENTIAL EVOLUTION
The search space exploration of NSGA-III is executed by the operation of crossover and mutation. The parent population seeks the optimal solution along their way until the offspring population, which is created by Simulated Binary Crossover, dominated the parent population. But this is not conducive to raise the convergence speed and optimization efficiency to some extent. In this paper, a selection of DE algorithm is applied as the acceleration operation to accelerate the searching process and enhance the diversity of the populations in the approach. In order to enhance the search efficiency, DE/best/1/bin is adopted. The individual which is the best of the three randomly selected different parent individuals is regarded as their donors. Hence, the individuals at the mth iteration in the operation are produced as follows: where F acc is the scaling factor, R acc is the crossover factor, Np is number of the operation that is equal to the number of the parent populations, D is the number of optimization parameters, iteration max represents the maximum iteration, r1, r2 and best identify three different random numbers uniformly distributed in [1,Np], the best is selected from the three individuals based on non-dominated operator. The offspring will replace its parents as the rank of the offspring is better than the parents, or based on the reference points the offspring's diversity is better than the parents. To our knowledge, the operation of DE acceleration has not been employed in the way described in our approach for the Multiobjective problem.

2) LÉVY FLIGHT STRATEGY
Lévy flight is random walk that is named after Paul Lévy, a French mathematician. Various studies have shown that VOLUME 8, 2020 flight behavior is more efficient than the random walk, because Lévy flight owns an inverse square power-law distribution of flight lengths. The symmetrical Lévy stable distribution can be defined as: (14) where y ∈ R, 0 ≤ λ ≤ 2 defines the index and determines the shape of the distribution. γ > 0 is a scaling factor that selects the scale unit of the distribution. By Lévy flight strategy, offspring of the population is calculated as: where α is the step size which is related to the scales of the problem of interest. In the proposed approach, α is random number for all dimensions of the individuals. Combining the Lévy flight strategy with the acceleration of DE, equation (13) will be transformed into: when the current number of iteration is considered in the stepsize, the equation which is adjusted along with the iteration will be converted to as follows: where u and v are normal distributed. That is, where the is the standard Gamma function, stepsize * (x m r1,j −x m r2,j ) is added to update the selected best individual to find a new individual. And the simulation studies in the next section will verify its efficiency in searching for a superior solution.

3) TREATMENT MECHANISM OF MULTITUDINOUS AND COMPLEX CONSTRAINS
The equality constraints and the inequality constrains of dependent variables in IENGS-M system are all non-linear functions of the control variables. Therefore, the multiobjective problem is a multitudinous and complex constraints optimization problem. How to deal with the equality and inequality constraints is an important difficulty to be solved. In this paper, a mixed constraints handling mechanism is proposed. The handling mechanism adopts the repair strategy and penalty terms to tackle the constraints. The repair strategy deals with the upper and lower bounder constraints of control variables. Within the repair strategy, the value of control variables can be satisfied through the resetting operation shown as formula (19) during the solution procedures. The repair strategy to the equality of electricity power flow at each bus is carried out through load flow calculation. The rest constraints of inequality contain transmission lines loading, gas well capacities, compression ratio. These constraints are added to the objective function in form of the quadratic penalty terms. Therefore, the objective function can be augmented as follows (20): where η, ψ and χ are the defined as penalty factors. fun vio (x) represents the violation function of the constraints x. In this paper, η, ψ and χ are set 50000.

4) DETAIL STEPS OF MNSGA-III FOR THE DAY-AHEAD SCHEDULING PROBLEM OF IENGS-M
According to the description in the previous section, the flowchart of the proposed method to solve the day-ahead scheduling problem of IENGS-M is expressed in Fig. 3.

IV. DECISION MAKING FOR A FINAL OPTIMAL SOLUTION
In general, the objectives in multi-objective problems usually can't be minimized simultaneously for the reason they conflict with each other. Therefore, one solution can't be said to be better than the other according to their fitness. When the final non-dominated solutions are obtained, it's vital to choose the best compromise solution in the decision making process. In this paper, the fuzzy membership approach is applied to find out the best compromise solution from the non-dominated solutions of the Pareto front. Due to the imprecise nature of the decision make's judgment, the ith objective function f i of the individual k is defined as follows: where f max  (22) where M is the number of objectives considered, N is the number of non-dominated solutions, φ i is the weight factor for the ith objective function. In this paper, we consider φ i is 1. The solution with the maximum membership function (µ k ) is selected as the best compromise solution.

A. SYSTEM DESCRIPTION
In this section, simulation studies are conducted to verify the performance of the proposed multi-objective optimization algorithm. Simulation studies are carried out on the modified IEEE 39-bus system and 15-node gas system, which is depicted in Fig. 4.
The two network are connected via three DHCs and two gas-fired units. DHCs are located at node 13, 14, 15 of gas network and bus 15, 16, 24 of electricity network. Gas-fired units are located at node 8, 9 of gas network and bus 8, 4 of electricity network. The detail parameter data can be found in [25]- [28]. The scheduling horizon is 24h. The day curves of electricity demand, gas demand, heating demand, cooling demand as well as wind speed and illumination intensity are summarized in Fig. 5 and Fig. 6.
In the electricity network, the power configurations of each microgrid are shown in Table 1. The predicted load demand of each time period is equally distributed to all PQ buses of the system. WT, PV, Battery don't produce atmospheric pollutants during power generation. Therefore, we just need to take the pollution emission of DG and natural gas consumed devices (gas-fired generators and DHCs) into consideration. The pollutant emission factors and the environmental values of different kinds of pollutants are summarized in Table 2.

B. ALGORITHMS COMPARISONS
In terms of algorithms, in order to demonstrate the validity and effectiveness of the MNSGA-III, 51 independent runs are carried out with comparison to the other two algorithms (NSGA-II and NSGA-III) for the coordinated operation of IENGS-M system, and their parameters are chosen as the same with each other.
Applying the proposed MNSGA-III, the multi-objective optimization problem is optimized with comparisons with the NSGA-II and the NSGA-III. Fig. 7 shows the best Pareto fronts obtained by MNSGA-III and the other two algorithms. As shown in this figure, it is obvious that MNSGA-III outperforms the NSGA-II and the NSGA-III in terms of searching for better converged and more evenly distributed non-dominated solutions in case of the coordinated operation of IENGS-M considering total operation cost and the pollution.
For the sake of brevity, the fuzzy function is applied to select the best compromise solution in these alternatives. The objective values associated with the best compromise solution acquired by different algorithms are present in Table 3. Judging from Table 3, it can be seen that best value of total operation cost by using MNSGA-III is 1.2844×10 4 $, with an average cost of 1.3527×10 4 $ and a worst value of 1.4271×10 4 $ which are less in comparison to the corresponding results obtained by the other two algorithms. The best value of pollution by using MNSGA-III is 131.918$, with an average cost of 135.0643$ and a worst value of 139.2679$ which are also less in comparison to the corresponding results obtained by the other two algorithms.
In order to more intuitively reflect the improved effect of the Pareto frontier obtained by proposed MNSGA-III, we compared the results of different algorithms with several important indicators which include SP-Metric, HV-Metric, I D -Metric, and H D -Metric. These indicators mainly consider the performance of multi-objective optimization algorithms from two aspects: convergence and diversity. Not only that, when we are analyzing these indicators, Wilcoxon signed rank test is executed at a significance level of 0.05 between the proposed MNSGA-III and those by other reported methods. In addition, we use some notations to clarify the comparison results. (1) The p-value of the Wilcoxon signed rank is the probability of a hypothesis of equal median of two-sided test. (2) The h-value is the result of the Wilcoxon signed rank. If h = 0, it indicates that the median of the difference between the proposed MNSGA-III and the compared algorithms is zero. Otherwise, there is a significant difference, if h = 1.
The results obtained by these metrics comparisons are presented in Table 4. As shown in the table, the proposed MNSGA-III can obtain more Pareto-optimal solutions than the NSGA-II and NSGA-III under the same number of iterations. Moreover, the SP-Metric and HV-Metric show that the solutions obtained by MNSGA-III have better quality of  convergence and diversity than those obtained by the two algorithms. Meanwhile, the smaller the I D -Metric and H D -Metric of MNSGA-III is, the better the distribution of the solutions acquired by MNSGA-III is. Not only that, the median differences between MNSGA-III and the other two compared algorithms are significant in the test case shown in the Table 4. As a consequence, it can be concluded that the improvement of the MNSGA-III can guarantee that it can find a better set of Pareto-optimal solutions compared to the NSGA-II and NSGA-III.

VI. CONCLUSION
This paper has proposed an effective MNSGA-III multiobjective algorithm to solve the complicated constrained optimization problem. In order to validate the applicability and performance of the MNSGA-III algorithm in real-world problems, the MNSGA-III is tested on the day-ahead scheduling problem of IENGS-M system.
Compared with the NSGA-II and NSGA-III, the performance of MNSGA-III on the indicator of SP-Metric, HV-Metric, I D -Metric, H D -Metric and the number of Paretooptimal solutions is better, which means that the MNSGA-III is more efficient in searching for the superior Pareto-optimal solutions in terms of convergence and diversity.
Furthermore, a decision making method based on fuzzy function approach has been applied to determine the best compromise solution for the IENGS-M. From the best compromise solutions of MNSGA-III and the other reported algorithms in this paper, it also can be concluded that the proposed algorithm is more efficient in solving the optimal operation problem of the IENGS-M.