Three-Layer Framework Integrating Optimal Placement of SCADA Measurements with Clustering-Based Electric Substations Selection for State Estimation of Medium-Voltage Distribution Networks

1. Introduction

1.1. Motivation of the Research

Due to the increasing need for distribution automation, the monitoring and control concepts integrated into the modern electric distribution networks (EDNs) have significantly changed from manual to digital in the past few years [1]. Such a challenge has led to the accelerated transition of the smart SCADA (Supervisory Control and Data Acquisition) system, a computer-based system to monitor and control the various components of the EDNs.

During the late 20th century, the EDNs started incorporating the SCADA system. This type of system comprises various components, such as sensors, communication units, and programmable logic controllers. The EDNs incorporate multiple problems, such as power interruptions caused by weather conditions and equipment failure. A SCADA system is a tool to manage these networks to improve their efficiency and reliability. The functions of a SCADA system are carried out by collecting and analyzing data from various sources such as power devices, circuit breakers, and remote terminal units (RTUs). These data are then processed and sent to a Control Centre. When an outage occurs, the system sends an alarm to the operators [2].

Although most electric distribution networks can integrate the SCADA system, the Distribution Network Operators (DNOs) cannot fully incorporate all of them due to the high investment. This economic aspect could be the main reason for partially integrating the SCADA system into the electric distribution network. For example, Figure 1 presents the status of the integration in the SCADA system in the electric distribution substations (Medium Voltage/Low Voltage) at 45 DNOs from the European Union Countries given in [3]. The data analysis highlights an integration degree of 30 % at most DNOs (93%). Only 3 of 45 DNOs (7%) have exceeded this value.

The DNOs can use the SCADA system to monitor and control various system elements, such as power transformers, electric lines, etc. With the help of data acquisition capabilities, it can also provide significant insight into the operation of the power system. SCADA systems are commonly used on the distribution side to automate the various components of the networks. They contain the following in their structure: a remote terminal unit (RTU), a master control centre, and a telecommunication network. All European DNOs use the SCADA system, which includes one or more functions: substation control, feeder control, or end-user load control [4], see Figure 2.

Substation control. The SCADA system monitors the status of a substation's various equipment. It then sends control signals to remote-control units. It also collects historical data about the facility and provides alarms in case of faults or electrical faults. The survey indicates that about 96% of DSOs use a SCADA system to control their substations. These operators mainly use it to maintain bus voltage, load balancing, circulating current, overload protection, and transformer fault detection.

Feeder control. Automated switching systems can improve the reliability of electric distribution networks using sophisticated algorithms and control systems. In a survey conducted by DSOs, over 80% of them use a supervisory control and data acquisition (SCADA) system for their feeder control. They also reported using a similar system for voltage control and variable-rate electricity. Automated switching systems can improve electric distribution networks' operational flexibility and reliability. More than 80% reported using SCADA systems for automatic feeder switching, feeder voltage control, and VAR control.

End-user load control. Automation on the end-user side enables various functions, such as remote load control and automatic meter reading. It helps manage the power consumption of commercial, industrial, and residential customers. This type of automation can also detect the presence of energy theft. It can automatically disconnect the service and reconnect if the DNO has the issue disappeared.

Of the total DNOs representing the base of the study presented in the study [3], 70% utilized a SCADA system for end-user load control. Among the DNOs that use this function, 62% utilized a similar system for remote load control and generating bills and 76% for automatic meter reading.

Following best practices is essential for DNOs to improve their SCADA system integration. These practices help it identify the various components of the system that need to be integrated and establish a clear understanding of the EDNs' requirements. In addition to conducting a comprehensive analysis of the current infrastructure, these practices can also help implement in the optimal locations (EDSs) for its integration. But, without the proper mathematical tools, DNOs often face issues regarding the state estimate process of the EDNs. Today, various algorithms and methods are available to satisfy these requirements.

1.2. Literature Review

Determination of the optimal number of electric distribution substations (EDSs) integrated into the SCADA system to perform on-site measurements used in the state estimation process of the EDNs represents a challenge both from a technical and economic point of view for DNOs. It involves developing a robust framework incorporating algorithms to model the EDN's conditions accurately. This paragraph focuses on the various methods used in this process.

