# Dynamic Environmental/Economic Scheduling for Microgrid Using Improved MOEA/D-M2M

Authors: Xin Li; Yanjun Fang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The environmental/economic dynamic scheduling for microgrids (MGs) is a complex multiobjective optimization problem, which usually has dynamic system parameters and constraints. In this paper, a biobjective optimization model of MG scheduling is established. And various types of microsources (like the conventional sources, various types of renewable sources, etc.), electricity markets, and dynamic constraints are considered. A recently proposed MOEA/D-M2M framework is improved (I-MOEA/D-M2M) to solve the real-world MG scheduling problems. In order to deal with the constraints, the processes of solutions sorting and selecting in the original MOEA/D-M2M are revised. In addition, a self-adaptive decomposition strategy and a modified allocation method of individuals are introduced to enhance the capability of dealing with uncertainties, as well as reduce unnecessary computational work in practice and meet the time requirements for the dynamic optimization tasks. Thereafter, the proposed I-MOEA/D-M2M is applied to the independent MG scheduling problems, taking into account the load demand variation and the electricity price changes. The simulation results by MATLAB show that the proposed method can achieve better distributed fronts in much less running time than the typical multiobjective evolutionary algorithms (MOEAs) like the improved strength Pareto evolutionary algorithm (SPEA2) and the nondominated sorting genetic algorithm II (NSGAII). Finally, I-MOEA/D-M2M is used to solve a 24-hour MG dynamic operation scheduling problem and obtains satisfactory results.


---

## Body

## 1. Introduction

The energy shortages and environmental pollution have given an impetus to the development of the microgrids (MGs), which are more flexible, energy saving, and environmentally friendly [1–4]. A microgrid system usually contains not only the conventional sources but also various types of renewable sources, such as photovoltaics, wind power, gas turbines, and microturbines [5, 6]. Usually, battery banks are also needed to balance the effects of the volatility of renewable energy. To reduce the cost and emission simultaneously, one of the key strategies is to optimize the output power of the distributed generators (DGs) and the distributed storage (DS) units by using multiobjective optimization approaches, known as MG environmental/economic scheduling. For the engineering application, there are mainly two kinds of MG operation scheduling modes, namely, static mode and dynamic mode. And the dynamic mode can adapt to the complicated changes in the real-world MG systems better than the static one. However, the multiobjective optimization problems (MOPs) of MG dynamic scheduling always have dynamic system parameters, time-varying constraints, and more uncertainties to deal with.Several studies have been made on the improvement of the nonlinear optimization methods, such as dual and quadratic programming [7], and their application on MOPs of scheduling for all kinds of power systems. However, these methods cannot handle nonconvex objective functions [8]. In recent years, many efforts have been made on the development of multiobjective evolutionary algorithms (MOEAs) in dealing with the MG scheduling problems. In [9], an adaptive modified particle swarm optimization (AMPSO) framework was proposed. The application on a day-ahead environmental/economic scheduling problem for MG system showed that the approach performed better in convergence and population diversity than the standard PSO did. In [10], a heuristics-based ant colony optimization (ACO) method was proposed and applied on a multiobjective scheduling problem for MG. The development of the basic ACO had enhanced the capability of handling complex constraints. Besides, MOEAs are widely used in other power systems scheduling problems, which have similarities with MG systems. In [11, 12], the strength Pareto evolutionary algorithm (SPEA) based approaches were presented and used to solve the environmental/economic power dispatch problem for traditional power system and hybrid power system, which includes wind and solar thermal energies. The results show that the methods could solve the MOPs efficiently, obtain well distributed results, and have good population diversities. In [13], the nondominated sorting genetic algorithm II (NSGAII) was improved to solve a biobjective dynamic scheduling problem for a hydrothermal power generation system, which could adequately track Pareto-optimal frontiers online. However, some necessary practical constraints considering time change have not been involved.It is evident from a growing number of researches on MG scheduling MOPs that these problems, especially the dynamic problems, offer serious challenges to the MOEAs for the following reasons:(1)
An MG system usually contains various types of components, which have diverse model functions of fuel cost, maintenance cost, emission, and so forth. Thus, the search space of the MOP is usually nonlinear, nonconvex, and nondifferentiable [14].(2)
In a real-world MG dynamic scheduling problem, there are always dynamic changes in system parameters (like the electricity prices) and constraints (such as the load demand), as a result of which the search space varies with time and the optimization algorithm may not always work well.(3)
In a dynamic procedure, the operators and the decision-makers usually cannot evaluate the quality of the results in a very short period of time without professional optimization knowledge.(4)
The dispatch sequences are dependent on their variation in previous time steps. That is to say, the selected solution will certainly be affected by the former results, and hence the errors between different optimization methods would be accumulated.(5)
In engineering problems, the intermediate trade-off solutions in the Pareto front are what the decision-makers really care about. However, in a dynamic scheduling MOP, since the search spaces using different MOEAs may not be the same every time and the true Pareto front is always unknown, it is difficult to compare the selected intermediate solutions directly.Considering the above challenges, the optimization approaches in handling the MOPs of MG dynamic scheduling should obtain the set of solutions as close to the true Pareto front as possible every time, which means the robustness of the algorithm performance is highly required under uncertainties. However, the algorithm robustness and the causes of failure in solving these MOPs have not been further studied. Meanwhile, the optimization methods need to be simple and self-adaptive to meet the frequent changes of the system parameters.In the remainder of this paper, the system models and objective functions are introduced in Section2. The MOEA/D-M2M algorithm is described and improved in Section 3. Thereafter, the improved algorithm is applied to several practical scenarios to study the effects of different parameters and the algorithm performance by comparing the results with the other two typical MOEAs in Section 4. Finally, the conclusions of the study are presented. The Appendix contains the parameters used in this paper.


## 2. Problem Statement

### 2.1. System Description

MG systems can have various structures. In this paper, the MG consists of three photovoltaic (PV) arrays, three wind turbines (WTs), two microturbines (MTs), two fuel cells (FCs), and a battery (BAT) bank. Assume that the output power of PVs and WTs is uncontrollable and has neither cost nor emission. On the other hand, the MTs and FCs generate power with natural gas as fuel and exhaust emissions, while the output power can be controlled by operators. Battery banks are electrochemical devices that store energy when the DGs generate more power than the load demand and supply the load when the power from DGs and the main grid fails to meet the load.


### 2.2. System Components Models

#### 2.2.1. Uncontrollable DGs

In this MG system, three 10 kW PVs and three 10 kW WTs are considered, of which the output power at timet can be calculated as follows: (1)PPVt=ISTCGAC⁡1+kTc-Tτ,PWTt=0,Vt<Vci,awVt3+bwVt2+cwVt+dw,Vci<Vt<Vr,Pr,Vr<Vt<Vco⁡,0,Vt>Vco⁡,where ISTC is the standard test condition (STC) coefficient; GAC⁡ is the incident irradiance; k is the temperature coefficient of power; Tc and Tτ are the cell temperature and the reference temperature, respectively; Pr is the rated power; Vci and Vco⁡ are the cut-in and cut-out wind speed, respectively; Vr and Vt are the rated and actual wind speed at time t, respectively; aw, bw, cw, and dw are the parameters depending on the wind turbine types. The values of all the parameters above can be found elsewhere [15].


#### 2.2.2. Controllable DGs

In this paper, two MTs (65 kW Capstone C65) and two FCs (40 kW IFC PC-29) are used as controllable DGs in the MG system. Their fuel cost models at timet can be expressed below, and the parameter settings of each DG are taken from [15]. Hence,(2)CG,itPit=Cng·1LHVng·Pitηit,where Pit is the output power of ith DG at time t; Cng is the natural gas price to supply the DGs; LHVng is the fuel lower heating rate; ηit is the efficiency of the DGs at time t.The efficiency curves of MT and FC can be found from manufacturers. Thus, the efficiency functions for the controllable DGs are obtained using curve-fitting methods, which are shown below:(3)ηit=ai·Pit653+bi·Pit652+ci·Pit65+di,i=1,2,ηit=ai·Pit+bi,i=3,4,where ηit is the efficiency at time t. If i=1 or 2, ηit is the efficiency of the MTs; otherwise, it is the efficiency of the FCs. The coefficient values in the models are different depending on the types of DGs, which are shown in Table 3.


#### 2.2.3. Distributed Storage Units

The state of charge (SOC) of the battery bank can be described by calculating the discharge/charge power of the battery, as expressed below:(4)SOC=SOCmax-P-+P+,where SOCmax is the maximum state of charge and P- and P+ are the discharge and charge power, respectively.


### 2.3. Optimization Problem

MG dynamic optimal scheduling is to optimize the output power of every distributed energy resource (DER) to meet the load demand as well as satisfy a series of dynamic constraints through the whole optimization process, while taking into account the economic and environmental effect of the MG system. In this section, the proposed objective functions and constraints are discussed.


