Distributed Optimal Energy Management for Energy ... - IEEE Xplore

IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS, VOL. 13, NO. 6, DECEMBER 2017

3081

Distributed Optimal Energy Management for Energy Internet Huaguang Zhang

, Fellow, IEEE, Yushuai Li, David Wenzhong Gao and Jianguo Zhou

Abstract—In this paper, a novel energy management framework for energy Internet with many energy bodies is presented, which features multicoupling of different energy forms, diversified energy roles, and peer-to-peer energy supply/demand, etc. The energy body as an integrated energy unit, which may have various functionalities and play multiple roles at the same time, is formulated for the system model development. Forecasting errors, confidence intervals, and penalty factor are also taken into account to model renewable energy resources to provide tradeoff between optimality and possibility. Furthermore, a novel distributed-consensus alternating direction method of multipliers (ADMM) algorithm, which contains a dynamic average consensus algorithm and distributed ADMM algorithm, is presented to solve the optimal energy management problem of energy Internet. The proposed algorithm can effectively handle the problems of power-heat-gascoupling, global constraint limits, and nonlinear objective function. With this effort, not only the optimal energy market clearing price but also the optimal energy outputs/demands can be obtained through only local communication and computation. Simulation results are presented to illustrate the effectiveness of the proposed distributed algorithm. Index Terms—Distributed algorithm, energy Internet, energy management, multiagent.

I. INTRODUCTION HE concept of “Energy Internet” is recently presented to meet the challenges of developing sustainable and environmentally friendly energy resources, hybrid energy utilization model, flexible energy management, and secure system control. Several literatures about energy Internet (EI), mainly focusing on system architecture [1], energy router [2], and voltage control [3], etc., have been proposed in recent years. Likewise, energy management problem (EMP), as one of fundamental issues in

T

Manuscript received April 9, 2017; accepted May 28, 2017. Date of publication June 9, 2017; date of current version December 1, 2017. This work was supported in part by the National Natural Science Foundation of China (NSFC) Key Program under Grant 61573094, and in part by the NSFC under Grant 61621004 Grant 61428301. Paper no. TII-17-0714. (Corresponding author: Huaguang Zhang.) H. Zhang, Y. Li, and J. Zhou are with the School of Information Science and Engineering, Northeastern University, Shenyang 110819, China (e-mail: [email protected]; [email protected]; jianguozhou.neu@ gmail.com). D. W. Gao is with the Department of Electrical and Computer Engineering, University of Denver, Denver, CO 80208 USA (e-mail: Wenzhong. [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TII.2017.2714199

, Senior Member, IEEE,

power system, is encouraged to be revisited to achieve the envisioned EI conception. EMP is typically formulated as an optimization problem. A number of approaches have been proposed to solve the above problem and can be classified into two categories, i.e., centralized approaches and distributed approaches. With regard to centralized approaches [4]–[7], they have the advantage in the aspect of obtaining the optimal solution. However, as the physical power system tends to become distributed, centralized approaches have to confront many challenges. Specifically, the centralized approaches, relying on a central controller, require high bandwidth communication capabilities to act on the system-wide gathered information, resulting in high susceptibility to single-point failures and modeling errors. Moreover, either the physical or communication network of the future power grid is likely to have a variable topology, which may undermine the efficacy of the centralized approach. In addition, the distributed energy resources owners are generally unwilling to disclose their own information to the external centralized controller to protect their privacy [8]. On the contrary, the distributed approaches, with better robustness [9], faster computation, and less communication [10], are promising options for the future power system. A number of literatures, concerned about solving the EMP in a distributed manner, have been proposed for smart grid or microgrid recently, which can be classified into two main categories. The one, which does not consider the demand response, is also called economic dispatch problem (EDP), whose object is minimizing the total operating cost while meeting some equality and inequality constraints. From this prospective, the λ-consensus algorithm has been first proposed in [11], which can solve the traditional EDP in a fully distributed manner. The work in [12] has presented a new distributed algorithm that does not rely on the projection scheme and also can achieve the underlying power flow control. In order to address nonquadratic cost function, a distributed bisection method and a distributed algorithm based on projected gradient method have been presented in [13] and [14], respectively. Moreover, the ramping rate limits [15], transmission losses [16], finite-time [17], nonconvex conditions [18], and frequency recovery [19] have been further taken into consideration and studied in the EDP. The other one aims at maximizing the social welfare through the coordination of the suppliers’ generations and customers’ demands [20]–[22]. Since the demand response is taken into consideration, the power system operation becomes more flexible and economic.

1551-3203 © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.

3082


Different from the studies above for the traditional system or smart grid, the EMP for EI is a relatively new but difficult problem, which will bring many challenges. On one hand, EI is a hybrid energy network with many new characteristics, which can be summarized in three aspects. First, EI integrates various kinds of energy networks. As the main energy medium, it is urgent for power, heat, and gas flows to be co-planned. Second, each unit in EI can be not only an energy consumer but also an energy supplier. The form of energy generation and consumption are tending to be diversified. Third, distributed energy suppliers or consumers within EI have peer-to-peer relationship, which is different from the traditional, vertical supply–demand relationship. The detailed features and requirements of EI will be discussed in Section II-A. Therefore, how to establish a system model, which can describe aforementioned features of EI in a better way, is one of the major challenges. On another hand, EI should meet the plug-and-play nature and topology variabilities, etc. Thus, it is desirable to develop a distributed method to solve the EMP of EI. More importantly, different from smart grid, there exists strong coupling among power, heat, and gas in the process of energy-generation, energy-conversion, and energyconsumption for EI. Meanwhile, there also exist many global constraint limits and nonlinear objective function. However, the existing distributed methods discussed above cannot deal with this kind of distributed nonlinear coupling optimization problem. Therefore, how to establish a fully distributed method to solve the EMP of EI is another major challenge. To address above challenges, this paper focuses on the static EMP of the future EI considering multienergy-networks and intermittency of renewable energy resources, where a novel distributed-consensus-ADMM algorithm is proposed to solve this problem in a fully distributed fashion. The major contributions of this paper are as follows: 1) The energy management framework is developed for the future EI featuring multiple couplings of different energy forms, diversified energy roles, and peer-to-peer energy supply/demand, etc. Subsequently, the energy body, seen as both energy supplier and customer, is presented for system model development. 2) Forecasting errors, confidence intervals, and penalty factor are also considered to model renewable energy resource. With this effort, the optimality and generation possibility of the renewable energy can be traded off by designing the penalty on curtailment of renewable generation. 3) By dividing the global computation process into the distributed participants, our implementation fashion can make each participant locally calculate its optimal operation. It can result in enhanced reliability, flexibility, and scalability, etc., and is more suitable for the future EI to integrate distributed multienergy resources. 4) The proposed algorithm can effectively overcome the weakness of conventional distributed ADMM algorithm, which assumes that part of global information is known via designing a dynamic average consensus algorithm for estimating the global information. Moreover, the optimality and convergence analysis are strictly proved.

5) We get the explicit iterative form of coupled variables by transforming and decoupling the max-min problem, and make the coupled inequality constraints be solved by locally calculating submin problem. With this effort, the strongly power-heat-gas-coupling problem can be effectively solved by implementing the proposed algorithm. The remainder of the paper is organized as follows. Section II identifies the main features and requirements, and builds the model of the future EMP of EI. Section III first introduces some basic graph theory and then the distributed-consensus-ADMM algorithm is presented. In Section IV, several case studies are presented to show the effectiveness of the proposed algorithm. Section V concludes the paper. II. ENERGY MANAGEMENT FRAMEWORK OF EI A. Structure and Features Analysis The anticipated structure of EI, shown in Fig. 1, consists of various kinds of distributed energy suppliers and/or customers which are referred as energy bodies (EBs). Fed by multienergycarriers, the energy resources of each EB can be divided into four classes, i.e., power-only devices which contain an equivalent distributed renewable generator (DRG), distributed fuel generator (DFG), and distributed power storage device (DPSD), heat-only devices which contain an equivalent distributed renewable heating device (DRHD), distributed fuel heating device (DFHD), and distributed heat storage device (DHSD), an equivalent distributed combined heat and power (DCHP) device and an equivalent distributed gas producer (DGP). Meanwhile, the energy loads (ELs) of each EB can be divided into three classes, i.e., power load, heat load, and gas load, in which every type of the load contains an equivalent controllable load and a must-run load. Note that each EB contains at least one type of energy resource or load which are connected to the corresponding energy bus. More importantly, the EB can be seen as small as a house or as large as a town based on the requirements of different scenarios. In order to achieve the anticipated EI, the major features and requirements of the proposed energy management framework of EI are summarized as follows: 1) As for EI, one of the features is to integrate multienergynetworks, such as electricity, transportation, nature gas, oil, and heat networks, etc. Since power, gas, and heat are the major energy that can be flexibly transferred and directly consumed, it is urgent that they can be comanaged and be the primary energy medium of energy-generation, energy-transmission, and energy-consumption, which brings many advantages. For example, since there are many kinds of alternative energy resources in EI, consumers can diversify their energy supply mix and determine the amount of energy utilization in an overall optimal manner. This will create significant incentives to increase the consumer participation and put full use of the function of “virtual storage” caused by demand dispatch to flatten out the variability of energy generation. Moreover, The co-planning benefits improvement of the energy efficiency and economy and decreasing carbon emission, etc. 2) In the anticipated EI, each EB may be equipped with its own energy generation, conversion, and storage devices and

