Distribution Locational Marginal Pricing for Combined Thermal and Electric Grid Operation

Distribution locational marginal prices (DLMPs) have been formulated for electric distribution grids to economically dispatch distributed energy resources (DERs) while addressing operational constraints of the electric grid. This paper proposes to extend this methodology to thermal grids, i.e, district heating or cooling systems, and specifically the combined operation of thermal and electric grids. To this end, thermal and electric grid models are formulated in a linear approximate model fashion. Then, the derivation and decomposition of DLMPs is formulated based on the combined optimal operation problem, assuming a linear state space model form for flexible loads (FLs). The ability of the DLMPs to reflect operational constraints in the thermal grid is demonstrated for a test case with 22 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, the concept of distribution locational marginal prices (DLMPs) was developed to economically dispatch DERs, while keeping within the operational constraints of the electric grid. If not operated appropriately, DERs may lead to increased losses and congestion in the distribution grid, whereas DLMPs have been demonstrated to efficiently aid with loss reduction [1] and congestion management [2], [3].
The concept for DLMPs is related to location marginal prices (LMPs), which are deployed by electric transmission grid operators, e.g., independent system operators (ISOs) in the US and UK [4]. LMPs refer to nodal electricity prices, which express the marginal cost of supplying load at a particular grid node. To this end, LMPs incorporate information on electric grid losses and congestion in addition to the marginal generation cost. For example, LMPs will increase at nodes which are served through a congested grid line, therefore serving as an incentive to reduce the load at these nodes to alleviate the congestion. Hence, LMPs are integrated into the day-ahead electricity market to add consideration for the grid operation limits. DLMPs are essentially an extension of LMPs to the distribution grid [4].
DLMPs are formulated based on the numerical optimization problem for the operation of the electric grid, which aims cost minimization for the DERs, while maintaining voltage drop, line loading and loss balance constraints of the electric grid. By obtaining the dual variables associated with these constraints, DLMPs can be formulated to essentially reflect the operational requirements of the electric grid into price signals to incentivize the contribution of DERs, e.g., by rescheduling loads. This has led to the envisioning of distribution system markets based on DLMPs to address the integration of DERs in [5].
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 shift of paradigms similar to the electric distribution grid. Particularly, combined heat and power (CHP) plants, heat pumps (HPs), solar-thermal plants have motivated research into bi-directional thermal grids [6]. 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. To this end, the concept of DLMPs which was proposed for the electric distribution grid can readily be extended for thermal grids, to aid the economic dispatch of thermal DERs in a similar fashion.
This paper formulates a methodology to obtain the DLMPs for the combined optimal operation of the thermal and electric distribution grids, such that the DLMPs reflect the operational constraints of the combined grid. To begin with, thermal and electric grid models are formulated in a linear approximate model fashion, which later facilitates the derivation of the DLMPs. Along with this, a linear FL model in state-space form is defined. Then, the combined optimal operation problem is formulated for the thermal and electric grids with interconnected FLs, which serves as the basis for the calculation and decomposition for the DLMPs. Lastly, the approach is applied for a test case with 22 FLs and congruent thermal and electric distribution grids based in Singapore, where the DLMPs are shown to represent the operational constraints of the combined grid.

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 diag(x) constructs a diagonal matrix with the entries of x. The Hadamard product is denoted by . Prices are in Singapore Dollar (SGD) which is denoted by S$.

II. MODELLING
A. Thermal grid model The following thermal grid model formulation assumes that 1) the thermal grid forms a connected tree graph, i.e., a radially connected grid, 2) a distributed pumping scheme is applied, i.e., distribution pumps are installed at the energy transfer station (ETS) of each FL rather than centrally, 3) thermal losses in the distribution branches are neglectable, 4) the layout of supply and return piping is symmetric and 5) thermal power is transmitted with constant supply and return temperatures, i.e., a constant absolute enthalpy difference between supply and return side. Note that points 4) and 5) render the proposed formulation more appropriate for district cooling systems, whereas for district heating systems, the model will require an extension to model dynamic supply and return temperatures.
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 vec- 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, which is formulated based on [7], [8] as: The matrix A n th ,b th ∈ R n th ,b th is the thermal branch to node incidence matrix with entries A n th ,b th n th ,b th = 1 if branch b th starts at node n th ; A n th ,b th n th ,b th = −1 if branch b th ends at node n th ; A n th ,b th n th ,b th = 0 otherwise. Note that A n th ,b th must be square and non-singular, which is satisfied if the thermal grid forms a connected tree graph. The matrix A n th ,f is the FL to thermal node incidence matrix with entries A n th ,f n th ,f = 1 if FL f is connected at node n th ; A n th ,f n th ,f = 0 otherwise. The scalars ρ wa and ∆h wa define the density of water and the absolute enthalpy difference experienced by the distribution water between the primary and secondary side. To this end, eqs. (2a) and (2b) are equivalent to formulating the nodal flow balance for the thermal grid nodes and the branch pressure head balance based on the Darcy-Weisbach equation: where F n th are the FLs connected at node n th , B n th ,1 /B n th ,2 are the branches starting/ending at node n th and n th,b th ,1 /n th,b th ,2 denotes the start/end node of branch with entries f f r b th describes the Darcy-Weisbach friction factor for branch b th depending on its absolute roughness ε b th , diameter d b th and the Reynolds number Re ref b th of the branch flow. Note that this formulation based on the Swamee-Jain Formula assumes that the Reynolds number stays in 4 · 10 3 ≤ Re ref b th ≤ 10 8 and the branch parameters stay in 10 −6 ≤ ε b th d b th ≤ 10 −2 [9]. Note also that minor friction losses, i.e., losses due to curvatures and junctions, are neglected in the model formulation. Lastly, the vector h ets ∈ R f describes the head loss across the ETS, i.e., heat exchanger, at full load for each FL f ∈ F.
A global approximation for the sensitivity matrices M v,p th , M h,p th , M p pm ,p th can be obtained based on the formulation of the thermal grid power flow problem in eq. (2). The global approximation essentially forms a linear expression as a tangent between the no-load point and the nominal operation point of the thermal grid for each non-linear model equation.
In fact, only the pressure head loss terms eqs. (2b) to (2d) require linearization whereas eqs. (2a) and (2e) are are simply translated into the matrix form. The derived global sensitivity matrices is provided as: 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- 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. (5) 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, input 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 input 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 the FL c f,t = c f,t (p th f,t , p f,t , q f,t ). 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 reactive power, is assumed. The detailed formulation of the state space model for such buildings can be obtained from [10].

