Coordinated Market Clearing for Combined Thermal and Electric Distribution Grid Operation

To economically dispatch distributed energy resources (DERs) while addressing operational concerns of the electric grid, distribution locational marginal price (DLMP)-based market frameworks have been formulated for the electric distribution grid. This paper proposes to extend this methodology to thermal grids, i.e. district heating or cooling systems, and presents a market-clearing mechanism based on the alternating direction method of multipliers (ADMM) for coordinating between thermal grid, electric grid and DER operation. The ability of the proposed mechanism to achieve the market equilibrium for a combined thermal and electric grid is demonstrated for a test case with 22 flexible loads (FLs).


I. INTRODUCTION
With the increasing integration of distributed energy resources (DERs) such as flexible loads (FLs) and distributed generators (DGs) into the electric distribution grid, local energy market topics are enjoying elevated research attention. If not operated appropriately, DERs may lead to increased losses and congestion in the electric distribution grid. Therefore, it is necessary to design appropriate market organization and pricing mechanisms that allow for the fair allocation of the operation cost of the electric grid as well as the cost-effective integration of DERs. In this scope, distribution locational marginal prices (DLMPs) have proved to aid with loss reduction [1] and congestion management [2], [3]. DLMPs essentially translate information on losses and congestion into incentive signals for the DERs to encourage operational contributions, e.g. the rescheduling of loads. The DLMP method naturally lends itself as a framework for local electric distribution grid markets, since it relates operational constraints into monetary objectives. DLMP-based market frameworks have been studied for the electric distribution grid in theoretical as well as real-life projects [4], [5], [6].
In district heating and cooling systems, i.e. thermal grids, thermal DERs in the form of distributed heating and cooling sources are leading to a paradigm shift similar to the electric distribution grid. Particularly, combined heat and power (CHP) plants, heat pumps (HPs) and solar-thermal plants have motivated research into bi-directional thermal grids [7]. On the demand side, buildings equipped with heating, ventilation and air-conditioning (HVAC) systems have already been proposed as FLs for the electric grid and can be retrofitted in a similar fashion to increase operational flexibility in the thermal grid. Our recent work [8] formulates a thermal grid model and methodology for deriving DLMPs for thermal grid operation.
The thermal grid DLMPs were shown to adequately represent congestion and pumping losses in the thermal grid similar to electric grid DLMPs. With the availability of the DLMP technique in both the thermal grid and the electric grid, there remains the issue of local market organization across thermal grid and electric grid boundaries.
The combined operation problem of the thermal grid, the electric grid and the DERs can in principle be formulated as a centralized economic dispatch problem, i.e. an optimization problem for the minimization of energy costs with respect to the operational constraints of the thermal grid, the electric grid and the DERs. However, the operation of each system is typically managed by a dedicated organization, i.e. thermal grid operator, electric grid operator and DER aggregators, where each entity may prefer to retain control and data of their respective systems. To this end, distributed optimization techniques such as alternating direction method of multipliers (ADMM) can be utilized to facilitate the energy trade between different entities while retaining organizational boundaries [1]. In general, ADMM serves as a framework for distributed optimization and is widely proposed for applications in power systems due to its good scalability and robustness [9].
In this work, we propose a decentralized local market organization for combined thermal and electric distribution grid operation, enabled by the ADMM solution methodology as a market-clearing mechanism. Along with this, linear model formulations for the thermal grid, the electric grid and FLs are presented. The proposed market-clearing mechanism is shown to achieve the market equilibrium for a test case with 22 FLs connected to thermal grid and electric grid, even in the presence of congestion.

A. Thermal grid model
The thermal grid is modelled with a linear approximate model as: The vectors h t ∈ R n th , v t ∈ R b th are the pressure head at thermal grid nodes n th ∈ N th and the branch volume flow at thermal grid branches b th ∈ B th for time step t ∈ T . The scalar p pm t denotes the total electric distribution pumping power demand for time step t. The reference point for each property is denoted by () ref , which in the following is chosen to be the nominal operation point of the thermal grid. The vector ∆p th t ∈ R f is the thermal power change at FLs f ∈ F for time step t. The matrices M h,p th ∈ R n th ,f , M v,p th ∈ R b th ,f , M p pm ,p th ∈ R 1,f are the sensitivity matrices for the change of the respective properties to the thermal power change, which in turn is defined as ∆p th t = p th t −p th,ref . The vector p th t , ∈ R f is the absolute thermal power demand and the vector p th,ref t is the thermal power demand reference which is chosen to be the nominal load, i.e. peak load, of the FLs f . The reference properties for the linear approximate model in eq. (1) are then obtained by solving the reference power flow problem of the thermal grid. The detailed model formulation is omitted here for the sake of brevity, but is discussed in detail in [8].