ZHANG et al.: DISTRIBUTED OPTIMAL ENERGY MANAGEMENT FOR ENERGY INTERNET

Fig. 1.

3083

Structure, features, and benefits of EI.

schedulable loads, so it plays a role of not only energy (i.e., power, heat, and gas) supplier but also energy consumer. For example, a EB may be a power supplier and a heat consumer at the same time. Meanwhile, motivated by the Internet, each EB is also both the information supplier and information consumer, and they only exchange information with its neighbors through local network. Moreover, all EBs are cooperative to achieve the maximization of social welfare while meeting a set of local and global constraints. Then, each EB can determine its corresponding energy-production or/and energy-consumption. It is worth noting that each EB may compensate its self-energy-demand, and even make additional profit by selling the part of overabundance energy to other EB. In this case, it can be seen as an energy-supplier. Otherwise, it is seen as an energy consumer. In such an energy management framework, EB can somehow weaken dependence on utility grid and promote network flattening; meanwhile, there will be more interactions and flexibility among EBs. More importantly, driven by profit, it is able to create sufficient incentives for each EB to make self-installation of renewable energy resources, resulting in increased penetration of renewable energy resources. 3) To integrate and coordinate the large number EBs with the plug-and-play feature, both the physical and communication topology of EI are subject to variability and uncertainty. Therefore, it is imperative to develop a distributed control strategy, with decision made by the EBs themselves rather than a central controller. This also means that each EB should have the function of local information processing and computing. The distributed implementation fashion benefits for enhancement of system reliability, flexibility for possible expansions of distributed energy resources, and scalability in terms of the computation and communications burden, etc. Moreover, in the proposed energy management framework, it achieves the

vertical energy supply mode to peer-to-peer mode. To be specific, EBs have peer-to-peer relationship. The renewable and conventional fuel devices are in peer-to-peer mode to accomplish their own function. In addition, the diversified energy resources, including electric power, gas, and heat, are also in the peer-to-peer mode to collectively achieve their optimization, management, and control. From the discussion above, the major differences of EI and the smart grid on EMP are summarized as follows. First, the smart grid mainly focuses on the optimization and management of electrical power on power system aspect, in which the electrical power is the major energy medium in the process of energygeneration, energy-conversion, and energy-consumption. However, EI focuses on integrating various kinds of energy resources to fulfill the coplanning of different energy networks, which contains multitype coupled energy mediums such as power, heat, and gas. In addition, each energy network in EI is in peer-topeer position to fulfil its functionalities. Second, each unit of smart grid can be only seen as an electrical power supplier or consumer in one scheduling horizon. On the contrast, the EB in EI, regarded as both energy supplier and consumer, can play multiple roles at the same time. Last but not least, both smart grid and EI are eager to develop distributed methods for the EMP. However, the strongly power-heat-gas coupling situation of EI, which exists in the objective function and constraint limits, vastly increases the solving difficulty and is more complicated than electricity. B. System Model 1) DRG and DRHD Models: For DRGs, the solar radiation and wind are the major energy sources with zero fuel cost. Thus, the direct cost of renewable generators is often the

3084


operating cost which is proportional to the produced power. It is worth noting that renewable generators cannot be regarded as dispatchable units due to the intermittency and randomness features. To take them into consideration in our EMP, the forecast technologies are used in this paper. The mean forecasting value, seen as the reference during each scheduling horizon, can be calculated from day-ahead forecast curve given by t+T e,† = prij,e,†t¯ dt¯ /(T ) (1) prij,t t

prij,e,†t¯

where is the forecasting power generation of jth DRG of ith EB at time t¯; T is the scheduling horizon. Note that real value may not be the forecasting value because of the forecast error. In this paper, it is assumed that the forecast error, denoted by e,† , obeys Gaussian distribution whose feasibility analysis Δprij,t has been discussed in [23]. Then, the corresponding probability density function can be modeled as 2 1 e,† r e,† re 2 =√ exp − Δp / 2δ . (2) f Δprij,t ij,t ij,t re 2πδij,t Furthermore, let the confidence level be 100(1 − ℘ij )%, then r e,†,down r e,†,up we can get the confidence intervals [pij,t , pij,t ] by using the method proposed in [24]. Therein, ℘ij denotes the corresponding significance level. In addition, we assume that the power generation of each DRG can be fully consumed by loads, then we can let e,† e,† e e e e = prij,t + Δprij,t → prij,t ≤ prij,t ≤ prij,t prij,t

(3)

r e,†,up e,†,down e,† e,† e e where prij,t = prij,t + prij,t , prij,t = prij,t + pij,t , and re pij,t is the power output of jth DRG of ith EB. Based on the above discussion, the cost function of DRG is modeled as

e e re prij,t − prij,t re (4) C pij,t = bij pij,t + εij exp γij r e e pij,t − prij,t

where bij > 0 and εij > 0 are the cost coefficients, γij < 0 is the penalty factor. The first item of the above equation denotes the direct operating cost while the second one denotes the penalty on curtailment of renewable generation. The purpose and significance of (4) are discussed in Remark 1. In addition, the modeling process and cost function of DRHD are designed e e and C(hrij,t ) denote the accordingly as for DRG. And let hrij,t corresponding heat generation and cost function, respectively. Remark 1: Compared with DFGs, the operating cost of each DRG is very small. As a result, from the viewpoint of optimization, it may tend to run at its corresponding upper e e . However, the higher value the prij,t is calbound, i.e., prij,t culated, the lower possibility the practical power generation has to meet this requirement. On the contrary, if the renewre , there is 100(1 − ℘ij )% posable generator runs at point pij,t sibility to believe that its practical capability can meet this power output, but that may decrease the optimality. Thus, the exponential penalty on curtailment of renewable generation is designed to provide tradeoff between the optimality and

Fig. 2.

Cost function model of DRG.

possibility. In this way, the DFG and DRG can be in the same form. Based on the aforementioned discussions, the concept of the cost function model of DRGs is described in Fig. 2. Further, a simple example is supplied to clearly illustrate the idea of (4). We consider a system composed of a DRG with direct operating cost CD (pr e ) = 0.1095pr e [25], a DFG with cost function C(pf e ) = 0.04(pf e )2 + 25pf e + 99 + 50 exp(0.01pf e ) [30(p.u.) ≤ pf e ≤ 150 (p.u.)] [13], and a mustrun load L = 220 (p.u.). The confidence level and the forecasting value are, respectively, set as 90% and 93.90(p.u.), s.t., 84.3 (p.u.) ≤ pr e ≤ 103.5 (p.u.). Before considering the penalty on curtailment of renewable generation, the goal is to minimize Goal = CD (pr e ) + C(pf e ) while meeting the supply–demand balance constraint, i.e., pr e + pf e = L, and the inequality constraints mentioned above. Since the direct operating cost is proportional to pr e , the incremental cost of the DRG is a constant, i.e., 0.1095, which is very small compared to the one of DFG. The calculation results show that pr e = 103.5 (p.u.) and pf e = 116.5 (p.u.). In this version, although it can obtain better optimality or less cost, its corresponding possibility is very small. In other words, only when the real value is over 103.5 (p.u.), the practical capability can meet this requirement. Note that the lower the value of pr e is, the higher the possibility becomes. Thus, to guarantee possibility, we can let pr e = 84.3 (p.u.), resulting in decreased optimization. Different from the two cases, the exponentialpenalty on curtailment of renewable generation, 103.5 − pr e set as 748 exp −1.1 , is further taken into con103.5 − 84.3 sideration. Then, the scheduling problem is recalculated with results of pr e = 125.19 (p.u.) and pf e = 94.81 (p.u.). From Fig. 2, it can be seen that, after considering the penalty on curtailment, the incremental cost of DRG exponentially increases as the increasing of pr e . It also means that the higher the value of pr e is, the greater the penalty becomes. As a consequence, the optimality and possibility can be traded off in a suitable position, but not the upper or lower bound to unilaterally ensure the optimization or possibility only. 2) DFG, DFHD, and DGP Models: To deal with and investigate the ramping rate limits of DFGs, its discrete form is always formulated into a knapsack problem [15]. For DFG, its


cost function is modeled as the following nonlinear function [4], [13], 2 e e e e = aij pfij,t + bij pfij,t + cij + εij exp ηij pfij,t C pfij,t (5) e,m in e e,m ax ≤ pfij,t ≤ pfij,t pfij,t e,ramp e,ramp e e − pfij,t ≤ pfij,t − pfij,t−1 ≤ pfij,t

(6) (7)

