Distributed Optimization of District Heating Networks Using Optimality Condition Decomposition

: The optimal operation of District Heating Networks (DHNs) is a challenging task. Current or future optimal dispatch energy management systems attempt to optimize objectives, such as monetary cost minimization, emission reduction, or social welfare maximization. Typically, this requires highly nonlinear models and has a substantial computational cost, especially for large DHNs. Consequently, it is difﬁcult to solve the resulting nonlinear programming problem in real time. In particular, as typical applications allow for no more than several minutes of computation time. However, a distributed optimization approach may provide real time performance. Thereby, the solution of the central optimization problem is obtained by solving a set of small-scale, coupled optimization problems in parallel. At runtime, information is exchanged between the small-scale problems during the iterative solution procedure. A well-known approach of this class of distributed optimization algorithms is Optimality Condition Decomposition (OCD). Important advantages of this approach are the low amount of information exchange needed between the small-scale problems and that it does not require the tuning of parameters, which can be challenging. However, the DHNs model equation structure brings along many difﬁculties that hamper the application of the OCD approach. Simulation results demonstrate the applicability range of the presented method.


Introduction
Optimal operation of District Heating Networks (DHNs) is a major task to lower emissions, operation costs, network losses or to maximize the social welfare in marketbased forms of DHN operation.Thereby, the optimal DHN operation regarded within this work is understood as a form of dispatch for producers and flexible loads, which is performed after a previously applied unit commitment (unit commitment is not within the scope of this work); see, e.g., [1] for operation and control hierarchies in DHN operation.The operational optimization problem is formulated as a Nonlinear Programming (NLP) problem, and an overview on the kinds of operational optimization problems used for DHNs is given in [2].
The most efficient form of DHN operation comprises variables mass flows [3], Variable Mass Flow Directions (VMFDs) [4] and variable temperatures [5].Thereby, the operation points of pumps, valves and the heat power injection and demand of all flexible network participants are set by an Energy Management System (EMS).Note that these flexible network participants can be producers, consumers, energy converters or thermal storage systems.When representing the physics of variable flows and variable temperatures in a mathematical model for operational optimization in an EMS, the computational burden becomes high for large real-world DHNs [6].In certain cases, the resulting calculation times can surpass the desired time interval used for the operational optimization.
Simulations performed by the authors of this work, using a rolling horizon approach, where the operational NLP problem was implemented in GAMS and solved by the IPOPT solver, for one of the 15 largest DHNs in Germany, showed calculation times above 15 min for prediction horizons larger than four time steps.Since prediction horizons in this range are low compared with the thermal transient propagation of temperature fronts throughout the DHN, the results of the optimizations can therefore most likely be increased for larger optimization horizons.
As time step intervals used in energy management systems for market and control frameworks for these kinds of systems are often within 15 min [7], further measures are needed to enable real time capability of the EMSs.The demand for further measures to reduce calculation times for future network operation becomes even more apparent in the light of an increasing number of flexible network participants needing to be coordinated in future DHNs [8].These result from the necessity to decarbonize future heat supplies, and thus, increasingly, small producers, such as renewable energy sources, heat pumps, powerto-heat devices or waste heat suppliers (such as data centers, refrigerated warehouses and purification plants), are being connected to DHNs [8], replacing the earlier centralized power plants.
A well-known measure to reduce calculation times is the parallelization of computations by decomposing the original central operational optimization problem into multiple coupled small-scale subproblems, which are then solved in a coordinated manner.Different approaches for parallelization are presented in [9,10].A large field of literature is provided for the distributed optimization of electric power networks or coupled electric power and district heating networks [11][12][13][14][15][16][17][18][19].Most of the distributed optimal operation approaches use the Alternating Direction Method of Multipliers (ADMM) [20] or the Optimality Condition Decomposition (OCD) [21] algorithms.
All the stated approaches decompose the operational optimization problems for the energy systems into small-scale subproblems, by dividing the energy network models into multiple zones, which are then optimized in parallel by multiple EMSs.However, the mentioned approaches all decompose the energy networks within electric power networks or at the boundary of electric power and district heating networks.Thus far, almost no work exists that decomposes the DHNs themselves into subproblems.Important advantages of the OCD approach are the low amount of information exchange needed between the small-scale problems and that it does not require the tuning of parameters, which can be challenging [16].These tuning procedures are needed for ADMM approaches [9].Further, the OCD approach enables real parallelization, while ADMM uses sequential calculations in its standard form [20].
Note that these sequential calculations also include computations by a central coordinating instance for ADMM approaches.Thus, three sequential calculations need to be performed for every iteration of the ADMM algorithm [20].Furthermore, a convergence criterion is provided in the literature for OCD [21], which enables to compare different forms of decomposition by convergence speed.This criterion can be used to compare different forms of network decomposition with each other [22].
Through the symmetrical flow conditions in DHNs in the supply and return networks, the coupling of the subproblems arising from a certain decomposition is inherently large in DHNs.This impedes the straight-forward application of the OCD approach, as strong couplings between the subproblems prevent fulfillment of the convergence condition.Thus, an adequate reformulation of the operational optimization model is needed, which represents the physical effects with similar precision but reduces the coupling of the subproblems.
To the best of the authors' knowledge, thus far, the field of distributed optimization of DHNs is barely researched.An approach for distributed optimization of DHNs based on the ADMM is presented in [23].As far as we know, the literature does not provide any approaches applying OCD to DHNs.Thus, we would like to bridge that gap and provide a distributed optimization approach that can be solved in parallel.Therefore, the contributions of this paper are:

