Energy Conversion and Management 174 (2018) 475–488 Contents lists available at ScienceDirect Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman Whole-day optimal operation of multiple combined heat and power systems by alternating direction method of multipliers and consensus theory T ⁎ Huynh Ngoc Trana, , Tatsuo Narikiyoa, Michihiro Kawanishia, Satoshi Kikuchib, Shozo Takabab a b Control System Laboratory, Toyota Technological Institute, Nagoya, Japan Higashifuji Technical Center, Toyota Motor Corporation, Shizuoka, Japan A R T I C LE I N FO A B S T R A C T Keywords: Combined heat and power system Fuel cell Fuel cell eﬃciency Distributed energy management system Alternating direction method of multipliers Multi-agent system This paper proposes a distributed energy management system-based algorithm to solve the optimal schedule of multiple combined heat and power systems in which exponential eﬃciencies of fuel cells are considered. To deal with the variation in fuel cell eﬃciencies, the proposed algorithm utilizes the alternating direction method of multipliers, consensus theory, and quadratic programming to solve the optimization problem with constant eﬃciencies. Then, the eﬃciency matching is checked iteratively to verify the obtained result. The optimization algorithm is constructed in a whole-day distributed form, such that all time slot variables can be solved simultaneously in a distributed way. With the help of this simultaneous solution method, operation of hot tanks and combined heat and power systems can be scheduled globally. The proposed algorithm is tested successfully with a test system having 4 combined heat and power systems, showing fast convergence in numerous simulations. 1. Introduction Recently, combined heat and power (CHP) systems or co-generation systems have been developed to increase the eﬃciency of energy systems. In CHP systems, the eﬃciency of energy conversion increases to over 80% as compared to an average of 30–35% for conventional fossil fuel ﬁred electricity generation systems [1]. In particular, Japan’s industry has developed CHP systems based on fuel cells (FCs) as one of the critical solutions to the high demand for clean energy [2]. The penetration of CHP systems into energy systems requires studies on system structure and energy management algorithm to optimally schedule CHP operation. Parameter sizing or structure design problems have been studied by various approaches in many papers for single integrated CHP system or fuel cell based systems. Ref. [3] proposes a soft-run strategy which takes battery size into consideration for real-time and multi-objective control algorithm design. In [4], a two-loop framework is used to solve a multiobjective optimization problem which considers both fuel economy and system durability. In [5], genetic algorithm is utilized to perform a multi-objective sizing optimization on a solar-hydrogen CHP system integrated with solar-thermal collectors. Multi-objective optimization approach is also utilized in [1] for designing integrated CHP systems for housing complexes in which the selection of technologies, the size of required units and operating modes of equipment are taking into ⁎ account. Ref. [6] proposes a new combined heat and power - heat pump integration system and utilizes genetic algorithm to optimize key design parameters, such as the prime mover capacity (PM), the outlet temperature from the heat pump and the decision value to run the PM. In [7], a new hub planning formulation is proposed to exploit assets of midsize/large CHPs. Optimal operation, planning, sizing and contingency operation of energy hub components are integrated and formulated as a single comprehensive mixed-integer linear programming problem. A probabilistic approach is proposed in [8] for optimal sizing of CHP systems in which the eﬀect of long-term uncertainty in energy demand is analytically investigated. As for the microgrid with penetration of CHP systems, [9] presents a power source sizing PSO-based strategy with integrated consideration of characteristics of distributed generations, energy storage and loads in an autonomy microgrid. As a basic problem for the parameter sizing/structure design problem, optimal scheduling or operational strategy for CHP system and microgrid with multiple CHP systems has been introduced in many papers. In [10], four CHP operational strategies are tested to choose their appropriate demand proﬁle. Ref. [11] proposes a strategy of sharing surplus heat and electricity produced by CHP plants in diﬀerent types of buildings, which could yield total energy cost savings of 1–9%. Dealing with uncertainties is one of important tasks of optimal scheduling problem, therefore it is taken into account in some papers. Ref. [12] proposes a robust optimization scheduling method to attenuate the Corresponding author. E-mail address: tranhuynhngoc@toyota-ti.ac.jp (H.N. Tran). https://doi.org/10.1016/j.enconman.2018.08.046 Received 11 March 2018; Received in revised form 9 August 2018; Accepted 12 August 2018 0196-8904/ © 2018 Elsevier Ltd. All rights reserved. Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Nomenclature η Fp ηBL ηex ηi,Fej ηi,Fh j μHT cp CkWhMJ G Fp Gi, j GiF, j GiBL ,j HTi, j m n P jB fuel cell eﬃciency of processing natural gas to hydrogen (assumed as 0.95) boiler eﬃciency (assumed as 0.99) hot water exchange eﬃciency (assumed as 0.95) electricity output eﬃciency of fuel cells at the ith building at time interval j heat recovery eﬃciency of fuel cells at the ith building at time interval j hot tank leakage ratio (assumed as 0.01) the speciﬁc heat of water (J/kg °C) conversion factor from kWh to MJ (3.6 MJ/kWh) processed gas energy of fuel cells at the ith building at time interval j (MJ) total gas energy consumed at the ith building in time interval j (MJ) natural gas energy consumed by fuel cells at the ith building at time interval j (MJ) gas energy consumed by the ith gas boiler at time interval j (MJ) hot tank level at the ending time of interval j at the ith building PiD, j PiF, j Pr G Pr je Tj Tpdem TpHT Tpin WiD, j WiBL ,j Wi,frj WiHT ,j Wito, j number of time intervals in a whole day number of buildings in the microgrid power bought from main grid to microgrid via substation at time interval j (kW) power demand at the ith building at time interval j (kW) power generated by fuel cells at the ith building at time interval j (kW) gas price (JPY/MJ) electricity price at time interval j (JPY/kWh) time length of time slot or time interval j (h) temperature of hot-water demand (assumed as 40 °C) temperature of hot water in hot tank (assumed as 70 °C) temperature of tap water (assumed as 20 °C) hot water demand at the ith building at time interval j (L) hot water from gas boiler at the ith building at time interval j (L) hot water sent from the ith building at time interval j to its exchanged building water discharged from hot tank at the ith building at time interval j (L) hot water received at the ith building at time interval j from its exchanged building from the demand side and large computational time. Therefore, in this paper, we introduce a distributed energy management system (EMS)based algorithm for scheduling a microgrid with penetrated multiple CHP systems. Our algorithm is based on the alternating direction method of multipliers (ADMM) [25], an analytical fast-convergence method which is suitable for large-scale optimization problems. There have been some papers applying ADMM approach for operational optimization of energy systems. Ref. [26] utilizes ADMM to develop a robust co-optimization scheduling model to study the coordinated optimal operation of the two energy systems. In [27], a day-ahead scheduling framework of integrated electricity and natural gas system is proposed at a distribution level based on the fast ADMM with restart algorithm considering demand-side response and uncertainties. As for the distributed EMS system, there have been some proposed algorithms [28,29] utilizing the ADMM approach and consensus theory for multiagent systems. These algorithms are proposed for application to optimal power ﬂow [29], demand response, or real-time pricing problems [28] in transmission power systems. Owing to the characteristics of such problems, their corresponding algorithms are based on time-slot sequential solving in which the whole-day schedule problem is divided into time interval schedule problems to be solved individually. Furthermore, the objective functions of these problems have a convex form, such that the ADMM approach can be applied easily. The microgrid with CHP systems presents some challenges for scheduling by these proposed algorithms. First, owing to the eﬃciency characteristics of the FC model, the relation between the FC input (energy of natural gas) and output (generated electricity and heat energy) is not linear. In this study, we assume an exponential eﬃciency curve for these relations based on a certain FC eﬃciency curve shape. Second, with the existence of energy storage (hot tank) in each CHP system, it is much better to solve the whole day/whole time problem than to sequentially solve time-interval problems. In this paper, we propose a distributed-EMSbased algorithm to overcome the two diﬃculties mentioned above. First, we introduce an iteration loop to change the nonconvex problem to a convex problem, such that the ADMM approach can be applied. The numerical results show that the added iteration loop has a good convergence characteristic. Second, we use the ADMM approach, consensus theory, and quadratic programming to solve our problem in a whole-day form, i.e., all time-interval problems are solved simultaneously. Furthermore, based on the case study results, we point out disturbance of uncertainty, and derive the day-ahead scheduling decision under two strategies including electrical load tracking and thermal load tracking. Ref. [13] presents a stochastic programming framework for conducting optimal 24-h scheduling of CHP-based microgrids considering demand response programs and uncertainties. A revision approach [2] is proposed to deal with the error in predicted demand for the optimal operation of fuel-cell-based residential energy systems. In [14], a colonial competitive algorithm is used in an optimal dispatching problem with uncertain electricity load. In terms of approach utilizing, various approaches including mix-integer linear/nonlinear programmings and meta-heuristic optimization have been developed for the optimization problem. Ref. [15] proposes a novel method based on information gap decision theory to evaluate a proﬁtable operation strategy for combined heat and power units in a liberalized electricity market. In [16], an optimal operation management system that integrates demand prediction into optimal scheduling by mixed-integer linear programming and real-time control for cogeneration units is proposed. Ref. [17] introduces a meta-heuristic optimization that could ﬁnd the global solution for a complex energy system. Ref. [18] builds a multiobjective optimization model for integrated CHP systems and compares its performance with the thermal-load-tracking control strategy. Ref. [19] introduces a multiobjective mixed-integer nonlinear programming model for the operational optimization of a large-scale combined cooling, heat, and power system. Through multiobjective optimization, Ref. [20] compares various system conﬁgurations that have on-site heat and electricity generation components. In [21], a solid oxide fuel cells based combined-cooling heating-and-power system design and operation optimization model has been developed using the Mixed Integer Non-linear Programming approach. A new isochronous governor control strategy for the CHP systems is proposed in [22] to provide zero-steady-state-error frequency regulation for a microgrid. Ref. [23] proposes a multiparty energy management framework with electricity and heat demand response for the CHP-Microgrid. Ref. [24] utilizes a coalitional game approach to propose a hybrid energy sharing framework with CHP system and photovoltaic (PV) prosumers. Although these mentioned studies cover a wide range of approaches for the structure design and optimal scheduling problems of CHP systems, all of their corresponding optimal algorithms/strategies are centralized, which may pose some diﬃculties for their application on the multiple CHP systems, owing to such factors as restricted information 476 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. eﬃciencies described in [30]. Owing to a long required startup time of FC units, we assume that all the FC units are always in the ON state. Each FC unit has a rating of 700 W output power and a minimum output power of 50 W. The number of FC units installed at each building can be diﬀerent, subject to the electricity and hot-water demands. This study considers only the economic dispatch of multiple CHP system, therefore transient and saturation features of FC are neglected. some features of the optimal scheduling of CHP systems at the residential level. The next section describes a microgrid model with penetrated residential CHP systems. Section 3 presents the problem formulation and a preliminary discussion on the solution method. Section 4 describes our proposed whole-day distributed ADMM approach, which is developed from a sequential distributed ADMM approach [28]. Section 5 demonstrates the whole proposed algorithm. Section 6 presents a 4-CHP test system and Section 7 shows the performance of the proposed algorithm on test system. The obtained results show that our proposed algorithm always converge in numerous simulations, which is better than performance of a particle swarm optimization approach on the same test system. In Section 8, conclusions drawn and future works are proposed. 2.2. Combined heat power model Fig. 2b depicts the CHP system [2] at one building and Fig. 3 describes the hot-water exchange between two CHP system of two buildings. In this study, we assume that hot water is exchanged between the demand sides but not between hot tanks. Based on the gas, power, and hot-water ﬂows in Figs. 2a and 3, the relations between electricity, gas, and hot-water energies in the CHP system i at time interval j are expressed as follows. 2. System modelling In this study, we consider a microgrid model, as shown in Fig. 1, with the following features: 2.2.1. Water demand The hot-water demand at building i is the sum of the hot-water volume from the mixer, boiler, and exchanged hot-water volume from/ to the building exchanging hot water with building i. • The grid is a lossless distribution network consisting of one substa- • • • tion and n buildings. The microgrid buys electrical power from the main grid via the substation. To maintain the electrical frequency in the microgrid, the microgrid has to be continuously connected to the main grid. Therefore, the bought power P jB at each time interval j is B required to be larger than a lower limitation of Pmin . Every building has its electricity and hot-water demand, which are predicted accurately. For simplicity’s sake, we do not consider other heating demands in this study. A CHP system is installed at each building to partially meet the building’s electricity demand and provide heat energy to the building’s hot-water system. The CHP model consists of FC units, mixers, a hot tank, and an auxiliary gas boiler, which is similar to that in [2]. A detailed description of the model is presented in the next subsections. To increase the usefulness of the CHP system, we consider the hotwater exchange between buildings in the microgrid model. It is assumed that each building exchanges hot water only with its nearest neighboring building. The optimization of the hot-water exchange structure is out of scope of this work. fr BL to WiD, j = WiHT , j Tp + Wi, j −Wi, j + Wi, j , where Tp = ηi,Fej = f Fe (PiF, j ) = ae−be e ηi,Fh j = f Fh (PiF, j ) = ah−bh e F Pimax ; −kh be = Fe ηmax −η0Fe PiF, j F Pimax , bh = 1−e−ke ; ae = be + η0Fe , Fh ηmax −η0Fh 1−e−kh , ah = bh + , and (3) is the hot-water volume from the 2.2.2. Hot-water exchange The hot-water exchange volume between buildings i and r is expressed by: ηex Wi,frj −Wrto, j = ηex Wrfr, j −Wito, j = 0, Wi,frj (4) Wito, j and are the hot-water volume sent from and received at where building i, respectively, and ηex (< 1) is the hot-water exchange eﬃciency. 2.2.3. Fuel cell and hot tank The natural gas energy G F , processed gas energy GpF , and FC output power P F at one time interval are related via the processing eﬃciency η Fp and FC electricity eﬃciency η Fe : A FC unit processes natural gas with an constant eﬃciency of η Fp and consumes the processed natural gas to provide both electricity and heat energy with eﬃciencies of η Fe and η Fh , respectively, to meet their respective demands. η Fe and η Fh are dependent on the FC operation mode and are not often given publicly. In [2], these two eﬃciencies are assumed to be proportional to the FC output power, beyond a minimum portion. However, based on the curve shape of these eﬃciencies in [2,10] and in an eﬀort to represent the FC features accurately, we assume these two eﬃciencies to be exponential functions of the FC output power, which are depicted in Fig. 2a and expressed by following formulas: −k e Tpdem − Tpin WiHT , j Tp mixer corresponding to the amount of water WiHT , j discharged by the hot tank. 2.1. Fuel cell model PiF, j TpHT − Tpin GiF, j = GiFp ,j = GiFp ,j η Fp , (5) PiF, j Tj CkWhMJ . ηi,Fej (6) The hot-water volume released from the FC to charge the hot tank is calculated from the FC heat energy output, as shown in (7). Fh F WiF, j (TpH T −Tpi n) cp = QiF, j = GiFp , j ηi, j = Pi, j Tj CkWhMJ ηi,Fh j η(Fe i, j ) The hot-tank current level is a result of the previous-interval level, as well as the hot-water charging and discharging, as follows: (1) HTi, j = (1−μHT ) HTi, j − 1 + WiF, j −WiHT ,j η0Fh , (2) where Fe ηmax = 0.465, η0Fe = 0.15, (7) Fh ke = 7; ηmax = 0.445; η0Fh = 0.25; kh = 6. Fig. 1. Microgrid system with multiple buildings. Fh Fe The values of ηmax and ηmax are based on the solid oxide fuel cell 477 (8) Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Fig. 2. (a) Assumption of FC eﬃciency. (b) Combined heat power system [2]. Fig. 3. Hot-water exchange between two buildings. where HTi, j − 1 is the hot-tank level at the ending time of the previous time interval, and HTimax is the maximum capacity (L) of the hot tank i. 2.2.4. Boiler The gas energy required by the boiler to heat a water volume of WiBL ,j at temperature Tpin to Tpdem is calculated by (9). GiBL ,j = ηex Wi,frj −Wrto, j = ηex Wrfr, j −Wito, j = 0 (16) fr D HT to WiBL , j = Wi, j −Wi, j Tp + Wi, j −Wi, j ⩾ 0 (17) fr D HT to GiBL , j = (Wi, j −Wi, j Tp + Wi, j −Wi, j ) (Tpdem−Tpin ) cp ηBL WiBL , j (Tpdem−Tpin ) cp ηBL (9) 0 ⩽ (1−μHT ) HTi, j − 1 + PiF, j (1−μHT ) HTi, j − 1 + PiF, j 2.2.5. Equality and inequality constraints of the CHP system The total gas energy of the CHP system is the sum of the gas energies of the boiler and the FC: F Gi, j = GiBL , j + Gi, j Tj CkWhMJ ηiFh ,j ηiFe , j (TpHT − Tpin) cp (18) −WiHT , j ⩽ HTimax , −WiHT , j = HTi0 , i = 1…(m−1) i=m (19) In (15)–(19), (10) PiF, j , WiHT ,j , of PiF, j Wi,frj , and Wito, j are scheduled variables, ηi,Fej and ηi,Fh are functions as shown in Fig. 2a, and WiD, j is the predicted j demand of the CHP system i at time interval j. From (3), we have: fr D HT to WiBL , j = Wi, j −Wi, j Tp + Wi, j −Wi, j Tj CkWhMJ ηiFh ,j ηiFe , j (TpHT − Tpin) cp BL ⩽ Gimax (11) 3. Problem formulation From (4) and (5), we have: GiF, j = The scheduling algorithm aims to minimize the total cost of the microgrid in a whole-day period: PiF, j Tj CkWhMJ ηi,Fej η Fp (12) m Minimizing F = Substituting (12) in (6), we have ∑ j=1 WiF, j = PiF, j Tj CkWhMJ ηi,Fh j ηi,Fej (TpHT −Tpin ) cp where Gi, j = (WiD, j −WiHT , j Tp + (13) (Tpdem−Tpin ) cp ηBL (14) As a result, the total gas energy consumed by the CHP system i at one time interval can be expressed by: fr to Gi, j = (WiD, j −WiHT , j Tp + Wi, j −Wi, j ) (Tpdem−Tpin ) cp ηBL + PiF Tj CkWhMJ ηi,Fej η Fp (Tp − Tp ) cp Wi,frj −Wito, j ) demη in BL (20) + PiF, j Tj Ck WhMJ Fp ηiFe ,j η is the energy of the gas consumed by the CHP system at the ith-building in the jth-time interval. The electricity price Pr je is given by the time-of-use price. The gas price Pr G is the same for all time intervals. In this study, we consider the short-term schedule/dispatch of CHP system, therefore the investment and maintenance costs does not aﬀect much to the result. Hence, these costs are neglected in the objective function. We will consider these costs in the next step of this study in which the parameter sizing problem is considered. From (9) and (11), we have fr D HT to GiBL , j = (Wi, j −Wi, j Tp + Wi, j −Wi, j ) n ⎛ e B ⎞ G ⎜Pr j P j Tj + Pr ∑ Gi, j ⎟ i=1 ⎝ ⎠ , (15) 3.1. Constraints Furthermore, the CHP variables need to satisfy the constraints of water exchange, boiler water volume, boiler limitation, and hot-tank level, as shown in (16)–(19), respectively: The microgrid power has to satisfy a lossless power balance between the supply and demand sides, i.e., 478 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. n P jB + ∑ n PiF, j = i ∑ PiD, j , 4. Whole-day distributed alternating direction method of multipliers j = 1…m. (21) i As mentioned above, the power bought from the main grid has a B lower limitation Pmin : We ﬁrst denote some vectors and matrixes to be used in the proposed algorithm. Let: B P jB ⩾ Pmin , P jF = [P1,Fj P2,Fj…PnF, j ]T , j = 1…m. (22) Wjfr = [W1,frj W2,frj…Wnfr, j ]T , Corresponding to the output range of the FC unit from 50 to 700 W, the CHP system at each building constrains its output power, depending on the number of FC units at each system. F F Pimax ⩾ PiF, j ⩾ Pimin , j = 1…m , Additionally, each CHP system needs to satisfy following constraints, which are derived from (16)–(19): (24) fr to D WiHT , j Tp−Wi, j + Wi, j ⩽ Wi, j , (25) fr to D WiHT , j Tp−Wi, j + Wi, j ⩾ Wi, j − 0 ⩽ (1−μHT ) HTi, j − 1 + PiF, j (1−μHT ) HTi, j − 1 + PiF, j BL ηBL Gimax , (Tpdem−Tpin ) cp Tj CkWhMJ ηiFh ,j ηiFe , j (TpHT − Tpin) cp Tj CkWhMJ ηiFh ,j ηiFe , j (TpHT − Tpin) cp (26) −WiHT , j ⩽ HTimax , −WiHT , j = HTi0 , i = 1…(m−1) i=m (27) 4.1. An approach of alternating direction method of multipliers 3.2. Preliminary discussion on problem-solving When the FC eﬃciencies are constant, the objective function F in Section 3 becomes convex. Therefore, it can be transformed into the following problem to ﬁt with the ADMM approach in [28,25]: Owing to the exponential characteristics of the FC eﬃciencies ηi,Fh j ηi,Fej , to to to T W to j = [W1, j W2, j…Wn, j ] B T 4n + 1, j = 1…m , is the vector of unPj ≜ [P jF W jHT Wjfr W to j Pj ] ∈ R known variables at time interval j, P ≜ [P1 P2…Pj…Pm ]T ∈ R(4n + 1) m is the vector of all unknown variables of the whole-day scheduling, P ≜ [P1 P2…Pj…Pm] ∈ R(4n + 1) × m is the matrix that consists of all elements of vector P but in a matrix form. From now on, when we mention elements of vector P , we will show them in the form of matrix P to avoid using excessive subscript levels, and similarly for elements of and matrix X, as mentioned below. vector X With these deﬁnitions of vector P and matrix P, variables fr to PiF, j , WiHT , j , Wi, j and Wi, j are indicated by Pi, j , Pi + n, j , Pi + 2n, j and Pi + 3n, j reB spectively. P j , which is the power purchased from the main grid at time interval j, is indicated by P4n + 1, j . ≜ [X1 X2 …Xj …Xm ]T ∈ R(4n + 1) m X and Let: Xj ∈ R 4n + 1, X ≜ [X1 X2 …Xj …Xm ] ∈ R(4n + 1) × m . and matrix X have the same structures as vector P and Vector X matrix P, respectively. We utilize these additional variables to satisfy the inequality constraints of the given problem, using an ADMM approach. Hence, the given problem can be separated into equality constraint and inequality constraint optimization problems, as shown in the following subsection. (23) ηex Wi,frj −Wrto, j = ηex Wrfr, j −Wito, j = 0, HT HT T W jHT = [W1,HT j W2, j …Wn, j ] PiF, j and the function of Gi, j with respect to in (20) is nonconvex. Therefore, the optimization problem is nonconvex, with both equality and inequality constraints, and it could not be directly solved by convex optimization methods. To overcome the nonconvex problem, we ﬁrst try to solve the optimization problem in which the FC eﬃciencies are constant values. Second, the FC eﬃciency matching is used to check whether the result is appropriate. The procedure is repeated until the result converges. This proposal is simple, but works eﬃciently on the given problem, as shown in Section 6, owing to the linearity of the optimization problem with constant FC eﬃciencies. As for the ﬁrst step in which FC eﬃciencies are constant values, we develop a whole-day distributed ADMM algorithm to solve the whole-day optimization problem in a distributed manner. Our approach is developed based on the sequential distributed ADMM algorithm presented in [28]. The sequential distributed ADMM algorithm divides the whole-day objective function into time-interval objective functions and then solves them sequentially. However, such an approach could not apply to our problem, in which the objective functions of all time intervals need to be solved at the same time to optimize the energy storage (hot tank) schedule. For instance, the hot-tank levels at a time interval, which is solved in a sequential algorithm, may limit the optimality of those in the next time intervals, owing to their relation in (19). In our approach, we form the X-update problem of the ADMM algorithm into a wholeday decentralized optimization problem, such that it can be solved in a decentralized manner, using quadratic programming. With this improvement, our approach is also valid for this class of optimization problem with additional types of energy storage, such as batteries or electric vehicles. In the next section, we describe the ﬁrst step of our proposed algorithm, which is the whole-day distributed ADMM algorithm, to solve the optimization problem with constant FC eﬃciencies. The whole algorithm is described in Section 5. m Minimizing FADMM = ∑ j=1 n ⎞ ⎛ e B G ⎜Pr j Tj P j + Pr ∑ Gi, j ⎟ + I1 (P ) + I2 (X ), i=1 ⎠ ⎝ (28) = 0, such that P−X where 0 if P ∈ Π1 I1 (P ) = ⎧ ⎨ ⎩∞ if P ∉ Π1 (29) Π1 = {P ∈ R(4n + 1) m : n P jB + ∑ n PiF, j = i=1 ∑ PiD, j , j = 1…m (30) i=1 ηex Wi,frj −Wrto, j = ηex Wrfr, j −Wito, j = 0, (i, r ) given} (31) Π1 is the set of P satisfying all equality constraints of the given problem. 0 if X ∈ Π2 I2 (X ) = ⎧ if X ∉ Π2 ∞ ⎨ ⎩ (32) ∈ R(4n + 1) m : Π2 = {X B X4n + 1, j ⩾ Pimin F Pimax ⩾ X i, j ⩾ (33) F Pimin Xi + n, j Tp−Xi + 2n, j + Xi + 3n, j ⩽ (34) WiD, j Xi + n, j Tp−Xi + 2n, j + Xi + 3n, j ⩾ WiD, j − 479 (35) BL ηBL Gimax (Tpdem−Tpin ) cp (36) Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. 0 ⩽ (1−μHT ) HTi, j − 1 + Xi, j Tj CkWhMJ η(Fh i, j ) η(Fe i, j ) (TpHT − Tpin ) cp i = 1…(m−1) HTi0 ⩽ (1−μHT ) HTi, j − 1 + Xi, j (44) can be re-written as −Xi + n, j ⩽ HTimax , k+1 k , uk ), = argminLρP (P, X ⎧P ⎪ n n ⎪ gj (P ) = ∑ Pi, j + P4n + 1, j− ∑ PiD, j = 0, j = 1…m ⎨ i=1 i=1 ⎪ ⎪ hir , j (P ) = ηex Pi + 2n, j−Pr + 3n, j = 0, (i, r ) given. ⎩ (37) Tj CkWhMJ η(Fh i, j ) η(Fe i, j ) (TpHT − Tpin ) cp −Xi + n, j ⩽ HTi0, i = m} (38) where satisfying all inequality constraints of the given Π2 is the set of X problem. Note that the second constraint in (27) involves two time interval variables, and thus it is transformed to last two inequality constraints in (38). With this transformation, all equality constraints of the given problem can be separated into time interval constraints. Based on the framework of the ADMM approach [28,30], we deﬁne the augmented Lagrange function as follows: m , η) = Lρ (P, X , u) = LρP (P, X (46) n k , uk ) s. t. g (P ) = 0, Minimize LρP (P, X j ρ ⎞ ⎛ e B G ⎜Pr j Tj P j + Pr ∑ Gi, j ⎟ + I1 (P ) + I2 (X ) + 2 ‖P j=1 ⎝ i=1 ⎠ ‖22 + +ηT (P−X ) −X (39) k be updated iteratively by: k k+1 k+1 uk + 1 = uk + P −X k ⩽ eprimal and ‖s k‖2 ⩽ edual , where k ‖r k‖2 = ‖P −X ‖2 , k −X k − 1)‖2 , ‖s k‖2 = ‖ρ (X eprimal = edual = ‖2 }, n ∊abs + ∊rel max{‖P ‖2 , ‖X k ∂gj (P ) ⎧ ∂LρP (P, X , u ) = λ , i = 1…n, j, k + 1 ∂P ⎪ ∂P4n + 1, j 4n + 1, j ⎪ k k ∂ g P ( ) ⎪ ∂LρP (P , X , u ) = λj, k + 1 j , i = 1…n, ∂Pi, j ∂Pi, j ⎪ ⎪ ∂LρP (P, X k , uk ) ∂gj (P ) ⎪ , i = 1…n, = λj, k + 1 ∂P i + n, j ⎪ ∂Pi + n, j (48) ⇔ ∂LρP (P, X k , uk ) ∂hir , j (P ) ⎨ , i = 1…n, = λjh,irk + 1 ∂P i + 2n, j ⎪ ∂Pi + 2n, j ⎪ ∂L (P, X k , uk ) ∂hir , j (P ) ⎪ ρP , i = 1…n = λjh,irk + 1 ∂P i + 3n, j ⎪ ∂Pi + 3n, j ⎪ g (P ) = 0 ⎪ j ⎪ hir , j (P ) = 0, (i, r ) given ⎩ (43) when ‖r k‖2 (48) λjh,irk + 1 where are the Lagrange scaled multipliers (at the iteration k + 1). Here, we see that (48) consists of msimilar, independent problems corresponding to m time intervals. It explains that which was mentioned in Section 3.2 about the constrains of variables P. In the next paragraph, without loss of generality, we show the solution of (48) for time interval j only. 1 The iteration is terminated k λjg, k + 1, where u = ρ η ∈ R(4n + 1) m is a scaled dual variable or scaled Lagrange , and u can multiplier. Following the procedure presented in [12], P, X (42) ⎧ ∂LρP (P, X , u ) = λ g ∂gj (P ) + ∑ λ hir ∂hir , j (P) j, k + 1 ∂P j, k + 1 ⎪ ∂P ∂P (i, r ) ⎨ ⎪ gj (P ) = 0, hir , j (P ) = 0, j = 1. .m ⎩ ρ ⎞ ⎛ e B G ⎜Pr j Tj P j + Pr ∑ Gi, j ⎟ + I1 (P ) + I2 (X ) + 2 ‖P −X j=1 ⎝ i=1 ⎠ Tη η + u‖22 − 2ρ (40) k + 1 = argminL (P k + 1, X , uk ) X ρ k n ∊abs + ∊rel ‖ρuk‖2 . By solving (41)–(43), the original optimization problem is separated into an equality constraint optimization problem of P and inequality . constraint optimization problem of X Note that constraints on variables P are global, equality but sepa rated into time interval constraints, whereas constraints on variables X are local, inequality but could not be decomposed into time interval constraints, owing to constraints (37) and (38). Therefore, we can adopt the average consensus protocol to solve variables P in a distributed manner, as the sequential distributed ADMM algorithm does. In con are easily divided into local variables at each agent trast, variables X (building), but need to be solved in a whole-day period. In the next subsections, detailed solution of P-update and X-update in (41) and (42), respectively, are given. k (49) Elaborating (49) give us (50)–(56) Pr je Tj + ρ (P4n + 1, j−X4kn + 1, j + u4kn + 1, j ) = λj, k + 1, Pr G Tj CkWhMJ η Fe η Fp −Pr GTp Pr G + ρ (Pi, j−Xik, j + uik, j ) = λj, k + 1, (Tpdem−Tpin ) cp ηBL (Tpdem−Tpin ) cp ηBL (50) i = 1…n, + ρ (Pi + n, j−Xik+ n, j + uik+ n, j ) = 0, (51) i = 1…n + ρ (Pi + 2n, j−Xi + 2n, j + ui + 2n, j ) = ηEx λjh,irk + 1, (52) i = 1…n (53) 4.2. P-update In this subsection, we present the solution of P in (41). From (41) and the deﬁnition of I1 (P ) , we have: k+1 k , uk ), = argminLρ (P, X ⎧ P ⎨ I1 (P ) = 0 ⎩ j = 1…m Because (47) is a convex optimization, it can be solved using the Karush-Kuhn-Turker condition, as follows: ∑ (41) hir , j (P ) = 0, (47) n k+1 k , uk ) P = argminLρ (P, X n ρ ⎞ ⎛ e 2 G ⎜Pr j Tj P1 + 4n, j + Pr ∑ Gi, j ⎟ + 2 ‖P−X + u‖2 . i=1 ⎠ ⎝ k+1 is the solution of the optimization From (45), we have that P problem where ρ > 0 is a scaled penalty parameter and η ∈ R(4n + 1) m is a Lagrange multiplier. Then, (39) can be rewritten as: , u) = Lρ (P, X m ∑ j=1 ∑ m (45) −Pr G (Tpdem−Tpin ) cp ηBL + ρ (Pi + 3n, j−Xi + 3n, j + ui + 3n, j ) = −λjh,irk + 1, i = 1…n (54) gj (P ) = (44) n ∑ i=1 480 n Pi, j + P4n + 1, j− ∑ PiD, j = 0, i=1 (55) Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. hir , j (P ) = ηex Pi + 2n, j−Pr + 3n, j = 0, (i, r ) given −1 ⎧ (1 + max (card (Ni ), card (Nl ))) , ⎪1− ∑ aik , l = i ail = ⎨ k ∈ Ni ⎪ 0, others ⎩ (56) From (50), (51) and (55), we have: Pr je Tj + nPr G Tj CkWhMJ ηFe ηFp n ⎞ ⎛ + ρ ⎜∑ (PiD, j −Xik, j + uik, j )−X4kn + 1, j + u4kn + 1, j ⎟ i = 1 ⎠ ⎝ Ni is the set of agents that communicate to agent i, which are called neighbor agents of agent i, card(Ni) is the number of neighbor agents of agent i, and Tc is the consensus calculating time. n+1 1 When Tc → ∞, as proved in [31], all x i, j → x∗, j = n + 1 ∑i = 1 Si, j . By that consensus, λj, k + 1 can be calculated at each agent by the following formula: = (n + 1) λj, k + 1 ⇒ λ j, k + 1 = Pr je Tj n+1 + nPr GTj CkWhMJ (n + 1) ηFe ηFp n + ρ ⎛ ∑ (PiD,j −Xik,j + uik,j)−X4kn+1,j n + 1 ⎜ i=1 ⎝ ⎞ + u4kn + 1, j ⎟ ⎠ λ j, k + 1 = (57) (TpHT −Tpin ) cp ρηBL , i = 1…n ⇒ + = (1 + ηEx ) Pr G (Tpdem−Tpin ) cp ηBL k+1 (1 + ηEx ηEx ) , i = 1…n, (62) , uk ) ,X (63) = argmin‖P ⎧X k+1 ⎨ I (X )=0 ⎩2 + uk‖22 , −X k+1 (64) Elaborating (64) and omitting superscripts for convenience, we obtain Eqs. (65)–(71): (59) = argmin‖P−X + u‖22 , X X4n + 1, j ⩾ F Pimax in (59) in a distributed way. It is easy to see that λjh,irk + 1 can be calculated at each building i and r if there exists a communication line between them. To calculate λj, k + 1 at each building and the substation in a distributed way, we apply an average consensus protocol, which is described in the next subsection. (65) B Pimin , (66) F Pimin (67) Xi + n, j Tp−Xi + 2n, j + Xi + 3n, j ⩽ WiD, j (68) Xi + n, j Tp−Xi + 2n, j + Xi + 3n, j ⩾ wi, j (69) ⩾ Xi, j ⩾ j j−t j ⎧ HTimax ⩾ αHT HTi0 + ∑t = 1 αHT (βit Xi, t −Xi + n, t ) ⎨ α j HTi0 + ∑ j α j − t (β Xi, t −Xi + n, t ) ⩾ 0, it t = 1 HT ⎩ HT 4.3. Average consensus j = 1…(m−1) (70) j j−t j ⎧ HTi0 ⩾ αHT HTi0 + ∑t = 1 αHT (βit Xi, t −Xi + n, t ) j ⎨ α j HTi0 + ∑ α j − t (β Xi, t −Xi + n, t ) ⩾ HTi0 it t = 1 HT ⎩ HT We note that in (57), the global term related to all agents (buildings and the substation) is a nonweighted summation term n ∑i = 1 (PiD, j −Xik, j + uik, j )−X4kn + 1, j + u4kn + 1, j . Such a summation term can be calculated in a distributed way, using average consensus theory. In this subsection, we apply an average consensus protocol to perform this calculation. This consensus protocol is also applied to other summation terms such as ‖r k‖22 , ‖s k‖22 , ‖P k‖22 , and ‖X k ‖22 in checking the stop criteria of the whole-day distributed ADMM algorithm. Let Si, j ≜ PiD, j −Xik, j + uik, j , i = 1…n , and Sn + 1, j ≜ −X4kn + 1, j + u4kn + 1, j . Then, (57) can be rewritten as: + + ρx∗, j. k+1 k+1 ⇒ ηBL (1 + ηEx ηEx ) n+1 (n + 1) ηFe ηFp = argminLρ (P ⎧X k+1 ⎨ I (X )=0 ⎩2 , i = 1…n Obviously, Pi + n, j is calculated locally in (58) without knowledge of other buildings. Moreover, if λj, k + 1 and λjh,irk + 1 are known, Pi, j, Pi + n, j, Pi + 2n, j, Pi + 3n, j and P4n + 1, j can also be calculated using only local information at each building i and the substation. Therefore, to solve variable P in a distributed manner, we calculate λj, k + 1 in (57) and λjh,irk + 1 λ j, k + 1 = nPr GTj CkWhMJ ) , we have: From (45) and deﬁnition of I2 (X (Tpdem − Tpin ) cp ρ (−ηEx Xi + 2n, j + ηEx ui + 2n, j + Xr + 3n, j −ur + 3n, j ) Pr je Tj + 4.4. X-update by quadratic programing (1 + ηEx ηEx ) λirh, j = ρ (−ηEx Xi + 2n, j + ηEx ui + 2n, j + Xr + 3n, j −ur + 3n, j ) λjh,irk + 1 n+1 (58) From (53), (54) and (56), we have: + (1 + ηEx ) Pr G Pr je Tj Hence, λj, k + 1 is calculated in a distributed way at each building and the substation, by applying the average consensus protocol to an extra local variable x i, j , as shown in (61) and (62). From (52), we have: Pi + n, j = Xik+ n, j −uik+ 2n, j + Pr G l ∈ Ni nPr GTj CkWhMJ (n + 1) ηFe ηFp ρ n+1 + ρ n+1 i=1 Si, j αHT = 1−μHT is the coefficient of all hot tanks, βit = Tt CkWhMJ ηFh (i, t ) ηFe (i, t )(TpHT − Tpin ) cp wi, j = WiD, j − (Tp , BL ηBL Gimax dem − Tpin ) cp . The optimization problem (65)–(71) could not be decomposed into time interval problems, owing to multiple-time-interval constraints (70) and (71). However, all constraints from (66)–(71) are local constraints at each building or the substation (each agent). Therefore, the X-update problem can be divided into local problems, as we show next. Denote: (60) n+1 ∑i = 1 Si, j directly, we apply multiagent Instead of calculating average consensus [31,32] to calculate this term in a distributed manner. Let Xid ≜ [Xi,1 Xi,2 …Xi, m …Xi + 3n,1 Xi + 3n,2 …Xi + 3n, m ]T ∈ R 4m , i = 1…n , which consists of all variables X of the i−th building in a whole-day period, Xnd+ 1 = [X 4n + 1,1 X 4n + 1,1…X 4n + 1, m ]T ∈ Rm , which consists of all variables X of the substation in a whole-day period, x i, j (0) = Si, j, x i, j (Tc + 1) = aii x i, j (Tc ) + ∑ ail xl, j (Tc ), (71) where n+1 ∑ j=m (61) where 481 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Pid ≜ [Pi,1 Pi,2…Pi, m…Pi + 3n,1 Pi + 3n,2…Pi + 3n, m ]T ∈ R 4m , i = 1…n , which consists of all variables P of the i-th building in a whole-day period, Pnd+ 1 = [P4n + 1,1 P4n + 1,1…P4n + 1, m ]T ∈ Rm , which consists of all variables P of the substation in a whole-day period, and uid ≜ [ui,1 ui,2…ui, m…ui + 3n,1 ui + 3n,2…ui + 3n, m ]T ∈ R 4m , i = 1…n , und+ 1 = [u4n + 1,1 u4n + 1,1…u4n + 1, m ]T ∈ Rm . average consensus protocol, which is similar to that described in Section 4.3. The latter is checked by ﬁrst transforming it to a logic value (0 or 1) at each agent and then applying the minimum consensus, as follows: Fh x i = ∧ (max‖EriFe , j (c + 1)‖ ⩽ e ηe , max‖Eri, j (c + 1)‖ ⩽ e ηh ), j = 1…m , x i (Tc + 1) = min(x i (Tc ), xl (Tc )), Based on the new variable assignment, the optimization problem (65)–(71) can be formed into whole-day period problems at each agent, as follows: where EriFe , j (c + 1) = d d d 2 d ⎧ Xi = argmin‖Pi −X + ui ‖2 ⎨ Ai Xid ⩽ bi , i = 1…n + 1, ⎩ Fe ηiFe , j (c + 1) − ηi, j (c ) EriFh , j (c + 1) = (72) ηiFe , j (c + 1) , Fh ηiFh , j (c + 1) − ηi, j (c ) ηiFh , j (c + 1) where and agents l are connected to agent i. It is easy to see that after limited times of calculation Tc , all x i = x∗ = minx i (0) , i.e., each agent can conﬁrm whether the stop criteria are satisﬁed at all agents. The stop criteria used in this study are set as follows: ρ = 0.06 ; ∊abs = ∊rel = 10−4; e ηe = e ηh = 0.01; Tc = 50 . I ⎡ Ai ⎤ Ai = ⎢ AiII ⎥ ∈ R(4m + 2m) × 4m , i = 1…n; An + 1 = Im ∈ Rm × m ⎥ ⎢ ⎢− AiII ⎥ ⎦ ⎣ B bi = [biI biII biIII ]T ∈ R(4m + 2m) , i = 1…n; bn + 1 = −Pmin 1m ∈ Rm ⎡ 1 ⎢− 1 I Ai = ⎢ 0 ⎢ 0 ⎣ 0 0 0 ⎤ 0 0 0 ⎥ 4m × 4m , Tp − 1 0 ⎥ ⊗ Im ∈ R Tp 1 − 1⎥ ⎦ 6. Case study We apply the proposed algorithm to a microgrid consisting of a substation and 4 buildings (agents) with CHP system at each building. Buildings 1 and 2 are residential ones and buildings 3 and 4 are complex (oﬃce and restaurant) ones. The given parameters of the tested microgrid are as follows: i = 1…n 0 α1 = [ αHT 0 … 0] ∈ R1 × m; j−1 j−2 0 αj = [ αHT … αHT αHT 0 … 0 ] ∈ R1 × m; 6.1. Power and water demands 0 Bi1 = [ αHT βi1 0 … 0] ∈ R1 × m; The power and hot-water demands of four buildings are shown in Fig. 5a and b, respectively. We estimated the power demand and hotwater demand curve shapes of the residential buildings based on the statistical data given in [33]. The power demand curve shape of a complex building is estimated from a report of the Japanese Agency for Natural Resources and Energy in 2011 [34]. Here, we assume that the power and hot-water demands are constant during each 1-h time interval. The ratio between the maximum power demand and the maximum hot-water demand is diﬀerent for each building. j−1 j−2 0 βi1 αHT βi2 … αHT βij 0 … 0 ⎤ ∈ R1 × m; Bij = ⎡ αHT ⎣ ⎦ AiII ⎡ B1 − α1 ⎤ ⋮ ⎢⋮ ⎥ ⎢ = Bj − αj 0m × 2m ⎥ ∈ Rm × 4m ⎢ ⎥ ⋮ ⎢⋮ ⎥ ⎢ ⎥ ⎣ Bm − αm ⎦ F ⎡ Pimax 1m ⎤ ⎢ − PF 1 ⎥ imin m ⎥ ∈ R 4m biI = ⎢ ⎢ {WiD, j }j ∈ {1.. m} ⎥ ⎥ ⎢ − ⎢ ⎦ ⎣ { wi, j }j ∈ {1.. m} ⎥ biII 6.2. Power limitations, hot-tank parameters, and gas-boiler capacity Table 1 lists the power limitations, hot-tank capacities and initial levels, and boiler capacities of the CHP systems at all buildings. Based on the power demand at each building, the numbers of FC units installed in each CHP system are 26, 18, 23, and 16 respectively. The hottank capacities are assigned as the same level as the water demand peak at each building, as shown in Table 1. Gas boilers with the capacities of 950 (L/h), 500 (L/h), 950 (L/h), and 500 (L/h) are installed at each CHP system, respectively, to ensure that the water demand at each building can be fulﬁlled without the FCs. ⎡ αHT HTi0 ⎤ ⎡ HTimax −αHT HTi0 ⎤ 2 2 ⎢ αHT ⎢ HTimax −αHT HTi0 ⎥ HTi0 ⎥ III m m =⎢ ⎥∈R . ⎥ ∈ R ; bi = ⎢ ⋮ ⋮ ⎥ ⎢ ⎥ ⎢ m m ⎢ ⎢ ⎦ ⎣ αHT HTi0−HTi0 ⎥ ⎦ ⎣ HTi0−αHT HTi0 ⎥ The optimization problem in (72) can be solved by quadratic programming, such as the interior-point method, which is embedded in Matlab/Simulink. 6.3. Energy price 5. Algorithm The electricity price is depicted in Fig. 6. It is a time-of-use two-level price [2], in which the higher level is applied from the 8-th time interval to the 22-nd time interval and the lower level is applied for the remaining time. The proposed algorithm is depicted in Fig. 4, in which two iteration loops are included. The ﬁrst, outer iteration loop (green) is the c-loop to deal with the variation of the FC eﬃciencies. The second, inner iteration loop (blue1) is exactly the whole-day distributed ADMM approach described in Section 4. The stop criteria of both the inner and outer loops can be checked in a distributed way. The former is checked by the 31.66 (JPY/kWh), Pr e (t ) = ⎧ ⎨ 13.10 (JPY/kWh), ⎩ 7 am ⩽ t ⩽ 10 pm, othertimes The gas price is 2.24 (JPY/MJ) [2], or 8.084 (JPY/kWh). At this gas price, a FC unit output of 1 kWh costs Pr Gas /(η Fpη Fe ) (JPY). Therefore, the cost range of 1 kWh from an FC unit is: 1 For interpretation of color in Fig. 4, the reader is referred to the web version of this article. 482 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Fig. 4. Algorithm ﬂowchart. Fig. 5. (a) Power demand vs. time slots. (b) Hot Water demand vs. time slots. Table 1 Power limitations, hot-tank and gas-boiler capacity at buildings. Buildings Type 1 Residential 2 Residential 3 Complex 4 Complex PDmax (kW) FC units 25.65 26 18.2 16.65 18 12.6 21.14 23 16.1 14.4 16 11.2 F Pmax (kW) F (kW) Pmin WDmax (L/h) HTmax (L) HT0 (L) BL Wmax (L/h) 18.3 (JPY) = 1.3 0.9 1.15 0.8 905.67 950 260 950 487.09 600 180 500 905.31 950 230 950 451.59 500 160 500 Fig. 6. Electricity price from main grid. Pr Gas Pr Gas Pr Gas ⩽ Fp Fe ⩽ Fp Fe = 34.45 (JPY), Fe η η η Fpηmax η ηmin Fe Fe = fFe (700 W) = 0.465 and ηmin = fFe (50 W) = 0.274 . It is easy where ηmax to see that when operating at a high percentage of output power, a FC could generate 1 kWh with a cost within the price range of the main grid. However, when an FC is utilized around its minimum output power, it costs more than the highest price of the main grid. Fig. 7. Microgrid networks. 6.4. Power, hot-water, and information networks Fig. 7 shows the microgrid power, hot-water, and information 483 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. utilized at maximum capacity, as shown in Fig. 11. Furthermore, hot water volume sent from building 1 to building 3 is larger than the volume building 1 received from building 3. The gas boilers are not needed all day. networks. In the hot-water network, building 1 exchanges hot water with building 3, and building 2 exchanges hot water with building 4. The information network is utilized for the distributed calculation (in consensus calculation). Consensus weightings ail and neighboring sets Ni can be easily calculated from the information network structure. Based on the convergence result of the consensus calculation presented in [28], we choose Tc = 50 (times) to reduce the computational cost while ensuring a good convergence in the consensus calculation. 7.1. Convergence features of the proposed method To provide information on the convergence feature of the proposed algorithm, we show the stop criteria parameters in Fig. 15, where 7. Test results n E (η Fe ) = The schedule results for total electricity are depicted in Fig. 8a. During the 1th –7th and 23th –24th time intervals, when the electricity price is low, all FC operate at the minimum output power to utilize the cheap power cost. During the 8th –22nd time intervals, when the electricity price is high, the total FC output power (Fig. 8a) shows its ﬂuctuation in a sawtooth shape. To maximize the FC eﬃciency, the optimal schedule tends to operate the FCs at high output in as many intervals as possible. However, owing to the limitation of hot-water demand and the constraint of the hot-tank level at the last time interval, the FCs cannot continuously operate at high output. Therefore, there must be time intervals during which the FCs operate at low output. These low outputs should be minimized to reduce the energy portion corresponding to low FC eﬃciency. Therefore, the FC power ﬂuctuates and forms a sawtooth shape, as mentioned above. Note that in this study FC capacities are at around 70% peak of power demand of their respective building. If FC capacities are designed at a lower percentage, FC utilization might be higher. However, design of optimal FC capacities is out of scope of this study. It is related to a structure, parameter sizing problem which is our future work. Fig. 8b shows the schedule results of the total hot-water variables. The total exchanged hot water in the sending direction is slightly larger than that in the receiving direction, owing to the heat losses in the exchange water pipes. These losses also lead to a slight diﬀerence between the total hot water from the mixers and the total hot-water demand during the 8th –23rd time intervals. The schedule results of buildings 1, 3, 2, and 4 are depicted in Figs. 9, 10, 12 and 13, respectively. Hot-tank schedules of buildings 1 and 3 are shown in Fig. 11 and hot-tank schedules of buildings 2 and 4 are shown in Fig. 14. Because electricity exchange is allowed in the microgrid, the FCs at buildings 1, 2, 3, and 4 sometimes operate at an output higher than power demand of their respective buildings, based on the optimal solution. Hot-water exchange occurs between buildings 1 and 3, and buildings 2 and 4. During the daytime, when the hot-water demands of buildings 1 and 2 are low, these buildings send hot water to buildings 3 and 4, respectively. On the contrary, at night, when the hotwater demands of buildings 1 and 2 are high, buildings 3 and 4 send hot water to buildings 1 and 2, respectively. Buildings 1 and 3 have the same hot-tank capacity, but building 3 requires a larger total hot-water volume than that of building 1, therefore hot-tank of building 3 is E (η Fh) = m ∑∑ i=1 j=1 n m ∑∑ i=1 ηi,Fej (c )−ηi,Fej (c−1) ηi,Fej (c ) Fh ηi,Fh j (c )−ηi, j (c−1) ηi,Fh j (c ) j=1 Maxe (η Fe ) = Max ηi,Fej (c )−ηi,Fej (c−1) Maxe (η Fh) = Max ηi,Fej (c ) Fh ηi,Fh j (c )−ηi, j (c−1) ηi,Fh j (c ) , , , , i = 1…n, j = 1…m , c = 1… 5, . The convergence of the c-loop is depicted in Fig. 15b, and the convergence of the whole-day distributed ADMM at the last c-loop iteration is depicted in Fig. 15a. It can be seen that the whole-day distributed ADMM approach converges at the last iteration of c-loop, as two conditions ‖r k‖2 ⩽ eprimal and ‖s k‖2 ⩽ edual are both satisﬁed. Furthermore, the general algorithm is terminated after ﬁve iterations of the c-loop, when the condition of matching eﬃciencies is satisﬁed. 7.2. Cost comparison To evaluate the proposed method strictly, a Particle swarm optimization (PSO) approach is also used to search for the optimal schedule of the test system. The approach utilized is based on the approach presented in [35,36] and was embedded into toolbox of Matlab. The detail algorithm and parameter setting of the approach is presented in [37]. To make inequality constraints (25)–(27) suitable for the PSO framework, we add hot-tank levels and hot water volume from boilers as new decision variables into the test system. With this modiﬁcation, all inequality constraints are transformed into lower and upper bounds of decision variables, which are available in the PSO framework, but the number of decision variables increases from 408 (in the proposed method) to 600. As for all equality constraints, square of their violations are added into the objective function (20) with a weighting factor of 1000. The algorithm is terminated if the relative change in the objective value over 20 iterations is less than a small tolerance (10−6 ). We conduct numerous simulations with the PSO approach, and see Fig. 8. (a) Total electricity schedule. (b) Total hot-water schedule. 484 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Fig. 9. (a) Electricity schedule, Bldg. 1. (b) Hot-water schedule, Bldg. 1. Fig. 10. (a) Electricity schedule, Bldg. 3. (b) Hot-water schedule, Bldg. 3. Fig. 11. (a) Hot-tank schedule, Bldg. 1. (b) Hot-tank schedule, Bldg. 3. Fig. 12. (a) Electricity schedule, Bldg. 2. (b) Hot-water schedule, Bldg. 2. function is corresponding to the 65-th simulation in which the total cost is 32,689 (JPY) and the total violation of equality constraints is 0.0313. It turns out that the optimal cost by PSO approach is larger than that of the proposed method, which is shown in Table 2. As the number of variables of the test system is quite large, PSO approach might have some diﬃculty to have a good convergence. Table 2 shows the cost of bought power from the main grid, the gas that all simulations do not converge to the same result. The convergence feature of these simulations is presented in Fig. 16. The PSO approach is applied with various values of swarm size from 50 to 600 in 120 simulations, as shown in Fig. 17. It can be seen that the total violation of equality constraints is quite large when the swarm size is less than 100, and reduced when the swarm size is in the range from 100 to 600. Among results of 120 simulations, the lowest objective 485 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Fig. 13. (a) Electricity schedule, Bldg. 4. (b) Hot-water schedule, Bldg. 4. Fig. 14. (a) Hot-tank schedule, Bldg. 2. (b) Hot-tank schedule, Bldg. 4. Fig. 15. (a) ADMM convergence in the terminal c-loop iteration. (b) Eﬃciency errors vs. c-loop iterations. Fig. 17. PSO approach - Swarm size vs. simulations. Table 2 Comparison of costs in cases of non-CHP and CHP scheduled by diﬀerent methods. Fig. 16. PSO approach - Result of 120 random initial scenarios. cost of FCs, and the gas cost of boilers, in three cases: non-CHP, CHP scheduled by the proposed method, and CHP scheduled by the PSO method. It is evident that the CHP schedule by the proposed method show signiﬁcant cost savings compared to the non-CHP and CHP-PSO 486 Unit (JPY) Cost of bought power Cost by FC Cost by gas boiler Total cost Total cost saved Non-CHP CHP-ADMM CHP-PSO 31,773 19,097 22,078 0 8638 10,611 3230 0 0 35,003 27,735 32,689 0% 20.76% 6.61% Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. Fig. 18. (a) Maximum eﬃciency errors vs. c-loop iterations in 250 scenarios. (b) Costs vs. 250 random initial scenarios. feature. To evaluate the proposed method, an PSO approach is also applied for searching for the optimal schedule of the CHP system. However, the PSO approach does not show a good convergence result as the proposed method. Owing to the structure of the whole-day optimization problem, the number of variables is quite large, and it might aﬀect the exploration and convergence features of the PSO method. case. 7.3. Convergence feature of random initial scenarios To verify the robustness of the proposed algorithm for the optimization problem, the algorithm is executed 250 times by randomly as u, i = 1…n, j = 1…m . signing the initial values of ηi,Fej , ηi,Fh j , P , X , and The convergences and costs of 250 random scenarios are shown in Fig. 18a and b, respectively. The maximum deviations in the bought electricity costs, FC gas costs, and total costs of these 250 initial scenarios are calculated as follows: dCostElectricity = dCostFC = max(ECostRand (j ) ) − min(ECostRand (j ) ) min(ECostRand (j ) ) max(FCostRand (j ) ) − min(FCostRand (j ) ) dCostTotal = min(FCostRand (j ) ) Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.enconman.2018.08.046. j = 1…250 References = 2.91%, max(TCostRand (j ) ) − min(TCostRand (j ) ) min(TCostRand (j ) ) = 2.47%, Appendix A. Supplementary material = 0.8%, j = 1…250 [1] Fuentes-Cortés LF, Ponce-Ortega JM, Nápoles-Rivera F, Serna-González M, ElHalwagi MM. Optimal design of integrated chp systems for housing complexes. Energy Convers Manage 2015;99:252–63. https://doi.org/10.1016/j.enconman. 2015.04.036. [2] Aki H, Wakui T, Yokoyama R. Development of an energy management system for optimal operation of fuel cell based residential energy systems. Int J Hydrogen Energy 2016;41(44):20314–25. https://doi.org/10.1016/j.ijhydene.2016.09.079. [3] Hu Z, Li J, Xu L, Song Z, Fang C, Ouyang M, et al. Multi-objective energy management optimization and parameter sizing for proton exchange membrane hybrid fuel cell vehicles. Energy Convers Manage 2016;129:108–21. https://doi.org/10. 1016/j.enconman.2016.09.082. [4] Xu L, Mueller CD, Li J, Ouyang M, Hu Z. Multi-objective component sizing based on optimal energy management strategy of fuel cell electric vehicles. Appl Energy 2015;157:664–74. https://doi.org/10.1016/j.apenergy.2015.02.017. [5] Assaf J, Shabani B. Multi-objective sizing optimisation of a solar-thermal system integrated with a solar-hydrogen combined heat and power system, using genetic algorithm. Energy Convers Manage 2018;164:518–32. https://doi.org/10.1016/j. enconman.2018.03.026. [6] Li H, Kang S, Lu L, Liu L, Zhang X, Zhang G. Optimal design and analysis of a new chp-hp integrated system. Energy Convers Manage 2017;146:217–27. https://doi. org/10.1016/j.enconman.2017.05.034. [7] Moradi S, Ghaﬀarpour R, Ranjbar AM, Mozaﬀari B. Optimal integrated sizing and planning of hubs with midsize/large chp units considering reliability of supply. Energy Convers Manage 2017;148:974–92. https://doi.org/10.1016/j.enconman. 2017.06.008. [8] Urbanucci L, Testi D. Optimal integrated sizing and operation of a chp system with Monte Carlo risk analysis for long-term uncertainty in energy demands. Energy Convers Manage 2018;157:307–16. https://doi.org/10.1016/j.enconman.2017.12. 008. [9] Liu Z, Chen Y, Zhuo R, Jia H. Energy storage capacity optimization for autonomy microgrid considering chp and ev scheduling. Appl Energy 2018;210:1113–25. https://doi.org/10.1016/j.apenergy.2017.07.002. [10] Aki H, et al. Operational strategies of networked fuel cells in residential homes. IEEE Trans Power Syst 2016;21(3):1405–14. https://doi.org/10.1109/TPWRS. 2006.879270. [11] Hirvonen J, Kayo G, Hasan A, Siren K. Local sharing of cogeneration energy through individually prioritized controls for increased on-site energy utilization. Appl Energy 2014;135:350–63. https://doi.org/10.1016/j.apenergy.2014.08.090. [12] Wang L, Li Q, Sun M, Wang G. Robust optimisation scheduling of cchp systems with multi-energy based on minimax regret criterion. IET Gener, Transm Distrib 2016;10(9):2194–201. https://doi.org/10.1049/iet-gtd.2015.1344. [13] Alipour M, Mohammadi-Ivatloo B, Zare K. Stochastic scheduling of renewable and chp-based microgrids. IEEE Trans Ind Inf 2015;11(5):1049–58. https://doi.org/10. 1109/TII.2015.2462296. [14] Karami H, Sanjari M, Gooi H, Gharehpetian G, Guerrero J. Stochastic analysis of j = 1…250 Based on Fig. 18a, it can be seen that the proposed algorithm always converges in all 250 random initial scenarios after only 14 or fewer iterations of the c-loop. Combining the cost results in Fig. 18b and the maximum cost deviations, these are convincing evidence that the proposed algorithm achieves good convergence on the optimization problem. 8. Conclusion and future work In this study, we address an optimization problem of a microgrid consisting of multiple CHP systems. A model of the CHP system is given in detail. Furthermore, we consider the exponential curve of the FC eﬃciency, based on some eﬃciency curve shapes. The optimization problem is a whole-day cost optimization with a nonconvex objective function, equality and inequality constraints. To render the objective function linear, an approximation is applied in which the FC eﬃciencies are assumed constant, and an iteration c-loop is introduced to update them, corresponding to the eﬃciency curves. Next, the proposed wholeday distributed ADMM is used to solve the problem of the linear objective function. The test result shows that the proposed method lead to a reasonable scheduling of CHP systems in a whole-day operation. With the combination of ADMM method and the consensus protocol, any buildings and the substation can schedule their variables without acquiring full information of the microgrid. With the use of ADMM method, all constraints of the optimization problem are satisﬁed. Therefore, the scheduling result is applicable. Because the ADMM method has been proved always converge for the convex optimization problem, the weakness of the proposed method is the iterative approach. However, the test result shows that the iterative approach converges very fast to the same cost result in many simulations with random initial variables. In the future, we will present an analytical result that could prove this convergence 487 Energy Conversion and Management 174 (2018) 475–488 H.N. Tran et al. [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] residential micro combined heat and power system. Energy Convers Manage 2017;138:190–8. https://doi.org/10.1016/j.enconman.2017.01.073. Aghaei J, Agelidis VG, Charwand M, Raeisi F, Ahmadi A, Nezhad AE, et al. Optimal robust unit commitment of chp plants in electricity markets using information gap decision theory. IEEE Trans Smart Grid 2017;8(5):2296–304. https://doi.org/10. 1109/TSG.2016.2521685. Wakui T, Sawada K, Kawayoshi H, Yokoyama R, Iitaka H, Aki H. Optimal operations management of residential energy supply networks with power and heat interchanges. Energy Build 2017;151:167–86. https://doi.org/10.1016/j.enbuild.2017. 06.041. Andrea LF, Stefano U. Meta-heuristic optimization for a high-detail smart management of complex energy systems. Energy Convers Manage 2018;160:341–53. https://doi.org/10.1016/j.enconman.2018.01.035. Ermanno LC, Davide B, Francesco D, Corrado S. Future distributed generation: an operational multi-objective optimization model for integrated small scale urban electrical, thermal and gas grids. Energy Convers Manage 2017;143:384–9. https:// doi.org/10.1016/j.enconman.2017.04.006. Zhu Q, Luo X, Zhang B, Chen Y, et al. Mathematical modelling and optimization of a large-scale combined cooling, heat, and power system that incorporates unit changeover and time-of-use electricity price. Energy Convers Manage 2017;133:385–98. https://doi.org/10.1016/j.enconman.2016.10.056. Delgado BM, Cao S, Hasan A, Sirén K. Multiobjective optimization for lifecycle cost, carbon dioxide emissions and exergy of residential heat and electricity prosumers. Energy Convers Manage 2017;154:455–69. https://doi.org/10.1016/j.enconman. 2017.11.037. Jing R, Wang M, Wang W, Brandon N, Li N, Chen J, et al. Economic and environmental multi-optimal design and dispatch of solid oxide fuel cell based cchp system. Energy Convers Manage 2017;154:365–79. https://doi.org/10.1016/j.enconman. 2017.11.035. Sun T, Lu J, Li Z, Lubkeman D, Lu N. Modeling combined heat and power systems for microgrid applications. IEEE Trans Smart Grid. https://doi.org/10.1109/TSG. 2017.2652723. Liu N, He L, Yu X, Ma L. Multiparty energy management for grid-connected microgrids with heat- and electricity-coupled demand response. IEEE Trans Ind Inf 2018;14(5):1887–97. https://doi.org/10.1109/TII.2017.2757443. Liu N, Wang J, Yu X, Ma L. Hybrid energy sharing for smart building cluster with [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] 488 chp system and pv prosumers: a coalitional game approach. IEEE Access 2018;6:34098–108. https://doi.org/10.1109/ACCESS.2018.2847023. Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found Trends Mach Learn 2011;3(1):1–122. He C, Wu L, Liu T, Shahidehpour M. Robust co-optimization scheduling of electricity and natural gas systems via admm. IEEE Trans Sustain Energy 2017;8(2):658–70. https://doi.org/10.1109/TSTE.2016.2615104. Chen J, Zhang W, Zhang Y, Bao G. Day-ahead scheduling of distribution level integrated electricity and natural gas system based on fast-admm with restart algorithm. IEEE Access 2018;6:17557–69. https://doi.org/10.1109/ACCESS.2018. 2818756. Nguyen DH, Narikiyo T, Kawanishi M. Optimal demand respond and real-time pricing by a sequential distributed consensus-based admm approach. IEEE Trans Smart Grid. https://doi.org/10.1109/TSG.2017.2676179. Erseghe T. Distributed optimal power ﬂow using admm. IEEE Trans Power Syst 2014;29(5):2370–80. https://doi.org/10.1109/TPWRS.2014.2306495. Kasuh T. Why does japan believe in domestic fuel cell? Adaptation to european market. In: European gas technology conference. Olfati-Saber R, Fax J, Murray RM. Consensus and cooperation in networked multiagent system. Proc IEEE 2007 2007;95(1):215–33. https://doi.org/10.1109/ JPROC.2006.887293. Xiao L, Boyd S, Lall S. A scheme for robust distributed sensor fusion based on average consensus. In: Proc. int. conf. information processing in sensor networks 2005; 2007. p. 63–70. https://doi.org/10.1109/IPSN.2005.1440896. Aki H, et al. Analysis of measured data on energy demand and activity patterns in residential dwellings in Japan. IEEJ Trans Electr Electron Eng 2018;13:157–67. https://doi.org/10.1002/tee.22509. METI J. A report of japanese agency for natural resources and energy; 2011. Available at < http://www.meti.go.jp/setsuden/20110513taisaku/16.pdf > . Mezura-Montes, Coello CA. Constraint-handling in nature-inspired numerical optimization: past, present and future. Swarm Evol Comput 2011:173–94. Pedersen ME. Good Parameters for Particle Swarm Optimization. Luxembourg: Hvass Laboratories; 2010. Mathworks. Particle swarm optimization algorithm; 2018. < https://www. mathworks.com/help/gads/particle-swarm-optimization-algorithm.html > .

1/--страниц