#### 2.3.1. The Proposed Objective Functions

One of the main goals of MG dynamic optimal scheduling is to minimize the economic cost at timet, such as fuel cost, maintenance cost, depreciation cost, and the electricity exchange cost with the main grid, which can be described below: (5)min⁡CtP=min⁡∑i=1NCG,itPit+OMitPit+CDP,itPit+CGrid,tPGrid,t,where N is the amount of all the DERs in the MG; P is the decision variable vector; CG,it(Pit) is the fuel cost of ith DER at time t; OMit(Pit) is the maintenance cost of ith DER at time t; CDP,it(Pit) is the depreciation cost of ith DER at time t; PGrid,t is the power exchanged with the grid at time t; CGrid,t(PGrid,t) is the cost of purchased power at time t if PGrid,t>0 or the income of sold power at time t if PGrid,t<0.The maintenance cost and depreciation cost at timet are described as follows, respectively. And all the parameter settings in the functions can be found in Table 3. Hence, (6)OMitPit=KOMiPit,CDP,itPit=ADCCiPr,i×8760×cfi×Pit,where KOMi, ADCCi, Pr,i, and cfi are the maintenance cost factor, the average depreciation cost, the maximum output power, and the capacity factor of the ith DER in the MG, respectively.In addition to the economic cost minimization, this paper mainly takes into account the carbon emission minimization as another optimization objective, as expressed below:(7)min⁡EtP=∑i=1MαiPit2+βiPit+γi,where M is the amount of all the controllable DGs; the values of emission factors αi, βi, and γi can be found in Table 3.


#### 2.3.2. Constraints Description

Power Balance Constraint. All the output power, including the power exchanged with the main grid, should meet the load demand (PL(t)) at time t, which can be expressed as(8)∑i=1MPit+PPVt+PWTt+PBTt-PLt=0.Rated Power Constraints. All the DGs have power generation limits, which can be described as(9)Pimin≤Pit≤Pimax,where Pimin is the minimum power output of ith DG, which is set to be zero in this paper, and Pimax is the maximum power output of the ith DG.State of Charge Constraints. The battery bank cannot be overcharged or overused, so the limits of the state of charge (SOC) of the battery bank are as follows:(10)SOCmin≤SOCt≤SOCmax,where SOCmin and SOCmax are the minimum and maximum state of charge, which are set to be 30 kW and 100 kW, respectively.Power Exchange Constraint. In this paper, the power exchanged between the MG system and the main grid has constraints as described below: (11)PGrid,-≤PGrid,t≤PGrid,+,where PGrid,- and PGrid,+ are the two bounds of exchange power. In this paper, PGrid,- and PGrid,+ are set to be −50 kW and 50 kW, respectively.Ramp Rate Constraint. The increase/decrease of output power of controllable DGs in unit time is called ramp rate, which reflects the performance of the DGs. This constraint can be expressed as (12)Rdown,i·Δt≤Pit-Pi,t-1≤Rup,i·Δt,where Rdown,i and Rup,i are ramp-down and ramp-up rate of the output power of the ith DG, respectively. The values of them are shown in Table 3. Δt is the time of the process.Charge/Discharge Rate Constraint. In this paper, the maximum allowable charge/discharge current is 20% of the battery bank capacity, as shown below:(13)PBAT,+≤0.2×Vsys×UBATΔt,PBAT,-≤0.2×Vsys×UBATΔt,where Vsys is the system voltage at the DC bus and UBAT is the battery bank capacity in AH. The value settings can be found in [15].


## 2.1. System Description

MG systems can have various structures. In this paper, the MG consists of three photovoltaic (PV) arrays, three wind turbines (WTs), two microturbines (MTs), two fuel cells (FCs), and a battery (BAT) bank. Assume that the output power of PVs and WTs is uncontrollable and has neither cost nor emission. On the other hand, the MTs and FCs generate power with natural gas as fuel and exhaust emissions, while the output power can be controlled by operators. Battery banks are electrochemical devices that store energy when the DGs generate more power than the load demand and supply the load when the power from DGs and the main grid fails to meet the load.


## 2.2. System Components Models

### 2.2.1. Uncontrollable DGs

In this MG system, three 10 kW PVs and three 10 kW WTs are considered, of which the output power at timet can be calculated as follows: (1)PPVt=ISTCGAC⁡1+kTc-Tτ,PWTt=0,Vt<Vci,awVt3+bwVt2+cwVt+dw,Vci<Vt<Vr,Pr,Vr<Vt<Vco⁡,0,Vt>Vco⁡,where ISTC is the standard test condition (STC) coefficient; GAC⁡ is the incident irradiance; k is the temperature coefficient of power; Tc and Tτ are the cell temperature and the reference temperature, respectively; Pr is the rated power; Vci and Vco⁡ are the cut-in and cut-out wind speed, respectively; Vr and Vt are the rated and actual wind speed at time t, respectively; aw, bw, cw, and dw are the parameters depending on the wind turbine types. The values of all the parameters above can be found elsewhere [15].


### 2.2.2. Controllable DGs

In this paper, two MTs (65 kW Capstone C65) and two FCs (40 kW IFC PC-29) are used as controllable DGs in the MG system. Their fuel cost models at timet can be expressed below, and the parameter settings of each DG are taken from [15]. Hence,(2)CG,itPit=Cng·1LHVng·Pitηit,where Pit is the output power of ith DG at time t; Cng is the natural gas price to supply the DGs; LHVng is the fuel lower heating rate; ηit is the efficiency of the DGs at time t.The efficiency curves of MT and FC can be found from manufacturers. Thus, the efficiency functions for the controllable DGs are obtained using curve-fitting methods, which are shown below:(3)ηit=ai·Pit653+bi·Pit652+ci·Pit65+di,i=1,2,ηit=ai·Pit+bi,i=3,4,where ηit is the efficiency at time t. If i=1 or 2, ηit is the efficiency of the MTs; otherwise, it is the efficiency of the FCs. The coefficient values in the models are different depending on the types of DGs, which are shown in Table 3.


### 2.2.3. Distributed Storage Units

The state of charge (SOC) of the battery bank can be described by calculating the discharge/charge power of the battery, as expressed below:(4)SOC=SOCmax-P-+P+,where SOCmax is the maximum state of charge and P- and P+ are the discharge and charge power, respectively.


## 2.2.1. Uncontrollable DGs

In this MG system, three 10 kW PVs and three 10 kW WTs are considered, of which the output power at timet can be calculated as follows: (1)PPVt=ISTCGAC⁡1+kTc-Tτ,PWTt=0,Vt<Vci,awVt3+bwVt2+cwVt+dw,Vci<Vt<Vr,Pr,Vr<Vt<Vco⁡,0,Vt>Vco⁡,where ISTC is the standard test condition (STC) coefficient; GAC⁡ is the incident irradiance; k is the temperature coefficient of power; Tc and Tτ are the cell temperature and the reference temperature, respectively; Pr is the rated power; Vci and Vco⁡ are the cut-in and cut-out wind speed, respectively; Vr and Vt are the rated and actual wind speed at time t, respectively; aw, bw, cw, and dw are the parameters depending on the wind turbine types. The values of all the parameters above can be found elsewhere [15].


## 2.2.2. Controllable DGs

In this paper, two MTs (65 kW Capstone C65) and two FCs (40 kW IFC PC-29) are used as controllable DGs in the MG system. Their fuel cost models at timet can be expressed below, and the parameter settings of each DG are taken from [15]. Hence,(2)CG,itPit=Cng·1LHVng·Pitηit,where Pit is the output power of ith DG at time t; Cng is the natural gas price to supply the DGs; LHVng is the fuel lower heating rate; ηit is the efficiency of the DGs at time t.The efficiency curves of MT and FC can be found from manufacturers. Thus, the efficiency functions for the controllable DGs are obtained using curve-fitting methods, which are shown below:(3)ηit=ai·Pit653+bi·Pit652+ci·Pit65+di,i=1,2,ηit=ai·Pit+bi,i=3,4,where ηit is the efficiency at time t. If i=1 or 2, ηit is the efficiency of the MTs; otherwise, it is the efficiency of the FCs. The coefficient values in the models are different depending on the types of DGs, which are shown in Table 3.


## 2.2.3. Distributed Storage Units

The state of charge (SOC) of the battery bank can be described by calculating the discharge/charge power of the battery, as expressed below:(4)SOC=SOCmax-P-+P+,where SOCmax is the maximum state of charge and P- and P+ are the discharge and charge power, respectively.


## 2.3. Optimization Problem

MG dynamic optimal scheduling is to optimize the output power of every distributed energy resource (DER) to meet the load demand as well as satisfy a series of dynamic constraints through the whole optimization process, while taking into account the economic and environmental effect of the MG system. In this section, the proposed objective functions and constraints are discussed.


### 2.3.1. The Proposed Objective Functions

