A Three-Stage Optimal Approach for Power System Economic Dispatch Considering Microgrids

Abstract

The inclusion of microgrids (MGs) in power systems, especially distribution-substation-level MGs, significantly affects power systems because of the large volumes of import and export power flows. Consequently, power dispatch has become complicated, and finding an optimal solution is difficult. In this study, a three-stage optimal power dispatch model is proposed to solve such dispatch problems. In the proposed model, the entire power system is divided into two parts, namely, the main power grid and MGs. The optimal power dispatch problem is resolved on the basis of multi-area concepts. In stage I, the main power system economic dispatch (ED) problem is solved by sensitive factors. In stage II, the optimal power dispatches of the local MGs are addressed via an improved direct search method. In stage III, the incremental linear models for the entire power system can be established on the basis of the solutions of the previous two stages and can be subjected to linear programming to determine the optimal reschedules from the original dispatch solutions. The proposed method is coded using Matlab and tested by utilizing an IEEE 14-bus test system to verify its feasibility and accuracy. Results demonstrated that the proposed approach can be used for the ED of power systems with MGs as virtual power plants.

1. Introduction

Energy efficiency, economic dispatch (ED), and greenhouse gas emissions are major issues associated with utilities. Although electrical energy has been used for over one hundred years since the first electric network was established in 1882 at the Pearl Street Station in New York by Thomson Edison [1], the development of power systems has evolved into generation, transmission, and distribution with a vertical system architecture. In a vertical system structure, a generation system is composed of major large centralized thermal power plants, minor large reservoir hydroelectric power plants, and a small portion of renewable energy generation. Carbon dioxide and other greenhouse gas emissions caused by thermal power plants exacerbate global warming because of the usage of carbon-based fuels in such plants. The electrical energy produced is delivered through transmission and distribution systems until it reaches customers. Electrical energy losses during transmission are considerable. Such losses not only affect the environment, but also reduce the overall system efficiency. Consequently, the traditional large centralized power system architecture is impractical in terms of high efficiency, energy savings, and environmental protection. Conversely, Distributed Energy Resources (DERs) composed of distributed generations (DGs) and energy storage devices [2] are incorporated in local high- or low-voltage distribution networks to satisfy the load demand. As such, DERs can reduce the losses caused by power delivery, and renewable energy-based wind turbine, solar power, or relatively clean units, such as fuel cells and microturbines, can be utilized to reduce greenhouse gas emissions and environmental impact.

DERs interconnected into an existing power grid to form an area grid, which can be operated in grid-tied or autonomous modes, are called MGs [3]. The power grid integrated with information and communication technologies is a smart grid. The development of low-carbon central generation and local energy networks has emerged as a new trend. Europe, the USA, Japan, and Canada have been actively involved in pilot system construction and related research. For instance, the European SmartGrids Technology Platform for future power grids began its work in 2005. This platform aims to formulate and promote a vision for the development of electricity networks in Europe in 2020 and beyond [4]. Future power systems will be composed of low-emission power plants, solar energy, wind turbines, biomass, ocean wave energy, hydroelectric, microturbines, gas turbines, fuel cells, and energy storage devices. As a result, the operation of the future electricity network will be shared between central and distributed generators. Thus, the control of distributed generators can be aggregated to form MGs or virtual power plants (VPPs) to facilitate their integration in physical systems and in markets. Overall, the planning, design, control, and operation of MGs are essential, especially the on-line management or optimal dispatch of various types of DG in an MG under grid-tied and autonomous operation modes. The optimal dispatch of future power grids integrated with MGs as VPPs for safety, economic, and energy loss reduction planning and operation is also relevant in power system research.