where aij , bij , cij , εij , and ηij represent nonnegative cost coefficients, which are derived from the energy emission of the e is the power output of jth DFG of ith thermal unit [4]; pfij,t f e,m in f e,m ax EB; pij,t and pij,t are the lower and upper bounds; and f e,ramp pij,t is the ramp rate limit. In addition, if the above DFG is a gas turbine, the amount of gas consumption, denoted by gaspij,t , cannot exceed the gas transmission maximum capability deax . For a gas turbine, the relationship between noted by gasp,m ij,t power generation and the amount of gas consumption can be approximately evaluated by the following equation 2 p fe fe (8) gasij,t = Θ μij pij,t + νij pij,t + θij where μij , νij , θij , ξij , and ωij are heat rate coefficients. Θ is the conversion ratio from MW to SCM/h, and the value of Θ is about 84. The model of DFG can be applied for DFHD. However, in different condition, the ramp rate limits of DFHDs are not taken e e and C(hfij,t ) denote the into consideration [4]. And let hfij,t corresponding heat generation and cost function, respectively. If the heating device is a gas boiler, the amount of gas consumption can be approximatively evaluated by e fe (9) gashij,t = Θ hfij,t /ηij,t fe ax where ηi,t is the efficiency. And we let gash,m represent the ij,t gas transmission maximum capability. Finally, the cost function of DGP is modeled as g 3 g 2 g g = aij gij,t + bij gij,t + dij gij,t + cij (10) C gij,t g ,m in g g ,m ax 0 ≤ gij,t ≤ gij,t ≤ gij,t

(11)

g where aij , bij , cij , and dij are nonnegative cost coefficients; gij,t g ,m in g ,m ax is the gas output of jth DGP of ith EB; and gij,t and gij,t are the corresponding lower and upper bounds. It is not difficult to verify that function (10) is convex in the region determined by (11). 3) DCHP Model: The cost function of DCHP is often modeled as the following convex function [5] 2 2 chp chp chp chp = a p h , h + b p + α C pchp ij ij ij ij,t ij,t ij,t ij,t ij,t chp chp + βij hchp (12) ij,t + σij pij,t hij,t + cij chp eij,m pchp ij,t + fij,m hij,t + zij,m ≥ 0, (m = 1, 2, 3, 4) chp,ramp chp chp,ramp ≤ pchp − pij,t ij,t − pij,t−1 ≤ pij,t

(13) (14)

3085

where aij , bij , αij , βij , σij , and cij are nonnegative cost coefficients; eij,m , fij,m , and zij,m are the coefficients of mth linear inequality constraint of DCHP caused by feasible operchp ating region; pchp ij,t and hij,t are the power and heat outputs of chp,ramp jth DCHP of ith EB; and pij,t is the ramp rate limit. If the DCHP unit is fed by nature gas, the total gas consumption can be approximatively evaluated by chp chp (15) pchp gaschp ij,t = Θ ij,t + hij,t /ηij,t chp chp,m ax where ηij,t is the total efficient of the DCHP. gasij,t represents its gas transmission maximum capability. e and 4) DPSD and DHSD Models: We denote pbij,t be SOCij,t as the exchanged power and the stored energy of jth e be DPSD of ith EB at time t, respectively. Therein, we let pbij,t positive for discharging and negative for charging actions. Each DPSD, which cannot be charged and discharged simultaneously, should meet the following dynamic constraints [26], [27]: ax ds,m ax e ≤ pbij,t ≤ pij,t − pch,m ij,t

(16)

in be m ax SOCm ij,t ≤ SOCij,t ≤ SOCij,t e e ch = SOCbij,t−1 − (ςijch δij,t−1 + SOCbij,t

(17) 1 ςijds

ds e δij,t−1 )pbij,t−1 T

(18) ch ds + δij,t−1 ≤1 δij,t−1

(19)

ax ds,m ax where pch,m and pij,t are the maximum charging and disij,t m in ax charging rates; SOCij,t and SOCm ij,t are the lower and upper bounds for allowed energy stored in the corresponding DPSD; ςijch and ςijds are charging and discharging coefficients; ch ds δij,t−1 , δij,t−1 ∈ {0, 1} denote the operation state of the corresponding DPSD determined by the previous scheduling horizon. ch ds = 1 or δij,t−1 = 1 represents charging or disTherein, δij,t−1 charging state. During each scheduling horizon, each DPSD shall be charged when the buying price is cheaper and vice versa. According to [21], the following cost function is used to capture the operations as follows e e = aij (pbij,t + bij )2 (20) C pbij,t

where aij and bij are the cost coefficients. In addition, the model above can also be applied for DHSD. And let hbe ij,t and C(hbe ) denote the corresponding heat output and cost function, ij,t respectively. 5) EL Models: During each scheduling horizon, it is assumed that the amount of equivalent must-run power, heat, and rl mrl mrl gas loads, denoted by pm ij,t , hij,t , and gij,t , are fixed. The utility function of EL demand is the sum of three quadratic functions of power, heat, and gas demands as follows cl cl (21) Ui,t = U pcl ij,t + U hij,t + U gij,t ⎧ l,m ax rl ⎪ 0 ≤ pcl − pm ⎪ ij,t ≤ pij,t ij,t ⎨ l,m ax rl (22) 0 ≤ hcl − hm ij,t ≤ hij,t ij,t ⎪ ⎪ ⎩ l,m ax cl mrl 0 ≤ gij,t ≤ gij,t − gij,t

3086


where ⎧ cl mrl 2 mrl U pij,t = apij pcl + bpij pcl ⎪ ij,t + pij,t ij,t + pij,t ⎪ ⎨ cl h mrl 2 mrl + bhij hcl U hcl ij,t = aij hij,t + hij,t ij,t + hij,t ⎪ ⎪ cl cl ⎩ cl mrl 2 mrl + gij,t + bgij gij,t + gij,t U gij,t = agij gij,t cl cl pcl ij,t , hij,t , and gij,t represent controllable power, heat, and gas load demands, respectively. apij , ahij , and agij are negative utility coefficients while bpij , bhij , and bgij are positive utility coefficients. ax l,m ax l,m ax pl,m are upper bounds of power, heat, and ij,t , hij,t , and gij,t gas loads, respectively. It is worth noting that some energy load demands can be satisfied by various forms of energy supplies. In other words, power, heat, and gas load demands can be transformed into each other in some circumstances. Inspired by Qiu et al. [34], the concept above can be mathematically expressed as ⎧ cl cl m ax ⎪ φm in ≤ pcl ij,t / pij,t + gij,t /Θ ≤ φij,g to p ⎪ ⎨ ij,g to p in cl cl cl m ax φm (23) ij,h to p ≤ pij,t / pij,t + hij,t ≤ φij,h to p ⎪ ⎪ ⎩ φm in ≤ hcl / hcl + g cl /Θ ≤ φm ax ij,t

ij,g to h

ij,t

ij,t

ij,g to h

in m in m in m ax m ax where φm ij,g to p , φij,h to p , and φij,g to h , and φij,g to p , φij,h to p m ax and φij,g to h express the corresponding lower and upper translating percentages.

C. Energy Management of EI In this paper, we mainly focus on the hourly EMP of EI to achieve the co-planning of power, heat, and gas. The objective function is to cooperatively maximize the social welfare shown in (24), while meeting all of the inequality constraints mentioned above and three global equality constraints shown in (25): max F =

n

(Wi,t )

(24)

i=1 n i=1

Δpi,t = 0;

n i=1

Δhi,t = 0;

n

Δgi,t = 0

(25)

i=1

where ⎧ prt Δpi,t + hrt Δhi,t + grt Δgi,t ⎪ ⎪ Wi,t = Ui,t − Ci,t (•) + ⎪ ⎪ re re ⎪ fe fe ⎪ + C h + C p + C h C (•) = C p ⎪ i,t ij,t ij,t ij,t ⎪ ⎪ ij,t ⎪ be be ⎨ chp chp g +C pij,t + C hij,t + C pij,t , hij,t + C gij,t ⎪ Δp = pr e + pf e + pbe + pchp − pm r l − pcl ⎪ i,t ⎪ ij,t ij,t ij,t ij,t ij,t ij,t ⎪ ⎪ chp ⎪ fe re be mrl ⎪ Δhi,t = hij,t + hij,t + hij,t + hij,t − hij,t − hcl ⎪ ij,t ⎪ ⎪ ⎩ g mrl cl Δgi,t = gij,t − gij,t − gij,t . Wi,t is the benefit function of ith EB. Δpi,t , Δhi,t , and Δgi,t are the power, heat, and gas supply–demand mismatches of ith EB, respectively. prt , hrt , and grt are the market power, heat, and gas clearing price at time t, respectively. Note that the ith EB serves as a power, heat, or gas supplier when Δpi,t , Δhi,t , or Δgi,t is positive and vice versa. Remark 2: The coupling relationships among power, heat, and gas of EI are embodied in three aspects: 1) For some energy conversion devices, the natural gas is the clean fuel for power

and heat generation; 2) For DCHPs, the power and heat are cooperatively generated; meanwhile, they also should run in the corresponding feasible operating region; 3) Power, heat, and gas demands can be transformed into each other in some situations. For convenience, let pij,t (hij,t ) represent the power (heat) generation of DRG (DRHD), DFG (DFHD), DPSD (DHSD), DCHP or controllable power (heat) load demand of ith EB, and let gij,t represent the gas supply of DGP or controllable gas load demand of ith EB. The number of decision variables of different participant may be different, e.g., DFG, DCHP, and EL have one, two, and three decision variables, respectively. Without loss of generality, we let eij,t = [pij,t , hij,t , gij,t ]T be a three-dimensional decision variable. Then, each participant of each EB can extend its variable to the form of r l dimensions mrl mrl T and Er t = eij,t . Furthermore, let q i,t = pm ij,t , hij,t , gij,t [prt , hrt , grt ]T . We rewrite V (eij,t ) as cost function or negative utility function of each participant, i.e., V (eij,t ) = Ci,t (•) or V (eij,t ) = Ui,t . Then, by some operations, the EM of EI can be rewritten as