One of the main goals of MG dynamic optimal scheduling is to minimize the economic cost at timet, such as fuel cost, maintenance cost, depreciation cost, and the electricity exchange cost with the main grid, which can be described below: (5)min⁡CtP=min⁡∑i=1NCG,itPit+OMitPit+CDP,itPit+CGrid,tPGrid,t,where N is the amount of all the DERs in the MG; P is the decision variable vector; CG,it(Pit) is the fuel cost of ith DER at time t; OMit(Pit) is the maintenance cost of ith DER at time t; CDP,it(Pit) is the depreciation cost of ith DER at time t; PGrid,t is the power exchanged with the grid at time t; CGrid,t(PGrid,t) is the cost of purchased power at time t if PGrid,t>0 or the income of sold power at time t if PGrid,t<0.The maintenance cost and depreciation cost at timet are described as follows, respectively. And all the parameter settings in the functions can be found in Table 3. Hence, (6)OMitPit=KOMiPit,CDP,itPit=ADCCiPr,i×8760×cfi×Pit,where KOMi, ADCCi, Pr,i, and cfi are the maintenance cost factor, the average depreciation cost, the maximum output power, and the capacity factor of the ith DER in the MG, respectively.In addition to the economic cost minimization, this paper mainly takes into account the carbon emission minimization as another optimization objective, as expressed below:(7)min⁡EtP=∑i=1MαiPit2+βiPit+γi,where M is the amount of all the controllable DGs; the values of emission factors αi, βi, and γi can be found in Table 3.


### 2.3.2. Constraints Description

Power Balance Constraint. All the output power, including the power exchanged with the main grid, should meet the load demand (PL(t)) at time t, which can be expressed as(8)∑i=1MPit+PPVt+PWTt+PBTt-PLt=0.Rated Power Constraints. All the DGs have power generation limits, which can be described as(9)Pimin≤Pit≤Pimax,where Pimin is the minimum power output of ith DG, which is set to be zero in this paper, and Pimax is the maximum power output of the ith DG.State of Charge Constraints. The battery bank cannot be overcharged or overused, so the limits of the state of charge (SOC) of the battery bank are as follows:(10)SOCmin≤SOCt≤SOCmax,where SOCmin and SOCmax are the minimum and maximum state of charge, which are set to be 30 kW and 100 kW, respectively.Power Exchange Constraint. In this paper, the power exchanged between the MG system and the main grid has constraints as described below: (11)PGrid,-≤PGrid,t≤PGrid,+,where PGrid,- and PGrid,+ are the two bounds of exchange power. In this paper, PGrid,- and PGrid,+ are set to be −50 kW and 50 kW, respectively.Ramp Rate Constraint. The increase/decrease of output power of controllable DGs in unit time is called ramp rate, which reflects the performance of the DGs. This constraint can be expressed as (12)Rdown,i·Δt≤Pit-Pi,t-1≤Rup,i·Δt,where Rdown,i and Rup,i are ramp-down and ramp-up rate of the output power of the ith DG, respectively. The values of them are shown in Table 3. Δt is the time of the process.Charge/Discharge Rate Constraint. In this paper, the maximum allowable charge/discharge current is 20% of the battery bank capacity, as shown below:(13)PBAT,+≤0.2×Vsys×UBATΔt,PBAT,-≤0.2×Vsys×UBATΔt,where Vsys is the system voltage at the DC bus and UBAT is the battery bank capacity in AH. The value settings can be found in [15].


## 2.3.1. The Proposed Objective Functions

One of the main goals of MG dynamic optimal scheduling is to minimize the economic cost at timet, such as fuel cost, maintenance cost, depreciation cost, and the electricity exchange cost with the main grid, which can be described below: (5)min⁡CtP=min⁡∑i=1NCG,itPit+OMitPit+CDP,itPit+CGrid,tPGrid,t,where N is the amount of all the DERs in the MG; P is the decision variable vector; CG,it(Pit) is the fuel cost of ith DER at time t; OMit(Pit) is the maintenance cost of ith DER at time t; CDP,it(Pit) is the depreciation cost of ith DER at time t; PGrid,t is the power exchanged with the grid at time t; CGrid,t(PGrid,t) is the cost of purchased power at time t if PGrid,t>0 or the income of sold power at time t if PGrid,t<0.The maintenance cost and depreciation cost at timet are described as follows, respectively. And all the parameter settings in the functions can be found in Table 3. Hence, (6)OMitPit=KOMiPit,CDP,itPit=ADCCiPr,i×8760×cfi×Pit,where KOMi, ADCCi, Pr,i, and cfi are the maintenance cost factor, the average depreciation cost, the maximum output power, and the capacity factor of the ith DER in the MG, respectively.In addition to the economic cost minimization, this paper mainly takes into account the carbon emission minimization as another optimization objective, as expressed below:(7)min⁡EtP=∑i=1MαiPit2+βiPit+γi,where M is the amount of all the controllable DGs; the values of emission factors αi, βi, and γi can be found in Table 3.


## 2.3.2. Constraints Description

Power Balance Constraint. All the output power, including the power exchanged with the main grid, should meet the load demand (PL(t)) at time t, which can be expressed as(8)∑i=1MPit+PPVt+PWTt+PBTt-PLt=0.Rated Power Constraints. All the DGs have power generation limits, which can be described as(9)Pimin≤Pit≤Pimax,where Pimin is the minimum power output of ith DG, which is set to be zero in this paper, and Pimax is the maximum power output of the ith DG.State of Charge Constraints. The battery bank cannot be overcharged or overused, so the limits of the state of charge (SOC) of the battery bank are as follows:(10)SOCmin≤SOCt≤SOCmax,where SOCmin and SOCmax are the minimum and maximum state of charge, which are set to be 30 kW and 100 kW, respectively.Power Exchange Constraint. In this paper, the power exchanged between the MG system and the main grid has constraints as described below: (11)PGrid,-≤PGrid,t≤PGrid,+,where PGrid,- and PGrid,+ are the two bounds of exchange power. In this paper, PGrid,- and PGrid,+ are set to be −50 kW and 50 kW, respectively.Ramp Rate Constraint. The increase/decrease of output power of controllable DGs in unit time is called ramp rate, which reflects the performance of the DGs. This constraint can be expressed as (12)Rdown,i·Δt≤Pit-Pi,t-1≤Rup,i·Δt,where Rdown,i and Rup,i are ramp-down and ramp-up rate of the output power of the ith DG, respectively. The values of them are shown in Table 3. Δt is the time of the process.Charge/Discharge Rate Constraint. In this paper, the maximum allowable charge/discharge current is 20% of the battery bank capacity, as shown below:(13)PBAT,+≤0.2×Vsys×UBATΔt,PBAT,-≤0.2×Vsys×UBATΔt,where Vsys is the system voltage at the DC bus and UBAT is the battery bank capacity in AH. The value settings can be found in [15].


## 3. The I-MOEA/D-M2M Algorithm

In this section, the original MOEA/D-M2M is briefly described, which is applied in this paper as the basic optimization framework. Then, the algorithm is modified and extended to improve the performance in dealing with real-world MOPs.


### 3.1. Description of MOEA/D-M2M

MOEA/D-M2M, introduced by Liu et al. [16] in 2014, is a variant of MOEA/D (MOEA based on decomposition). However, unlike MOEA/D using aggregation methods, MOEA/D-M2M decomposes a MOP into a set of simple sub-MOPs and solves them in one single run, which can be called “multiple to multiple (M2M).” This algorithm framework is designed for some MOPs which render it difficult for other MOEAs to achieve a good diversity of population. The procedure of the algorithm can be found in [16] and is reproduced in Figure 1.Figure 1
Flowchart of MOEA/D-M2M.The major advantages of MOEA/D-M2M include the following [16, 17]:(1)
MOEA/D-M2M transforms a MOP into several sub-MOPs, which is still an equivalent optimization process, whereas the original MOEA/D makes an approximate transformation.(2)
Compared with SPEA2, NSGAII, MOEA/D, and so forth, MOEA/D-M2M “protects” every subpopulation by decomposition and solves every subproblem using multiobjective optimization approaches, which balances the diversity and convergence simultaneously at each generation.(3)
The optimization process requires much less manual operation. The operators only need to set the direction vectors at the very beginning instead of adjusting various complex parameters, which means the algorithm is more adaptive for the MG dynamic scheduling problem with dynamic system parameters and constraints.(4)
The algorithm can be customized with elements from other MOEAs, according to different real-world MOPs. The Pareto-based MOEAs can be applied in the subproblems to handle specific difficulties with some prior knowledge, which can improve the accuracy and efficiency in the whole optimization procedure.


### 3.2. Improvement of MOEA/D-M2M