ED problems in traditional power systems should be solved in terms of the total minimum operation costs by adjusting the output of each generator to satisfy the load demand and related operational constraints. Since 1920, several power engineering experts and scholars have conducted relevant research. ED problems are based on circuit and optimal theories. The primary task aims to obtain real-time system states, and the power flow algorithm is applied to determine the exact bus voltages and line flows of power systems via traditional AC power flow method [5,6,7,8]. Although solutions generated by such methods are precise, they require time-consuming processes and iterations and exhibit a convergence problem. Newton-Raphson method and Gauss-Seidel method are not applicable for the real-time calculation of power flows. The methods of network sensitivity factors, such as generation shift distribution factor [9], generalized generation shift distribution factor [10], Z-bus distribution factor [11], and Jacobian-based distribution factor (JBDF) [12], have been proposed to improve the computation speed and to eliminate the need for iterations. The application of these approaches overcomes the weakness of Newton and Gauss-based power flow algorithms in ED. Common optimal algorithms in ED and unit commitment are based on Lagrange multiplier-based methods [13,14] and meta-heuristic approaches [15]. Soroudi and Rabiee [16] proposed a new model to solve the multi-area dynamic ED problem by considering the uncertainties in wind power generation, energy prices, and system demands. Optimality condition decomposition is also utilized to verify the proposed approach for the real-time operation of practical power systems. Niu et al. [17] presented an improved quantum-behaved particle swarm optimization (QPSO) to solve the ED problems with valve-point effects and multiple fuel options. This method connects the basic QPSO [18] with a selective probability operator to preserve the fast convergence characteristics and to reduce suffering from premature convergence. Cheng and Zhang [19] combined a dynamic economic dispatch model and an improved particle swarm optimization (PSO) approach to solve the ED problem considering the uncertain nature of wind generation. Li et al. [20] proposed a Quasi-Newton method to update the multiplier, together with an initialization strategy for the multipliers and the approximation matrix. It can be used in dynamic economic dispatch to large-scale integration of renewable generation in power systems. Zaman et al. [21] presented a self-adaptive differential evolution and a real-coded genetic algorithm to solve the dynamic dispatch problem, simulation results demonstrated that the approach is with high solution quality and reliability. Liu et al. [22] introduced an inexact mix-integer two-stage programming model to address low-carbon energy systems in terms of costs, energy security, and carbon dioxide emission. In a previous study [23], fruit fly algorithm is applied to solve an optimal dispatch problem in power systems with a VPP. This algorithm is a nonlinear mixed-integer programming with inter-temporal constraints, and the proposed approach can manage the optimal dispatch strategy among DGs in a VPP and a trading plan with electricity market. Although, these algorithms provide solution quality, convergence speed, and robustness depending on their applications. However, all the related studies mentioned above seldom consider the distribution substation-level MGs in ED problems in modern power systems. This is vital because of the inclusion of MGs in power systems will significantly affect the power dispatch due to the large magnitude of the import and export power flows. Therefore, in this study, an innovative three-stage approach is proposed for the ED of modern power systems comprising MGs as VPPs. The ED model that is derived by the network sensitivity factors in main power grid is with high efficiency and without iteration problem after load changes; besides, the improved direct search method (IDSM), which is with easy and fast characteristics is proposed for optimal dispatch in MGs. Finally, the power exchanges between the main power grid and MGs are modeled by linear incremental model for solving the entire system ED problem by linear programming method. Consequently, this new approach can be utilized for smart grid applications. This paper is divided into five sections. Section 1 introduces the background and objectives of this study. Section 2 describes the traditional Lagrange multiplier-based ED and optimal dispatch in MGs and then derives the proposed three-stage ED model. Section 3 explains the solution procedure of the proposed model. Section 4 discusses the simulation results. Section 5 presents our conclusions.

2. Problem Description and Formula Derivation

MGs play important roles as VPPs in power system operations. Therefore, power systems will be composed of the main power grid and several MGs. The generation dispatch corresponding to load demand will become highly complicated in the near future. In this section, the ED problem is divided into four parts, namely, system structure, ED model of the main power grid, ED model of the MG, and integrated ED model of the power system.

2.1. System Structure

Figure 1 illustrates the system structure of a power system that consists of the main power grid and N distribution substation-level MGs. In terms of transmission systems, the distribution networks with DERs and loads can be regarded as MGs at the primary side of the main transformer as the point of common coupling. According to a multi-area ED model, power systems can be partitioned into the main power grid and a MG in each area. The MG in each area is considered a VPP, which is composed of diesel engine generators, fuel cells, microturbine generators, wind turbine generators, photovoltaic systems, and loads.

2.2. ED of the Main Power Grid

The ED in the main power grid can be described as follows:

$$\min \mathrm{J}={\displaystyle \sum _{i=1}^{NG}{f}_i(P_{Gi})},$$

(1)

subject to:

$${\displaystyle \sum _{i=1}^{NG}{P}_{Gi}-P_D-P_L=0},$$

(2)

$$P_{Gi}^{\min }\le P_{Gi}\le P_{Gi}^{\max } \forall i=1,...,NG,$$

(3)

$$-P_m^{\max }\le P_m\le P_m^{\max } \forall m=1,...,NL.$$

(4)

where ƒ_i and P_Gi represent the cost function and power generation of the ith thermal generator unit, respectively; P_D denotes the total power demand in the main power grid; P_L is the transmission loss; and P_m is the real power flow in the ith transmission line. Using Lagrange multiplier and adding equality constraints, we can rewrite the optimal problem of ED as follows [24]:

$$L={\displaystyle \sum _{i=1}^{NG}{f}_i(P_{Gi})+\lambda (P_D+P_L-{\displaystyle \sum _{i=1}^{NG}{P}_{Gi}})}$$

(5)

On the basis of Equation (5), we let partial L to $P_{Gi}$ be 0, as shown in Equation (6):