min F =

mi n

V (eij,t )

i=1 j =1

⎧ n mi ⎨ Y e =Q ij ij,t s.t. i=1 j =1 ⎩ Φ (eij,t ) ≤ 0, eij,t ∈ R3 ⇔ eij,t ∈ X ij

(26)

where Q = ni = 1 q i,t represents the sum of must-run energy demands; mi is the number of participates of ith EB; X ij , determined by Φ (eij,t ), is the local closed convex set for eij,t . Y ij = −I if eij,t represents the controllable energy load demand; otherwise, Y ij = I. Φ (eij,t ) represents the local inequality constraints related to eij,t . Besides eij,t , Er t is also an important variable which is expected to be calculated in a distributed fashion. Although, it is not included in problem (26) caused by the three global equality constraints (25), we can also find its optimal solution by implementing the proposed algorithm, which will be discussed in Remark 4.

D. 5-Min Power Scheduling Adjustment Note that the variations of power loads and the fluctuations of DRGs are always seen as 5-min or 15-min, while the heat and gas loads variations are always seen as being hourly [28], [29], [34]. Thus, on the basis of the hourly EMP of EI, an adjustment strategy for power-only devices and controllable power loads is further developed for adapting to the shorter power scheduling problem, i.e., 5-min considered in this paper. Then, the 5-min power loads variations and fluctuations of DRGs can be accommodated by dispatching the power generations of DFGs, DPSDs, and controllable power loads. In this paper, the goal is to make the 5-min power scheduled result of each power participant follow its corresponding hourly scheduled result as close as possible. Meanwhile, in general, the larger the capacity upper is, the heavier task the participant has to smooth out the variations and fluctuations. This concept above can be


mathematically expressed as n 1 e e 2 be 2 min Obj = (pfij,t− −pfij,t ) +ij (pbe ij,t− −pij,t ) f e,m ax p ij,t i=1

1 cl 2 + l,m ax (pcl (27) ij,t− − pij,t ) pij,t

3087

where Γ (λ) = max −V (eij,t ) − λT Y ij eij,t , κ∗ = Q/ ei j , t ∈X i j n m , λ = [λ1 , λ2 , λ3 ]T is the Lagrange dual variable, and i i=1 λ1 , λ2 , and λ3 represent the power, heat, and gas dual variables, respectively. Note that λ is global variable, which is not suitable for distributed implementation. To solve (29) in a distributed fashion, a set of auxiliary variables are added into (29) which is rewritten as

where ij = min{ p ch,1m a x , p ds,1m a x }; the subscript “t− ” is utiij,t

ij,t

lized to represent the th 5-min interval of tth hourly interval. e For example, pfij,t− is the th 5-min power output of jth DFG of ith EB during tth hour, and the similar definitions for pbe ij,t− , chp re mrl cl pij,t− , pij,t− , pij,t− , and pij,tm . In addition, problem (27) should satisfy the 5-min power supply–demand balance constraint given by n

chp e e (prij,t− + pfij,t− + pbe ij,t− + pij,t−

i=1 rl cl − pm ij,t− − pij,t− ) = 0.

(28)

Since DCHPs generate not only power but also heat, they are not employed to participate in the 5-min power scheduling. chp Thus, the value of pchp ij,t− is the same as pij,t determined by its corresponding hourly scheduled result. Moreover, all the DFGs, DPSDs, and controllable power loads also should meet their 5-min timescale operation constraints which are similar to the corresponding hourly ones.

min

s.t. λij = ϕij, ij¯ ; λij¯ = ϕij, ij¯ where Γ (λij ) = max

In this paper, we consider that problem (26) is solvable, i.e., there exist a set of values {eij,t | i = 1, . . . , n; j = 1, . . . , mi } satisfying all the equality and inequality constraints. Moreover, it is worth noting that (26) is a convex optimization problem and Φ (eij,t ) in (26) is local and affine function. Thus, the refined Salter condition is satisfied [30, pp. 226-227]. That also means (26) has a zero duality gap. As a consequence, the primal problem and its dual problem hold the same optimal solution, which is the saddle point of the Lagrangian [31]. Further, in order to solve problem (26), we can transform it into its dual problem given by min

mi n Γ (λ) + λT κ∗ i=1 j =1

(29)

ei j , t ∈X i j

¯ ∈ Nij ; ∀ij

∀ij ∈ V (30)

−V (eij,t ) − λTij Y ij eij,t ; λij is lo-

L {λij } , ϕij, ij¯ , μ = +

mi n

Γ (λij ) + λTij κ∗

i=1 j =1

μTij, ij¯ λij − ϕij, ij¯ + μTij¯ ,ij λij¯ − ϕij, ij¯

¯ ∈N i j ij

c λij − ϕij, ij¯ 2 + λij¯ − ϕij, ij¯ 2 2 2 2

(31)

where μ := μij, ij¯ , μij¯ ,ij , μij, ij¯ , and μij¯ ,ij are corresponding Lagrange multipliers. We first assume that Γ (λij ) can be calculated, then the following distributed algorithm is proposed to obtain the optimal λij of problem (30), i.e., v ij (k) = v ij (k − 1) + c

λij (k) − λij¯ (k)

(32)

¯ ∈N i j ij

κij (k + 1) =

B. Main Algorithm

cal estimated dual variable. Equations (29) and (30) are equivalent optimization problems in essence. The reason is that, for each pair of neighbors, we add an auxiliary variable and let λij = ϕij, ij¯ = λij¯ , such that, over connected graph, all λij must be equal to the same value which is the λ in (29). In this way, (30) holds a dissociable structure. Then, the quadratically augmented Lagrangian function of (30) is given by

+

Consider an EI with n EBs, where each EB has mi participants. A graph G = (V, E) is used to model the EI, where V = {vij | i = 1, . . . , n; j = 1, . . . , mi } is a set of nodes representing participants and E ⊂ V × V is a set of edges. Therein, ¯ node can the edge (vij , vij¯ ) denotes that ijth node and ijth exchange information between each other as needed. The set of neighbors of ijth node is denoted by Nij with cardinality |Nij |. In this paper, we assume that graph G is connected.

Γ (λij ) + λTij κ∗

i=1 j =1

III. DISTRIBUTED ALGORITHM A. Graph Theory

mi n

mi n

wij, ij¯ κij¯ (k)

(33)

i=1 j =1

λij (k + 1) = arg min Γ (λij ) + λTij κij (k + 1) + v Tij (k) λij λi j

+c

λij − λij (k) + λij¯ (k) /22 2

(34)

¯ ∈N i j ij

where c is a constant for the algorithm; v ij (k) = ¯ ∈ 2 ij¯ ∈N i j μij, ij¯ (k); wij, ij¯ = 2/ |Nij | + Nij¯ + Λ if ij ¯ = Nij , wij, ij¯ = 1 − ij¯ ∈N i j 2/ |Nij | + Nij¯ + Λ if ij ij, otherwise, wij, ij¯ = 0 [21]. Therein, Λ is a small number. The initial values are set as μij, ij¯ (−1) = μij¯ ,ij (−1) = v ij (−1) = λij (0) = 0 and κij (0) =q i,t /mi . In this way, μij, ij¯ (k) = −μij¯ ,ij (k) is always satisfied, and then we can get the following

3088


fact which will be used for the convergence analysis, mi n

v ij (k) =

mi n

μij, ij¯ (k) − μij¯ ,ij (k) = 0.

¯ ∈N i j i=1 j =1 ij

i=1 j =1

Theorem 1: We assume the communication graph G is connected and Γ (λij ) can be calculated. Then for all ij ∈ V, staring form the corresponding initial conditions mentioned above, the algorithm composed of (32), (33), and (34) can make λij (k) → λ∗ as k → ∞. (The proof of Theorem 1 is presented in Appendix A.) However, in Theorem 1, it is very difficult to directly calculate λij (k + 1) of (34) in fact, because Γ (λij ) is a maximum subproblem determined by eij,t . To take the explicit form of λij (k + 1) into account, (34) is rewritten as λij (k + 1) = arg min max −V (eij,t ) λi j ei j , t ∈X i j

− λTij Y ij eij,t + λTij κij (k + 1) + v Tij (k) λij 2 λij − λij (k) + λij¯ (k) /2 2 +c ¯ ∈N i j ij

= arg

max min −V (eij,t )

ei j , t ∈X i j

where Θ =

λi j

λij (k) + λij¯ (k) − v ij (k)−Y ij eij,t +

ij ∈N i j κij (k+1)/c . The reason for obtaining the above equation is that the objective function is concave for eij,t when λij is given; meanwhile, the one is convex for λij when eij,t is given. Thus, the min-max problem and its corresponding max-min problem can be transformed into each other. Then, based on the final form of (35), the calculation of λij and eij,t can be decoupled as follows: λij (k) + λij¯ (k) − v ij (k) λij (k + 1) = ¯

¯ ∈N i j ij

− Y ij eij,t (k + 1) + κij (k + 1) /c /(2 |Nij |) eij,t (k + 1) = arg

(36)

min V (eij,t ) + c Θ 22

ei j , t ∈X i j

2 λij (k) + λij¯ (k) /2 2 . / (4 |Nij |) − c