The original MOEA/D-M2M solved the modified ZDT and DTLZ test instances [16, 18] well, in which there are no constraints and both the search space and the objective space are 0,1n. However, a real-world MOP, such as the MG dynamic scheduling problem, with some equality and inequality constraints can be generally formulated as follows:(14)min⁡Ftx=f1,tx,…,fm,txsubject  togk,tx=0,k=1,…,Khl,tx≤0,l=1,…,L,where f1,t(x),…,fm,t(x) are objective functions at time t; x is the decision vector at time t; gk,t(x) and hl,t(x) are equality and inequality constraints at time t, respectively. The real-world MOPs as described in (14) often have various constraints, dynamic variable bounds with the change of time, and uncertain objective spaces. Thus, the algorithm needs to be modified to meet the practical requirements. Also, some other changes are needed to improve the performance of MOEA/D-M2M.


#### 3.2.1. Constraints Handling

The constraint-handling method, presented by Deb et al. [19], is introduced in this paper during the process of sorting and selection. According to this method, a solution i has a better nondomination rank than solution j, if(1)
solutioni is feasible while solution j is not;(2)
both of them are feasible and solutioni dominates solution j;(3)
both of them are infeasible and solutioni violates the overall constraint less.This constrained-domination principle guarantees that all the feasible solutions have priority to be selected, and the computational complexity of the original algorithm does not change [19].


#### 3.2.2. Self-Adaptive Decomposition Strategy

In [16], the objective space and the Pareto fronts of the test MOPs are always 0,1n; hence, the starting point of the direction vectors can be fixed as the original point. However, in the real-world MOPs with constraints, usually the bounds of the Pareto fronts are unknown and may vary with time. Therefore, to achieve a more precise result and reduce unnecessary computational work in practice, it is important to revise the decomposition strategy.In this paper, the self-adaptive decomposition strategy is formulated based on the change of starting point in each generation. In the decomposition procedure at the first generation, the minimum values for each objectivef1,1,min,…,fm,1,min, which are calculated using all the feasible solutions in the current population, form the first starting point O1(f1,1,min,…,fm,1,min) for the decision variables. For the jth generation, the selected minimum values f1,j,min,…,fm,j,min for each objective in the current population are compared with the associated values of Oj-1 at the (j-1)th generation, respectively, and the smaller ones are picked to form Oj(f1,j,min′,…,fm,j,min′). In this way, the algorithm can focus on the effective search space and avoid being trapped into the useless or infeasible regions, which will improve the efficiency in dealing with the constrained real-world MOPs with uncertain Pareto fronts.


#### 3.2.3. Allocation of Individuals

When allocating individuals after decomposition, the remainingS-Pk solutions are randomly chosen from the combined population Q if Pk<S, according to Figure 1, whereas it is better to select the solutions which have small acute angles with the direction vector in the objective space, since they can conduct to a better front for the subregion Ωk. Thus, the method to allocate the remaining individuals in this situation can be modified as the following steps.Step 1.
Choose two subregions nearest toΩk. If the amount of feasible individuals in the two subregions is sufficient, select the S-Pk ones and add them to Pk; otherwise, select all the feasible individuals and then move to Step 2.Step 2.
Choose two more subregions nearest toΩk, and select the rest of the needed individuals from them. If the quantity of feasible individuals is still insufficient, repeat this step. If all the feasible individuals out of Ωk are selected and the total number is still less than S, move to Step 3.Step 3.
Select the rest of the individuals from the infeasible ones by comparing their overall constraint and add them toPk.


## 3.1. Description of MOEA/D-M2M

MOEA/D-M2M, introduced by Liu et al. [16] in 2014, is a variant of MOEA/D (MOEA based on decomposition). However, unlike MOEA/D using aggregation methods, MOEA/D-M2M decomposes a MOP into a set of simple sub-MOPs and solves them in one single run, which can be called “multiple to multiple (M2M).” This algorithm framework is designed for some MOPs which render it difficult for other MOEAs to achieve a good diversity of population. The procedure of the algorithm can be found in [16] and is reproduced in Figure 1.Figure 1
Flowchart of MOEA/D-M2M.The major advantages of MOEA/D-M2M include the following [16, 17]:(1)
MOEA/D-M2M transforms a MOP into several sub-MOPs, which is still an equivalent optimization process, whereas the original MOEA/D makes an approximate transformation.(2)
Compared with SPEA2, NSGAII, MOEA/D, and so forth, MOEA/D-M2M “protects” every subpopulation by decomposition and solves every subproblem using multiobjective optimization approaches, which balances the diversity and convergence simultaneously at each generation.(3)
The optimization process requires much less manual operation. The operators only need to set the direction vectors at the very beginning instead of adjusting various complex parameters, which means the algorithm is more adaptive for the MG dynamic scheduling problem with dynamic system parameters and constraints.(4)
The algorithm can be customized with elements from other MOEAs, according to different real-world MOPs. The Pareto-based MOEAs can be applied in the subproblems to handle specific difficulties with some prior knowledge, which can improve the accuracy and efficiency in the whole optimization procedure.


## 3.2. Improvement of MOEA/D-M2M

The original MOEA/D-M2M solved the modified ZDT and DTLZ test instances [16, 18] well, in which there are no constraints and both the search space and the objective space are 0,1n. However, a real-world MOP, such as the MG dynamic scheduling problem, with some equality and inequality constraints can be generally formulated as follows:(14)min⁡Ftx=f1,tx,…,fm,txsubject  togk,tx=0,k=1,…,Khl,tx≤0,l=1,…,L,where f1,t(x),…,fm,t(x) are objective functions at time t; x is the decision vector at time t; gk,t(x) and hl,t(x) are equality and inequality constraints at time t, respectively. The real-world MOPs as described in (14) often have various constraints, dynamic variable bounds with the change of time, and uncertain objective spaces. Thus, the algorithm needs to be modified to meet the practical requirements. Also, some other changes are needed to improve the performance of MOEA/D-M2M.


### 3.2.1. Constraints Handling

The constraint-handling method, presented by Deb et al. [19], is introduced in this paper during the process of sorting and selection. According to this method, a solution i has a better nondomination rank than solution j, if(1)
solutioni is feasible while solution j is not;(2)
both of them are feasible and solutioni dominates solution j;(3)
both of them are infeasible and solutioni violates the overall constraint less.This constrained-domination principle guarantees that all the feasible solutions have priority to be selected, and the computational complexity of the original algorithm does not change [19].


### 3.2.2. Self-Adaptive Decomposition Strategy

In [16], the objective space and the Pareto fronts of the test MOPs are always 0,1n; hence, the starting point of the direction vectors can be fixed as the original point. However, in the real-world MOPs with constraints, usually the bounds of the Pareto fronts are unknown and may vary with time. Therefore, to achieve a more precise result and reduce unnecessary computational work in practice, it is important to revise the decomposition strategy.In this paper, the self-adaptive decomposition strategy is formulated based on the change of starting point in each generation. In the decomposition procedure at the first generation, the minimum values for each objectivef1,1,min,…,fm,1,min, which are calculated using all the feasible solutions in the current population, form the first starting point O1(f1,1,min,…,fm,1,min) for the decision variables. For the jth generation, the selected minimum values f1,j,min,…,fm,j,min for each objective in the current population are compared with the associated values of Oj-1 at the (j-1)th generation, respectively, and the smaller ones are picked to form Oj(f1,j,min′,…,fm,j,min′). In this way, the algorithm can focus on the effective search space and avoid being trapped into the useless or infeasible regions, which will improve the efficiency in dealing with the constrained real-world MOPs with uncertain Pareto fronts.


### 3.2.3. Allocation of Individuals

When allocating individuals after decomposition, the remainingS-Pk solutions are randomly chosen from the combined population Q if Pk<S, according to Figure 1, whereas it is better to select the solutions which have small acute angles with the direction vector in the objective space, since they can conduct to a better front for the subregion Ωk. Thus, the method to allocate the remaining individuals in this situation can be modified as the following steps.Step 1.
Choose two subregions nearest toΩk. If the amount of feasible individuals in the two subregions is sufficient, select the S-Pk ones and add them to Pk; otherwise, select all the feasible individuals and then move to Step 2.Step 2.
Choose two more subregions nearest toΩk, and select the rest of the needed individuals from them. If the quantity of feasible individuals is still insufficient, repeat this step. If all the feasible individuals out of Ωk are selected and the total number is still less than S, move to Step 3.Step 3.
Select the rest of the individuals from the infeasible ones by comparing their overall constraint and add them toPk.


## 3.2.1. Constraints Handling

The constraint-handling method, presented by Deb et al. [19], is introduced in this paper during the process of sorting and selection. According to this method, a solution i has a better nondomination rank than solution j, if(1)
solutioni is feasible while solution j is not;(2)
both of them are feasible and solutioni dominates solution j;(3)
both of them are infeasible and solutioni violates the overall constraint less.This constrained-domination principle guarantees that all the feasible solutions have priority to be selected, and the computational complexity of the original algorithm does not change [19].


## 3.2.2. Self-Adaptive Decomposition Strategy