•
A new approach for the distributed optimal operation of DHNs based on OCD.
• In order to enable the application of OCD, the inherently strong coupling of small scale subproblems needs to be reduced by an adequate model reformulation.Thereby, the relevant physical effects are further accurately described.

Notation
Column vectors x ∈ R n are represented with boldface lowercase letters, matrices A ∈ R n×n with boldface capital letters.For simplicity, we represent a vector composed of two stacked column vectors as x = [x 1 , x 2 ].The notation ∇ x f is defined as the partial derivative of f w.r.t.x.We denote x as the transpose of vector x.

Structure
The rest of this paper is organized as follows: Section 2 explains the basic principles and properties of the OCD algorithm.Further, Section 3 provides the operational optimization model of the DHN, the reformulation needed for OCD and the final OCD approach for DHN operation.A proof of principle is given based on the simulation results in Section 4. Finally, a discussion is found in Section 5.

Optimality Condition Decomposition
In this section, the OCD algorithm first presented in [21] is briefly reviewed.The decomposition approach of OCD is presented here without loss of generality for the case of two subproblems/operational zones, marked with the superscripts za and zb below.An exemplary central optimization problem used to illustrate the approach is given as [10]: with the equality constraints divided into two parts, h za and h za .We assume the global objective function is the sum of the objective functions of the operational zones, as these represent cost, and we aim to minimize the sum of all costs.Note that the optimization problem cannot be solved separately by the two operational zones, since the optimization problems are coupled through the coupling constraints (1b) and (1c) (also called complicating constraints).For reasons of simplicity, we neglect constraints that depend on one group of variables only, since they can be readily included [10,21].Then, the Lagrangian of (1) is defined as: with the Nabla operator ∇ and the search directions of the central problem given by ∆ cent,za = ∆x cent,za , ∆λ cent,za and ∆ cent,zb = ∆x cent,zb , ∆λ cent,zb .Further, the KKT matrices are defined as [10]: The OCD approach decomposes the optimization problem (1) in a special manner into the following two subproblems, for Zone a: and Zone b: where the operator indicates that these (dual) variables are constant and considered as parameters in the respective subproblem.The complicating constraint of the other operational zone is included in the objective function of the regarded operational zone.All variables and constraints have been assigned to a certain operational zone, while, for the objective function, each operational zone only considers its own objective function.Further, (5b) and (6b) represent complicating constraints.Complicating constraints are not only considered within their proper operational zone by the approach but also as soft constraints within the objectives of the other operational zones related through the complicating constraint, see (5) and (6).In either case, the respective parameters are obtained from the (dual) variables of the other operational zone from the last Newton iteration v − 1.This indicates the basic OCD algorithm for every Newton iteration ν: Perform one Newton iteration of the two subproblems ( 5) and (6).

2.
Exchange the parameters between the subproblems/operational zones, here x za , x zb , λ za and λ zb .

