Time-Series Temperature-Dependent Power Flow Considering Unbalanced Thermoelectric Equivalent Circuits for Underground LV and MV Cables

—This manuscript proposes a time-series temperature-dependent power flow method for unbalanced distribution networks consisting of underground cables. A thermal circuit model for unbalanced three-phase multi-core cables is developed to estimate the conductor temperature and resistance of Medium and Low Voltage distribution networks. More specifically, a novel approach is proposed to model and estimate the parameters of the three-phase thermal circuit of 3/4-core cables, using the results of Finite Element Method and Particle Swarm Optimization. The proposed approach is generic and can be accurately applied to any kind of 3- or 4-core cables buried in homogeneous or non-homogeneous soil. Furthermore, it is applicable in cases where one or more adjacent cables exist. Using the proposed approach, the conductor temperature of each phase can be individually and precisely calculated even in networks with highly unbalanced loads. The proposed approach is expected to be an important tool for simulating the steady state of unbalanced distribution networks and estimating the conductor temperatures. The proposed thermal circuit is validated using two 4-core LV and one 3-core MV cables buried in different depths in homogeneous or non-homogeneous soil. Time-series power flow for a whole year is performed in a 25-bus unbalanced LV network consisting of multicore underground cables.


I. INTRODUCTION
OWER flow is a fundamental tool for several power system applications such as planning, state estimation, real time control and optimization, contingency and stability analysis, etc. Conventional power flow is based on the assumption of a fixed conductor resistance, i.e., the line resistances have pre-specified and invariable values, without considering their dependency on the conductor temperature. Therefore, conventional power flow approaches present inaccuracies, due to the neglect of continuous variation of conductor resistances.
Nowadays and 41% of all distribution lines between 1 kV to 100 kV, in Europe, are underground [1]. Looking into the future, more and more EU member states are moving towards underground electricity distribution networks [1]. In addition, a high uptake of Distributed Generators (DGs) and Electric Vehicles (EVs) is also observed, which is possibly going to force distribution networks to their limits [2]. As a result, accurate power flow modeling inclusive of the thermal characteristics of the unbalanced underground LV and MV distribution grid is necessary for the planning, analysis, operation, control, and supervision of these networks.

A. Challenges of modeling underground MV/LV cables
Overhead lines present low thermal inertia reaching their steady state temperature in a few minutes [3], [4]. As a result, the conductor temperature can be estimated based on the conductor current at the present time instant. On the contrary, underground cables have a much higher thermal inertia, due to the large specific heat constants of the insulation, shield, screen, sheath, and jacket of the cable as well as of the soil [4]. As a result, the historical thermal condition of the cable and soil have to be considered for the estimation of the conductor temperature, which increases the complexity of the calculation process.
Another challenge is the modeling of non-homogeneous soils, which is usually the case in installations where a backfill material exists with different thermal resistivity and specific heat constant than those of native soil [5]. IEC and IEEE standards provide formulas for the modelling of soil, but they are only applicable to a limited number of installation geometries [6], [7]. In practice, several layers of different thermal resistivities and specific heat constants may be present between the cable and the ground surface. In such cases, standard methods cannot be applied [5].
The complexity of modeling underground MV/LV cables is even further increased in cases where multiple cables are buried next to the examined cable, due to the mutual heating between the cables [8]. More specifically, the heat produced by a buried cable can drastically affect the temperature of an adjacent cable.
Moreover, MV and LV cables are often 3-core or 4-core cables, which are asymmetrically loaded in unbalanced distribution networks. The unbalances are expected to increase in the near future due to the connection of single phase DGs and EVs [9]. As a result, the conductors of the cable exhibit different temperatures, which should be individually calculated considering the load unbalances and not P simplistically assuming that the three phases are equally loaded, as typically occurs in the existing literature so far [6], [7], [13].