In [16], the objective space and the Pareto fronts of the test MOPs are always 0,1n; hence, the starting point of the direction vectors can be fixed as the original point. However, in the real-world MOPs with constraints, usually the bounds of the Pareto fronts are unknown and may vary with time. Therefore, to achieve a more precise result and reduce unnecessary computational work in practice, it is important to revise the decomposition strategy.In this paper, the self-adaptive decomposition strategy is formulated based on the change of starting point in each generation. In the decomposition procedure at the first generation, the minimum values for each objectivef1,1,min,…,fm,1,min, which are calculated using all the feasible solutions in the current population, form the first starting point O1(f1,1,min,…,fm,1,min) for the decision variables. For the jth generation, the selected minimum values f1,j,min,…,fm,j,min for each objective in the current population are compared with the associated values of Oj-1 at the (j-1)th generation, respectively, and the smaller ones are picked to form Oj(f1,j,min′,…,fm,j,min′). In this way, the algorithm can focus on the effective search space and avoid being trapped into the useless or infeasible regions, which will improve the efficiency in dealing with the constrained real-world MOPs with uncertain Pareto fronts.


## 3.2.3. Allocation of Individuals

When allocating individuals after decomposition, the remainingS-Pk solutions are randomly chosen from the combined population Q if Pk<S, according to Figure 1, whereas it is better to select the solutions which have small acute angles with the direction vector in the objective space, since they can conduct to a better front for the subregion Ωk. Thus, the method to allocate the remaining individuals in this situation can be modified as the following steps.Step 1.
Choose two subregions nearest toΩk. If the amount of feasible individuals in the two subregions is sufficient, select the S-Pk ones and add them to Pk; otherwise, select all the feasible individuals and then move to Step 2.Step 2.
Choose two more subregions nearest toΩk, and select the rest of the needed individuals from them. If the quantity of feasible individuals is still insufficient, repeat this step. If all the feasible individuals out of Ωk are selected and the total number is still less than S, move to Step 3.Step 3.
Select the rest of the individuals from the infeasible ones by comparing their overall constraint and add them toPk.


## 4. Simulation Results and Discussion

In this section, the MG optimization model described in the previous section is applied in three scenarios, considering the load demand variation, the electricity price change, and a 24-hour dynamic optimization process. In each scenario, three MOEAs including I-MOEA/D-M2M, NSGAII, and SPEA2 [20] are used and their optimization results are compared. The population size is 100 and the maximum generation number is 500 for each algorithm. For SPEA2, the archive size is 100. And, for I-MOEA/D-M2M, K=S=10. The simulated binary crossover [21, 22] and polynomial mutation methods are applied in the three algorithms with the same control parameters in [19]. And both NSGAII and SPEA2 use the same constraint-handling strategy as I-MOEA/D-M2M does. The reference point of the hypervolume values is (30,10) based on plenty of experimental results. All the experiments are accomplished using MATLAB 2013, on a PC with an Intel Core i3 (3.30 GHz) processor under Windows 7 using 6 GB of RAM.


### 4.1. Scenario One

In Scenario One, the load demand variation is considered as the primary factor to affect the optimization results. In this case, the other system parameters, like the environmental parameters, the electricity price, and so forth, are fixed, which are shown in Table4. The ramp rate constraint is not applied here since the optimization procedures are independent.Figures2(a)–2(i) show the optimization results under different load demands using I-MOEA/D-M2M, SPEA2, and NSGAII, respectively. The red points represent the population in the last generation, and the black ones are the feasible nondominated solutions. The comparisons of hypervolume values (h-values), feasible extreme solutions, and average CPU time are presented in Table 1.Table 1
Theh-values, extreme solutions, and CPU time using the three algorithms under different load demands.
Load demand (kW)
I-MOEA/D-M2M
SPEA2
NSGAII
50
h-value
Best
278.5705
278.3354
276.4289
Mean
278.5012
274.0246
272.5123
Extreme solutions
For obj11
(2.1137,  0.7416)
(2.1145,  0.7451)
(2.4410, 0.2894)
For obj2
(4.7214,  0.0012)
(4.6002,  0.0013)
(7.5578, 0.0015)
100
h-value
Best
261.0048
260.8247
257.4370
Mean
260.9123
257.2458
252.4999
Extreme solutions
For obj1
(3.8618,  0.7431)
(3.8788,  0.7358)
(4.1293, 0.3158)
For obj2
(6.9167,  0.0013)
(6.9822,  0.0015)
(8.2485, 0.0018)
150
h-value
Best
204.6214
201.2211
195.5587
Mean
202.2548
196.6741
187.5855
Extreme solutions
For obj1
(9.1871,  0.9048)
(9.3121,  0.8863)
(9.2865, 0.8563)
For obj2
(13.4595,  0.1249)
(13.2509,  0.1386)
(13.2488, 0.2360)
Average CPU time (s)
49.1254
292.5186
78.3404
1The terms obj1 and obj2 represent to the two optimization objectives, respectively.Figure 2
Columns from left to right:(1) the optimization results using I-MOEA/D-M2M; (2) the optimization results using SPEA2; (3) the optimization results using NSGAII. Rows: (1) the optimization results under 50 kW load demand; (2) the optimization results under 100 kW load demand; (3) the optimization results under 150 kW load demand.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)It is clear from Figures2(a), 2(d), and 2(g) that the results achieved by I-MOEA/D-M2M are always well distributed. According to Table 1, the best and mean h-values have very small difference when using I-MOEA/D-M2M, which means the algorithm can converge to the Pareto fronts in less iteration. The extreme solutions for each objective are the best compared with those using the other two algorithms. Therefore, this algorithm is more robust and can provide more diverse trade-off solutions for the decision-makers under various load demands. And since it needs much less CPU time, I-MOEA/D-M2M is more suitable for MG dynamic scheduling problems.For SPEA2, when the load is 50 kW, the population is distributed uniformly. The besth-value is similar to that obtained by I-MOEA/D-M2M, but the mean value is much smaller. The minimum cost and minimum emission do not have much difference compared with those obtained by I-MOEA/D-M2M. For instance, the minimum cost is 2.1147$, which is only 0.038% more than the I-MOEA/D-M2M optimized solution. With the load increasing, SPEA2 loses the population diversity gradually (shown in Figures 2(b), 2(e), and 2(h)). It can be seen from Figure 2(h) that part of the population gets trapped into the infeasible region using SPEA2 when the load is 150 kW. On the other hand, the h-values and the minimum cost and minimum emission are smaller than those using I-MOEA/D-M2M, and the difference between them increases with the load. Also, it should be noted that SPEA2 costs much more time (292.5186 s) than the other two methods do. Although it performs better than NSGAII within 500 generations, SPEA2 is not a good choice in dealing with this MG dynamic MOP.For NSGAII, the change is similar to that using SPEA2, while theh-values and the obtained extreme solutions are worse. It can be seen from Figures 2(c), 2(f), and 2(i) that the amount of the feasible nondominated solutions reduces with the load increasing, and some of the solutions always aggregate in a small region. For example, when the load demand is 100 kW, 45% of the population have similar costs in the objective space, which are between 4.1624$ and 4.2843$. In this paper, this region is called Region A. To analyze this, all the feasible solutions in Figure 2(f) are shown in Figure 3 by parallel coordinates. It is obvious that only the third and the fourth variables (i.e., output of FC1 and FC2) of the red solutions have different values, while others remain unchanged as their boundary values. Hence, it is easier to find nondominated solutions in Region A since only two variables need to be changed. As a result, NSGAII may get trapped into the regions like this, as convergence takes precedence over population diversity in its search strategy. The iterative processes of NSGAII and I-MOEA/D-M2M are shown in Figure 4, which describe how the population varies with generations.Figure 3
The feasible solutions obtained by NSGAII under 100 kW load demand.Figure 4
The iterative processes of NSGAII and I-MOEA/D-M2M under 100 kW load demand. Rows:(1) the optimization results after 600, 1000, and 2000 generations using NSGAII; (2) the optimization results after 100, 300, and 500 generations using I-MOEA/D-M2M.
(a)
(b)
(c)
(d)
(e)
(f)It can be seen from Figure4(a) that, at the 600th generation, more than half of the population stays in Region A by using NSGAII, and only a few solutions have the cost values higher than 4.2843$ in the objective space. In this paper, the region is called Region B, where the cost value is higher than 4.2843$. After 400 generations, the majority of the population extends to Region B and is distributed uniformly, and there are only a few solutions which have the cost values lower than 4.1624$ in the objective space. The region is called Region C, where the cost value is lower than 4.1624$. When the generation reaches 2000, there are more solutions in Region C and the whole front is similar to that obtained by I-MOEA/D-M2M. Although NSGAII can also achieve a good result finally, it costs more running time. Furthermore, when the load changes, it is difficult to decide how many generations the algorithm needs to get satisfying results, like the case in Figure 2(i). With the load demand constraints changing, the feasible search space varies frequently and uncertainly, and NSGAII cannot guarantee to find the Pareto front every time, since it may get trapped into some regions like Region A. SPEA2 cannot deal with these problems either (shown in Figure 2(h)) because they both adopt the convergence first and diversity second selection strategy [16], although SPEA2 performs better than NSGAII.In contrast, I-MOEA/D-M2M can get a well distributed front within just 100 generations due to its decomposition strategy, and the results at 300th generation are almost the same as those at 500th generation, which means the algorithm is much faster in handling these MOPs. By using the self-adaptive decomposition strategy, I-MOEA/D-M2M can avoid searching the useless regions and pay more attention to population diversity. So, whatever the load demand is, I-MOEA/D-M2M can always quickly obtain the feasible trade-off solutions.