controller or a leader. To be specific, each node only needs to exchange the information of λij (k) and κij (k) with its neighbors. And the updating of variables v ij (k), κij (k + 1), eij,t (k + 1), and λij (k + 1) are implemented via local calculation based on their own and neighbors information. Theorem 2: Let G be a connected graph. If the problem of (26) is feasible, then the Algorithm 1 is stable, and each eij,t can converge to its optimal point, i.e., eij,t (k) → e∗ij,t ,

+ c |Nij | λij − Θ/ (2 |Nij |) 22 − c Θ 22 / (4 |Nij |) λij (k) + λij¯ (k) /22 , +c (35) 2 ¯ ∈N i j ij

Algorithm 1: Distributed-Consensus-ADMM Algorithm. Initialize: Set v ij (−1) = 0, λ(0) = 0, eij,t any admissible values and κij (0) =q i,t /mi for each participant. Iteration: (k ≥ 0) 1) Transmit λij (k) and κij (k) to its neighbors 2) Update v ij (k) based on (32) 3) Update κij (k + 1) based on (33) 4) Update eij,t (k + 1) by solving problem (37) 4) Update λij (k + 1) based on (36) 5) Let k = k + 1, turn to step 1.

(37)

¯ ∈N i j ij

Thanks to (36) and (37), our main algorithm, referred as distributed-consensus-ADMM algorithm, is further proposed and summarized as in Algorithm 1. Remark 3: The distributed-consensus-ADMM algorithm is a fully distributed algorithm, which only requires the local communication and calculation among nodes without a central

k → ∞.

The proof of Theorem 2 is presented in Appendix B. The proposed distributed-consensus-ADMM algorithm, which can guarantee improved reliability, robustness, flexibility, and scalability, and provide better plug-and-play functionality, is a promising option for applying in the future EI. First, the work of each node is based only on the local sharing of information among neighbors to find the optimal operating point, resulting in alleviated single-node congestion. Moreover, the EI can work well even under single and even several communication link failure(s), as long as the communication network maintains connection. Thus, the EI becomes more robust against singlepoint failure by implementing the proposed algorithm. Second, a node can join in the EI by establishing communication links with its local neighbors only. Consequently, it is more flexible for possible integration and expansions of different energy resources. Third, regardless of the scale of the EI, each node only performs a local optimization. It also means that the communication and computation burden can be divided among the distributed processors, thus being able to encourage better scalability. Last but not least, distributed implementation fashion can provide desired plug-and-play capability to face the topology variability of the future EI. It is worth noting that problem (27), with all of its equality and inequality constraints, does not have complex coupling relationship among variables, which can be seen as a simplified version of problem (26). Therefore, the proposed distributedconsensus-ADMM algorithm can also be employed to solve the 5-min power scheduling problem. Remark 4: In our algorithm, λ has important physical significance, that is λ is the same as negative energy market clearing price, i.e., λ = −Er t

(38)


Fig. 3.

3089

Energy mismatch in case 1.

Fig. 5. Energy outputs/demands in case 1: (a) power outputs/demands, (b) heat outputs/demands, (c) gas outputs/demands.

distributed ADMM algorithm. In the end, not only the optimal energy outputs/demands but also the optimal energy market clearing price can be calculated by decoupling the max-min problem. IV. SIMULATION RESULTS

Fig. 4. Lagrange dual variable in case 1: (a) power dual variable, (b) heat dual variable, (c) gas dual variable.

which implies that each participant can estimate the global energy market clearing price in a distributed fashion because they can distributively calculate system λ. The proof of (38) is shown in Appendix C. Remark 5: In our algorithm, the updating of λ(k + 1) in (34) and (36) does not include κ∗ = Q/ ni = 1 mi , which implies this paper does not require the assumption that each participates can access κ∗ . However, the existing distributed ADMM algorithms, e.g., [32] and [33], etc., need this requirement unless κ∗ = 0. Note that, in the EI, κ∗ cannot directly be known by each participant because it is a global constant during each scheduling horizon. In this paper, a dynamic average consensus algorithm, i.e., (33), is used to estimate the κ∗ in a fully distributed fashion. Then, κij (k) is employed to design our

The configuration of the test system with five EBs is shown in Fig. A1 in Appendix D. The system cost/utility parameters and constraints are listed in Table A1 in Appendix D. Energy translating percentages and the heat rate coefficients can be found in [34] and [35], respectively. We unify the energy scale, i.e., 1 p.u.= 1 MW for power (or heat), 1 p.u.= 84 SCM/h for gas, and 1 p.u.=1 $/MWh for price. A. Convergence and Profit Analysis In this case study, the focus is on showing the convergence of the proposed algorithm and analyzing the profits for EBs. The must-run energy loads for EB1 to EB5 are initialized with [150 (p.u.), 124 (p.u.), 50 (p.u.)], [105 (p.u.), 150 (p.u.), 60 (p.u.)], [85 (p.u.), 135 (p.u.), 80 (p.u.)], [100 (p.u.), 90 (p.u.), 50 (p.u.)], and [50 (p.u.), 140 (p.u.), 0 (p.u.)], respectively. Running the proposed distributed-consensus-ADMM algorithm at t0 , the power, heat, and gas mismatches, Lagrange dual variables, and outputs/demands are obtained as shown in Figs. 3–5,

3090


TABLE I RESULTS OF ENERGY OUTPUTS AND DEMANDS Num EB1

Δp1 , t 0 EB2

Δp2 , t 0 EB3

Δp3 , t 0 EB4

Δp4 , t 0 EB5

Δp5 , t 0

Fig. 6.

pi j , t 0

hi j, t 0

gi j , t 0

Profits

N.11 N.12 N.13 N.14 N.15 N.16 N.17 N.18 = 135.6014

48.7267 0.0000 99.9955 0.0000 81.9282 0.0000 205.9113 175.0625 150.9603 0 0.0000 0.0000 0.0000 140.4004 0.0000 93.8493 Δ h 1 , t 0 = 285.3122

0.0000 18 323 0.0000 0.0000 0.0000 151.1224 697.1795 0.0000 0.0000 Δ g 1 , t 0 = 496.0571

N.21 N.22 N.23 N.24 N.25 N.26 N.27 = −256.3394

79.8126 0.0000 40.0000 0.0000 97.5582 124.7394 368.7102 234.9917 0.0000 149.9000 0.0000 94.9338 0.0000 −50.9579 Δ h 2 , t 0 = −66.3764

0.0000 16 007 0.0000 0.0000 640.0000 0.0000 0.0000 0.0000 Δ g 2 , t 0 = −700.0000

N.31 N.32 N.33 N.34 N.35 N.36 = −64.4906

43.0000 0.0000 86.3168 0.0000 139.1637 145.6758 247.9711 214.2604 0.0000 −66.9453 0.0000 115.3434 Δ h 3 , t 0 = −155.1865

0.0000 7032 0.0000 0.0000 100.0000 0.0000 0.0000 Δ g 3 , t 0 = −180.0000

N.41 N.42 N.43 N.44 N.45 N.46 = 199.4447

57.3022 0.0000 63.8648 0.0000 74.8884 0.0000 137.9981 168.0150 34.6088 169.6302 0.0000 0.0000 Δ h 4 , t 0 = −91.6152

0.0000 15 329 0.0000 0.0000 0.0000 77.3927 731.8679 Δ g 4 , t 0 = 604.4752

N.51 N.52 N.53 N.54 = −14.2161

59.7181 107.4697 131.4040 0.0000 Δh5 , t 0 =

0.0000 965 0.0000 220.5323 0.0000 Δ g 5 , t 0 = −220.5323

0.0000 155.1540 104.8087 117.5203 27.8656

Energy mismatch in case 2.

respectively. From Fig. 3, It can be observed that the nestimated n Δp , power, heat, and gas mismatches, i.e., i,t i=1 i=1 Δhi,t , and ni=1 Δgi,t in (25), converge to zero after about 75 iterations. This results also mean that the three global equality constraints are satisfied to further verify the power heat and gas supply–demand balance subplot. More importantly, the estimated power, heat, and gas dual variables of each participant can converge to three different values, i.e., λ1 = −32.6887 (p.u.), λ2 = −23.6611 (p.u.), and λ3 = −15.2825 (p.u.), respectively. Meanwhile, the final power, heat, and gas outputs/demands of

Fig. 7. Lagrange dual variable in case 2: (a) power dual variable, (b) heat dual variable, (c) gas dual variable.

each participant, shown in Table I, are within their operational ranges, which implies all the set of inequality constraints are satisfied. Therefore, the optimization goal is fulfilled. According to the results above, we can also get that the cooperative power, heat, and gas market clearing prices are 32.6887, 23.6611, and 15.2825 (p.u.), respectively. In addition, the final profit of each EB is also summarized in Table I. Each EB may play different roles based on its and neighbors’ state. Let us use EB4 as an illustration. Note that Δp4,t 0 > 0 and Δg4,t 0 > 0, so it serves as a power and gas supplier, and these parts of excess self-generated power and heat will be sold to its neighbors. Meanwhile, Δh4,t 0 < 0, so it also serves as a heat consumer and will buy deficit heat from its neighbors. B. Plug and Play Test In this case study, the focus is on testing the plug and play performance of the proposed algorithm. All the initial conditions are the same as in case 1. On the basis, at t1 , the DFG of EB1, DCHP of EB2, and EL of EB5 are removed from the system at the same time, and the variables related to them are set to zero. From Figs. 6–8, it can be seen that the total power, heat, and gas