B. Literature review
Temperature-dependent power flow, known also as ElectroThermal Coordination (ETC), has been introduced by several research groups to more accurately model distribution and transmission networks and exploit their maximum power transfer capacity by considering the thermal inertia of conductors and cables.
ETC for cable based transmission grids was introduced in [10], where the authors divided the cable and its environment into a large number of zones, and expressed them into a Thermoelectric Equivalent Circuit (TEE). The TEE interacts with a power flow tool in such a way that the current sources of TEE are updated based on the results of power flow, while the parameters of power flow, i.e., the cable resistances, are updated according to the results of TEE. The authors in [11] modelled underground cables, when real-time measurements of the temperature at the surface of the cable is available. The TEE of [10], [11] was improved by the authors in [12], by reducing the number of thermal zones of the soil, and therefore, the number of components of TEE, without sacrificing the accuracy of the model. In [8], the mutual heating effects caused by neighboring cables were introduced into the TEE, increasing its accuracy in case of multiple buried cables.
Authors in [13] implement Optimal Power Flow (OPF) to minimize the DG curtailment due to the thermal violation of underground cables, transformer and overhead conductors of the network. The power flow is combined with a TEE in a similar sense as in [10] to couple the current calculated from power flow with the conductor temperature calculated by TEE. The TEE is constructed using the approach of [14], in which the cable and soil were divided into a total of 4 zones using formulas from IEC standard [6] and the Van Wormer coeficients [15].
Finally, authors in [4] propose a single-phase time-series Newton-Raphson (NR) power flow method for transmission network considering the TEE and the different thermal inertia constants of transformer, overhead lines and underground cables.

C. Technical contribution of this paper
All the aforementioned works present considerable limitations when implemented in unbalanced LV and MV distribution networks. Firstly, they use the mathematical formulas of IEC standard for calculating the thermal resistance and capacitance of the soil in the TEE. However, these formulas are only applicable to a limited number of installation geometries. When the medium surrounding the cable is composed of materials with different thermal resistivities and specific heat constants, which is usually the case in installations with a non-homogeneous soil or backfill material, standard methods cannot be applied [5]. Furthermore, the aforementioned works assume either singlecore or 3-core cables with balanced three-phase loading. In reality, a large portion of MV and LV cables are 3-core or 4core cables with a highly unbalanced loading in each conductor, and thus, the temperatures of the cores can significantly deviate from each other.
In this manuscript, we propose a temperature-dependent power flow method for underground distribution networks with the following distinct features: • We propose a TEE model for the thermal simulation of the underground cables with dynamic conductor states, the parameters of which are estimated and validated based on the results of finite element method (FEM) software. Therefore, it is generic and applicable to every installation geometry and type of soil. • Our proposed model also considers the mutual heating caused by adjacent cable installations.
• Our model accounts for the unbalanced loading of each phase in multicore cables, calculating individually and accurately the temperature of each conductor.
• The power flow is based on the method of [16], [17], which is applicable to AC and hybrid AC/DC networks operated in both islanded and grid-connected mode as well as in networks with different voltage levels [18].

D. Paper structure
Section II of this manuscript describes shortly the power flow solver employed. Section III explains the proposed approach for modeling the thermal behaviour of unbalanced underground multicore cables. Section IV validates the proposed TEE model of Section III. Section V presents a detailed step-by-step solution procedure of the proposed power flow approach, which accounts for the thermal condition of the cables, while section VI present a time-series power flow study over a whole year using the proposed power flow approach. Finally, section VII concludes the paper.

II. CONVENTIONAL POWER FLOW ALGORITHM
The full configuration of an unbalanced LV network is presented in Fig. 1, including the neutral conductor and the grounding resistances. Let us define the current and voltage vectors of node i as follows: where and denote the load current and voltage (in complex form) of node i at phase r = {a, b, c} and conductor y = {a, b, c, n, g}, respectively, as shown in Fig. 1.
For a network with m nodes, the load vectors can be expressed as a function of the voltage vectors according to Eq. (3): where coresponds to the slack node. Subsequently, for the power flow solution, we remove the first five rows of Eq. (3) that correspond to the slack node, and Eq. (4) is obtained.
Then, by transferring all the voltage variables of the lefthand side of Eq. (4) to the right-hand side, Eq. (5) is derived, where and are the modified current and admittance matrices.
As a last step, we define the final matrices ′ and . The first one consists of the first five columns of , while the second one consists of the remaining columns such that = [ ′ ]. Eq. (6) is then derived from Eq. (5) by subtracting the product ′ from both equation sides.
Using Eq. (6), we finally derive Eq. (7), which is iteratively solved until a certain preset tolerance is reached. In Eq. (7), k denotes the iteration number, and the vector contains the voltages of all nodes except the slack.