### 4.2. Scenario Two

In Scenario Two, the effect of electricity price change is studied. The environmental parameters are held constant as shown in Table4, and the load demand is 150 kW. The ramp rate constraint is still not applied here as Scenario One. The electricity prices settings are shown in Table 5. Figure 5 shows the optimization results in three cases using I-MOEA/D-M2M, SPEA2, and NSGAII, respectively.Figure 5
Columns from left to right:(1) the optimization results using I-MOEA/D-M2M; (2) the optimization results using SPEA2; (3) the optimization results using NSGAII. Rows: (1) the optimization results in Case One; (2) the optimization results in Case Two; (3) the optimization results in Case Three.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)It is evident that I-MOEA/D-M2M performs well in all the three cases, the solutions obtained are distributed uniformly in the objective space, and it can be seen from Figures5(a), 5(d), and 5(g) that the fronts are quite different in the three MOPs. For SPEA2, when the electricity purchase price drops to 0.12$/kWh, the quality of the results is better than that in Case One since all the nondominated solutions are feasible, and the front is similar to that using I-MOEA/D-M2M, whereas the amount of the feasible nondominated solutions reduces when the electricity sale price increases to 0.13$/kWh, and nearly half of the solutions are infeasible. For NSGAII, no matter how the electricity prices change, it only finds part of the Pareto front within 500 generations, according to Figures 5(c), 5(f), and 5(i). And there are a large number of infeasible or dominated solutions left. In conclusion, I-MOEA/D-M2M can adapt to the MG MOP with electricity prices varying, while the quality of results using SPEA2 is based on the combination of electricity purchase and sale prices, and NSGAII may need far more than 500 generations to get similar Pareto fronts to I-MOEA/D-M2M.


### 4.3. Scenario Three

With the above independent optimization runs and analyses of the algorithms performance, it is time to demonstrate the use of the proposed optimization methodology to a complex 24-hour MG dynamic operation scheduling problem involving economic and environmental considerations. The load demand in a typical day is shown in Figure6, and the environmental data for wind speed, irradiance, and air temperature can be found in [15]. Besides, the time-varying electricity prices are considered, shown in Table 6.Figure 6
Variation of load demand in 24 hours.The pseudoweight vector method [23, 24] is introduced in this paper to select the solutions based on the decision-makers’ preference. After the optimization process in each hour, the pseudoweight ωi for the ith objective for every solution is calculated using the following equation:(15)ωi=fimax-fi/fimax-fimin∑j=1Mfjmax-fj/fjmax-fjmin,where fimax and fimin are the maximum and minimum values for the ith objective in the obtained results, respectively; fi is the value of the ith objective function for every solution; M is the number of the optimization objectives. Thereafter, the following decision function can be applied:(16)Min⁡dn=ωn-ω∗,where dn and ωn are the decision function value and the pseudoweight vector for the nth solution in the final population and ω∗ is the weight vector set by the decision-maker. The solution with the minimum decision function value will be selected, which is satisfied by the decision-maker. In this paper, three cases are studied, of which the ω∗ values are (1,0), (0,1), and (0.5,0.5), respectively. The optimization results by I-MOEA/D-M2M are shown in Figure 7.Figure 7
Optimization results using I-MOEA/D-M2M in the three cases of the 24-hour MG dynamic operation scheduling problem.
(a)
(b)
(c)Figure7(a) shows the optimization results for the 24-hour MG dynamic operation scheduling problem when ω∗=(1,0), which means the solution with the minimum cost in the objective space is selected in each hour. The power output curves of PVs and WTs are not shown here since they are used first in each case and can be calculated by the forecast environmental data. It can be seen that, for the two MTs, MT1 is used first and MT2 starts to work only if MT1 outputs full power under the relative constraints. For the two FCs, similarly, FC1 always generates more power than FC2, while sometimes FC2 starts running even though FC1 does not reach the full power under its constraints. The differences can be explained by computing the average prices of power generated by the DGs using (2). For the BAT, to avoid frequent charging and discharging, it will be charged when the state is 30 kW (SOCmin) and will not discharge until it reaches 100 kW (SOCmax). Since there is no need to provide more power than the load demand by the other DGs in the MG in this paper, the BAT is mainly charged using the free power obtained by uncontrollable DGs. As a result, the output curve of BAT does not change much in the three cases. For the power exchanged with the main grid, the value depends on the electricity price and the average prices of output power of the controllable DGs. For instance, it can be calculated by (5) that no matter how much power the DGs output the average prices are always higher than the electricity purchase price in valley period (00:00–07:00). Thus, it is more economical to buy power from the grid before using the components in the MG.Figure7(b) shows the optimization results when ω∗=(1,0), which means minimizing the emission level is the most concerned objective. Interestingly, the output curves of the two MTs have the same varying tendency under different load demands, while MT2 supplies a little more power (less than 5 kW) than MT1 does every time. For the FCs, the situation is similar. This is quite different from Case One but can also be explained by (7). For the power exchanged with the main grid, the value is always the maximum since it is assumed to produce no harmful emission in this paper.Figure7(c) shows the optimization results when ω∗=(0.5,0.5), which consider trade-offs between cost and emission. It can be seen that the curves of the power output by the DGs and the power exchanged with the main grid change irregularly with time, which cannot be derived mathematically by the MG models in Section 2. However, in practical engineering problems, the intermediate solutions like those in Case Three are the ones the decision-makers really care about. Therefore, the extreme solutions obtained in Case One and Case Two are of important reference significance for evaluating the performance of the optimization method. With the above simulations, it is clear that I-MOEA/D-M2M does better in searching the extreme solutions than the other two algorithms. Table 2 shows the total cost and emission obtained by the three algorithms with different ω∗ values. It can be seen that the results using I-MOEA/D-M2M have got the minimum total cost in Case One and total emission in Case Two compared with those using the other algorithms. Particularly, in Case Two, the minimum emission solution using I-MOEA/D-M2M incurs a total emission of 19.5747 kg, which is 13.47% and 27.57% better than the SPEA2 and NSGAII optimized solutions, respectively. In Case Three, the cost and emission levels by I-MOEA/D-M2M are lower than those by NSGAII. But the emission by I-MOEA/D-M2M is more than that by SPEA2. According to Case One and Case Two, the former result is still the first choice compared with the latter, since I-MOEA/D-M2M can approximate the true Pareto front better than SPEA2 and reflect the decision-makers’ preference more precisely.Table 2
The total cost and emission obtained by the three algorithms in different cases.
Case One
Case Two
Case Three
Cost ($)
Emission (kg)
Cost ($)
Emission (kg)
Cost ($)
Emission (kg)
Results
I-MOEA/D-M2M
323.8993
45.0830
370.3897
19.5747
339.7723
24.6929
SPEA2
337.6289
45.1581
373.3720
22.6218
343.6473
23.9759
NSGAII
354.7322
48.3517
388.6219
27.0249
359.3728
31.3728Table 3
Parameter settings of the MG system components.
Parameters
MT
FC
BAT
PV
WT
i
1
2
3
4
5
6–8
9–11
K
O
M
i ($/kWh)
7.58 × 10−3
9.35 × 10−3
16.12 × 10−3
17.74 × 10−3
7.3 × 10−3
0
0
A
D
C
C
i ($)
342.82
342.82
1222.2
1222.2
22.84
1259.8
623.41
c
f
i (%)
54.99
54.99
36.73
36.73
32.67
29.34
22.13
a
i
0.752
0.0692
−0.0023
−0.0026
—
—
—
b
i
−0.3095
−0.2904
0.6735
0.5821
—
—
—
c
i
0.4174
0.3752
—
—
—
—
—
d
i
0.1068
0.1090
—
—
—
—
—
α
i (kg/kW)
4.09 × 10−4
2.54 × 10−4
5.33 × 10−4
4.87 × 10−4
—
—
—
β
i (kg/kW)
−5.55 × 10−4
−6.05 × 10−4
−3.55 × 10−4
−4.23 × 10−4
—
—
—
γ
i (kg/kW)
6.49 × 10−4
5.64 × 10−4
3.38 × 10−4
2.46 × 10−4
—
—
—
Ramp rate (kW/min)
0.33
0.33
0.25
0.25
—
—
—Table 4
Parameter settings of the environmental factors and the electricity price.
Wind speed (m/s)
Irradiance (W/m2)
Air temperature (°C)
Electricity sale price ($/kWh)
Electricity purchase price ($/kWh)
5.5
400
21
0.15
0.11Table 5
Parameter settings of the electricity prices.
Case One
Case Two
Case Three
Electricity sale price ($/kWh)
0.11
0.11
0.13
Electricity purchase price ($/kWh)
0.15
0.12
0.15Table 6
The electricity price settings in 24 hours.
Period
Time (h)
Price ($/kW)
Sale
Purchase
Peak period
10:00–15:0018:00–21:00
0.11
0.15
Normal period
07:00–10:0015:00–18:0021:00–23:00
0.06
0.08
Valley period
00:00–07:0023:00-24:00
0.02
0.03