The representative load profiles are estimated using Adaptive Neuro-Fuzzy Inference Systems in reference [5] to analyze historical data. Different characteristics, such as the current value of the load, the previous 24-hour load, and the weekly system model factors, have been used to generate the best fit for a particular week or day. The proposed algorithm has been tested on a medium voltage EDN, leading to a MAPE value of around 4%. Xie et al. proposed in [6] a control and prevention method for operating a distribution network based on ultra-short-term load forecasting. The DNO then estimated the steady state of the network, and the risk assessment was performed. Finally, the risk assessment results were sorted and screened before identifying the most critical incidents. Benato et al. proposed in [7] a general expression used to estimate the various factors that affect the power demand of a distribution network during a restoration process. It also considers the relationship between the frequency behavior and the distributed generation.

Reference [8] proposed a decentralized load estimating method integrated as support in the real-time optimization of a distribution network's voltage and load characteristics. The process is based on local data and limited information from its neighbor areas, enabling the selection and partitioning of the load profiles. It aims to limit the communication needs while ensuring the accuracy of the estimates. Kong et al. [9] developed a dynamic estimation method used to monitor the operation of an active distributed power network. It considered the various features of the operation scenario of the distribution network and its distribution generation. The method contains two phases associated with an adaptive estimation and an integrated multi-model algorithm. The first phase is to improve and incorporate the two standard estimators into a framework. The second phase identifies the operation modes and produces joint estimates.

Wang et al. developed a novel method that can estimate the baseline load for residential electric loads in different conditions [10]. It uses a deep learning approach to select the appropriate customer baseline load based on the historical data of the customers and non-participants. This method can also benefit from the concurrent load data collected from the participants in the Demand Response program and non-participants. Reference [11] considered the various aspects of electric load, such as the control mode and the vehicle's capability to run on a network, producing a comprehensive estimate of the total energy consumption of plug-in electric vehicles. It then considers the aggregate demand of different load categories at the feeder-head. The authors from reference [12] developed an optimal prediction interval estimation method used for various information regarding the feeders. It can be applied using reinforcement learning, which enables the integration of two online tasks.

The hybrid approach has been proposed in [13] to provide a short-term load forecasting solution. It uses a weighted least squares state estimation, neural network, and adaptive Neuro-Fuzzy Inference system. Thus, the Decision Maker can determine the optimal ranges for the membership functions used by the fuzzy system. The proposed framework in [14] for the baseline load estimation used the Fuzzy C-Means features and the fuzzy membership matrix. Correlation theory allows [15] to develop fuzzy load models using an algorithm for a supply of HV/MV EDS used as a reference. The algorithm's starting point was a statistical analysis of the data collected from the MV/LV EDSs. The decision maker estimated the hourly loads in each MV/LV EDS based on the simple regression models established, which had HV/MV EDS as the reference.

Rong proposed in [16] a different method with three steps for loop distribution networks. He considered aspects regarding network reduction, state estimation, and load forecasting. The analysis results revealed that the average error has been around 3%. A hybrid approach proposed in [17] led to a lower value than 5% for MAPE. Ding et al. built neural network-based models used in the load forecasting process of the MV/LV substations but with a relatively large MAPE (10% approximately) [18]. The benchmark process utilized the multiple linear regression model. All approaches used the SCADA database to estimate the load in the EDSs from the EDNs. An unsupervised learning technique has been proposed in [19] for the profiling process to estimate the electric distribution system load profiles. The errors for an EDN with 34 EDSs had the values of about 5%.

All approaches supposed that the data comes from the SCADA system with all EDSs from the EDN integrated, which implies a 100% monitoring degree. This assumption contradicts the results presented in the study from [3], where the integration degree in the SCADA system of the MV/LV EDSs at the level of most DNOs is approximately 30%.

1.3. Original Contributions

In this context, the main contribution is designing and developing a three-layer framework for optimal placement of SCADA measurements with clustering-based EDS selection to estimate the state of the medium voltage distribution networks. Each layer integrates original approaches to solving the assumed objectives:

• The first layer allows the determination of the classes of the EDSs with similar features from the viewpoint of requested loads based on the K-means clustering algorithm.
• The second layer identifies the "candidate" classes and the pilot EDSs (representing the optimal solution) with the SCADA measurements placed. The optimal placement corresponds to the minimization of the estimation errors obtained using the multiple linear regression models between the EDSs from the classes not included in the set of the "candidate" classes and the pilot EDSs.
• The third layer allows the state estimation of the EDNs based on the load values measured in the pilot EDEs (with the SCADA system implemented) and the other EDSs obtained through the regression models. Also, the layer contains a module that verifies that it satisfies all the technical constraints, having integrated the functions to implement the strategies for optimal operation of the EDNs.

The results obtained in the case of a real urban EDN with 39 EDSs demonstrated that the proposed framework could have lower integration costs (the optimal solution corresponds to 12 of 39 EDSs integrated into the SCADA system) and substantially improve the state estimation of the EDN.

1.4. Paper Structure

The rest of the paper includes the following sections: Section 2 presents the characteristics associated with the proposed framework containing details on the mathematical tools used to fulfil the objectives of each layer, Section 3 consists of the case study where the tests on urban EDNs with 39 EDSs highlight the performances of the proposed approach, and Section 4 integrates the conclusions.

2. Materials and Methods

This section presents the information about the algorithms, methods and performance metrics included in the proposed three-layer iterative framework for the optimal placement of SCADA measurements with clustering-based EDS selection to estimate the state of the MV EDNs. Figure 3 illustrates the flow chart of the proposed framework.

The three-layer framework improves the solution by continuously updating the subsets of the decision variables. The details presented below regarding each layer include the developed mathematical tool.

2.1. Layer 1

Step 1. The information collected from the electric distribution substations (the active and reactive power profiles) following the measurement campaigns carried out by the DNOs represents input data. These profiles are recorded in the database and uploaded in the matrices [P] and [Q] having the sizes TxN_EDS, where T represents the period used in the measurement campaigns to record the data with various sampling steps (usually 60 minutes) and NEDS correspond with the number of the EDSs from the analysed EDN.