3.
Check if the stopping criterion is fulfilled.If it is not met yet, then calculate the next Newton iteration v + 1 in step 1. Different stopping criteria are stated in the given literature [16,21].Computing these distributed subproblems with Newton's method is equal to solving the following set of equations within every Newton iteration [16]: with the distributed search directions ∆ dist,za = ∆x dist,za , ∆λ dist,za and ∆ dist,zb = ∆x dist,zb , ∆λ dist,zb as well as the approximated KKT matrix KKT.By comparing (7) with (3), it is understood that KKT is obtained from KKT by setting the matrices on the secondary diagonals to 0. Now, convergence of the distributed solution towards the central solution of (1) can be shown if, at the provided second-order KKT point y * = x za , x zb , λ za , λ zb , the following holds [21]: Assumption 1.The second-order derivatives of f za , f zb , h za and h zb are Lipschitz continuous.
Note that Assumption 2 is fulfilled if the gradients of the constraints are linearly independent.This is not fulfilled if the system has, e.g., identical, parallel lines between two nodes.Last but not least, the condition for convergence given below needs to hold at y * [21]: with ρ( ) defining the spectral radius, which equals the absolute value of the maximal eigenvalue of the given matrix.Further, KKT * represents the KKT matrix at the optimal solution.If Assumptions 1 and 2 are fulfilled, and Condition 1 is fulfilled as well, then iteratively solving ( 5) and ( 6) will converge towards y * with at least linear rate ρ ocd, * for all starting points y 0 sufficiently close to y * .The literature often refers to ρ ocd, * as the coupling factor.As can be understood from Equation (8), the coupling factor is smaller for decompositions possessing less coupling, since then less non zero entries on the secondary diagonals of the KKT block matrix KKT za,zb and KKT zb,za are neglected in KKT * .In the context of optimal dispatch of energy systems, this refers to less tie lines or pipelines connecting the decoupled operational zones [16].
Within the following remarks further important aspects of the OCD approach w.r.t. the present work are described.
Remark 1.The OCD approach, which was presented here for the simple case of an optimization problem with complicating constraints, can be simply expanded to the case of additional non complicating equality and inequality constraints, which only include variables of the respective operational zone.Furthermore, complicating inequality constraints, such as flow limits over border edges, can be readily included in the same manner as complicating equality constraints [24].
Remark 2. There are indications that symmetrically decomposed problems are more likely to fulfill Condition 1. Symmetrical decomposition means that a complicating constraint in both operational zones is present, which has an identical part and incorporates the same (parameterized) variables.The first indication is that the majority of the published applications of OCD to energy system dispatch problems decompose the central problems in a symmetrical form.For example, in optimal power flow problems, the power flow over a tie line P flow i,j is part of the power flow constraints, of both border nodes in the two adjacent operational zones and given by: This power flow is dependent on the voltage angle θ and voltage amplitude V of the bordering nodes i and j, which are thus part of both complicating constraints [16,25].The parameters B bus and G bus represent the susceptance and conductance of the nodal admittance matrix used to model the electric power network for alternating current power flow calculations.Furthermore, it becomes clear by comparing the structure of the KKT and KKT matrices that the distributed search directions ∆ dist in (7) will most likely be closer to the central search direction ∆ cent in (3) if all operational zones take information of the other (bordering) operational zones into account, while calculating ∆ dist .

Remark 3.
The decomposition of a central optimization problem into many coupled but smaller subproblems achieves a considerable reduction of the computation time.In [21], it is shown that even solving the small scale subproblems sequentially on only one processor, the computation time is less than computing the solution of the central problem directly.Note that this computation time can be lowered by orders of magnitude by further parallelizing the subproblems and solving them in a coordinated manner.This is due to the bad scaling of the computation time with the number of optimization variables in nonlinear programming; it is more efficient to solve many small scale problems sequentially.

Optimal Dispatch of District Heating Networks Using Optimality Condition Decomposition
This section first introduces the central optimization problem used in this work to optimize DHN operation.Based on this, the model reformulation needed to reduce the coupling of the subproblems when applying OCD is presented next.Finally, the application of OCD to the resulting problem is provided including the complicating constraints and the resulting objective functions arising from the OCD approach for each zone are provided.

Central Optimization Problem for Optimal Dispatch of District Heating Networks
The nonlinear programming problem used for operational optimization of DHNs used within this paper is given in the following.This comprises the objective function, the equality constraints, as well as box constraints limiting the permissible range of all variables.

Objective Function
The objective function f , which can, in general, represent emissions, operation costs, network losses or renewable energy sources curtailment, which should by minimized.A further possibility that is used in this work is to maximize the social welfare in a marketbased form of DHN operation, such as those provided in [26].