$$\frac{\partial L}{\partial P_{Gi}}=0$$

(6)

Using Equation (6), we derive Equation (7):

$$\frac{\partial f_i}{\partial P_{Gi}}=\lambda (1-\frac{\partial P_L}{\partial P_{Gi}}),$$

(7)

we let:

$$\frac{\partial L}{\partial \lambda }=0,$$

(8)

equality constraints, where the generation is equal to the total demands plus losses, can be obtained, as shown in Equation (2). In Equation (7), $\partial P_L/\partial P_{Gi}$ is the incremental transmission loss (ITL), and Equation (7) is rearranged as:

$$(\frac{1}{1-IT{L}_i})(\frac{\partial f_i}{\partial P_{Gi}})=\lambda ,$$

(9)

we rewritten Equation (9) as:

$$P{f}_i(\frac{\partial f_i}{\partial P_{Gi}})=\lambda ,$$

(10)

where Pf_i is equal to $1/(1-IT{L}_i)$ corresponding to the penalty factor. Although ED problems can be solved by artificial intelligence-based optimal algorithms, the Lagrange multiplier-based methods remain the most reliable, most accurate, and fastest in practical system applications. On the basis of Equations (9) and (10), we observe that the solution kernel of ED depends on the rapid calculation of ITL_i, Pf_i, and P_L, and the first task is to derive the mathematical model of P_L to rapidly compute ITL_i and Pf_i. The formula derivation is described in the following sections.

2.3. ED in a Microgrid

The DERs included in MGs are shown in Figure 1. In terms of the multi-area ED concept, the large thermal power plants, transmission systems, and distribution substations in case of the Taipower system regarded as equivalent loads can be formed as the main power grid, and partial distribution stations whose distribution networks with DERs and loads are developed as MGs. Consequently, the primary task of the integrated ED problem in modern power systems is the optimal dispatch in MGs. Thus, the first step is to establish the cost function in an MG.

Similar to large thermal units, the fuel cost of a diesel engine generator is the quadratic polynomial function of its actual power output. However, the fuel cost of a fuel cell generation system is the linear function of its actual power output and efficiency. Therefore, the aggregate cost function of the MG is combined with nonlinear and linear functions. This problem can be solved by direct or indirect search algorithms, where the direct search algorithm called nongradient or zeroth-order method only considers objective function and the indirect search algorithm should consider the partial derivatives of objective function. In this study, an IDSM is proposed to solve the optimal dispatch problem in MGs effectively. In general, a direct search method (DSM) [25] is unsuitable for large-scale systems. The first step of this search approach is to set the incremental cost (IC) value (ΔP). The unit with the minimum IC value is then selected to increase the generation until the IC of each unit is equal and thus satisfies the power balance equation. The weakness of this search method is the start point, that is, it usually starts from 0. As such, the searching time becomes lengthy. The proposed IDSM can improve the limitation of the DSM. The first step is to obtain the ICs of all units, including fuel cell and microturbine with linear costs, as shown in Equations (11) and (12):

$$f_{FCi}={\beta }_i{P}_{FCi}, \forall i=1,...,NFC$$

(11)

$$f_{MTi}={\gamma }_i{P}_{MTi}, \forall i=1,...,NMT$$

(12)

where $f_{FCi}$ denotes the fuel cost of the ith fuel cell; ${\beta }_i$ and $P_{FCi}$ represent the cost coefficient and power generation of the ith fuel cell unit, respectively. $f_{MTi}$, ${\gamma }_i$, and $P_{MTi}$ are the fuel cost, cost coefficient, and power generation of the ith microturbine unit, respectively. NFC and NMT represent the unit number of fuel cell and microturbine in an MG, respectively. IC can be derived as follows:

$$I{C}_{FCi}=\frac{\partial f_{FCi}}{\partial P_{FCi}}={\beta }_i,$$

(13)

$$I{C}_{MTi}=\frac{\partial f_{MTi}}{\partial P_{MTi}}={\gamma }_i.$$

(14)

The ICs of fuel cell and microturbine are constant, as shown in Equations (13) and (14). Therefore, they can be sequentially arranged as follows:

$$I{C}_A=\left[I{C}_1,I{C}_2,...,I{C}_N\right],$$

(15)

where $I{C}_1<I{C}_2<...<I{C}_N$, and N is the summation of ${\mathrm{N}}_{FC}$ and ${\mathrm{N}}_{MT}$ in MG-A. The fuel cost of diesel engine generator can be expressed as follows:

$$f_{DGi}=a_{DGi}+b_{DGi}{P}_{DGi}+c_{DGi}{P}_{DGi}^2 \forall i=1,...,NDG,$$

(16)

where $f_{DGi}$ denotes the fuel cost of the ith diesel engine unit; $a_{DGi}$, $b_{DGi}$, and $c_{DGi}$ are the cost coefficients. NDG represents the unit number of the diesel engine unit in an MG. The IC can be derived as follows:

$$I{C}_{DGi}=\frac{\partial f_{DGi}}{\partial P_{DGi}}=b_{DGi}+2c_{DGi}{P}_{DGi},$$

(17)

substituting the $I{C}_1$ of Equation (15) in Equation (17) results in:

$$P_{DGi}^0=\frac{I{C}_1-b_{DGi}}{2c_{DGi}},$$

(18)

the initial generation, $P_{DGi}^0$, is started, as shown in Equation (18). P_DGi is gradually adjusted according to the proposed IDSM method until this procedure satisfies the power balance equation, as shown in Equation (19):

$${\displaystyle \sum _{i=1}^{NDG}{P}_{DGi}+{\displaystyle \sum _{i=1}^{NFC}{P}_{FCi}}+{\displaystyle \sum _{i=1}^{NMT}{P}_{MTi}}=P_{DA}+P_{LA}-P_{PV}-P_{WT}},$$

(19)

where $P_{DA}$ is the total power demand in MG-A, $P_{DA}$ is the real power loss in MG-A, $P_{PV}$ is the power generation of photovoltaic in MG-A, and $P_{WT}$ is the power generation of wind turbine in MG-A.

2.4. ED Integrated with Main Power Grid and MGs

In this section, the integrated model of the main power grid and MGs is derived on the basis of Section 2.2 and Section 2.3. According to the multi-area dispatch concept, the completed ED problem involves three stages:

• Stage I: ED in main power grid
• Stage II: ED in each MG
• Stage III: integrated ED of main power grid and MGs

The ED in stage III is based on those in stages I and II. The detailed derivation in each stage to establish the linear increment model is described as follows.

2.4.1. Thermal Unit in Main Power Grid

The linear cost model of the thermal power unit in the main power grid is expressed as Equation (20). The coefficient of each unit depends on the power plant type:

$$\Delta f_i=\frac{\partial f_i}{\partial P_{Gi}}\Delta P_{Gi}=(b_i+2c_i{P}_{Gi}^0)\Delta P_{Gi}$$

(20)

2.4.2. NonRenewable Unit in MGs

The linear cost models of diesel engine, fuel cell, and microturbine units in an MG are expressed as Equations (21)–(23), respectively:

$$\Delta f_{DGi}=\frac{\partial f_{DGi}}{\partial P_{DGi}}\Delta P_{DGi}=(b_{DGi}+2c_{D\mathrm{G}i}{P}_{DGi}^0)\Delta P_{DGi}$$

(21)

$$\Delta f_{F\mathrm{C}i}=\frac{\partial f_{FCi}}{\partial P_{FCi}}\Delta P_{FCi}={\beta }_i\Delta P_{FCi}$$

(22)

$$\Delta f_{MTi}=\frac{\partial f_{MTi}}{\partial P_{MTi}}\Delta P_{MTi}={\gamma }_i\Delta P_{MTi}$$

(23)

2.4.3. Linear Programing Model of ED

The linearized dispatch model is established in stage III on the basis of the linear cost model described in Section 2.4.1 and Section 2.4.2 and is solved by the linear programming method. The objective function is expressed as follows:

$$\min J={\displaystyle \sum _{i=1}^{NG}\Delta f_i+{\displaystyle \sum _{i=1}^{NDG}\Delta f_{DGi}}+{\displaystyle \sum _{i=1}^{NFC}\Delta f_{FCi}}+{\displaystyle \sum _{i=1}^{NMT}\Delta f_{MTi}}},$$

(24)

subject to:

$${\displaystyle \sum _{i=1}^{NG}\Delta {\mathrm{P}}_{\mathrm{G}i}-\Delta P_L-P_{eA}}=0,$$

(25)

$${\displaystyle \sum _{i=1}^{NDG}\Delta {\mathrm{P}}_{DGi}}+{\displaystyle \sum _{i=1}^{NFC}\Delta {\mathrm{P}}_{FCi}}+{\displaystyle \sum _{i=1}^{NMT}\Delta {\mathrm{P}}_{MTi}}-\Delta P_{LA}+\Delta P_{eA}-\Delta P_{DA}=0,$$

(26)

$$\Delta P_{Gi}^{\min }\le \Delta P_{Gi}\le \Delta P_{Gi}^{\max },$$

(27)

$$\Delta P_{DGi}^{\min }\le \Delta P_{DGi}\le \Delta P_{DGi}^{\max },$$

(28)

$$\Delta P_{FCi}^{\min }\le \Delta P_{FCi}\le \Delta P_{FCi}^{\max },$$

(29)

$$\Delta P_{MTi}^{\min }\le \Delta P_{MTi}\le \Delta P_{MTi}^{\max },$$