3091

Fig. 8. Energy outputs/demands in case 2: (a) power outputs/demands, (b) heat outputs/demands, (c) gas outputs/demands.

Fig. 9.

Average number of iterations required for convergence.

mismatches converge to zero; meanwhile, the estimated power, heat, and gas dual variables converge to their corresponding optimal values. Of course, the remaining participates have to adjust their power, heat, and gas outputs/demands to compensate for the amount of power, heat, and gas previously supplied or consumed by the disconnected DFG, DCHP, and EL, and finally converge to a new solution. Further, at t2 , they are all plugged again back to the system. The simulation results are also shown in Figs. 6–8. It can be seen that the system again converges to

Fig. 10.

5-min power dispatching results.

3092


the new solutions responding to the new topological change. In addition, the final convergence solutions after t2 are the same as the one prior to the disconnection. This implies that the proposed algorithm provides good plug and play capability. C. Scalability Analysis This case study focuses on analyzing the computational scalability of the proposed distributed-consensus-ADMM algorithm with respect to the number of nodes (or participants). The number of nodes is set as [10, 30, 50, . . . , 970] with total 50 groups. Meanwhile the percentage of DRG, DFG, DPSD, DRHD, DFHD, DHSD, DCHP, DGP, and EL are set as 10%, 10%, 10%, 10%,10%, 10%,10%, 10%, and 20% in each example. The parameters of each type of participant are randomly selected and the corresponding ones as shown in Table A1. The stopping criterion of the proposed algorithm is based on calculating the solution accuracy of the following equation MSE =

n mi 1 eij,t (k) − e∗ij 2 / e∗ij 2 N i=1 j =1

(39)

where N is the number of nodes. The algorithm is set to stop once MSE below the preset accuracy denoted as acc. Next, we will analyze the trend of the average number of iterations, expressing the average computational complexity, with the growth of the number of nodes under the same acc. For each group of nodes, it will be randomly tested with 100 times. By employing the same data fitting method in [36], the fitting curve and functions with acc = 10−4 are shown in Fig. 9. Therein, R2 ≤ 1 represents the coefficient of determination, and the higher value of R2 is, the better fitting effect becomes. From Fig. 9, it can be observed that the average number of iterations for the proposed algorithm is approximately logarithmic (but not exponential) growth with the increase of the number of nodes. Thus, we can verify that the proposed algorithm can provide good scalability on the part of computations burden. In addition, for each participant, the number of communication times directly comes from the number of computation iterations. Therefore, it can be concluded that the proposed algorithm also exhibits better scalable behavior in terms of communications burden. D. 5-Min Power Scheduling Adjustment This case study focuses on analyzing the 5-min power dispatching problem on the basis of 1-h dispatching results in case 1. The 5-min power must-run loads variations and fluctuations of DRGs are shown in Fig. 10(a) and (b), while the corresponding dispatching results of DFGs, DPSDs, and controllable loads during each scheduling horizon are shown in Fig. 10(c)–(e). It can be seen that we can not only adjust the power generations of DFGs and DPSDs but also the power demand of controllable loads to smooth out the renewable energy fluctuations and loads variations. Note that, for each DFG, DPSD or controllable power load, the 5-min dispatching results can track its hourly dispatching result in a better way, to further

ensure the effectiveness of maximum social welfare in hourly scheduling. In addition, the participant, which holds larger capacity, is able to undertake task for accommodating the power variations and fluctuations. For example, the changes of the ax 5-min dispatching curve of participant N.15 with pl,m 15,t = 1000 ax (p.u.) is larger than that of participant N.53 with pl,m 53,t = 750 (p.u.) as seen in Fig. 10(e). From the discussion mentioned above, the designed 5-min power scheduling adjustment strategy cannot only meet the requirement of hourly scheduling goal but also smooth out 5-min loads variations and fluctuations of DRGs. V. CONCLUSION In this paper, an innovative energy management framework has been proposed for the future EI. Along with the EB, the system model has be established, which can better reflect the features and requirements of EI in a more precise way. Furthermore, the EMP of EI has been finally formulated as a distributed nonlinear coupling optimization problem, which can be effectively solved by the proposed distributed-consensus-ADMM algorithm. In addition, the proposed algorithm is a fully distributed algorithm, where each participant requires only the information from its neighbors to estimate the optimal energy market clearing price and calculate the optimal energy outputs/demands. The effectiveness of the proposed algorithm has been demonstrated by several simulation results. APPENDIX A A. Proof of Theorem 1 According to the definition of wij, ij¯ , (33) makes up an average-consensus protocol with κij (0) =q i,t /mi . Then, after a sufficient long time K, we can get κij (K) =

mi n i

=

j

n

κij (0)/

n

(mi q i,t /mi ) /

i

mi

i n

mi =κ∗

(40)

i

which implies that if we construct a variable κij for each participant, then, starting from κij (0) =q i,t /mi , κij (k) can converge to κ∗ by implementing (33). Note that each participant can easily access q i,t /mi . Now, for (31), ADMM results in the following updating iterative process 1) Step 1: Updating of local estimate λij (k + 1) = arg min L {λij } , ϕij, ij¯ (k) , μ(k) λi j

= arg min Γ (λij ) + λTij κij (k + 1) λi j

+2

μTij, ij¯ (k) λij

¯ ∈N i j ij

+c

λij − ϕij, ij¯ (k)2 . 2

¯ ∈N i j ij

(41)


The reason for obtaining the above equation is that participant ¯ is the neighbor of participant ij; meanwhile, participant ij is ij ¯ Moreover, we use κij (k + 1) computed also the neighbor of ij. by (33) to replace κ∗ , since each participant cannot directly access κ∗ . 2) Step 1: Updating of auxiliary

ϕij, ij¯ (k + 1) = arg min L {λij (k + 1)} , ϕij, ij¯ ϕi j , i¯j

, μ(k)

Finally, by plugging (44) and v ij (k) = 2 (41), λij (k + 1) is simplified to

2 c + λij (k + 1) − ϕij, ij¯ 2 2 2 + λij¯ (k + 1) − ϕij, ij¯ 1 = (μij, ij¯ (k) + uij¯ ,ij (k)) 2c 1 + (λij (k + 1) + λij¯ (k + 1)). 2

(42)

λi j

+c

λij − λij (k) + λij¯ (k) /22 . 2

It can be seen that (46) and (47) are the same as (32) and (34), respectively.

In Theorem 1, it has been verified that λij can converge a common value, i.e., λ∗ . The rest is to verify that eij,t (k + 1) can converge to its corresponding optimal point. To begin with, according to the KKT (Karush–Kuhn–Tucker) condition of (37), we can get

k →∞

+ c(λij (k + 1) − ϕij, ij¯ (k + 1)) μij¯ ,ij (k + 1) = μij¯ ,ij (k) + c(λij¯ (k + 1) − ϕij, ij¯ (k + 1)).

(43)

The next is to verify that (33) and (41)–(43) can be simplified to (32)–(34). Note that, for (43), if the initial values of μij, ij¯ and μij¯ ,ij are set to zero, then μij, ij¯ (k + 1) = −μij¯ ,ij (k + 1) is always satisfied for all k ≥ 0. Thus, (42) can be further simplified to 1 (λij (k + 1) + λij¯ (k + 1)). 2

(44)

/(2 |Nij |) /∂eij,t (k + 1) = 0,

Furthermore, let v ij (k) = 2

¯ ∈N i j ij

= 0,

μij, ij¯ (k − 1) + c

¯ ∈N i j ij

= v ij (k − 1) + c

mi n

¯ ∈N i j ij

¯ ∈N i j ij

=c

(46)

(49)

Y ij eij,t (k + 1) −

mi n

κij (k + 1)

i=1 j =1

mi n

2λij (k + 1) − λij (k) − λij¯ (k)

¯ ∈N i j i=1 j =1 ij

(45)

λij (k) − λij¯ (k)

λij (k) − λij¯ (k) .

eij,t ∈ X ij .

Moreover, upon summing (36), we can get

+

ing to (45), we can get

(48)

k →∞

i=1 j =1

μij, ij¯ (k). Then, accord-

eij,t ∈ X ij ,

T ∗ ⇒ lim ∂V (eij,t (k + 1)) + Y ij λ /∂eij,t (k + 1)

By substituting (44) in (43), μij, ij¯ (k + 1) can be rewritten as the following form μij, ij¯ (k + 1) = μij, ij¯ (k) c + [λij (k + 1) − λij¯ (k + 1)]. 2

¯ ∈N i j ij

− v ij (k) − Y ij eij,t (k + 1) + κij (k + 1) /c

μij, ij¯ (k + 1) = μij, ij¯ (k)

(47)

T λij (k) + λij¯ (k) lim ∂V (eij,t (k + 1)) + Y ij

3) Step 3: updating of Lagrange multiplier

v ij (k) = 2

μij, ij¯ (k) in

Proof of Theorem 2

2

ϕij, ij¯ (k + 1) =

¯ ∈N i j ij

¯ ∈N i j ij

− (μTij, ij¯ (k) + μTij¯ ,ij (k))ϕij, ij¯

ϕi j , i¯j

λij (k + 1) = arg min Γ (λij )+λTij κij (k + 1) + v Tij (k) λij