## 4.1. Scenario One

In Scenario One, the load demand variation is considered as the primary factor to affect the optimization results. In this case, the other system parameters, like the environmental parameters, the electricity price, and so forth, are fixed, which are shown in Table4. The ramp rate constraint is not applied here since the optimization procedures are independent.Figures2(a)–2(i) show the optimization results under different load demands using I-MOEA/D-M2M, SPEA2, and NSGAII, respectively. The red points represent the population in the last generation, and the black ones are the feasible nondominated solutions. The comparisons of hypervolume values (h-values), feasible extreme solutions, and average CPU time are presented in Table 1.Table 1
Theh-values, extreme solutions, and CPU time using the three algorithms under different load demands.
Load demand (kW)
I-MOEA/D-M2M
SPEA2
NSGAII
50
h-value
Best
278.5705
278.3354
276.4289
Mean
278.5012
274.0246
272.5123
Extreme solutions
For obj11
(2.1137,  0.7416)
(2.1145,  0.7451)
(2.4410, 0.2894)
For obj2
(4.7214,  0.0012)
(4.6002,  0.0013)
(7.5578, 0.0015)
100
h-value
Best
261.0048
260.8247
257.4370
Mean
260.9123
257.2458
252.4999
Extreme solutions
For obj1
(3.8618,  0.7431)
(3.8788,  0.7358)
(4.1293, 0.3158)
For obj2
(6.9167,  0.0013)
(6.9822,  0.0015)
(8.2485, 0.0018)
150
h-value
Best
204.6214
201.2211
195.5587
Mean
202.2548
196.6741
187.5855
Extreme solutions
For obj1
(9.1871,  0.9048)
(9.3121,  0.8863)
(9.2865, 0.8563)
For obj2
(13.4595,  0.1249)
(13.2509,  0.1386)
(13.2488, 0.2360)
Average CPU time (s)
49.1254
292.5186
78.3404
1The terms obj1 and obj2 represent to the two optimization objectives, respectively.Figure 2
Columns from left to right:(1) the optimization results using I-MOEA/D-M2M; (2) the optimization results using SPEA2; (3) the optimization results using NSGAII. Rows: (1) the optimization results under 50 kW load demand; (2) the optimization results under 100 kW load demand; (3) the optimization results under 150 kW load demand.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)It is clear from Figures2(a), 2(d), and 2(g) that the results achieved by I-MOEA/D-M2M are always well distributed. According to Table 1, the best and mean h-values have very small difference when using I-MOEA/D-M2M, which means the algorithm can converge to the Pareto fronts in less iteration. The extreme solutions for each objective are the best compared with those using the other two algorithms. Therefore, this algorithm is more robust and can provide more diverse trade-off solutions for the decision-makers under various load demands. And since it needs much less CPU time, I-MOEA/D-M2M is more suitable for MG dynamic scheduling problems.For SPEA2, when the load is 50 kW, the population is distributed uniformly. The besth-value is similar to that obtained by I-MOEA/D-M2M, but the mean value is much smaller. The minimum cost and minimum emission do not have much difference compared with those obtained by I-MOEA/D-M2M. For instance, the minimum cost is 2.1147$, which is only 0.038% more than the I-MOEA/D-M2M optimized solution. With the load increasing, SPEA2 loses the population diversity gradually (shown in Figures 2(b), 2(e), and 2(h)). It can be seen from Figure 2(h) that part of the population gets trapped into the infeasible region using SPEA2 when the load is 150 kW. On the other hand, the h-values and the minimum cost and minimum emission are smaller than those using I-MOEA/D-M2M, and the difference between them increases with the load. Also, it should be noted that SPEA2 costs much more time (292.5186 s) than the other two methods do. Although it performs better than NSGAII within 500 generations, SPEA2 is not a good choice in dealing with this MG dynamic MOP.For NSGAII, the change is similar to that using SPEA2, while theh-values and the obtained extreme solutions are worse. It can be seen from Figures 2(c), 2(f), and 2(i) that the amount of the feasible nondominated solutions reduces with the load increasing, and some of the solutions always aggregate in a small region. For example, when the load demand is 100 kW, 45% of the population have similar costs in the objective space, which are between 4.1624$ and 4.2843$. In this paper, this region is called Region A. To analyze this, all the feasible solutions in Figure 2(f) are shown in Figure 3 by parallel coordinates. It is obvious that only the third and the fourth variables (i.e., output of FC1 and FC2) of the red solutions have different values, while others remain unchanged as their boundary values. Hence, it is easier to find nondominated solutions in Region A since only two variables need to be changed. As a result, NSGAII may get trapped into the regions like this, as convergence takes precedence over population diversity in its search strategy. The iterative processes of NSGAII and I-MOEA/D-M2M are shown in Figure 4, which describe how the population varies with generations.Figure 3
The feasible solutions obtained by NSGAII under 100 kW load demand.Figure 4
The iterative processes of NSGAII and I-MOEA/D-M2M under 100 kW load demand. Rows:(1) the optimization results after 600, 1000, and 2000 generations using NSGAII; (2) the optimization results after 100, 300, and 500 generations using I-MOEA/D-M2M.
(a)
(b)
(c)
(d)
(e)
(f)It can be seen from Figure4(a) that, at the 600th generation, more than half of the population stays in Region A by using NSGAII, and only a few solutions have the cost values higher than 4.2843$ in the objective space. In this paper, the region is called Region B, where the cost value is higher than 4.2843$. After 400 generations, the majority of the population extends to Region B and is distributed uniformly, and there are only a few solutions which have the cost values lower than 4.1624$ in the objective space. The region is called Region C, where the cost value is lower than 4.1624$. When the generation reaches 2000, there are more solutions in Region C and the whole front is similar to that obtained by I-MOEA/D-M2M. Although NSGAII can also achieve a good result finally, it costs more running time. Furthermore, when the load changes, it is difficult to decide how many generations the algorithm needs to get satisfying results, like the case in Figure 2(i). With the load demand constraints changing, the feasible search space varies frequently and uncertainly, and NSGAII cannot guarantee to find the Pareto front every time, since it may get trapped into some regions like Region A. SPEA2 cannot deal with these problems either (shown in Figure 2(h)) because they both adopt the convergence first and diversity second selection strategy [16], although SPEA2 performs better than NSGAII.In contrast, I-MOEA/D-M2M can get a well distributed front within just 100 generations due to its decomposition strategy, and the results at 300th generation are almost the same as those at 500th generation, which means the algorithm is much faster in handling these MOPs. By using the self-adaptive decomposition strategy, I-MOEA/D-M2M can avoid searching the useless regions and pay more attention to population diversity. So, whatever the load demand is, I-MOEA/D-M2M can always quickly obtain the feasible trade-off solutions.


## 4.2. Scenario Two

In Scenario Two, the effect of electricity price change is studied. The environmental parameters are held constant as shown in Table4, and the load demand is 150 kW. The ramp rate constraint is still not applied here as Scenario One. The electricity prices settings are shown in Table 5. Figure 5 shows the optimization results in three cases using I-MOEA/D-M2M, SPEA2, and NSGAII, respectively.Figure 5
Columns from left to right:(1) the optimization results using I-MOEA/D-M2M; (2) the optimization results using SPEA2; (3) the optimization results using NSGAII. Rows: (1) the optimization results in Case One; (2) the optimization results in Case Two; (3) the optimization results in Case Three.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)It is evident that I-MOEA/D-M2M performs well in all the three cases, the solutions obtained are distributed uniformly in the objective space, and it can be seen from Figures5(a), 5(d), and 5(g) that the fronts are quite different in the three MOPs. For SPEA2, when the electricity purchase price drops to 0.12$/kWh, the quality of the results is better than that in Case One since all the nondominated solutions are feasible, and the front is similar to that using I-MOEA/D-M2M, whereas the amount of the feasible nondominated solutions reduces when the electricity sale price increases to 0.13$/kWh, and nearly half of the solutions are infeasible. For NSGAII, no matter how the electricity prices change, it only finds part of the Pareto front within 500 generations, according to Figures 5(c), 5(f), and 5(i). And there are a large number of infeasible or dominated solutions left. In conclusion, I-MOEA/D-M2M can adapt to the MG MOP with electricity prices varying, while the quality of results using SPEA2 is based on the combination of electricity purchase and sale prices, and NSGAII may need far more than 500 generations to get similar Pareto fronts to I-MOEA/D-M2M.