(30)

where $P_{eA}$ is the power interchange between the main power grid and MG-A, and positive $P_{eA}$ means that power is fed into MG-A. These equations can be expressed as a standard linear model as follows:

$$\min J{=\mathrm{C}}^{\mathrm{T}}x,$$

(31)

subject to:

$$\{\begin{array}{c}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1N}{x}_N=b_1\\ ⋮⋮⋮⋮ \\ a_{N1}{x}_1+a_{N2}{x}_2+\cdots +a_{NN}{x}_N=b_N\end{array},$$

(32)

$$\{\begin{array}{c}{x}_1^{\min }\le x_1\le x_1^{\max }\\ ⋮ \\ x_N^{\min }\le x_N\le x_N^{\max }\end{array},$$

(33)

each extreme point is gradually searched in accordance with the linear programming method. The optimal power generation of each unit is derived, as shown in Equations (34)–(37). The optimal ED solution of power systems considering MGs as VPPs can be obtained:

$${\mathrm{P}}_{\mathrm{G}i}=P_{Gi}^0+\Delta {\mathrm{P}}_{\mathrm{G}i},$$

(34)

$${\mathrm{P}}_{DGi}=P_{DGi}^0+\Delta {\mathrm{P}}_{DGi},$$

(35)

$${\mathrm{P}}_{FCi}=P_{FCi}^0+\Delta {\mathrm{P}}_{FCi},$$

(36)

$${\mathrm{P}}_{MTi}=P_{MTi}^0+\Delta {\mathrm{P}}_{MTi}.$$

(37)

3. Solution Procedure

In this section, the detailed solution algorithm and model of the main power grid and MGs are explained. The proposed approach can be applied in multi-area complex power systems. The solution procedure is described as follows.

3.1. Stage I

The JBDF [12]-based Lagrange multiplier method is proposed to solve the ED in the main power grid. A power flow program should be executed to solve the ED problem in each iteration by the traditional Lagrange multiplier method. However, this procedure is time consuming and may result in divergence condition. The key point is to obtain the ITL and penalty factor (Pf) efficiently after the load demand and generation are changed. The JBDF-based formula is derived to obtain ITL and Pf in the ED procedure instead of running a power flow program and to overcome the weakness of ITL and Pf derivation after load demand changes.

3.1.1. Models of Transmission Line Loss and Line Flow

P_m and Q_m represent the line m flows of real and reactive powers, and the transmission line loss P_L is expressed as follows:

$${\mathrm{P}}_{\mathrm{L}}={\displaystyle \sum _{m=1}^{NL}{R}_m\cdot (P_m^2+Q_m^2)} \forall m=1,...,NL,$$

(38)

where R_m denotes the resistance of line m, and NL is the number of transmission lines. The real and reactive power line flows can be expressed as [12]:

$${\mathrm{P}}_m=P_m^0+{\displaystyle \sum _{i=1}^{NG}{F}_p(m,i)}\cdot \Delta {\mathrm{P}}_{Gi}+{\displaystyle \sum _{i=1}^{NG}{K}_p(m,i)}\cdot \Delta {\mathrm{Q}}_{Gi}$$

(39)

$${\mathrm{Q}}_m={\mathrm{Q}}_m^0+{\displaystyle \sum _{i=1}^{NG}{F}_q(m,i)}\cdot \Delta {\mathrm{P}}_{Gi}+{\displaystyle \sum _{i=1}^{NG}{K}_q(m,i)}\cdot \Delta {\mathrm{Q}}_{Gi},$$

(40)

in the ED procedure, the real power is generally only adjustable using generators corresponding to load demand changes, and the reactive power is kept constant. Consequently, Equations (39) and (40) are simplified as Equations (41) and (42), respectively. $\Delta {\mathrm{P}}_{Gi}$ is the change in the real power of the ith generator from the base case:

$${\mathrm{P}}_m=P_m^0+{\displaystyle \sum _{i=1}^{NG}{F}_p(m,i)}\cdot \Delta {\mathrm{P}}_{\mathrm{G}i}$$

(41)

$${\mathrm{Q}}_m={\mathrm{Q}}_m^0+{\displaystyle \sum _{i=1}^{NG}{F}_q(m,i)}\cdot \Delta {\mathrm{P}}_{\mathrm{G}i}$$

(42)

$$\Delta {\mathrm{P}}_{Gi}={\mathrm{P}}_{Gi}-P_{Gi}^0$$

(43)

where F_p and F_q are the distribution factors of JBDF, which are expressed as Equations (44) and (45), respectively:

$$F_p(m,i)=\frac{\partial \left|{\mathrm{V}}_p\right|}{\partial P_i}\cdot \frac{\partial P_m}{\partial \left|{\mathrm{V}}_p\right|}+\frac{\partial \left|{\mathrm{V}}_q\right|}{\partial P_i}\cdot \frac{\partial P_m}{\partial \left|{\mathrm{V}}_q\right|}+\frac{\partial {\delta }_p}{\partial P_i}\cdot \frac{\partial P_m}{\partial {\delta }_p}+\frac{\partial {\delta }_q}{\partial P_i}\cdot \frac{\partial P_m}{\partial {\delta }_q}$$

(44)

$$F_q(m,i)=\frac{\partial \left|{\mathrm{V}}_p\right|}{\partial P_i}\cdot \frac{\partial Q_m}{\partial \left|{\mathrm{V}}_p\right|}+\frac{\partial \left|{\mathrm{V}}_q\right|}{\partial P_i}\cdot \frac{\partial Q_m}{\partial \left|{\mathrm{V}}_q\right|}+\frac{\partial {\delta }_p}{\partial P_i}\cdot \frac{\partial Q_m}{\partial {\delta }_p}+\frac{\partial {\delta }_q}{\partial P_i}\cdot \frac{\partial Q_m}{\partial {\delta }_q}$$

(45)

where $\partial \left|{\mathrm{V}}_p\right|/\partial P_i$, $\partial \left|{\mathrm{V}}_q\right|/\partial P_i$, $\partial {\delta }_p/\partial P_i$ and $\partial {\delta }_q/\partial P_i$ represent the magnitude and phase angle of the bus voltage sensitive to real power injections. These terms can be calculated in the inverse Jacobian matrix of the base case power flow solution.

3.1.2. Formula of ITL and Penalty Factor

The partial derivatives of Equation (38) and the JBDF-based ITL can be derived as Equation (46). Pf can be calculated by $1/(1-IT{L}_i)$:

$$IT{L}_i=\frac{\partial P_L}{\partial P_{Gi}}=2{\displaystyle \sum _{i=1}^{NL}\left[P_m\cdot F_p(m,i)+Q_m\cdot F_q(m,i)\right]},$$

(46)

3.1.3. Model of JBDF-Based ED

Based on Equations (2) and (7), the model of JBDF-based ED conditions is rewritten as follows:

$$F_i=\frac{\partial L}{\partial P_{Gi}}=b_i+2c_i\cdot P_{Gi}+\lambda (1-IT{L}_i)=0,$$

(47)

$$F_{NG+1}=\frac{\partial L}{\partial \lambda }={\displaystyle \sum _{i=1}^{NG}{P}_{Gi}-P_D-P_L}=0,$$

(48)

extending Equation (48) to all units results in (NG + 1) simultaneous nonlinear equations in a matrix form as follows:

$$\left[\begin{array}{c}0-F_1\\ ⋮\\ 0-F_{NG+1}\end{array}\right]=\left[\begin{array}{cc}{J}_1& J_2\\ J_3& J_4\end{array}\right]\left[\begin{array}{c}\Delta P_{G1}\\ ⋮\\ \Delta P_{GNG}\\ \Delta \lambda \end{array}\right],$$

(49)

Equation (49) can be solved by the Newton-based algorithm, where J₁, J₂, J₃, and J₄ are the sub-matrices of Jacobian matrix, as explained in the following section.

The diagonal and off-diagonal elements of J₁ are expressed as follows:

$$\frac{\partial F_i}{\partial P_{Gi}}=2c_i-\lambda (\frac{\partial IT{L}_i}{\partial P_{Gi}})=0,$$

(50)

$$\frac{\partial F_i}{\partial P_{Gj}}=-\lambda (\frac{\partial IT{L}_i}{\partial P_{Gj}}), j\ne i,$$

(51)

the partial derivatives of ITL are shown in Equations (52) and (53):

$$\frac{\partial IT{L}_i}{\partial P_{Gi}}=2R_m{\displaystyle \sum _{i=1}^{NL}\left[F_p^2(m,i)+F_q^2(m,i)\right]},$$

(52)

$$\frac{\partial IT{L}_i}{\partial P_{Gj}}=2R_m\left[F_p(m,i)\cdot F_p(m,j)+F_q(m,i)\cdot F_q(m,j)\right], j\ne i,$$

(53)

the elements of J₂ are:

$$\frac{\partial F_i}{\partial \lambda }=1-IT{L}_i,$$

(54)

the elements of J₃ are:

$$\frac{\partial F_{NG+1}}{\partial P_{Gi}}=1-IT{L}_i,$$

(55)

the elements of J₄ are:

$$\frac{\partial F_{NG+1}}{\partial \lambda }=0,$$

(56)

the model of the JBDF-based ED in the main power grid can be derived according to these equations. The solution flow chart is shown in Figure 2. The solution is obtained when the convergence conditions, $\Delta P$ and $\Delta \lambda$, are less than a specified accuracy, i.e., $\epsilon$. On the contrary, it is divergence, the iterations will reach to maximum iteration number and then stop the iteration procedure. It is worth noting that the Newton-based algorithm is sensitive to the initial guess value. In general, the physical meaning of $\lambda$ is the system incremental cost of delivered power; therefore, the initial value could be guessed in a reasonable range according to the average cost per megawatt for avoiding the unreasonable solution or divergence.

3.2. Stage II

The DERs are divided into nonrenewable and renewable units. The cost functions of nonrenewable units and the prediction of the available power of renewable units are essential for optimal dispatch in an MG. The characteristics of the DERs are described as follows.

For the nonrenewable units, the fuel cost of a diesel engine generator is expressed as Equation (16), and the total fuel cost of diesel engine units in MG-A are shown in Equation (57). Molten carbonate fuel cells and solid oxide fuel cells are used in this study. The fuel cost of the fuel cell unit is shown in Equation (11), and the total fuel cost of fuel cell units in MG-A is shown in Equation (58). Similar to the cost function of the fuel cell unit, the fuel cost of the fuel cell unit is shown in Equation (12), and the total fuel cost of fuel cell units in MG-A is shown in Equation (59):

$$F_{DG(A)}={\displaystyle \sum _{i=1}^{NDG}{a}_{DGi}+b_{DGi}{P}_{DGi}+c_{DGi}{P}_{DGi}^2}$$

(57)

$$F_{FC(A)}={\displaystyle \sum _{i=1}^{NFC}{\beta }_i{P}_{FCi}}$$

(58)

$$F_{MT(A)}={\displaystyle \sum _{i=1}^{NMT}{\gamma }_i{P}_{MTi}}$$

(59)

The power output of the wind turbine unit is expressed as Equation (60), where C_p is the power coefficient, $\nu$ is the air density, $V_w$ is the wind velocity, and A is the rotor swept area. The total power generation of wind power units in MG-A is shown in Equation (61). The power output of the photovoltaic module can be expressed as Equation (62), with the measured solar irradiation (W/m²) and ambient temperature (°C), where P_stc is the maximum power at the standard test condition, G_ing is the incident irradiance, G_stc is the irradiance at the standard test condition (1000 W/m²), k is the temperature coefficient of power, T_c is the cell temperature, and T_r is the reference temperature. The total power generation of photovoltaic units in MG-A is shown in Equation (63). Therefore, the renewable units are regarded as negative power loads in optimal dispatch in MGs:

$$P_{WT}=\frac{1}{2}{\mathrm{C}}_p\nu {\mathrm{V}}_w^3\mathrm{A}$$

(60)

$$P_{WT(A)}=\frac{1}{2}{\displaystyle \sum _{i=1}^{NWT}{\mathrm{C}}_{pi}\nu {\mathrm{V}}_w^3\mathrm{A}}$$

(61)

$$P_{PV}={\mathrm{P}}_{STC}\frac{G_{ING}}{G_{STC}}\left[1+K_i(T_c-T_r)\right]$$

(62)

$$P_{PV(A)}={\displaystyle \sum _{i=1}^{NPV}{\mathrm{P}}_{PV}}$$

(63)

On the basis of the cost function and power generation of each unit, we can express the objective function for the optimal dispatch problem in MG-A as follows:

$$\min F_{T(A)}={\displaystyle \sum _{i=1}^{NDG}{f}_{DGi}+{\displaystyle \sum _{i=1}^{NFC}{f}_{FCi}}+{\displaystyle \sum _{i=1}^{NMT}{f}_{MTi}}},$$

(64)

subject to:

$${\displaystyle \sum _{i=1}^{NDG}{P}_{DGi}}+{\displaystyle \sum _{i=1}^{NFC}{\mathrm{P}}_{FCi}}+{\displaystyle \sum _{i=1}^{NMT}{\mathrm{P}}_{MTi}}=\Delta P_{D(A)}+\Delta P_{L(A)}-\Delta P_{WT(A)}-\Delta P_{PV(A)},$$

(65)

$$P_{DGi}^{\min }\le P_{DGi}\le P_{DGi}^{\max },$$

(66)

$$P_{FCi}^{\min }\le P_{FCi}\le P_{FCi}^{\max },$$

(67)

$$P_{MTi}^{\min }\le P_{MTi}\le P_{MTi}^{\max }.$$

(68)

Finally, the IDSM is proposed to solve the optimal dispatch in an MG, and the solution flow chart is shown in Figure 3.

3.3. Stage III

In stage III, the linear incremental model of the entire power system is established on the basis of the optimal dispatch solutions of the two previous stages. Considering the power interchange of MG-A, P_eA, into the main power grid, we can obtain the linearly optimized model of the main power grid and MG-A. The linearly optimized model in stage III can be derived by increasing N number of MGs in the power system, as shown in Equation (69):

$$\min J={\displaystyle \sum _{i=1}^{NG}\Delta f_i+{\displaystyle \sum _{j=A}^N{\displaystyle \sum _{i=1}^{NDG}\Delta f_{DG(j,i)}}}+{\displaystyle \sum _{j=A}^N{\displaystyle \sum _{i=1}^{NFC}\Delta f_{FC(j,i)}}}}+{\displaystyle \sum _{j=A}^N{\displaystyle \sum _{i=1}^{NMT}\Delta f_{MT(j,i)}}}$$

(69)

According to the linear model and satisfying the related constraints, the optimal dispatch generations are solved by linear programming. The optimal solution of the entire power system with MGs as VPPs can also be found.

4. Numerical Results and Discussion

An effective ED program is coded in the Matlab environment according to the proposed three-stage model and solution algorithms. The proposed approach is first used to evaluate the ED problem of the original IEEE 14-bus test system and to verify the accuracy and effectiveness of the proposed approach, and the related parameters for power flow and ED are obtained from previous studies [25,26]. In this scenario, the MG is assumed to be operated autonomously and therefore not considered in ED in the main power grid for the comparison of the traditional Lagrange multiplier method and the proposed approach based on the known results. The numerical results (

Table 1. Comparison of the simulation results of the traditional Lagrange multiplier method and the proposed approach for the ED of the original IEEE 14-bus test system.

Bus No.	Traditional Lagrange Multiplier (A)	Proposed Approach (B)	Error %
	Generation (MW)	Generation (MW)	|1 − B/A| × 100
1	160.450	160.460	0.0062
2	68.803	68.810	0.0102
6	39.567	39.550	0.0430
Total cost ($/h)	1137.7	1137.7	0.0000
Real power loss (MW)	9.82	9.82	0.0000
λ	405.4473	405.4590	0.0029
CPU execution time (s)	0.188	0.047	-

) demonstrate that the developed ED program is valid. The ED problem of the main power grid considering MGs is then solved by using the three-stage ED program.

The IEEE 14-bus test system and an MG (Figure 4 and Figure 5) are modified as a sample system to verify the feasibility and accuracy of the proposed approach. The equivalent lumped daily load demand curves of the main power grid are illustrated in Figure 6. The equivalent daily load demand curves that correspond to different load types at each bus are shown in Figure 7. The purchased/sold electricity cost of MG-A are depicted in Figure 8. This cost comprises the peak-load (0.2392 $USD/kW·h), half-peak-load (0.0894 $USD/kW·h), and off-peak-load (0.0447 $USD/kW·h) pricing from the main power grid during the summer season. On the basis of the three-stage ED algorithm and constraints, we show the numerical results of the optimal dispatch of each generator in the main power grid in Figure 9, and the optimal power interchange between the main power grid and MG-A is shown in Figure 10. The numerical results of the generation cost per kWh and the total generation cost of the sample system are illustrated in Figure 11a,b, respectively. The numerical results also demonstrate that the bus voltage, current flow, and power generation limits of each unit are consistent with the constraints. The executing time of the proposed approach is around 25% of the traditional Lagrange multiplier method in the IEEE 14-bus test system with a MG; furthermore, there is no divergence problem for power flow solution of the proposed approach because of no iteration after load demand changes. Therefore, the performance of both computation time and reliability of the proposed approach is superior to those of the traditional Lagrange multiplier method. Hence, the proposed approach exhibits high efficiency and rapidly solves the ED of power systems considering MGs as VPPs.

5. Conclusions

In this study, effective and fast three-stage ED model and approach are proposed to solve the problem of modern power systems considering MGs as VPPs. Using the traditional Lagrange multiplier method, we derive the JBDF-based ED algorithm for stage I dispatch according to the network sensitivity factors in the main power grid, and the IDSM for stage II is developed for the optimal dispatch in an MG. The linear incremental model and solution algorithm for stage III are obtained and integrated with stages I and II. For the algorithms proposed in this paper, the power systems are divided into the main power grid and area MGs. The ED problem is solved on the basis of the concept of multi-area ED by using the proposed three-stage approach. Besides, the traditional methods of ED required iterations in each solution procedure; nevertheless the proposed approach only need to execute one iteration procedure by Newton algorithm as base case and then without any iteration after load demand changes. It is absolutely verified that not only the computing speed but also the convergence are superior to other methods that need iterations after load demand changes in real-time dispatch. The numerical results of the test system demonstrate that the proposed approach exhibit high accuracy and efficiency. Therefore, this approach is helpful for optimal dispatch because of the large amounts of DERs interconnected into distribution systems. The proposed approach can also be utilized to solve the area power dispatch problems because of the distribution systems formed as MGs, which are regarded as VPPs for future smart grid applications.