A) Introduction of the proposed thermoelectric-equivalent circuit
In this section, we propose a TEE circuit model for modeling a 4-core underground cable as illustraded in Fig. 2. Ιt consists of the examined 4-core underground cable, burried near to an adjacent thermal body. The soil is modeled through a set of thermal resistances and capacitances between the surface of the investigated cable and the surface of the ground ( , 1 , 2 , 3 , 4 , 6 ), a set of thermal resistances and capacitances between the surface of the investigated cable and the surface of the adjacent cable ( , 1 , 2 , 3 , 4 , 6 ) and a set of thermal resistances and capacitances between the surface of the adjacent cable and the surface of the ground ( , 1 , 2 , 3 , 4 , 6 ), as depicted in Fig 2. Inside the cable, four current sources exist ( , , , ), at the position of the conductors, which denote their Joule losses due to phase current flows. The heat path between the conductors is denoted with a set of thermal resistance and capacitance ( , ). Finally, the heat path between the conductor and the surface of the cable is denoted with a set of thermal resistances and capacitance ( , 1 , 2 ), which represent the various layers of the cable. It is noted that the parameters , , , 1 , 2 are the same in all conductors, due to the geometrical symmetry of the cable, which implies that the paths of the heat of the four conductors inside the cable will have equal thermal properties.
In order to calculate the parameters of our TEE circuit, we model a 4-core underground cable in a FEM software. The thermal resistances of the TEE circuit of Fig. 2 are calculated based only on the results of FEM simulation, while the capacitances are calculated using the results of finite element method (FEM) simulation and particle swarm optimization (PSO), as explained below.

B) Calculation of thermal resistances
For the calculation of thermal resistances the following simulatios are executed in FEM: Simulation 1: Calculation of , , . We load the four conductors with a current of equal magnitude, so that the heat flowing through the thermal resistances is zero. More specifically, the four conductors are equally loaded with a current ∠0 o , ∠90 o , ∠−90 o , ∠180 o , respectively, so that they generate an equal heat . Furthermore, the adjacent cable is loaded with a suitable current so that the temperatures at the surfaces of the two cables is equal and no heat is flowing through the resistances . The total thermal resistance from the conductor to the surface of the cable (3 ) is calculated by Eq. (8): where is the temperature of the conductor, while _ is the temperature at the surface of the investigated cable, both calculated and measured in the FEM software. The total thermal resistance from the surface of the investigated cable to the surface of the ground (6 ) where _ is the temperature at the surface of the ground, while 4 denotes the total heat produced by the equally loaded four conductors of the cable.
The total thermal resistance from the surface of the adjacent cable to the surface of the ground (6 ) is calculated by Eq.
where _ is the temperature at the surface of the adjacent cable, while denotes the total heat produced by the adjacent cable, both measured by the FEM simulation software.
Simulation 2: Calculation of . We load the conductors , , and with ∠0 o , 2 ∠180 o and ∠0 o , respectively. Moreover, the adjacent cable is loaded with a suitable current so that the temperatures at the surfaces of the two cables are equal and no heat is flowing through the resistances . The heat produced by the conductor of phase flows through the thermal resistances towards the conductors and as well as through the resistances towards the surface of the cable according to Eq. (11).
where , , denote the heat produced by the conductor , the temperatures of conductors , , , respectively, all measured in the FEM simulation software, while is known from Simulation 1 above. is easily calculated from Eq. (11).
where , _ , _ , _ are all measured in the FEM simulation software, while has been calculated in Simulation 1. Therefore, is easily calculated from Eq. (12).