## 4.3. Scenario Three

With the above independent optimization runs and analyses of the algorithms performance, it is time to demonstrate the use of the proposed optimization methodology to a complex 24-hour MG dynamic operation scheduling problem involving economic and environmental considerations. The load demand in a typical day is shown in Figure6, and the environmental data for wind speed, irradiance, and air temperature can be found in [15]. Besides, the time-varying electricity prices are considered, shown in Table 6.Figure 6
Variation of load demand in 24 hours.The pseudoweight vector method [23, 24] is introduced in this paper to select the solutions based on the decision-makers’ preference. After the optimization process in each hour, the pseudoweight ωi for the ith objective for every solution is calculated using the following equation:(15)ωi=fimax-fi/fimax-fimin∑j=1Mfjmax-fj/fjmax-fjmin,where fimax and fimin are the maximum and minimum values for the ith objective in the obtained results, respectively; fi is the value of the ith objective function for every solution; M is the number of the optimization objectives. Thereafter, the following decision function can be applied:(16)Min⁡dn=ωn-ω∗,where dn and ωn are the decision function value and the pseudoweight vector for the nth solution in the final population and ω∗ is the weight vector set by the decision-maker. The solution with the minimum decision function value will be selected, which is satisfied by the decision-maker. In this paper, three cases are studied, of which the ω∗ values are (1,0), (0,1), and (0.5,0.5), respectively. The optimization results by I-MOEA/D-M2M are shown in Figure 7.Figure 7
Optimization results using I-MOEA/D-M2M in the three cases of the 24-hour MG dynamic operation scheduling problem.
(a)
(b)
(c)Figure7(a) shows the optimization results for the 24-hour MG dynamic operation scheduling problem when ω∗=(1,0), which means the solution with the minimum cost in the objective space is selected in each hour. The power output curves of PVs and WTs are not shown here since they are used first in each case and can be calculated by the forecast environmental data. It can be seen that, for the two MTs, MT1 is used first and MT2 starts to work only if MT1 outputs full power under the relative constraints. For the two FCs, similarly, FC1 always generates more power than FC2, while sometimes FC2 starts running even though FC1 does not reach the full power under its constraints. The differences can be explained by computing the average prices of power generated by the DGs using (2). For the BAT, to avoid frequent charging and discharging, it will be charged when the state is 30 kW (SOCmin) and will not discharge until it reaches 100 kW (SOCmax). Since there is no need to provide more power than the load demand by the other DGs in the MG in this paper, the BAT is mainly charged using the free power obtained by uncontrollable DGs. As a result, the output curve of BAT does not change much in the three cases. For the power exchanged with the main grid, the value depends on the electricity price and the average prices of output power of the controllable DGs. For instance, it can be calculated by (5) that no matter how much power the DGs output the average prices are always higher than the electricity purchase price in valley period (00:00–07:00). Thus, it is more economical to buy power from the grid before using the components in the MG.Figure7(b) shows the optimization results when ω∗=(1,0), which means minimizing the emission level is the most concerned objective. Interestingly, the output curves of the two MTs have the same varying tendency under different load demands, while MT2 supplies a little more power (less than 5 kW) than MT1 does every time. For the FCs, the situation is similar. This is quite different from Case One but can also be explained by (7). For the power exchanged with the main grid, the value is always the maximum since it is assumed to produce no harmful emission in this paper.Figure7(c) shows the optimization results when ω∗=(0.5,0.5), which consider trade-offs between cost and emission. It can be seen that the curves of the power output by the DGs and the power exchanged with the main grid change irregularly with time, which cannot be derived mathematically by the MG models in Section 2. However, in practical engineering problems, the intermediate solutions like those in Case Three are the ones the decision-makers really care about. Therefore, the extreme solutions obtained in Case One and Case Two are of important reference significance for evaluating the performance of the optimization method. With the above simulations, it is clear that I-MOEA/D-M2M does better in searching the extreme solutions than the other two algorithms. Table 2 shows the total cost and emission obtained by the three algorithms with different ω∗ values. It can be seen that the results using I-MOEA/D-M2M have got the minimum total cost in Case One and total emission in Case Two compared with those using the other algorithms. Particularly, in Case Two, the minimum emission solution using I-MOEA/D-M2M incurs a total emission of 19.5747 kg, which is 13.47% and 27.57% better than the SPEA2 and NSGAII optimized solutions, respectively. In Case Three, the cost and emission levels by I-MOEA/D-M2M are lower than those by NSGAII. But the emission by I-MOEA/D-M2M is more than that by SPEA2. According to Case One and Case Two, the former result is still the first choice compared with the latter, since I-MOEA/D-M2M can approximate the true Pareto front better than SPEA2 and reflect the decision-makers’ preference more precisely.Table 2
The total cost and emission obtained by the three algorithms in different cases.
Case One
Case Two
Case Three
Cost ($)
Emission (kg)
Cost ($)
Emission (kg)
Cost ($)
Emission (kg)
Results
I-MOEA/D-M2M
323.8993
45.0830
370.3897
19.5747
339.7723
24.6929
SPEA2
337.6289
45.1581
373.3720
22.6218
343.6473
23.9759
NSGAII
354.7322
48.3517
388.6219
27.0249
359.3728
31.3728Table 3
Parameter settings of the MG system components.
Parameters
MT
FC
BAT
PV
WT
i
1
2
3
4
5
6–8
9–11
K
O
M
i ($/kWh)
7.58 × 10−3
9.35 × 10−3
16.12 × 10−3
17.74 × 10−3
7.3 × 10−3
0
0
A
D
C
C
i ($)
342.82
342.82
1222.2
1222.2
22.84
1259.8
623.41
c
f
i (%)
54.99
54.99
36.73
36.73
32.67
29.34
22.13
a
i
0.752
0.0692
−0.0023
−0.0026
—
—
—
b
i
−0.3095
−0.2904
0.6735
0.5821
—
—
—
c
i
0.4174
0.3752
—
—
—
—
—
d
i
0.1068
0.1090
—
—
—
—
—
α
i (kg/kW)
4.09 × 10−4
2.54 × 10−4
5.33 × 10−4
4.87 × 10−4
—
—
—
β
i (kg/kW)
−5.55 × 10−4
−6.05 × 10−4
−3.55 × 10−4
−4.23 × 10−4
—
—
—
γ
i (kg/kW)
6.49 × 10−4
5.64 × 10−4
3.38 × 10−4
2.46 × 10−4
—
—
—
Ramp rate (kW/min)
0.33
0.33
0.25
0.25
—
—
—Table 4
Parameter settings of the environmental factors and the electricity price.
Wind speed (m/s)
Irradiance (W/m2)
Air temperature (°C)
Electricity sale price ($/kWh)
Electricity purchase price ($/kWh)
5.5
400
21
0.15
0.11Table 5
Parameter settings of the electricity prices.
Case One
Case Two
Case Three
Electricity sale price ($/kWh)
0.11
0.11
0.13
Electricity purchase price ($/kWh)
0.15
0.12
0.15Table 6
The electricity price settings in 24 hours.
Period
Time (h)
Price ($/kW)
Sale
Purchase
Peak period
10:00–15:0018:00–21:00
0.11
0.15
Normal period
07:00–10:0015:00–18:0021:00–23:00
0.06
0.08
Valley period
00:00–07:0023:00-24:00
0.02
0.03


## 5. Conclusion

In this paper, an improved MOEA/D-M2M algorithm is proposed for solving the dynamic MG environmental/economic scheduling problem. In this I-MOEA/D-M2M framework, a constraint-handling method is suggested during the process of sorting and selection. Meanwhile, the self-adaptive decomposition and the individuals’ allocation strategies are employed to enhance the capability of dealing with uncertainties and make the algorithm more efficient. And then I-MOEA/D-M2M is applied to a series of independent static MG environmental/economic scheduling problems. The results are compared with those obtained by SPEA2 and NSGAII in terms of hypervolume values, extreme solutions, and CPU time. It is evident that I-MOEA/D-M2M can achieve a well distributed front every time with much less time, while the other two algorithms cannot guarantee the quality of the solutions set and take more time. The optimization processes in solving the two problems are studied and the shortcomings of the two typical MOEAs are pointed out, which also illustrates that the I-MOEA/D-M2M framework can adapt to the search space changes well. Finally, the feasibility of the I-MOEA/D-M2M is verified by the simulations of a complex 24-hour MG dynamic operation scheduling problem. Further studies are still necessary on more realistic MG models and the effects of other dynamic parameters.



---

*Source: 2167153-2016-05-17.xml*