Equality Constraints
This section briefly describes the DHN models used within this work.In particular, the hydraulic and thermal models of nodes, pipelines, consumers and producers are provided.The DHN network graph is spanned by the node edge incidence matrix, defined as (cf.[27]) where i represents a node from the set of all network nodes S i , and e is a respective edge from all edges S e .

Hydraulic Model
The models of all typical hydraulic elements are given below.For all components representing DHN edges e, a general stationary relation of the mass flow ṁ and the pressure difference ∆p can be defined for all regarded time steps k ∈ S k as: where β and µ are the variable component coefficients and p i represents the nodal pressure at node i.For pipelines, the values of β and µ are calculated as given in [4], and the Prandtl-Colebrook equation is simplified using the Haaland approximation [28].Thereby, turbulent flow conditions are assumed [4].For consumers, we assume β = 0 and µ is to be within the box constraints, which are defined by the parameters of the valve and the heat exchanger within the consumer.Further, producers are assumed to contain a pump, a heat exchanger and a valve.The pumps can vary the value of the box constrained β, while the rest of the producers are modeled as consumers, thus, with a constrained value µ.Note that ( 11) is continuously differentiable as the absolute value function of the mass flow ṁ is replaced by a continuous differentiable approximation given by: wherein ∆ε > 0 is a sufficiently small parameter.In case of a predefined positive flow direction (positive flow direction equals the edge orientation), when the e is not element of the set of all edges with VMFDs, S vmfd e , (11) can be written as The law of conservation of mass is represented by ( 14) for all mass flows entering and leaving a node i: Thermal Model For all nodes S vmfd i ⊂ S i where, at least on one connected edge, a flow direction change can take place, the node temperature T i is calculated by the heat balance equation, considering all possible flow conditions as given in [4]: The right-hand side of Equation ( 15) contains the sum of all temperatures entering this node and leaving the respective edge T out e , or T in e in the case of negative flow, weighted by their mass flows ṁ.This has to equal the node temperature T i , which is multiplied by the mass flows leaving the node.All possible cases, based on the varying mass flow directions in the pipelines, are distinguished by combining the incoming and leaving node edge incidence matrix elements A + i,e and A − i,e and the approximated maximum operator.The approximation is analogous to the absolute value function approximation in (12) and is defined as: The elements of the incoming and leaving edge node incidence matrices A + and A − of the DHN are defined as: For nodes without edges with VMFDs entering the node, the model Equation ( 15) is simplified to the form, e.g., found in [29].

Consumer and Producer
The heat power supply or demand Φ of producers and consumers is defined by the following equation throughout the existing literature, such as those given in [30]: Pipelines are modeled as in [31] (Node Method Version 1b); thereby, we neglect the steal core model, as future DHN will not possess steal cores due to lower temperature operation values [8], and we precalculate the values y, z (in [31], y and z are written as n and m, which are already in use for other symbols in this paper), R and S before every optimization based on the mass flow values from the last optimization.
This precalculation of the respective parameters enables reducing the obtained mixed integer nonlinear programming to a nonlinear programming problem [32].Further, the length of stay is calculated based on the formula presented in [32] and a predefined flow direction is assumed for model simplicity reasons.Note that pipelines crossing zone borders will be differently modeled in Section 3.3.

Model Reformulation for OCD
In the following, the necessity of the reformulation of the model presented in Section 3.1.2is provided.Before starting, the following definition is provided.For the spatial form of decomposition applied in this work, all network elements will be assigned to a specific operational zone.Thereby, the following special network elements are defined: Definition 1.

•
Border edge: An edge e crossing the border of two neighboring operational zones, e.g., Zone a and Zone b.Thereby, this border edge will be defined as a border edge in both neighboring operational zones a and b. • Border node: A node i connected to a border edge e.