C) Calculation of thermal capacitances
As shown in the TEE of Fig. 2, there is a total of 23 unknown thermal capacitances. Το calculate their values, we perform time-series simulations in FEM software by setting specific current values for the conductors of the cable in specified time periods as shown in Tables I, II, III, and IV. A time-varying waveform for the conductor temperatures is obtained such as the one shown in Figure 8, 9 , 10, and 11. The thermal capacitances are estimated so that the response of the proposed TEE circuit of Fig. 2 fits the response of the FEM simulation software for the same conductor current values. The estimates are realized using PSO, which is a heuristic optimization approach with global optimum searchability [19]. Although PSO presents an increased execution time, it is executed offline only once for every cable installation and not in every load variation.
The thermal circuit of Fig. 2 is mathematically expressed in a form of state-space equations, as follows: where ( ) = [ 1 ( ) 2 ( ) … 34 ( )] is a vector including the temperatures of each thermal capacitance. = where , , , denote the heat produced by the four conductors of the investigated cable, is the heat produced by the adjacent cable, while _ is the temperature at the surface of the ground. The vector is the input of the thermal circuit of Fig. 2, which determines its response. The thermal resistances and capacitances are included in the matrices , .
Assuming that the states 1 ( ), 2 ( ), 3 ( ), 4 ( ) of the vector denote the temperatures of conductors , , , , respectively, the following equation expresses the objective function that is minimized with PSO.
where is the total simulation time by the FEM software. ( ) for = { , , , } is the temperature of conductor of the examined cable at time instant calculated by the FEM software. The optimal variables are the 22 unknown thermal capacitances. It is reminded that in Fig. 2, some nodes inside the cable have equal capacitance values, due to the geometrical symmetry of the cable. In practice, Eq. (14) searches for the opimal thermal capacitance values in order to minimize the deviation between the conductor temperatures obtained from FEM simulation and the proposed TEE circuit of Fig. 2.

IV. VALIDATION OF THE PROPOSED TEE CIRCUIT MODEL
The proposed TEE circuit model in Section III is validated here, using four different underground cable installations. The validation is performed by comparing the results of the proposed TEE circuit model against FEM simulation.

A) Investigated installations
Four different underground cable installations are simulated in the FEM software as explained below: • Installation 1: A 4-core LV cable is buried at a depth of 0.5 m underground inside backfill material as presented in Fig. 3. The thermal resistivity and specific heat constant of the backfill are 0.45 Km/W and 850 J/(kgK), while for the native soil they are 1 Km/W and 1500 J/(kgK). The density of the backfill and native soil are 1700 kg/ 3 and 2600 kg/ 3 , respectively [11]. The 4-core cable was taken from [21, page 26] with a nominal cross-sectional area 35 2 . All conductor cores have the same cross-sectional area. An adjacent cable is installed at a distance of 25 cm. • Installation 2: A 4-core LV cable is buried at a depth of 1.5 m underground. The ground is composed of two layers of soils having different thermal properties as presented in Fig.  4. The thermal resistivity, specific heat constant, and density of the top layer are 0.45 Km/W, 850 J/(kgK), and 1700 kg/ 3 , respectively, while for the bottom layer they are 1 Km/W, 1500 J/(kgK), and 2600 kg/ 3 [11]. The 4-core cable is the same as in Installation 1. All conductor cores have the same cross-sectional area. An adjacent cable is installed at a distance of 1m. • Installation 3: The 3-core MV cable of Fig. 6 is installed in homogeneous soil with a thermal resistivity, specific heat constant and density value equal to 1 Km/W, 1500 J/(kgK) and 2600 kg/ 3 [11], respectively. The cable is buried in 2 m as presented in Fig. 5. • Installation 4: Similar to Installation 1 but the 4-core LV cable has a reduced neutral conductor. More specifically, the phases have a cross-sectional area 35 2 , while the neutral 16 2 (see page 26 of ref. [21]).