= arg min

3093

mi n

v ij (k),

eij,t ∈ X ij .

i=1 j =1

It is worth noting that, as k → ∞, λij (k + 1) = λij (k) = λij¯ (k) = λ∗ , κij (k + 1) = κ∗ and the sum of v ij (k) is equal to zero mentioned above, then we can obtain lim

k →∞

mi n

Y ij eij,t (k + 1)−Q = 0,

eij,t ∈ X ij . (50)

i=1 j =1

Note that (49) and (50) are the optimality conditions of (26). Thus, we can conclude that eij,t (k + 1) → e∗ij,t as k → ∞.

3094

Fig. A1.


Test system with five EBs.

Proof of (38) In the interest of clarity, we analyze the case where the local inequality constraints are not taken into account. For a given energy selling (or bidding) price denoted by Er ij,t . The estimated profit function of EB i is

energy market clearing price, is determined by the conditions ⎧ ∂V (eij,t ) ⎪ − Er T Y ij = 0 ⎪ ⎨ ∂e ij,t . (52) mi n ⎪ ⎪ ⎩ Y ij eij,t = Q i=1 j =1

Wi,t =

mi j =1

Wij,t =

mi

(−V (eij,t ) + Er Tij,t (Y ij eij,t −q i,t )).

j =1

The corresponding maximum profit can be calculated by ∂Wij,t ∂V (eij,t ) =− + Er Tij,t Y ij = 0. ∂eij,t ∂eij,t

(51)

To obtain maximum profit, each participant tends to adjust its energy outputs (or energy demands) until the marginal cost (or utility) is equal to Er ij,t . Stimulated by energy supplier-demand mismatches, each participant will adjust its selling (or bidding) to meet the supplier-demand balance. For example, a supplier is willing to increase its selling price and produce more energy for deficit market energy demand. We consider that the market is comparatively competitive. Thus, the market is cleared when all the selling and bidding price are equal to the same value, i.e., Er ij,t = Ert for all participates. And Ert , i.e., so-called

Now, we consider the Lagrangian function of (26), i.e., ⎛ ⎞ mi mi n n V (eij,t ) + λT ⎝ Y ij eij,t − Q⎠ . L(F ) = i=1 j =1

i=1 j =1

(53) The optimal operating point is determined by the following KKT conditions ⎧ ∂V (eij,t ) ⎪ + λT Y ij = 0 ⎪ ⎨ ∂e ij,t . (54) mi n ⎪ ⎪ ⎩ Y ij eij,t = Q i=1 j =1

Compared (52) with (54), it can be observed that λ = −Er t . That indicates the calculation of λ and Er t is equivalent. Configuration and Parameters of The Text System See Fig. A1 and Table A1.


3095

TABLE A1 PARAMETERS OF THE TEST SYSTEM ap

bp

ah

bh

ag

bg

pl , m a x

hl , m a x

gl,m ax

− 0.0163 − 0.0150 − 0.0120 − 0.0175 − 0.0083

42.50 46.75 40.80 37.40 35.70

− 0.0100 − 0.0057 − 0.0075 − 0.0150 − 0.0107

25.50 28.05 28.90 31.45 28.90

− 0.01670 − 0.00625 − 0.01420 − 0.03500 − 0.02770

22.0 29.7 20.9 24.2 27.5

1000 700 550 600 750

800 600 500 1000 550

500 700 180 500 600

a

b

α

β

σ

c

p chp, ramp

N.14

0.0345

14.5

0.030

4.2

0.031

230

45

N.23

0.0435

36

0.027

0.6

0.011

124

40

N.33

0.0289

24.5

0.022

9.6

0.022

250

35

N.44

0.0350

8.9

0.025

5.6

0.041

167

35

N.52

0.0443

28

0.012

2.1

0.013

146

40

DPSD N11 N31 N41

a 0.028 0.025 0.021

b 535 600 721

ςch 0.98 0.97 0.95

ςds 0.98 0.97 0.95

−p c h , m a x − 90 − 100 − 120

pd s , m a x 90 100 120

S O C m in 35 52 59

SOC m ax 350 520 590

S O C 0b e 240 95 260

DHSD N27 N35

a 0.013 0.010

b 961 1250

ςch 0.98 0.98

ςds 0.98 0.98

−H c h , m a x -150 -200

H ds,m ax 150 200

S O C m in 62 75

SOC m ax 620 750

S O C 0b e 560 270

DFG N.13 N.22 N.32 N.42

a 0.040 0.032 0.023 0.054

b 25 32 28 27

c 99 150 110 50

ε 50 45 45 28

η 0.010 0.011 0.008 0.013

pf e , m in 30 40 30 50

pf e , m a x 150 180 160 130

pf e , r a m p 45 35 30 45

DFH N.18 N.26

a 0.027 0.031

b 18 17

c 60 30

ε 35 30

η 0.008 0.010

hf e , m in 50 40

DGP N.16 N.46

a 2 × 10 −6 4 × 10 −6

b 0.006 0.004

c 50 80

d 4 3

g g , m in 100 80

gg ,m ax 1500 1400

DRG N.12 N.21 N.43 N.51

b 0.11 0.15 0.11 0.10

ε 748 1300 1215 1550

γ -1.1 -1.2 -1.5 -0.5

pr e 84.3 76.2 71.3 56.4

pr e 103.5 94.5 88.4 72.4

DRH N.17 N.25 N.36 N.54

EL N.15 N.24 N.34 N.45 N.53 DCHP

REFERENCES [1] A. Q. Huang, M. L. Crow, G. T. Heydt, J. P. Zheng, and S. J. Dale, “The future renewable electric energy delivery and management (FREEDM) system: The energy internet,” Proc. IEEE, vol. 99, no. 1, pp. 133–148, Nov. 2011. [2] Q. Sun, Y. Zhang, H. He, D. Ma, and H. Zhang, “A novel energy functionbased stability evaluation and nonlinear control approach for energy internet,” IEEE Trans. Smart Grid, vol. 8, no. 3, pp. 1195–1210, May 2017. [3] Q. Sun, R. Han, H. Zhang, J. Zhou, and J. M. Guerrero, “A multiagent-based consensus algorithm for distributed coordinated control of distributed generators in the energy internet,” IEEE Trans. Smart Grid, vol. 6, no. 6, pp. 3006–3019, Nov. 2015. [4] T. Niknam et al., “A new multi-objective reserve constrained combined heat and power dynamic economic emission dispatch,” Energy, vol. 41, no. 1, pp. 530–545, 2012.

feasible operating region p chp +0.178h chp ≤ 247 −p chp +1.78h chp ≤ 105.7 −p chp −0.17h chp ≤ −98.8 −h chp ≤ 0 p chp +0.115h chp ≤ 125.8 −p chp +1.158h chp ≤ 46.89 −p chp −0.053h chp ≤ −44 −h chp ≤ 0 p chp +0.1375h chp ≤ 185 −p chp +1.51h chp ≤ 79.35 −p chp −0.13h chp ≤ −60 −h chp ≤ 0 p chp +0.125h chp ≤ 159 −p chp +1.205h chp ≤ 64.46 −p chp −0.122h chp ≤ −55 −h chp ≤ 0 p chp +0.202h chp ≤ 225 −p chp +2.4h chp ≤ 264.9 −p chp −0.059h chp ≤ −70 −h chp ≤ 0

hf e , m a x 250 180

b 0.12 0.13 0.15 0.22

ε 534 1760 756 570

γ -1.3 -0.9 -1.5 -1.0

hr e 133.2 149.9 112.2 109.0

re

h 148.2 166.4 127.0 124.5

[5] E. Abdollahi, H. Wang, and R. Lahdelma, “An optimization method for multi-area combined heat and power production with power transmission network,” Appl. Energy, vol. 168, pp. 248–256, 2016. [6] M. R. Benam, S. S. Madani, S. M. Alavi, and M. Ehsan, “Optimal configuration of the CHP system using stochastic programming,” IEEE Trans. Power Del., vol. 30, no. 3, pp. 1048–1056, Jun. 2015. [7] T. Fang and R. Lahdelma, “Optimization of combined heat and power production with heat storage based on sliding time window method,” Appl. Energy, vol. 162, pp. 723–732, 2016. [8] D. Papadaskalopoulos and G. Strbac., “Nonlinear and randomized pricing for distributed management of flexible loads,” IEEE Trans. Smart Grid, vol. 7, no. 2, pp. 1137–1146, Mar. 2016. [9] N. Rahbari-Asr, M. Y. Chow, J. Chen, and R. Deng, “Distributed real-time pricing control for large scale unidirectional V2G with multiple energy suppliers,” IEEE Trans. Ind. Informat., vol. 12, no. 5, pp. 1953–1962, Oct. 2016.

3096