•
Border loop: Border loops are loops l that have edges e in multiple operational zones.
Condition 1 is easily violated when OCD is applied to DHNs.This comes from the fact that decomposing DHN optimization model constraints as provided in Section 3.1.2into subproblems directly results in multiple complicating constraints in every zone.This becomes clear with regard to the small example DHN presented in Figure 1.A naive decomposition into two operational areas Zone a and Zone b applied in this example, yields four complicating constraints for every zone in every time step.This is due to the complicating constraints arising from ( 14) and ( 15) for both nodes in every zone.Based on the model presented in Section 3.1.2,the pipeline model equations would also bring along a complicating constraint.However, this can be prevented by using the Assumption 4, which enables replacing the thermal pipeline model by an identity function mapping the input temperature to the output temperature.From this example, it can be understood that the inherent strong coupling of DHNs results from the topology of DHNs, including the supply and a return network.In general, a stronger coupling of the subproblems, which is increased by every complicating constraint, leads to higher values of ρ ocd, * ; see the paragraph on OCD in Section 2. Thus, the more complicating constraints exist, the more likely it is that Condition 1 is not met.
In contrast to DHNs, applying OCD to electric power networks only leads to a single complicating constraint in every zone for every border node [16], which is why OCD has been applied successfully for distributed optimization of electric power systems.From now on, we will restrict ourselves to the case of two zones, to reduce the notational burden and facilitate understanding of the provided methodology.
To reduce the coupling of the given problem, the hydraulic node Equation ( 14) of the border nodes in the return network, i 3 and i 4 in Figure 1 are omitted.Thereby, two complicating constraints of the form (21)-introduced and explained later on-are eliminated.However, by neglecting the hydraulic border nodes of the return networks, the model becomes incomplete, as the relation between mass flows and differential pressure of supply and return network on the border edges is not incorporated in the model anymore.
This gap is bridged by explicitly introducing the equality ∆p 1 = ∆p 3 for Figure 1, which can be notated as p 1 − p 2 = p 4 − p 3 or sorted by zones as p 1 + p 3 = p 4 + p 2 into the model.This also implies ṁ1 = ṁ3 , based on (11).Using the loop edge incidence matrix B [27], this can be notated as: where the border loops are given in S bl l , and the border edges of Zone a S be,za e represent a subset of the edges of Zone a given in S za e .The same notation is used for edges of Zone b, equivalently.If more then four nodes are found in the regarded border loop l, the relevant border nodes are found through |B l,e ||A i,e | = 1 by checking this condition for all border nodes i ∈ S bi i and edges e, which are no border edges (all border edges are given in e ∈ S be e ) e ∈ S e , e / ∈ S be e of the respective zone.For example, the left side of ( 19) can be simplified to ∑ i∈S bi,za i p i,k if only one point of coupling/border edge (in the return and supply network) to another zone exists as in Figure 1.
A further key aspect represents the artificial extension of the optimization problem used to enable a symmetrical form of decomposition in the following.Applying OCD to a central optimization problem with (19) as a complicating constraint would lead to an unsymmetrical decomposition, as this single equation is assigned either to Zone a or Zone b.Furthermore, representing the exact same Equation ( 19) twice in the central problem would lead to an undesired rank deficiency of the KKT matrix as explained in Assumption 2. This can be overcome by replacing (19) A comparison of these Equations (20) with the initial form (19) shows that the pressure potentials p i of the nodes in the bordering zone have been replaced by the constant predefined pressure value p pre .As p pre is equal in both equations of (20), every solution fulfilling (20) will also satisfy Equation (19), since (20) is a special case of (19).The benefits of this new formulation are that the two equations in (20) are naturally assigned to their respective zones without creating further complicating constraints and thus reducing the coupling of the subproblems once more.Additionally, this procedure yields symmetrically decomposed subproblems.Further, the directions of the border edges are slightly varied as explained in the following remark.Remark 4. The decomposition only becomes fully symmetrical in the resulting border node Equations ( 21), ( 26), ( 29) and (31) outlaid below, by changing the edge directions of the border edges.As shown in Figure 1, edge e 1 is directed away from node i 1 and i 2 , leading to A i 1 ,e 1 = −1 and A i 2 ,e 1 = −1.Analogously, e 3 is directed towards i 3 and i 4 .Visually speaking, every zone calculates its solution with its own border edge directions, where the supply network border edge is directed away from the border node of the respective zone and the return network edge is directed to the proper border node of this zone.This is a deviation from the standard notation of the directed graph spanning up the DHN, where every edge is directed away from the start node and towards the end node.Still, the resulting system of equations is equivalent, as the differential pressure over the border edges ∆p e is defined using absolute values of A i,e ; see (23).