B) Validation results
The validation of the proposed TEE circuit model is carried out by comparing the dynamic response of the proposed TEE circuit against the FEM software. The conductors of the cables of the aforementioned 4 installations are loaded over time, as shown in Tables I-IV. In most time periods the three-phases are unbalanced loaded, which is usually the case in MV and LV networks.
The conductor temperatures calculated by FEM software and the proposed TEE circuit for the four installations are depicted in Figs. [8][9][10][11]. As shown, the proposed TEE circuit is in good agreement with the simulations from the FEM software over the whole investigated time period, confirming the accuracy of the proposed approach. More specifically for Installation 1, the maximum difference between the two models is 0.5 o C at the 2 nd hour. In Installation 2, the maximum difference is 0.4 o C at the 7 th hour, in Installation 3 it is 0.6 o C at the 45 th hour, while in Installation 4 it is 0.6 o C at the 3 rd hour.
Another important thing to note is the large deviation between the temperature of the conductors, in periods where they are unbalanced loaded. For instance, in Installation 1, the temperature of the conductors , , , at the 7 th hour is 80.52 o C, 90.73 o C, 79.57 o C, 71.1 o C, respectively. As a result, the assumption of a balanced loading and equal temperature between the conductors, which is typically considered so far in the thermal modeling of multi-core cables, could be significantly misleading in power system analysis.
It is pointed out that in Fig. 11, the neutral conductor has similar temperature with phase (74.44 o C and 74.24 o C, respectively at 3 nd hour), although the current through the neutral is significantly lower. This occurs due to the reduced cross-sectional area of the neutral in Installation 4.     (Table I) and underground Installation 1. Fig. 9. Conductor temperature for input current (Table II) and underground Installation 2.  (Table III) and underground Installation 3. Fig. 11. Conductor temperature for input current (Table IV) and underground Installation 4.

V. PROPOSED TEMPERATURE-DEPENDENT POWER FLOW ALGORITHM
Proceeding from the validation of the proposed TEE circuit model, we present the proposed temperature-dependent power flow algorithm in this section as follows: Step 1: Firstly, the TEE circuit parameters of Fig. 2 is calculated offline for each cable installation of the network, using the PSO based optimization proposed in Section III.
Step 2: After obtaining the parameters of the TEE in Step 1, we formulate Eq. (13), to mathematically express the TEE circuit.
Step 3: We discretize Eq. (13) so that it is expressed in the following form: In this manuscript, the discretization was executed in MATLAB with the c2dm command, because of MATLAB's fast and efficient inbuilt algorithms.
Step 4: Then, the iterative process begins by calculating the three-phase power flow using the algorithm described in Section II.
Step 5: From the currents obtained via the power flow solution in Step 4, we calculate the temperature of each conductor through Eq. (14). The vector is updated through the Joule's law, e.g.,

= 2
where I is the current calculated by the power flow in Step 4, while is the AC resistance of conductor .
Step 6: From the conductor temperatures calculated in Step 5, we update the conductor resistances, as follows: where is the resistance value at θ o C, 20 is the resistance value at 20 o C, is the temperature coefficient of resistivity, is the conductor temperature at θ o C, and 20 is the conductor temperature at 20 o C. This iterative process is updated from Step 4 until convergence.

VI. TIME-SERIES TEMPERATURE-DEPENDENT POWER FLOW SIMULATION
Time-series temperature-dependent power flows are executed here for a whole year in a 25-bus unbalanced LV network, using the proposed algorithm. The examined 25-Bus network is depicted in Fig. 12 [3]. Every bus (except bus 1) supplies three single-phase residential loads, one in each phase. Thus, the network supplies in total 72 single-phase customers. The network was assumed to operate in grid-connected mode although the proposed power flow algorithm can be applied in islanded operation as well [3].