A. Combined optimal operation problem
The combined optimal operation problem of the thermal and electric grid addresses the economic dispatch of FLs subject to the operational constraints of the thermal grid, the electric grid and the FLs: 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 are the Lagrangian multipliers associated the respective inequality and equality constraints. 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) and (7c) to (7e)) and the electric grid (in eqs. (5) and (7f) to (7j)) are coupled only through the FL operation (in eqs. (6) and (7k)).

B. DLMP derivation and decomposition
In order to obtain the distribution locational marginal prices (DLMPs), the Lagrangian function of eq. (7) is formulated for Note that for the loss and pumping components, the Lagrangian multipliers become λ p ls

C. Application of DLMPs
Considering an independent system operator which operates the combined thermal and electric distribution system and allocates prices for the electric and thermal consumption, the market procedures can be cast as follows: 1) FLs submit their instantaneous bids, energy requirements as well as their dispatch capabilities to the system operator; 2) The system operator forecasts its underlying electric/thermal demand and solves the combined optimal operation problem eq. (7) in a day-ahead manner; 3) The system operator obtains the respective DLMPs for electric and thermal consumption and then passes them to the FLs; 4) The FLs receive the price signals and solve the local FL operation problem. This results in an optimal dispatch which implicitly respects the operational constraints of the thermal and electric grids. Note that such a market organization enables the individual market participants like FLs to make bids purely based on the locally available information. This essentially enables the price-based control for the system operator over FLs that act as price takers for both electric and thermal demand.

IV. RESULTS AND DISCUSSION
The presented methodology is demonstrated for a test case based on a neighbourhood in Singapore, which was originally developed as part of the City Energy Analyst (CEA) [11]. The test case consists of thermal and electric grids with identical layout according to fig. 1 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 following, three scenarios are distinguished for 1) unconstrained operation, 2) constrained branch flow and 3) constrained pressure head.
Scenario 1: Unconstrained operation: In this scenario, the optimal dispatch of the FLs does not incur any thermal or electric grid constraint violations. Figure 2 depicts the thermal grid DLMPs of FLs "16" and "17". The observed thermal grid DLMPs consist only of energy and pumping components, where the variations in energy component and pumping component correspond to changes in c ref t . Thermal and electric demand are scheduled such that the apparent price peak at 09:30 is avoided.
Scenario 2: Constrained branch flow: For this scenario, an artificial branch volume flow constraint was introduced at the branch connecting nodes "A" and "B" (see fig. 1). As depicted in fig. 3  observed at FL "16", as it is adjacent to a different section of the thermal grid which does not connect to the source node through branch "A"-"B". Even with this minor price change, FL "17" is incentivized to reschedule its demand.
Scenario 3: Constrained pressure head: Scenario 3 considers an artificial node pressure head constraint at node "B". Figure 4 highlights that this introduces head support DLMP components for both FLs "16" and "17". This demonstrates that the node pressure head is a function of all upstream nodes' pressure head, as described by eq. (3) and as such interdependent between different grid sections. FL "17" in this case again reschedules its demand, apparently to avoid a slightly increased price peak at around 12:00.

V. CONCLUSION
This paper presented an extension for distribution locational marginal prices (DLMPs) to combined thermal and electric grids. Based on linear thermal and electric grid models, the derivation and decomposition of DLMPs was formulated for the combined optimal operation problem. The approach was  demonstrated for a test case based in Singapore. Future work will focus on the impact of thermal grid injections and model extensions to consider thermal losses in the thermal grid. The presented models and test case are implemented in [12].