Distributed Optimal Dispatch of DHNs Based on OCD
The following shows how optimal operation problems for DHNs can be solved in parallel by applying the OCD approach to the reformulated central optimization problem provided above within this Section 3.
The presented approach is based on the following assumptions, used for simplicity reasons: Assumption 3. Pipelines with VMFDs are the only component types for which edges can become border edges.These pipelines are not part of a mesh in the return or supply network.Furthermore, no control paths reach across zone borders.
This is a valid assumption, as pipelines are the most common components in DHNs; thus, the limitations in finding possible decompositions for this assumption are manageable.Furthermore, as discussed after the upcoming assumption, the edges over which the approach decomposes the DHN can also be virtual pipelines.Assumption 4. The pipelines representing the border edges are regarded as negligibly short, and thus transmission delay, temperature losses and pressure potential differences due to height differences ∆h (leading to β > 0 in the hydraulic pipeline model) are neglected here, as ∆h = 0 is assumed.
If this assumption would lead to unwanted deviations, due to a non negligible length of the pipeline, two artificial nodes and pipelines, one in the return and one in the supply network, could be introduced into the model in one of the zones at the end of the current border edge.Thereby, all parameters of the old border edge pipelines are kept the same, and the newly inserted pipelines, which are set as the new border edges, have the same parameters except for the pipeline length L, which is chosen as negligibly short.
For notational brevity, only the optimization problem of Zone a is fully stated, while the NLP of Zone b is directly evolvable thereof.All equations for every area, except the ones replaced by the complicating constraints and the new objectives stated here, stay the same as given in the previous subsections of Section 3.
The new hydraulic node equation, replacing (14) for supply network border nodes in Zone a, is given as (the extensive superscripts of h stand for: district heating network, hydraulic, supply network, border node and Zone a): The main difference here to ( 14) is that border edges e ∈ S be e and "normal" edges e ∈ S e , e / ∈ S be e are distinguished.Border edge mass flows can then be expressed in dependence of pressure potentials p as stated below.Thereby, taking into account Assumption 4, which entails ∆h = 0 and thus β e,k = 0 in Equation ( 11 The dependency of the pressure coefficient on the border edge mass flow µ e,k ( ṁe,k ) is eliminated by using a predefined pressure coefficient µ pre e,k .Alternatively, this pressure coefficient can be precalculated, based on the border edge mass flow of the previous optimization and ( 11), thus, using a similar procedure as effectively applied for the thermal pipeline model.Further, the differential pressure ∆p e,k over the considered border edge e can be replaced by the difference of the pressure potential of the two connected border nodes i, from which one is in Zone a and the other in Zone b.Thus, the differential pressure over the border edge connected a border node i of Zone a can be written as: By inserting ( 23) in ( 22) and exchanging the variable µ e,k ( ṁe,k ) by the parameter µ pre e,k , the mass flow on the border edge e itself is defined as: For the network shown in Figure 1, Equation ( 23) is given by ∆p e 1 ,k = p i 1 ,k − p i 2 ,k for Edge e 1 in Zone a. Based on this, the border edge mass flow ṁe 1 ,k for Zone a stemming from Equation (24) results in: The thermal node equation, replacing (15) for border nodes, is formulated as follows (note that the border nodes are written separately to simplify the derivation of ( 31)): Similarly to (21), the border mass flows are also defined by ( 24) in (26) above.With Assumption 4, supposing negligibly short pipelines between zones, the thermal pipeline model reduces to an identity function mapping the input temperature to the output temper-ature here, T i,k = T i ,k (however, more detailed models, such as a static pipeline model or a symmetric dynamic thermal pipeline model, could also be considered in the thermal node Equation ( 26)).For Node i 1 in Zone a in the DHN in Figure 1, Equation ( 26) is written as The resulting objective function for Zone a f za , is defined as given below: This objective function f za is composed of two parts.The first describes the decomposed objective of a central EMS f (x).This comprises the objective f (x) resulting from a modification term δ • I with a small parameter δ and the identity matrix I added to the Hessian of the Lagrangian of the centralized and the distributed optimization problems to enhance convergence properties [33] (p.574).The objective resulting from this modification is defined here as f rest,za for Zone a.The second part is constituted by the sum of the products of the Lagrange multipliers λ and the complicating constraints h coming from Zone b.
As indicated by the superscripts, the Lagrange multipliers come from Zone b of the last OCD iteration from the complicating constraints h hydr,sn,bi,zb i,k and h therm,bi,zb i,k of the respective border nodes.The complicating constraints itself are utilized in the slight modified form given below in ( 29)- (31) with the additional obj superscript.In comparison to (22), ( 24) and ( 26), all variables become OCD parameters and vice versa.Furthermore, all sets that were based on Zone a now become dependent on Zone b and the other way around.This leads to the following definitions: where the border edge mass flow is defined similarly as in ( 24) by: The thermal node equation is formulated as follows for border nodes: Therein, the border mass flows are defined by (30).
The overall subproblem for Zone a then comprises the objective function ( 28), the complicating constraints ( 21) and ( 26) for the respective border nodes of Zone a, the equality constraints listed in Section 3.1.2for all other network components in Zone a, as well as box constraints limiting the solution set.

Results
Two examples demonstrating the applicability range of the described method to the distributed optimal operation of DHNs are considered.The eight-node network, as depicted in Figure 2, extends the four-node network introduced in Section 3.2 and depicted in Figure 1, by modeling the thermal properties of the pipes dynamically.Each example network was solved in both a centralized and a distributed manner using OCD.We found that distributed and central optimization problems converged to identical solutions in both considered examples.However, the OCD convergence properties for the eight-node network is highly dependent on the used parameter set.The central problem was solved using the Interior Point Optimizer (IPOPT) [34].The steps of the OCD procedure were implemented as optimization problems with an iteration limit of one, calling IPOPT once per step and zone.

Parameterization
The the values and bounds of the variables used in the following examples are given in Table 1.The Lagrangian modification parameter δ ensures a non-singular Hessian of the Lagrangian by adding δI to the aforementioned.In addition to the parameters given in Table 1, the eight-node network considers identical pipes of diameter 100 mm with a surface roughness of 0.4 mm.The pipe length is set to 50 m.
All parameters are equal for all times steps k.Therefore, it is sufficient to consider a single time step in the four-node network; however, the dynamic elements in the eight-node network warrant the consideration of multiple time steps.The following considers six time steps of 5 min each.

Dimensionality Reduction
The four-node network problem optimizes a total of 28 variables, the eight-node network a total of 432.Thus, the results must be presented in aggregated form.This subsection introduces the methods used to do so.

Maximum Infeasibility
Given the optimization variables x and the the active constraints h(x) = 0, the maximum infeasibility at x is defined as the ∞ -norm of h(x).Note, that for any v ∈ R N , its ∞ -norm is given by v ∞ = max{v 1 , . . ., v N }.By definition a solution to any optimization problem must satisfy its constraints corresponding to a maximum infeasibility of zero.In practice, a value of zero is difficult to obtain numerically; thus, the following will be content with small values.An initial point may, and in the following does, violate the constraints yielding a non-zero maximum infeasibility.

Representative Pressures
Note that (20) allows the definition of p i 3 ,k and p i 4 ,k in terms of p i 1 ,k and p i 2 ,k and the combination of ( 14) and (11) gives an injective relationship between the nodal pressures p i,k and the mass flow rates ṁe,k .Therefore, p i 1 ,k and p i 2 ,k are sufficient for the study of any feasible hydraulic state of either example network without loss of generality.As shown in Figure 3, intermediate states may be infeasible in some steps of the iteration process., produces values on a similar order of magnitude independent of network size.The normalization constants correspond to the ranges set for the respective variables as defined in Table 1b, ensuring that individual variables will not typically dominate the value of the MSE and guaranteeing unit-less results.Note that this definition does neither assume nor require feasibility.

Initialization
The four-node network was solved with all optimization variables initially set to zero.This flat initialization ignores any constraints placed on the variables.Convergence to the central solution was reached after 28 iterations.If the hydraulic equations contained in the model are solved beforehand in a so called hydraulic pre-calculation the obtained starting point can reduce the number of iterations needed.In the four-node network example this reduces the number of iterations to 25.The eight-node network requires a hydraulic pre-calculation as it does not converge using flat initialization, which is known in the literature, see, e.g., [4] for similar observations.

Four Node Network with Flat Initialization
Figure 4a shows the nodal pressures p 1,1 and p 2,1 in the four-node network, and therefore it is in a hydraulic state.The dashed line indicates the central solution.As depicted, the distributed solution converges to the central solution in 32 iterations.IPOPT took 28 iterations to solve the central problem.Figure 4b shows the thermal powers of the consumers and producers, demonstrating that convergence is not isolated to the hydraulic variables.Figure 4c aggregates all relevant variables according to (32).
As shown in Figure 3, the initial point is infeasible.Throughout the OCD steps, the maximum infeasibility tends to decrease.

Computation Time
An Intel i7-10700T CPU was used at a clock speed of 4.5 GHz to test the run time of the given examples.Table 2 summarizes the obtained results.The time given is the total time spend evaluating objective functions, constraints, their derivatives and completing optimization steps in IPOPT but excluding any boilerplate code for setting up, progress monitoring and switching zones.
In the OCD approach the zones where optimized sequentially to avoid zones influencing each others performance by accessing the same cash but the time given in Table 2 assumes parallel execution.The computation time depends on the specifics of the machine used for testing, however the ratio of the run times for the OCD and central approaches as given in Table 2 stays ruffly constant across differing hardware and is therefore more meaningful than absolute computation times.Knowledge of the central solution was used to terminate the OCD iteration as soon as the MSE fell below 3 × 10 −10 .As shown in Table 2 the OCD solution is significantly slower as the central solution.However the advantage of the central solution is smaller for the larger example considered.This matches the results obtained in [22], where OCD becomes much faster as the central solution for larger networks.Hence, this leads to the assumption that OCD will be faster for a sufficiently large network.

Other Solutions
With other parameter sets, the decomposed problem was observed to converge to other solutions.For example, increasing the marginal bid price of the consumer on e 4 from 5 to 6 produces the results shown in Figure 6 and caused the OCD approach to converge to a solution distinct from that of the central approach.The MSE was found to plateau after 100 iterations at a value of 0.08.The final objective value of the distributed solution was 0.35 % better than in the central case; however, while the central solution solves all constraints to an infeasibility of at least 1 × 10 −13 the distributed solution shows maximum infeasibilites in the order of 1 × 10 −4 .
The results of further increasing the marginal bid price of the consumer e 4 to 12 are shown in Figure 7 and causes the OCD approach to fail as no convergence is observed with even the best set of values reached through out the OCD procedure violating the constraints by an infeasibility of at least 1 × 10 −2 .After iteration 130 the MSE oscillated around a value of 0.075 corresponding to a 4.9 % better objective value compared to the central case.
While the OCD approach converged to the central solution for the prevailing amount of parameters tested in the four-node network, the distributed solution of the eight-node network only converged toward the central solution in a minority of the considered parameter sets.A potential explanation for this result is that the eight-node network case, which in contrast to the four-node network includes the thermal pipeline models, Condition 1 might be violated regularly.

Discussion
The presented approach shows how OCD can be applied to the optimal operation of DHNs.As explained, several difficulties need to be overcome to enable the application of OCD to DHNs.This includes the strong coupling of the subproblems resulting from the symmetric flow conditions in supply and return networks in DHNs.While the concept shows good convergence properties for the four-node network for the majority of the tested parameter sets, the results indicate that further work is needed to enhance the convergence of the distributed solution towards the central solution for the regarded dynamic eight-node network and thus likely also for networks with even more components.Hence, the presented ideas to circumvent the strong coupling of the subproblems appearing from the application of the OCD method to DHNs may need to be further extended for large-scale DHNs.

Figure 1 .
Figure 1.Example: four-node network.DHN with four nodes decomposed into Zones a and b.

8 Figure 2 .
Figure 2. Example: eight-node network.DHN with eight nodes and dynamic elements decomposed into Zones a and b.

Figure 3 .
Figure 3. Maximum infeasibility per zone when solving the four-node network using flat initialization.

Figure 4 .
Figure 4. Four-node network from Figure 1 solved from an initial point with all variables set to zero.(a) Pressures at nodes one and two; (b) Thermal power output/consumption of consumers and producers; (c) Mean Square Error as defined in (32).

4. 5 .Figure 5 .
Figure5shows that the eight-node network example converges in 126 iterations.IPOPT required 69 iterations to solve the central problem.

Figure 6 .Figure 7 .
Figure 6.Eight-node network with the marginal bid price c e 4 increased to 6. Initialization based on Hydraulic Pre-Calculation.(a) MSE as defined in (32); (b) Maximum infeasability per zone.

Table 1 .
The parameters used for simulation.The values are equal for all time steps k ∈ S k .
4.2.3.Mean Square ErrorThe Mean Square Error (MSE) was used as a metric for measuring the difference of a central and distributed solution unifying all relevant variables in a single metric according to the definition where ∆ denotes the difference of the central and distributed solutions in the corresponding variable.Any edges that are neither consumers nor producers do not define a thermal power Φ e,k and are therefor ignored.Division by the number of variables summed in n, with n = |S k | |S e | + 2|S i | + S exch e

Table 2 .
Computation time and iteration counts for the considered examples.