B. Electric grid model
The electric grid is modelled with a linear approximate model as: The vectors u t ∈ R n el , |s f t | 2 , |s t t | 2 ∈ R b el are the voltage magnitude at electric grid nodes n el ∈ N el and the squared branch power flow in "from" and "to" direction at electric grid branches b el ∈ B el for time step t ∈ T . The scalars p ls t , q ls t denote the total active and reactive loss for time step t. The reference point for each property is denoted by () ref , which in the following is chosen to be the nominal operation point of the electric grid. The vectors ∆p t , ∆q t ∈ R f are the active and reactive power change at FLs f ∈ F for time step t. The ma- ,f are the sensitivity matrices for the change of the respective properties to the active and reactive power change, which in turn are defined as are the active and reactive power demand reference which is chosen to be the nominal load, i.e. peak load, of the FLs f . The reference properties for the linear approximate model eq. (2) are then obtained by solving the reference power flow problem of the electric grid, e.g. through a fixed-point solution methodology, and a global approximation for the sensitivity matrices can be obtained as a function of the reference properties. The detailed formulation is omitted here for the sake of brevity, but can be obtained from [3].

C. Flexible load model
The FL model is expressed in state space form as: The vectors x f,t , c f,t and d f,t are the state, control and disturbance vectors for FL f at time step t. The matrices A f , C f are the state and output matrix, and are the control and feed-through matrices, on the control and disturbance vectors respectively. Note that control vector is a function of the thermal, active and reactive power dispatch of . Therefore the FL model serves as the interconnection between the thermal and electric grids.
For the presented test case (section IV), air-conditioned buildings serve as FLs and a fixed power factor, i.e. a fixed ratio between active and apparent power, is assumed. The detailed formulation of the state space model for such buildings can be obtained from [10].

III. MARKET ORGANIZATION
As depicted in figure (fig.) 1, the considered participants in the local thermal energy market and the local electric energy market on district level are 1) thermal grid operator, 2) electric grid operator, and 3) flexible load aggregator. Note that for notational brevity, the following formulations are for a single flexible load aggregator, although there may in fact be several aggregators within a district. Also note that while not explicitly considered in this paper, the flexible load aggregator can be transformed into a DER aggregator to take into account the dispatch of other DER types, e.g. generators, fixed loads or energy storage systems.
The flexible load aggregator is responsible for dispatching the flexible loads with respect to their operational constraints and correspondingly trading thermal power and active / reactive power with the thermal grid operator and the electric grid operator respectively. The thermal grid operator is responsible for maintaining the operational constraints of the thermal grid, while trading thermal power with the flexible load aggregator and obtaining any active / reactive generation from the wholesale electricity market. The electric grid operator is responsible for maintaining the operational constraints of the electric grid, while trading active / reactive power with the flexible load aggregator and obtaining active / reactive generation from the wholesale electricity market. Note that wholesale electricity market clearing is assumed to be cleared independently in advance to the local market clearing. Thus, the scalar c ref t is the cleared wholesale electric energy price.
The local electric and thermal energy markets are cleared by utilizing the alternating direction method of multipliers (ADMM). ADMM serves as a framework for distributed optimization and can be interpreted as a price coordination mechanism to achieve the Walrasian tâtonnement process [11], which only requires the information exchange of proxy variables without revealing the sensitive local information like system models and constraints. To this end, the market participants exchange negotiation variables with the local thermal and electric energy markets.
The negotiation variables are depicted with dashed arrows in fig. 1. The vectors p th,f l t and p th,th t denote the thermal power requested by the the flexible load aggregator and the thermal power supply offered by the thermal grid operator. The vectors p th,ex t , π th,f l t and π th,th t are the resulting thermal power exchange and nodal prices for thermal power demand and supply. The vectors p f l t / q f l t and p el t / q el t denote the active / reactive power requested by the the flexible load aggregator and the active / reactive power supply offered by the electric grid operator. The vectors p ex t / q ex t , π p,f l t / π q,f l t and π p,el t / π q,el t are the resulting active / reactive power exchange and the nodal prices for active / reactive power demand and supply.
The final cleared power and price schedule are denoted with solid arrows in fig. 1. The vectors p th t and π th t denote the final cleared thermal power exchange and the associated nodal price for thermal energy trade between the thermal grid operator and the flexible load aggregator. The vectors p t , q t as well as π p t , π q t denote the final cleared active / reactive power exchange as well as the associated nodal price for active / reactive energy trade between the electric grid operator and the flexible load aggregator for time step t.