[10] W. Zeng and M. Y. Chow, “A reputation-based secure distributed control methodology in D-NCS,” IEEE Trans. Ind. Electron., vol. 61, no. 11, pp. 6294–6303, Nov. 2014. [11] S. Yang, S. Tan, and J. Xu, “Consensus based approach for economic dispatch problem in a smart grid,” IEEE Trans. Powers Syst., vol. 28, no. 4, pp. 4416–4426, Nov. 2013. [12] G. Chen, F. L. Lewis, E. N. Feng, and Y. Song,“ Distributed optimal active power control of multiple generation systems,” IEEE Trans. Ind. Electron., vol. 62, no. 11, pp. 7079–7090, Nov. 2015. [13] H. Xing, Y. Mou, M. Fu, and Z. Lin, “Distributed bisection method for economic power dispatch in smart grid,” IEEE Trans. Powers Syst., vol. 30, no. 6, pp. 3024–3035, Nov. 2015. [14] F. Guo, C. Wen, J. Mao, and Y.-D. Song, “Distributed economic dispatch for smart grids with random wind power,” IEEE Trans. Smart Grid, vol. 7, no. 3, pp. 1572–1583, May 2016. [15] Y. Xu, W. Zhang, and W. Liu, “Distributed dynamic programming-based approach for economic dispatch in smart grids,” IEEE Trans. Ind. Informat., vol. 11, no. 1, pp. 166–175, Feb. 2015. [16] G. Binetti, A. Davoudi, F. L. Lewis, D. Naso, and B. Turchiano, “Distributed consensus-based economic dispatch with transmission losses,” IEEE Trans. Powers Syst., vol. 29, no. 4, pp. 1711–1720, Jul. 2014. [17] G. Chen, J. Ren, and E. N. Feng, “Distributed finite-time economic dispatch of a network of energy resources, ”IEEE Trans. Smart Grid, vol. 8, no. 2, pp. 822–832, Mar. 2017, doi: 10.1109/TSG.2016.2516017. [18] G. Binetti, A. Davoudi, D. Naso, B. Turchiano, and F. L. Lewis, “A distributed auction-based algorithm for the nonconvex economic dispatch problem,” IEEE Trans. Ind. Informat., vol. 10, no. 2, pp. 1124–1132, May 2014. [19] Z. Wang, W. Wu, and B. Zhang, “A fully distributed power dispatch method for fast frequency recovery and minimal generation cost in autonomous microgrids, ”IEEE Trans. Smart Grid, vol. 7, no. 1, pp. 19–31, Jan. 2016. [20] N. Rahbari-Asr, U. Ojha, Z. Zhang, and M.-Y. Chow, “Incremental welfare consensus algorithm for cooperative distributed generation/demand response in smart grid,” IEEE Trans. Smart Grid, vol. 5, no. 6, pp. 2836–2845, Nov. 2014. [21] Y. Xu and Z. Li, “Distributed optimal resource management based on the consensus algorithm in a microgrid,” IEEE Trans. Ind. Electron., vol. 62, no. 4, pp. 2584–2592, Apr. 2015. [22] W. Zhang, Y. Xu, W. Liu, C. Zang, and H. Yu, “Distributed online optimal energy management for smart grids,” IEEE Trans. Ind. Informat., vol. 11, no. 3, pp. 717–727, Jun. 2015. [23] J. Wu et al., “Optimal economic dispatch model based on risk management for wind-integrated power system, ”IET Gener. Transm. Distrib., vol. 9, no. 15, pp. 2152–2158, 2015. [24] C. Wan, Z. Xu, P. Pinson, Z. Y. Dong, and K. P. Wong, “Probabilistic forecasting of wind power generation using extreme learning machine, ”IEEE Trans. Powers Syst., vol. 29, no. 29, pp. 1033–1044, May 2014. [25] N. Nikmehr and S. N. Ravadanegh, “Optimal power dispatch of multimicrogrids at future smart distribution grids,” IEEE Trans. Smart Grid, vol. 6, no. 4, pp. 1648–1657, Jul. 2015. [26] Z. Yang, R. Wu, J. Yang, K. Long, and P. You, “Economical operation of microgrid with various devices via distributed optimization,” IEEE Trans. Smart Grid, vol. 7, no. 2, pp. 857–867, Mar. 2016. [27] Y. Zhang, N. Rahbari-Asr, J. Duan, and M.-Y. Chow, “Day-ahead smart grid cooperative distributed energy scheduling with renewable and storage integration,” IEEE Trans. Sustain. Energy, vol. 7, no. 4, pp. 1739–1747, Oct. 2016. [28] Z. Bao, Q. Zhou, Z. Yang, Q. Yang, L. Xu, and T. Wu, “A multi time-scale and multi energy-type coordinated microgrid scheduling solution - Part I: Model and methodology,” IEEE Trans. Powers Syst., vol. 30, no. 5, pp. 2257–2266, Sep. 2015. [29] Z. Bao, Q. Zhou, Z. Yang, Q. Yang, L. Xu, and T. Wu, “A multi time-scale and multi energy-type coordinated microgrid scheduling solution - Part II: Optimization algorithm and case studies,” IEEE Trans. Powers Syst., vol. 30, no. 5, pp. 2267–2277, Sep. 2015. [30] S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U. K.: Cambridge Univ. Press, 2004. [31] N. Rahbari-Asr, Y. Zhang, and M. Y. Chow, “Consensus-based distributed scheduling for cooperative operation of distributed energy resources and storage devices in smart grids,” IET Gener. Transm. Distrib, vol. 10, no. 5, pp. 1268–1277, 2016.

[32] A. Teixeira, E. Ghadimi, I. Shames, H. Sandberg, and M. Johansson, “The ADMM algorithm for distributed quadratic problems: Parameter selection and constraint preconditioning,” IEEE Trans. Signal Process., vol. 64, no. 2, pp. 290–305, Jan. 2016. [33] T. H. Chang, M. Hong, and X. Wang, “Multi-agent distributed optimization via inexact consensus ADMM,” IEEE Trans. Signal Process., vol. 63, no. 2, pp. 482–497, Jan. 2015. [34] J. Qiu et al., “Multi-stage flexible expansion co-planning under uncertainties in a combined electricity and gas market,” IEEE Trans. Power Syst., vol. 30, no. 4, pp. 2119–2129, Jul. 2015. [35] A. Shabanpour-Haghighi and A. R. Seifi, “An integrated steady-state operation assessment of electrical, natural gas, and district heating networks,” IEEE Trans. Power Syst., vol. 31, no. 5, pp. 3636–3647, Sep. 2016, doi: 10.1109/TPWRS.2015.2486819. [36] N. Rahbari-Asr and M. Y. Chow, “Cooperative distributed demand management for community charging of PHEV/PEVs based on KKT conditions and consensus networks,” IEEE Trans. Ind. Informat., vol. 10, no. 3, pp. 1907–1916, Aug. 2014.

Huaguang Zhang (M’03–SM’04–F’14) received the B.S. and M.S. degrees in control engineering from Northeast Dianli University of China, Jilin City, China, in 1982 and 1985, respectively. He received the Ph.D. degree in thermal power engineering and automation from Southeast University, Nanjing, China, in 1991. He joined the Department of Automatic Control, Northeastern University, Shenyang, China, in 1992, as a Postdoctoral Fellow for two years. Since 1994, he has been a Professor and the Head of the Institute of Electric Automation, School of Information Science and Engineering, Northeastern University, Shenyang, China. His main research interests include fuzzy control, stochastic system control, neural networks based control, nonlinear control, and their applications. He has authored and coauthored more than 280 journal and conference papers, six monographs, and co-invented 90 patents. Dr. Zhang is the E-letter Chair of IEEE CIS Society, the former Chair of the Adaptive Dynamic Programming & Reinforcement Learning Technical Committee on IEEE Computational Intelligence Society. He is an Associate Editor of Automatica, the IEEE TRANSACTIONS ON NEURAL NETWORKS, the IEEE TRANSACTIONS ON CYBERNETICS, and Neurocomputing, respectively. He was an Associate Editor of the IEEE TRANSACTIONS ON FUZZY SYSTEMS (2008–2013). He was awarded the Outstanding Youth Science Foundation Award from the National Natural Science Foundation Committee of China in 2003. He was named the Cheung Kong Scholar by the Education Ministry of China in 2005. He received the IEEE Transactions on Neural Networks 2012 Outstanding Paper Award.

Yushuai Li received the B. S. degree in electrical engineering and automation from Northeastern University, Shenyang, China, in 2014. He is currently working toward the Ph. D. degree in control theory and control engineering with the School of Information Science and Engineering, Northeastern University, China. His current research interests include distributed optimization and its applications in microgrids, smart grid, and energy Internet.


David Wenzhong Gao (S’00–M’02–SM’03) received the M.S. and Ph.D. degrees in electrical and computer engineering, specializing in electric power engineering, from Georgia Institute of Technology, Atlanta, GA, USA, in 1999 and 2002, respectively. He is currenly with the Department of Electrical and Computer Engineering, University of Denver, Denver, CO, USA. His current teaching and research interests include renewable energy and distributed generation, microgrid, smart grid, power system protection, power electronics applications in power systems, power system modeling and simulation, and hybrid electric propulsion systems. He is an Editor of the IEEE TRANSACTIONS ON SUSTAINABLE ENERGY. He is the General Chair for the 48th North American Power Symposium and the IEEE Symposium on Power Electronics and Machines in Wind Applications.

3097

Jianguo Zhou was born in Yunnan, China, in 1987. He received the B.S. degree in automation, and the M.S. degree in control theory and control engineering from the Northeastern University, Shenyang, China, in 2011 and 2013, respectively, where he is currently working toward the Ph.D. degree in control theory and control engineering. His current research interests include power electronics, hierarchical and distributed control, distributed optimization and multiagent systems with application to microgrids and energy Internet.