A) Investigated Network
For the load, we used real measured data from [22], which were normalized, as shown in Fig. 13. Then each load connected to phases A, B, C is multiplied with a base load = 1700 + 1275, = 2700 + 2025, = 2200 + 1650, respectively, to form a highly unbalanced LV network. All loads were assumed to have a power factor equal to 0.8. A single-phase photovoltaic (PV) is connected to phase A of buses 13, 19, 25. Each PV has a base power of 5 kW multiplied by the normalized solar irradiance of Fig. 14. The solar irradiance is measured in the region of Sindos (a suburb of Thessaloniki, Greece) [23], for the whole year of 2018.
The distance between the buses was assumed 50 m, which is a typical distance in LV networks. The buses are connected with underground LV cables. Installation 4 (see section IV.A) was assumed for all the cables of the network without the presence of the adjacent cable. A grounding resistance (Zgri in Fig. 1) of 25 Ohm was assumed for all the buses of the network, which is a typical value for LV networks.

B) Power Flow Results and Analysis
The temperature of each conductor (three phases and neutral) for the cable connecting the buses 1-2 is depicted in Fig. 15. The smaller sub-figure on the top-left depicts the conductor temperature for the entire year. As expected, the temperature is proportional to the loading condition, which takes its maximum value in the summer period. The maximum conductor temperatures are observed on July 3. The corresponding conductor currents for this day are shown in Fig. 16. Fig. 13. Normalized residential load data [22]. As shown, the maximum peak temperatures of phases , , , and neutral are 87.61 o C, 92.62 o C, 89.65 o C and 82.25 o C, respectively. A temperature deviation of around 5 o C is observed between the conductors of phase and , which highlights the importance of considering the unbalanced loading in the thermal modeling of multicore cables. Moreover, the phase has slightly exceeded the maximum operating temperature limit of 90 o C referred in [21, page 18] for XLPE insulated cables, while the other phases remain within the limit. This indicates the excessive electro-thermal stress on conductor phase and, therefore, the importance of modelling multicore underground cables appropriately.
Finally, it should be noted from Fig. 16 that the conductor loading on July 3 is around 300 A, which is well above the static current rating of 158 A referred in [21, page 24]. Therefore, the proposed algorithm can constitute an accurate analysis tool for studying and utilizing the dynamic line rating of LV and MV multicore cables.
The resistance of each conductor of the cable connecting the buses 1-2 is depicted in Fig. 17 during the whole year. As shown, it significantly changes with the variation of conductor temperature. It is clarified that in this example, we used a resistivity temperature coefficient equal to 3.39•10 -3 K -1 (see α in Eq. (15)), which is a typical value for copper. The neutral conductor has a higher resistance due to its reduced crosssectional area. From this example, it is indicated that the proposed temperature-dependent power flow approach is more accurate than the conventional power flow, which usually assumes a constant resistance throughout the year.
Finally, the total power losses of the network over the whole year are calculated at 19.02 MWh using the proposed temperature-dependent power flow approach. The conventional power flow calculates it at 20.754 MWh assuming a constant resistance at 60 o C, resulting in a deviation of about 9 %.

C) Computation Time
The proposed approach needs on average 0.2 sec to solve the power flow for one time instant. In this example, we executed the power flow with one hour resolution for a whole year, i.e., 8760 time instants. Thus, the total computation time of the algorithm over the whole year is around 30 minutes. A PC with a 64-bit Intel Core i7, 3.4GHz CPU and 16GB RAM was used for all simulations.

VII. CONCLUSION
This paper proposes a time-series temperature-dependent three-phase power flow method for unbalanced distribution networks consisting of underground cables. A novel approach is proposed to estimate the parameters of the three-phase thermal circuit of each multicore cable, using the results of Finite Element Method (FEM) and Particle Swarm Optimization (PSO). The proposed approach is generic and can be accurately applied to any kind of 3-core or 4-core cables buried in homogeneous or non-homogeneous soil. It is also applicable in cases, where one or more adjacent heating bodies exist. Using the proposed approach, the conductor temperature of each phase can be separately and precisely calculated even in networks with highly unbalanced loads. The proposed approach can constitute an accurate tool for power flow analysis of unbalanced distribution networks and estimation of the conductor temperatures.