—Participation of distribution systems in frequency regulation has become an important factor for grid operation af- ter the integration of distributed energy resources (DERs). While available reserves from a single distribution system may not be sufﬁcient for frequency regulation, aggregated reserves from several distribution systems can provide frequency regulation at grid-scale. This paper proposes a cooperative game theory-based approach for secondary frequency regulation in distribution systems consisting of distributed energy resources (DERs). A two- stage strategy is proposed to effectively and precisely determine a change in active power set points of DERs following a command from system operators. In the ﬁrst stage, the value or worth of each DER and their coalitions are determined using initial rates of change of frequency (ROCOF) for each DER/coalition. In the second stage, the Shapley value, one of the solution concepts of cooperative game theory, is used to determine changes in active power set points of DERs. The proposed method is implemented on several distribution systems including a modiﬁed IEEE 13- node system and modiﬁed 33-node distribution system. Index Terms —Cooperative game theory, distributed energy resources, frequency regulation, Shapley Value. NOMENCLATURE N set of players of a cooperative game V;W characteristic functions S a coalition that is subset of N 2Npossible set of coalitions payoff vector fjg singleton or unit set of player j Snfjg coalition set without player j j Shapley value of player j n number of players Hi inertia constant of ith generator fi frequency of ith generator fc frequency of equivalent center of inertia fn rated system frequency k Shapley value of kth DER DFk distribution factor of kth DER Pd total disturbance power Pref reference power command received from the distribution system operator (DSO) Pref kchange in active power set-point of kth DERI. I NTRODUCTION The deregulation of electric utilities and large-scale inte- gration of distributed energy resources (DERs), such as pho- tovoltaic systems, wind turbines, battery storage, and electric vehicles, have contributed to frequency regulation challenges in modern power systems. When power imbalance between the generation and load plus losses occurs in power systems (e.g., due to sudden changes in system operating conditions, islanding of microgrids, large load changes, and generating unit tripping), frequency regulating schemes come into play to compensate for frequency deviations. Frequency control in power systems can be categorized into primary, secondary, and tertiary control. The primary frequency control is a local primary line of action which acts immediately (within millisec- onds) following a disturbance. If the primary frequency control cannot completely compensate for frequency deviations, the secondary frequency control adjusts the output power of partic- ipating resources (generating units, DERs, demand response, etc.) to restore the frequency to its nominal value. However, the resulting operating points of participating resources due to the secondary frequency control may not be economically optimal. In this context, tertiary frequency control schemes come into action to minimize cost of generations. Conventional generating units in isolated power systems are usually equipped with an integral controller to change refer- ence settings of turbine governors for secondary frequency regulation [1]. Also, conventional power systems consist of slow centralized secondary frequency controllers to restore the system frequency when the deviation is larger than a speciﬁed tolerance [2]. However, the secondary frequency regulation applied for conventional power systems may not be directly applicable for distribution systems consisting of DERs. Therefore, a reliable and quick methodology that can determine changes in active power set points of DERs is pivotal for the participation of distribution systems in the secondary frequency regulation. Various methods have been presented in the literature for the secondary frequency regulation and control in transmis- sion and distribution systems. A fully decentralized integral controller has been proposed in [3] for a secondary frequency control. A time-decoupled distributed secondary controller has been proposed in [4], which can be implemented for both frequency and voltage restoration. A hybrid control strategy that can be used for primary and secondary frequency regu- lation has been proposed in [5] for dynamic demand control. A Lagrange multiplier-based approach has been proposed in [6] for under-frequency load shedding, which can be used for frequency regulation. A distributed secondary frequency control has been proposed in [7] to coordinate between generating units, storage devices, and other energy sources. A distribution-level phasor measurement unit (D-PMU)-based secondary frequency control and response for islanded mi- crogrids has been proposed in [2]. A distributed secondary voltage/frequency control scheme that ensures reactive power sharing in addition to voltage and frequency regulation has been proposed in [8]. The distributed control scheme proposed in [8] can ensure power sharing and frequency regulation, but there is a necessity of communication between DER inverters, which can impair the ﬂexibility of microgrids [9]. Also, developed and employed secondary frequency regulation and control methods for transmission systems and microgrids cannot be directly applied at the distribution system level to participate in frequency regulation; enabling distribution sys- tems with high penetration of DERs to participate in secondary frequency regulation at the grid-scale is still a challenge. Cooperative game theory-based approaches have been suc- cessfully applied in various ﬁelds of power systems. A co- operative game theory-based approach has been implemented in [10] for loss reduction allocation of distributed generations using Shapley values. A cooperative game theory-based ap- proach has been proposed in [11] for under frequency load shedding control. A cooperative game theory-based approach for computing participation factors of distributed slack buses has been proposed in [12]. For secondary frequency regulation, the cooperative game theoretic approaches based on the Shap- ley value can ensure that the total disturbance power is fairly distributed among different DERs taking into account their marginal contributions. Therefore, this paper investigates the cooperative game theoretic approach for secondary frequency regulation in distribution systems. This paper proposes a cooperative game theory-based sec- ondary frequency regulation scheme that determines active power set points of DERs in a distribution system. The proposed approach determines active power set points of par- ticipating DERs based on their Shapley values. The proposed approach is implemented in two-stages. In the ﬁrst stage, char- acteristic functions of DERs and their coalitions are computed using initial rate of change of frequency (ROCOF). In the second stage, the disturbance or reference power is distributed among different participating DERs using Shapley values. The proposed secondary frequency regulation is implemented on the modiﬁed IEEE 13-node system and the modiﬁed 33-node distribution system. The remainder of the paper is organized as follows. Section II describes cooperative game theory including the Shapley value. Section III presents the formulation of the cooperative game model. Section IV explains the proposed approach for secondary frequency regulation. Section V describes casestudies on the modiﬁed IEEE 13-node and the modiﬁed 33- node distribution systems. Section VI provides concluding remarks. II. C OOPERATIVE GAME THEORY AND SHAPLEY VALUE In game theory, a game can be categorized into (a) cooper- ative game and (b) non-cooperative game. In non-cooperative games, there is no coalition between players while in cooper- ative game theory there is a cooperation or coalition between players. In cooperative games, each player can form alliances with other players to maximize its incentives. Since coali- tions among players are formed to increase their individual incentives, a coalition must always result in equal or greater incentives than individual player’s incentives. In this paper, we use cooperative games with the goal of maximizing beneﬁts of the grid and DER owners; therefore, moving forward we only discuss and use cooperative games. In each cooperative game, there are three components as follows [13]: 1) A ﬁnite set of players, denoted by N. 2) A set of coalitions that a player can form. 3) Preference of each player over all possible coalitions. Value or worth of each coalition in the cooperative game is represented using a characteristic function. In other words, a characteristic function is the total utility of all members of the coalition. The characteristic function can be represented asV(S), whereSis a coalition. The characteristic function is a real-valued function (i.e., V(S) : 2N!R) with an empty set having zero value (i.e., V( ) = 0 ). The total payoff or incentive is distributed among the players using solution concepts including the Shapley value, the Nucleolus, and Nash-bargaining solution. A. The Core of a Cooperative Game In game theory, the core refers to the set of feasible allocations that cannot be further improved through any other coalitions. Generally, outcomes of a cooperative game are expressed as n-tuples of utility: =f i:i2Ng, called payoff vectors that are measured in some common monetary unit [14]. The core is the set of imputations under which all sets of coalitions have values less than or equal to the sum of its members’ payoffs. Thus, is core if and only if, :eSV(S);8SN (1) :eN=V(N) (2) whereeSdenotes the n-vector having eS i= 1 ifi2Sand eS i= 0ifi2NS. Equations (1) and (2) denote, respectively, stability and efﬁciency criteria. B. The Shapley Value The Shapley value is one of the solution concepts of cooper- ative game theory. The Shapley value assigns a unique payoff vector that is efﬁcient, symmetric, and satisﬁes monotonicity. The Shapley value allocates the payoffs in such a way that is fair for cooperative solutions. The Shapley value of a cooperative game is given as follows [15]. G GGGGrid G Distributed energy resource (DER) Communication between DERs and substationSubstationDist. System Operator (DSO) Communication between DSO and substationFig. 1. Layout of the proposed cooperative game theoretic approach j(V) =X S22N;j2S(jSj1)!(njSj)! n![V(S)V(Snfjg)](3) wheren=jNjis the total number of players. The Shapley value satisﬁes the following axioms: 1)Efﬁciency: The efﬁciency axiom states that the sum of Shapley values of all players is equal to the worth of grand coalition so that the total gain is distributed among the players, i.e.,P j2N j(V) =V(N). 2)Individual Rationality: This axiom states that the Shap- ley value of each player should be greater than or equal to its individual worth, i.e., j(V)V(fjg);8j2N. 3)Symmetry: Ifjandkare such that V(S[fjg) =V(S[ fkg)for every coalition Snot containing jandk, then j(V) = k(V). 4)Dummy Axiom: Ifjis such that V(S) =V(S[fjg) for every coalition Snot containing j, then j(V) = 0 . 5)Additivity: IfVandWare characteristic functions, then (V+W) = (V) + (W). III. C OOPERATIVE GAME MODEL A cooperative game model is formulated to determine the Shapley value of each distributed energy resource (DER) participating in secondary frequency regulation. Fig. 1 shows the layout of the proposed cooperative game theory-based approach for secondary frequency regulation of an active dis- tribution system. In this paper, the task of secondary frequency regulation is regarded as a game and participating DERs being players of the game. Since participating DERs act in a cooperative manner to compensate frequency deviation after receiving a command from system operators and restore the frequency to its nominal value, the game becomes a cooper- ative game. As explained in Section II, a cooperative game should have characteristic functions in addition to a player set. Therefore, initial rate of change of frequency (ROCOF) is utilized to determine characteristic functions, which serves as a worth or value of each DER and its coalitions. The DER cooperative model formulation for the proposed approach can be explained as follows.1) Collect system data including substation data, line data, load data, etc., which serve as input to the cooperative game model. 2) Generate the list of all possible coalitions of DERs. For example, if three DERs ( DER 1,DER 2, andDER 3) are participating in secondary frequency regulation, the set of all possible coalitions, denoted by 2N, is as follows, where denotes an empty set. 2N=f ;fDER 1g;fDER 2g;fDER 3g; fDER 1;DER 2g;fDER 1;DER 3g; fDER 2;DER 3g;fDER 1;DER 2;DER 3gg. 3) For each DER and its possible coalitions, perform tran- sient analysis to determine the initial ROCOF. The formulated cooperative game model, which uses the procedure explained in the above steps, is essential for sec- ondary frequency regulation. IV. T HEPROPOSED COOPERATIVE GAME THEORETIC APPROACH In this section, an approach for the determination of dis- turbance power is explained and an analytical expression for change in active power set point of each participating DER is developed using a cooperative game theoretic approach. When a disturbance occurs in power systems, it gets re- ﬂected in the system frequency. System operators can calculate disturbance power, Pd, based on the initial ROCOF referred to the equivalent center of inertia (COI). ROCOF can be determined using different methods. Here we explain it using the swing equation of the ith machine with inertia constant Hias follows [16], [17]. (2Hi=fn)(dfi=dt) = pi (4) Pd=MX i=1pi= 2MX i=1Hi(dfc=dt)=fn (5) wherefnis the rated system frequency, fiis the frequency of theithmachine, piis the power deviation of the ithmachine, Mis the total number of machines, and fcis the frequency of the equivalent COI. The frequency of the equivalent COI describes the average system frequency at the time of electromechanical transients when individual machine frequencies are not the same, which can be calculated using (6). fc=MX i=1Hifi.MX i=1Hi (6) A distribution system may not be capable of contributing the total disturbance power, Pd, obtained using (5). System operators decide on how to distribute power imbalance among available resources such as power pants, control areas, and participating active distribution systems. If the system operator commands a particular distribution system to contribute active power equal to Pref, then it can be distributed among n Start Read system data Select a coalition of DERs from the list of possible coalitions Compute characteristic function of the coalition using initial ROCOF StopCompute the disturbance power P d using initial ROCOF after disturbanceEnumerate all possible coalitions of DERs participating in frequency regulation Compute ΔPref of each participating DER All coalitions evaluated? YesNo Compute Shapley values corresponding to each participating DERCompute distribution factor (DF) of each participating DERConnect only the DERs in that particular coalitionFig. 2. Flow chart of the proposed approach participating DERs based on the distribution factor, DFk, deﬁned for the kthDER as follows. DFk= k.nX k=1 k (7) where kis the Shapley value of the kthDER. The proposed framework is for grid connected sources. If a distribution system is self-sufﬁcient and is being operated in an islanded mode, then its controllers have to determine the reference power, Pref, locally. Now, since the sum of distribution factors of all participating DERs should be equal to unity, the change in active power set point of the kthDER is calculated as follows. Pref k=DFkPref (8) The proposed approach or the solution algorithm to de- termine the change in active power set points of DERs for secondary frequency regulation can be summarized as follows. 1) Provide system data related to lines, loads, transformers, and generators. 2) Perform secondary frequency regulation based on a command received from system operators. In case of an islanded distribution system, calculate power imbalance locally. Local calculation of frequency deviation requires local frequency measurement devices; however, deter- mining types and functionalities of local measurement devices is out of the scope of this work. G 646 645 632 633 634650 692 675 611 684 652671 680 G GDER 2 DER 1DER 3Fig. 3. The modiﬁed IEEE 13-node distribution system 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 172223 24 1819 20 21 26 27 28 29 30 31 32S/S 25 33G GGG DER 1 DER 2DER 4DER 3 Fig. 4. The modiﬁed 33-node distribution system 3) Enumerate all possible coalitions of participating DERs and compute their characteristic functions. 4) Compute Shapley value of each DER based on (3) and the respective change in active power set point using (7) and (8). The ﬂow chart of the proposed cooperative game theoretic approach is shown in Fig. 2. V. C ASE STUDIES AND DISCUSSIONS This section presents the implementation of the proposed cooperative game theoretic approach on a modiﬁed IEEE 13- node system and modiﬁed 33-node distribution system. Since we do not have ROCOF measurement data, a synchronous generator is added at the substation node to emulate the power grid on which the substation is connected and initial ROCOF is calculated based on the swing equation assuming constant inertia. Note that in practical applications, the ROCOF mea- sured using frequency measuring devices must be used. A. The Modiﬁed IEEE 13-node System The IEEE 13-node system is a 4:16kV distribution test system characterized by having overhead and underground lines, transformers, a voltage regulator, shunt capacitor banks, and unbalanced loading with constant current, power, and impedance models. The total real and reactive loads of this system are, respectively, 3577 kW and 1725 kV Ar. For the detailed data of the IEEE 13-node system, the readers are referred to reference [18]. In this paper, this system has been modiﬁed by including three DERs of capacities 200 kW, 300 kW, and 250 kW, respectively, at phase 1 of node 652, phase 2 of node 645, and phase 1 of node 675 as shown in Figure 3. Also, a synchronous generator of inertia constant (H) equal to1:01s is added at the substation node to emulate the calculation of the ROCOF. TABLE I CHARACTERISTIC FUNCTION OF POSSIBLE COALITIONS IN TERMS OF INITIAL ROCOF FOR THE MODIFIED IEEE 13- NODE SYSTEM Possible DER coalitions Initial ROCOF (Hz/s) DER 1 0:0618 DER 2 0:0891 DER 3 0:0786 DER 1;DER 2 0:1515 DER 1;DER 3 0:1388 DER 2;DER 3 0:1683 DER 1;DER 2;DER 3 0:2287 Since all three DERs of the modiﬁed IEEE 13-node system are allowed to participate in secondary frequency regulation, there are 7possible coalitions excluding an empty set. For each set of possible coalitions listed in Table I, the initial ROCOF is determined after participation of DERs in each set of coalitions. If only DER 1is allowed to participate in secondary frequency regulation, initial ROCOF is 0:0618 Hz/s. If only DER 2is allowed to participate in the secondary frequency regulation, the initial ROCOF is 0:0891 Hz/s. If DER 1and DER 2are allowed to participate in the secondary frequency regulation, the initial ROCOF is 0:1515 Hz/s. Similarly, the initial ROCOF for all other sets of coalitions can also be obtained, which are tabulated in Table I. It can be seen from the table that the initial ROCOF of a coalition of two or more DERs is not always equal to the sum of initial ROCOFs of individual DERs. This difference will impact the marginal contribution and hence the Shapley values of each player (here, DER). The value of initial ROCOF obtained for each set of possible coalitions shown in the table serves as the characteristic function of the game.TABLE II DISTRIBUTION FACTORS AND CHANGE IN DER S ETPOINTS FOR A REFERENCE POWER OF +500 KWFOR THE IEEE 13- NODE SYSTEM Distributed Energy Distribution Change in active power Resources Factors set points ( Pref) DER 1 0:2675 134 kW DER 2 0:3916 196 kW DER 3 0:3409 170 kW Based on the characteristic functions as shown in Table I, Shapley values of each DER can be computed using (3) and based on the Shapley values of each DER, the change in active power set point of each DER for a particular reference power, Pref, can be obtained using (7) and (8). For a command ref- erence power of +500 kW received from the system operator, the change in active power sets points of DER 1, DER 2, and DER 3in case of the modiﬁed IEEE 13-node system are 134 kW,196kW, and 170kW, respectively. B. The Modiﬁed 33-node System The33-node distribution test system is a 12:66kV radial distribution system with 33nodes and 32branches [19]. The total active and reactive power loads on this system are 3715 kW and 2300 kV Ar, respectively. In this paper, the 33-node system is modiﬁed by placing four DERs of capacity 300 kW,200 kW,400 kW, and 200 kW at nodes 7,14,24, and 32, respectively, as shown in Fig. 4. In addition to this, a synchronous generator of inertia constant (H) equal to 1:01s is added at the substation node. TABLE III CHARACTERISTIC FUNCTION OF POSSIBLE COALITIONS IN TERMS OF INITIAL ROCOF FOR THE MODIFIED 33-NODE SYSTEM Possible DER coalitions Initial ROCOF (Hz/s) DER 1 0:0928 DER 2 0:0626 DER 3 0:1236 DER 4 0:0632 DER 1;DER 2 0:1551 DER 1;DER 3 0:2160 DER 1;DER 4 0:1557 DER 2;DER 3 0:1860 DER 2;DER 4 0:1255 DER 3;DER 4 0:1864 DER 1;DER 2;DER 3 0:2781 DER 1;DER 2;DER 4 0:2178 DER 1;DER 3;DER 4 0:2786 DER 2;DER 3;DER 4 0:2485 DER 1;DER 2;DER 3;DER 4 0:3405 Since all four DERs of the modiﬁed 33-node system are allowed to participate in secondary frequency regulation, there are16possible coalitions including empty set and single set of DERs. For each set of possible coalitions listed in Table III, the transient analysis is performed to determine initial ROCOF of after participation of DERs in each set of coalitions. If only DER 1is allowed to participate in secondary frequency regulation, initial ROCOF is 0:0928 Hz/s. If only DER 2is allowed to participate in secondary frequency regulation, initial ROCOF is 0:0626 Hz/s. If DER 1and DER 2are allowed to participate in secondary frequency regulation, initial ROCOF is0:1551 Hz/s. In this way, the initial ROCOF for all sets of possible coalitions are obtained, which are tabulated in Table III. It can be seen from the table that the initial ROCOF of a coalition of two or more DERs is not always equal to the sum of initial ROCOFs of the individual DERs. This difference will impact the marginal contribution and hence the Shapley values of each player (here, DER). The value of initial ROCOF obtained for each set of possible coalitions in the table serves as the characteristic function of the game. TABLE IV DISTRIBUTION FACTORS AND CHANGE IN ACTIVE POWER SETPOINTS OF DER FOR A REFERENCE POWER OF +500 KWFOR THE 33-NODE SYSTEM Distributed Energy Distribution Change in active power Resources Factors set points ( Pref) DER 1 0:2712 136 kW DER 2 0:1828 91kW DER 3 0:3617 181 kW DER 4 0:1843 92kW Based on the characteristic functions as shown in Table III, Shapley values of each DER can be computed using (3) and based on the Shapley values of each DER, the change in active power set point of each DER for a particular reference power, Pref, can be obtained using (7) and (8). For a command ref- erence power of +500 kW received from the system operator, the change in active power sets points of DER 1, DER 2, DER 3, and DER 4in case of the modiﬁed 33-node system are 136kW, 91kW,181kW, and 92kW, respectively. VI. C ONCLUSION In this paper, a cooperative game theoretic two-stage ap- proach for secondary frequency regulation of active distri- bution systems has been proposed. The distribution system performs secondary frequency regulation based on a command received from system operators or measured locally. In the ﬁrst stage of the proposed approach, the characteristic functions of the game were determined based on the initial ROCOF computed after performing transient analysis. In the second stage, the characteristic functions were used to compute Shap- ley values of each player (here, DER) and Shapley values were used to determine active power set points of each DER. The proposed approach was implemented on the modiﬁed IEEE 13-node and the modiﬁed 33-node distribution test systems. The case studies exhibit the effectiveness of the proposed approach for secondary frequency regulation of an active distribution system. ACKNOWLEDGEMENT This paper is based upon work supported by the U.S. De- partment of Energy’s Ofﬁce of Energy Efﬁciency and Renew- able Energy (EERE) under the Solar Energy Technologies Of- ﬁce Award Number DE-EE0009022. This paper was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government norany agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.