A. Combined optimal operation problem
The combined optimal operation problem addresses the economic dispatch of FLs subject to the operational constraints of the thermal grid, the electric grid and the FLs. The ADMMbased market-clearing algorithm (alg.) is applied as a distributed optimization solution methodology for the combined optimal operation problem. To form the basis for deriving the individual sub-problems of the market clearing alg., first, the combined optimal operation is formulated as follows: The scalar c ref t is the electric energy price at the reference node and η ch is the coefficient of performance (COP) of the district heating or cooling plant. Note that this formulation translates to assuming 1) a constant-efficiency model for the performance of the heat pump at the district heating or cooling plant and 2) that electric power required for district heating or cooling is drawn at the source node of the electric grid. The vectors u − , u + ∈ R n el , |s f,+ | 2 , |s t,+ | 2 ∈ R b el describe the electric grid voltage and branch loading limits, whereas h − ∈ R n th , v + ∈ R b th describe the thermal grid head and volume flow limits. The scalars p src t , q src t , p th,src t describe the total active, reactive power demand and the electric power demand of the thermal grid. The vectors y − f,t , y + f,t describe the operational constraints of the FLs f , which may be timedependent to reflect set-back periods. Note that this formulation of the combined operation assumes that the district heating or cooling plant of the thermal grid is co-located with the electric grid source node, such that the thermal grid source node is not subjected to any electric grid constraints. Therefore, the operation of the thermal grid (in eqs. (1), (4c) and (4d)) and the electric grid (in eqs. (2), (4e) and (4f)) are coupled only through the FL operation (in eqs. (3) and (4g)).

B. Market clearing
For the ADMM-based market clearing, individual subproblems are formulated for 1) thermal grid operator, 2) electric grid operator and 3) flexible load aggregator. Each subproblem consists of a subset of the objective and constraints from eq. (4) which are augmented by ADMM Lagrangian terms. In the following, ρ > 0 denotes the ADMM penalty factor, which serves as the convergence tuning parameter for the market clearing alg. and remains constant during ADMM iterations.
1) Thermal grid operator: The optimal operation problem of the thermal grid operator is formulated as: is terminated once the residuals reach the desired termination threshold given by > 0. The ADMM-based decomposition method scales reasonably well with the increasing size of complicating variables, as it was shown through several examples to be well suited for a wide variety of large-scale distributed optimization problems [11].

IV. RESULTS AND DISCUSSION
The presented methodology is demonstrated for a test case for a neighbourhood in Singapore, which was developed as part of [12]. The test case consists of thermal and electric grids with identical layout according to fig. 2 and 22 commercial buildings modelled as FLs according to [10], where the thermal grid serves as a district cooling system (DCS). The source node electricity price c ref t is derived for one day of the Universal Singapore Energy Price (USEP). In the presented scenario, an artificial branch volume flow constraint is defined in thermal grid at the branch between nodes "A" and "B" (see fig. 2). Figures 3 and 4 depict the dispatch and price schedules for thermal and reactive power which are obtained from centralized operation, thermal grid operator, electric grid operator and flexible load aggregator. In this context, the centralized  The results for reactive power coincide with the results for thermal and active power and are omitted for brevity. Note that the dispatch schedules for thermal and active power exhibit a similar shape, because the demand is mainly driven by air conditioning, with thermal power serving the cooling demand and electric power the pumping and ventilation demand. The mid-day dip in the dispatch schedules is due to congestion, which highlights the ability of alg. 1 to determine a marketclearing solution despite the presence of congestion. Figure 5 depicts the primal residuals, where the alg. 1 converged after approx. 180 iterations for = 10 −6 . The ADMM parameters were not tuned for convergence speed according to [11], although this is strongly suggested for practical applications.
V. CONCLUSION This work presented a local market organization for combined thermal and electric distribution grid operation along with an ADMM solution methodology as a market-clearing mechanism. The proposed market-clearing mechanism was shown to achieve the market equilibrium in the presence of congestion. The presented models and test case are implemented in [13] and available open source 1 .

ACKNOWLEDGMENT
This work was financially supported by the Singapore National Research Foundation under its Campus for Research Excellence And Technological Enterprise (CREATE) programme.

NOMENCLATURE
Let R be the domain of real numbers. Non-bold letters x, X denote scalars R 1 , bold lowercase letters x denote vectors R n and bold uppercase letters X denote matrices R n,m . Bold numbers 0 and 1 denote vectors or matrices of zeros and ones of appropriate sizes. The transpose of a vector or matrix is denoted by () and the p-norm of a vector is denoted by ||x|| p .