$$\left[P\right]=\left[\begin{array}{cccccc}{P}_{1,1}& P_{1,2}& \dots & P_{1,n}& \dots & P_{1,N_{EDS}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ P_{t,1}& P_{t,2}& \dots & P_{t,n}& \dots & P_{t,N_{EDS}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ P_{T,1}& P_{T,2}& & P_{T,n}& & P_{T,N_{EDS}}\end{array}\right]$$

$$\left[Q\right]=\left[\begin{array}{cccccc}{Q}_{1,1}& Q_{1,2}& \dots & Q_{1,n}& \dots & Q_{1,N_{EDS}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ Q_{t,1}& Q_{t,2}& \dots & Q_{t,n}& \dots & Q_{t,N_{EDS}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ Q_{T,1}& Q_{T,2}& & Q_{T,n}& & Q_{T,N_{EDS}}\end{array}\right]$$

(1)

Step 2. The correlation coefficients are calculated between each two EDSs (n and m, n ≠ m, n, m ∈ {N_EDS}) from the analysed EDN and recorded in the matrix [CC], having the size (N_EDSxN_EDS).

$$C{C}_{P_n,P_m}={\displaystyle \frac{Cov\left(P_n,P_m\right)}{{\sigma }_{P_m}\cdot {\sigma }_{P_n}}}\hspace{1em}\left(\forall \right)\hspace{0.17em}n,m\in \left\{N_{EDS}\right\}\hspace{0.17em};\hspace{1em}n\ne m$$

(2)

where Cov(P_m, P_n) represents the covariance between the active power profiles from the set of the EDSs, generically noted n and m; σ_Pn and σ_Pm - the standard deviation of active powers associated with the EDSs m and n recorded in the time interval T.

$${\sigma }_n=\sqrt{{\displaystyle \frac{1}{T}}{\displaystyle \sum _{t=1}^T{\left(P_{t,n}-{\overline{P}}_n\right)}^2}};\hspace{1em}{\overline{P}}_n={\displaystyle \frac{1}{T}}{\displaystyle \sum _{t=1}^n{P}_{t,n}}\hspace{1em}n\in \left\{N_{EDS}\right\}$$

(3)

Step 3. The classes corresponding to the EDSs with similar features of the correlation coefficients are identified based on a clustering process. The optimal partition in the classes will be determined using the K-Means clustering algorithm [20] and an internal test based on a performance indicator, namely the Silhouette Coefficient (SC) [21]. The clustering process includes the following phases:

• The NEDS vectors associated with the columns of the matrix [CC] should be integrated into K classes:

$$\left[CC\right]=\left[\begin{array}{cccccc}C{C}_1& C{C}_2& \cdots & C{C}_n& \cdots & C{C}_{N_{\mathrm{EDS}}}\end{array}\right]$$

(4)

$$C{C}_n={\left[\begin{array}{cccc}C{C}_{n,1}& C{C}_{n,2}& \cdots & C{C}_{n,N_{EDS}}\end{array}\right]}_t\hspace{0.17em},\hspace{1em}n=1,\mathrm{...},N_{EDS}$$

(5)

2.

Determination the maximum number in which the NEDS electric distribution substations can be distributed using the relation [22]:

$$K_{\max }=\sqrt{N_{EDS}}$$

(6)

3.

The vectors CC_n, n = 1, …, N_EDS, will be randomly assigned in the K classes of the EDSs.

4.

The centroids C_k, k = 1, …, K and K = 2, …, K_max, representing vectors with the size (N_EDSx1) are calculated.

$$C_k={\displaystyle \frac{1}{N_{EDS}}}{\displaystyle \sum _{j=1}^{n_k}C{C}_j^{(k)}}$$

(7)

where: n_k – the number of the EDSs from the classes k, k = 1, …, K and K = 2, …,K_max; CC_j(k) – the vector associated with an EDS, generically noted j inside the class k from partition K, K = 2, …, K_max.

5.

The repartition of the EDSs in in one of the K classes, K = 2, …, K_max, will be based on the minimization of an objective function OF having the following expression:

$$\min (OF)=\min \left({\displaystyle \sum _{k=1}^K{\displaystyle \sum _{j=1}^{n_k}{‖C{C}_j^{(k)}-C_k‖}^2}}\right)$$

(8)

6.

The positions of the C_k centroids are re-adjusting through their recalculation using relation (7). In the case when all vectors CC_n, n = 1, …, N_EDS, are considered and re-labelled, Step 5 is repeated.

7.

The silhouette coefficient for each partition K = 2, …, K_max, will be calculated using the formula [23]:

$${c_s}^{(K)}={\displaystyle \frac{1}{N_{EDS}}}{\displaystyle \sum _{i=1}^{N_{EDS}}\left(\frac{x_i-y_i}{\max \left\{x_i,y_i\right\}}\right)}\hspace{0.17em},\hspace{1em}K=2,\mathrm{...},K_{\max }$$

(9)

$$x_i={\displaystyle \frac{{\displaystyle \sum _{\begin{array}{l}j\in X\\ j\ne i\end{array}}dist(v_i,v_j)}}{n_X-1}}\hspace{0.33em};\hspace{1em}{y}_i=\underset{X\ne Y}{\min }{\displaystyle \frac{{\displaystyle \sum _{\begin{array}{l}s\in Y\\ j\ne i\end{array}}dist(v_i,v_j)}}{n_Y}}$$

(10)

where: n_X and n_Y – the number of the vectors CC belonging to the EDSs from the classes X and Y associated with the partition K, K = 2, …, K_max; dist(v_i, v_j) – the distance calculated between the vectors v_i and v_j from the matrix [CC].

8.

The value of the silhouette coefficient SC(k) for each partition K = 2, …, K_max is recorded in the vector [SC] and the maximum value is identified:

$$\left[c_s^{\max },p^{\max }\right]=\max \left\{\left[\mathrm{SC}\right]\right\}$$

(11)

where c_s^max represents the maximum value of silhouette coefficient; p_max is the position from the vector [SC] where the maximum value has been identified.

9.

Determination of the optimal partition containing K_opt classes:

$${\mathrm{K}}_{\mathrm{opt}}{=\mathrm{p}}^{\max }+1$$

(12)

2.2. Layer 2

Step 1. Determining the “candidate” classes (notated with k_c, k_c ∈ {K_opt}) for placement of SCADA measurements based on a ranking established according to correlation coefficients, in reverse order beginning with those which have the highest values. The Decision Maker will impose a threshold (CC_threshold) for an average value of CC to choose the “candidate” classes [24].

Step 2. Identifying the set of the pilot EDSs which will be integrated in the SCADA system belonging to each “candidate” class.

$$\left\{N_{pilot}^{EDS}\right\}=\left\{n_k\in k_c|C{C}_k>C{C}_{threshold}\hspace{0.17em}\right\}$$

(13)

where n_k represents the number of the EDSs from each “candidate” class k_c, , k_c ∈ {K_opt}.

The pilot EDSs will represent the optimal places where be performed the SCADA measurements such that the load flow calculations to be carried out in the medium voltage distribution networks.

Step 3. Building the accurate, robust, and interpretable multiple linear regression models to estimate the active powers in the EDSs from the others classes that were not included in the set of the “candidate” classes. The models will use the pilot EDSs as regressors identified at Step 2.

$$\begin{array}{l}{P}_{t,m}^{(e)}={\beta }_{0,m}+{\displaystyle \sum _{p=1}^{n_k}{\beta }_{p,m}{P}_{t,p}^{(r)}},\hspace{1em}t=1,\mathrm{...},T\\ m\in \left\{k_{nc}\right\};\hspace{0.17em}\hspace{0.17em}p\in \left\{N_{pilot}^{EDS}\right\};\hspace{0.33em}\left\{N_{EDS}\right\}=\left\{N_{pilot}^{EDS}\right\}\cup \left\{N_{k_{nc}}\right\}\end{array}$$

(14)

where: {N_knc} – the set of the EDSs non-included in the “candidate” classes; {N_pilot^EDS} – the set of the pilot EDSs with SCADA measurements integrated; {N_EDS} – the set of the EDSs from the analysed EDN; P_t,m^(e) - the value of the active power estimated in the EDS "m" from the set {N_knc}, in [kW]; P_t,p^(r) - the value of the active power measured in the pilot EDS “p” from the set {N_pilot^EDS}, in [kW]; β_0,m – the constant coefficient of the linear regression model for each EDS from the set {N_knc}; β_p,m – the coefficients corresponding with each regressor associated with the pilot EDS “p” from the set {N_pilot^EDS} in the linear regression model for each EDS from the set {N_knc}.

Step 4. Analyzing the performance of the regression models based on the following metrics: percentage error (PE), average percentage error (APE), and mean absolute percentage error (MAPE):

$$P{E}_{t,m}={\displaystyle \frac{\left|P_{t,m}^{(r)}-P_{t,m}^{(e)}\right|}{P_{t,m}^{(r)}}}\cdot 100,\hspace{0.17em}\hspace{0.33em}t=1,\mathrm{...},T,\hspace{0.33em}\hspace{0.17em}m\in \left\{N_{k_{nc}}\right\}\hspace{0.17em},\hspace{1em}\left[\%\right]$$

(15)

$$AP{E}_m={\displaystyle \frac{1}{T}}{\displaystyle \sum _{t=1}^{T}P{E}_{t,m}},\hspace{0.17em}\hspace{0.33em}m\in \left\{N_{k_{nc}}\right\}\hspace{0.17em},\hspace{1em}\left[\%\right]$$

(16)

$$MAPE={\displaystyle \frac{1}{N_{EDS}}}{\displaystyle \sum _{n=1}^{N_{EDS}}AP{E}_n},\hspace{1em}\left[\%\right]$$

(17)

2.3. Layer 3

Step 1. Determination of the hourly reactive powers using the active powers measured in the pilot EDSs and estimated in the EDSs non-integrated in the “candidate” classes and the regression models established in Step 2 from Layer 1.

$$Q_{t,m}^{(e)}=f(P_{t,m}^{(e)},{\beta }_m)\hspace{0.17em},\hspace{1em}m\in \left\{k_{nc}\right\}\hspace{0.17em},\hspace{1em}t=1,\dots ,T$$

(18)

$$Q_{t,p}^{(e)}=f(P_{t,p}^{(r)},{\beta }_p)\hspace{0.17em},\hspace{1em}p\in \left\{N_{pilot}^{EDS}\right\},\hspace{1em}t=1,\dots ,T$$

(19)

Step 2. Power Flow Calculation using the Newton-Raphson method to estimate the state of the EDN. The state variables recorded refer to the voltages on the MV side of the HV/MV EDS, the active and reactive power flows, the power/energy losses, and the active/reactive powers injected in the slack bus (the MV bus of the HV/MV EDS).

Step 3. Verifying the technical constraints regarding the voltages and loading of the lines. If there are violations of some limits, then the DNO can apply technical measures to bring the state values to the admissible limits provided in the performance standards or imposed by the manufacturers.

3. Results

A real urban MV EDN has been used to test the proposed framework. The END, supplied from an HV/MV (110/20 kV) EDS, contains three MV feeders with 39 MV/LV (20/0.4 kV) EDSs. Figure 4 presents the topology of the test MV EDN.

Table A1 from Appendix A includes each branch's length between two EDSs associated with each feeder. The EDN has the cross-section of the first branch (EDS – 1, EDS -14, and EDS – 29), with a size of 185 [mm²] (r₀ = 0.157 [Ω/km] and x₀ = 0.112 [Ω/km]), which is different from the other branches, equal to 150 [mm²] ((r₀ = 0.194 [Ω/km] and x₀ = 0.115 [Ω/km]). Also, the branch length is in the range [0.1 – 0.5] km, with a total length for the second feeder of 6.07 km higher than the other two feeders (3.36 km for the first feeder and 3.54 km for the third feeder, respectively). Figure 5 presents the synthesis of this information.

The rated power of the transformers, denoted with S_r, from the EDSs supplied by the three feeders is between 400 and 1000 kVA, see Table A2 from Appendix A. The total installed powers (calculated as the sum of the rated powers of the distribution transformers supplied by each feeder) are 7.15 MVA for the first feeder, 12.21 MVA for the second feeder, and 9.52 MVA for the third feeder, resulting in a total power of the EDN of 28.88 MVA, see Figure 6.

Step 1 of the first layer from the proposed methodology considers the building of a database with the active and reactive power profiles recorded in the matrix form, as seen in the relationships (1), based on a measurement campaign performed by the DNO during a week in the most loaded period.

Figure 7 and Figure 8 present the boxplots of the active and reactive powers associated with the measurement campaign.

A boxplot is a method that shows the locality, skewness, and spread groups of data through their respective quartiles. In addition to the box, there are also whiskers, which are lines that extend from the box to indicate variability outside the lower and upper limits of the range. These are referred to as box-and-whiskers and box-and-whiskers, respectively [25,26].

Table A3 and Table A4, from Appendix A, show the quartiles associated with the boxplot representation and two other statistical indicators (mean and standard deviation). The analysis of the information highlights that there are 14 EDSs with a high variation of the powers (9, 10, 12, 14, 17, 22, 27, 29, 30, 31, 33, 34, 38, and 39), 7 EDSs with slight variations (1, 2, 3, 4, 23, 24, and 25), and 18 EDSs with normal variations.

Step 2 of the first layer involves the calculation of the correlation coefficient between each of the two EDSs from the analysed EDN, which is recorded in the matrix [CC] and has a size [39x39]. Figure 9 presents the heatmap chart associated with matrix [CC] containing the values between -0.1 and 0.99. The vast majority of values are above 0.6. Still, there are also values close to 0 or even negative, which means that as one variable increases, the other decreases proportionally. Application of the K-means algorithm-based clustering process included in the last step of the first layer led to the classes corresponding to the EDSs with similar features of the correlation coefficients.

The maximum number of the clusters (named classes in the following) corresponded to a partition in 6 clusters. The internal test used the Silhouette Coefficient, and the algorithm identified the optimal number equal to 5 (for which the silhouette coefficient has the highest value, SC = 0.7168), see Figure 10. The representation of the silhouette coefficient for the optimal partition, K_opt = 5, is given in Figure 11.

Figure 12 highlights each class's percentage and the number of EDS it integrates. The analysis of the data indicates that the representative class is C5, containing 12 EDSs (30.77%), followed by classes C1 (10 EDSs, 25.64%), C2 (10 EDSs, 25.64%), C3 (5 EDSs, 12.82%), and C4 (2 EDSs, 5.13%).

The allocation of the EDSs in each class can be observed in Figure 13.

Two (Feeder 1 and Feeder 2) from three feeders have the EDSs associated with all classes, and one feeder (Feeder 3) contains the EDSs allocated to only three classes (C1, C2, and C5). Table 1 presents the statistical indicators corresponding to the quartiles (Q0 – Q4), mean (M), and standard deviation (SD) for the obtained classes following the clustering process.

The class with the values of the highest average correlation coefficient is C5, the range is [0.75, 0.82], and the smallest values are associated with C4 (0.4). The trend of the mean is the same, with the highest value recorded for class C5 (M = 0.79), followed by C1 (M = 0.74), C2 (M = 0.62), C3 (M = 0.62), and C4 (M = 0.4). The highest variation of the correlation coefficient belongs to class C2 (between 0.57 and 0.68), followed by C3 (between 0.59 and 0.64), C1 (between 0.71 and 0.75), and C5 (between 0.77 and 0.81), see Figure 14. The class C4 (having only 2 EDSs) was not taken into account, being unrepresentative for the analysis.

Step 1 of the second layer aims to set the threshold for CC at 0.75 to choose the candidate classes. This value is widely used in the literature for the minimum threshold value. The only class that meets this condition for the average correlation coefficient is C5. Thus, the EDSs from this class, identified with blue in Figure 13, will represent the optimal locations to integrate the SCADA system.

In the third step, the algorithm determines the multiple linear regression models to estimate the powers in the EDSs from the other classes not included in the set of the “candidate” classes. The models use as regressors the EDSs identified in the previous step.

Table 2 and Figure 15 present the APEs obtained based on the previously determined regression models for each EDS without integration into the SCADA system.

The APEs are equal to 0 for the EDSs from the candidate class (in our case, class C5), identified in bold in the table, because the system uploads the measurements from the SCADA system during the analysis period for these EDSs.

The results obtained with the proposed hybrid clustering-multiple linear regression (C-MLR) method have been compared with two other methods from the literature to demonstrate the performance of the proposed approach: the simple linear regression (SLR) method and the typical load profiles (TLPs) method. Estimating the powers from the EDSs based on the SLR method proposed in [15] contains the linear regression models between the reference, considered the MV bus of the supply HV/MV EDS, and the MV/LV EDSs from the EDN. Regarding the second method [19], each EDS has assigned a typical active and reactive profile depending on the energy consumption category obtained due to a clustering process.

The results indicate a substantial improvement in the estimation process of the APE values at the EDSs level with the proposed method compared to the SLR and TLPs methods. Table 4 and Figure 16, containing the statistical indicators associated with the quartiles Q0 – Q4, confirms this conclusion with a decrease of the maximum value from 30.9 % (calculated for SLR method) and 14.4% (calculated for TLPs method) to 6.3% and a median value (quartile Q2) from 11.8 % (calculated for SLR method) and 9.0 % (calculated for TLPs method) to 2.8 %. The SLR method led to the highest variation of the APEs in the range [4.1%, 30.9%], followed by the TLPs method with a range [3.1%, 25.4%]. The proposed method has the least variation, in the range [0.0%, 6.3%].

Regarding the MAPEs, Figure 17 presents the values of MAPEs calculated for the three methods. The value of MAPE decreases from 12.8% in the case of the SLR method to 9.9% in the case of the TLPs method and 2.6% (3 % approximately) in the case of the proposed C-MLR method (10% approximately).

In the following, Layer 3 of the methodology will use the results obtained with the proposed C-MLR method to estimate the state of the EDN associated with a day.

Because the voltages calculated at the level of the EDSs with the Newton-Raphson method based on the estimated powers are very close to the real data, being located between the allowable limits [+10%, -10%] compared with the nominal voltage, only the state variables obtained from the power flow calculations using the real and estimated active and reactive powers from the MV/LV EDSs associated with the injected active and reactive powers (P_inj and Q_inj) from the slack bus (the 20 kV bus of the HV/MV EDS), the requested active and reactive powers (P_req and Q_req) at the level of the MV/LV EDSs from the EDN, the active and reactive power losses (ΔP and ΔQ), and capacitive reactive power of the cables (Q_cap) are analysed. Table A5 and Table A6 from Appendix A present the hourly values obtained for these state variables. Table 4 contain the hourly percentage errors calculated for each state variable.

The data analysis highlights that higher errors than 1% are associated with 4 hours (3, 7, 8, and 18) for all state variables (less _Qcap). However, the maximum value does not exceed 3% in the case of the active power loss at hour 7. Exceeding the 1% threshold is sporadically present, as other variables are, but they are tiny.

Table 5 and Figure 18 present the values corresponding to statistical indicators, which offer a clear image of the accuracy of the state estimation of the EDN. Each state variable has one or a maximum of two outliers, but the values cannot be considered as high as long as they are below a threshold accepted by the DNOs.

However, all state variables have MAPEs below 1%, as seen in Figure 19, which means an outstanding performance of the estimation process.

Finally, the proposed framework can offer the DNOs the state estimate synthesis based on the aggregated values quantified as energy amounts for the analysed period. Table 6 presents the values of the state variables obtained in both cases (real and estimated) and the errors.

The information presented in the table confirm again the accuracy of the estimation process where the state variables have the errors below 0.15%.

4. Conclusions and Discussions

A state estimation process analyses an electric network's operating conditions. It involves gathering all known information about the network to determine the most likely state. This process is commonly utilized in transmission networks to improve observability and optimize operating regimes. Although the estimation state process is frequently used in transmission networks, it has not been widely adopted in the distribution sector due to the lack of proper supervision and monitoring. Despite this, it can represent an essential EDN management component for DNOs.

In this context, determining the optimal number of EDSs integrated into a SCADA system to perform on-site measurements used in the state estimation process of the EDNs represents a challenge both from a technical and economic point of view for DNOs. It involves developing a robust framework incorporating algorithms to model the EDN's conditions accurately.

Thus, a three-layer framework for optimal placement of SCADA measurements with clustering-based EDS selection to estimate the state of the EDNs has been proposed. This framework has undergone rigorous testing to ensure its reliability and accuracy. The first layer includes the determination of the classes of the EDSs with similar features from the viewpoint of requested loads based on the K-means clustering algorithm. The second layer identifies the "candidate" classes and the pilot EDSs (representing the optimal solution) with the SCADA measurements placed. The optimal placement process aims to minimise the estimation errors based on the multiple linear regression models between the EDSs from the classes not included in the set of the "candidate" classes and the pilot EDSs. The third layer allows the state estimation of the EDNs based on the load values measured in the pilot EDEs (with the SCADA system implemented) and the other EDSs obtained through the regression models. The framework was tested in the case of a real urban EDN with 39 EDSs, where accuracy was ensured for on-site measurements in 12 of the 39 EDSs (representing 31% of EDSs integrated into the SCADA system). This solution is similar to the conclusions of the study from [3], where the integration degree in the SCADA system of the MV/LV EDSs at the level of most European DNOs is approximately 30%.

The power estimation process at the level of the EDSs for the analysed period indicated a decrease of MAPE from 12.8% in the case of the SLR method [15] to 9.9% in the case of the TLPs method [19] and 2.6% (3 % approximately) in the case of the proposed C-MLR method (10% approximately). This improvement in accuracy can significantly enhance the power estimation process, leading to more reliable network operations. Reading the state variables determined by a power flow calculation, the obtained MAPEs have been below 1%, which means an outstanding performance of the estimation process. Finally, the aggregated values quantified as energy amounts for the analysed period confirmed the accuracy of the estimation process, where the errors were below 0.15%.

The framework has been tested only in the EDNs without the distributed generation sources in this research stage. However, the authors are not stopping here. They are now working on an improved version that includes power injections from distributed generation sources. The degree of uncertainty at the level of the requested/injected powers is also being considered.