Interval Load Flow for Uncertainty Consideration in Power Systems Analysis

Abstract

Modern power systems must deal with a greater degree of uncertainty in power flow calculation due to variations in load and generation introduced by new technologies. This scenario poses new challenges to power system operators which require new tools for an accurate assessment of the system state. This paper presents an interval load flow (ILF) approach for dealing with uncertainty in power system analysis. A probabilistic load flow (PLF), based on Monte Carlo Simulation (MCS), was also implemented for comparative purposes. The ILF and PLF are used to estimate the network states. Both methods were implemented in Python^® using the IEEE 34-bus, IEEE 69-bus and 192-bus Brazilian distribution system. The results with the proposed ILF on the aforementioned benchmark test systems proved to be compatible with that of the MCS, evidencing the robustness and applicability of the proposed approach.

1. Introduction

Electric power systems (EPS) play a key role in modern societies since they enable the use of technology and provide electricity to homes and industries. Nonetheless, their proper functioning can only be guaranteed if they count with the right set of tools for their operation and planning. Within these studies, the most common analysis is the power flow or load flow calculation, which can be carried out by several well-known techniques [1]. Nonetheless, the main drawback of these techniques is that their results are as accurate as their input data [2]. Inaccuracies in power flow data may be due to measurement errors, inexact forecasts, or assuming the load with some limits or unscheduled interruptions [3].

The current political and economic situation of the electricity industry is increasingly encouraging the use of distributed energy resources (DER) which include small-scale generation units that are located on the consumer’s side of the meter, storage devices and demand respond. Therefore, the role of customers is progressively changing and adapting, becoming more active. Within this context, the load models adopted for power flow calculation must also adapt to characterize new uncertainties [4]. Interval analysis can be seen as an effective tool to take into account this type of uncertainties within power flow studies. In this case, the load is modeled through a set of values ranging from the smallest to the largest possible one; and therefore, nodal voltages, losses and power flows are calculated in an interval way [5,6,7,8,9]. This approach is known as Interval Load Flow (ILF).

To calculate an ILF it is necessary to use mathematical tools that guarantee its convergence. This may be done through the Newton Intervalar method [10] as well as with the Krawczyk method [11], which finds the solution of the interval system within a given tolerance.

In [8], an optimal distribution planning model is proposed with the addition of uncertainties via ILF, where the interval results are used as merit functions for an optimization algorithm. In addition to the ILF, interval mathematics are also used to add uncertainty to other problems, such as state estimation and short circuit calculations [12].

Another approach for adding uncertainties is through the Probabilistic Load Flow (PLF) [2,3], where loads are modeled through their probability distribution functions. Several works have proposed mathematical methods to solve this problem. The Monte Carlo Simulation (MCS) approach is the most widely used method for solving the PLF. This method consists of the solution of several deterministic load flows with random variations in their parameters, based on their probability density functions; thus, if the sample of results is large enough, a highly accurate result can be obtained [13]. Probabilistic uncertainties can also be evaluated with other methods such as Polynomial Chaos Expansion (PCE), which is able to work in a non-intrusive way and exhibits a high performance in stochastic modeling of random variables and processes [14].

ILF has proven to be an effective tool when dealing with uncertainties. In [15], the authors implemented an ILF analysis considering interval output of wind farms with the aim of providing accurate ranges of the power grid variables which can be used by the system operator to guarantee security of the power system. In [16], an affine arithmetic-based power flow algorithm is proposed to take into account regional control of unscheduled power fluctuations. The authors in [17] propose an ILF based on multi-stage affine arithmetic which is applied to an unbalanced distribution network to address the influence of distributed generation and load uncertainty. In [18], a hybrid probabilistic and interval load low is proposed by means of a clustering-based analytical method. In this case, the uncertainties of random and interval variables are simultaneously considered. The authors in [19] propose an ILF by means of the optimizing-scenarios method. In this approach, the interval uncertainties are considered as variables with varying bounds, and the objective function is set to calculate these unknown variables. In [20] the authors proposed an ILF based on the Taylor inclusion function, while in [21] an ILF approach with correlated uncertain power injections is presented. Other variants of the ILF are also proposed in [22,23].

The main objective of this paper is to employ the ILF and analyze its performance in calculating the states of the electrical network for benchmark distribution systems comparing the result with that of a conventional MCS where the same networks are analyzed with their load modeled by a function of uniform probability density over the same range of the interval load. Therefore, the main contributions of this work can be described as follows: (i) an exemplified and straightforward calculation methodology of the ILF is presented, (ii) a comparison of the proposed ILF with MCS is performed and, (iii) a sensibility analysis is carried out regarding the input variables of the ILF and its convergence.

The structure of this paper is as follows: In Section 2, basic concepts of Interval Arithmetic are presented. Section 3 presents, in detail, the methodology for the ILF. The main results of the ILF and its comparison with MCS are presented in Section 5. Finally, the conclusions of the work are presented in Section 6.

2. Interval Arithmetic

Interval arithmetic is based on operations using real intervals. So, a given variable X can be defined as $X=[x_1,x_2]=\{\overline{x}\in \mathit{R},x_1\le \overline{x}\le x_2\}$. Within this interval, the midpoint, diameter and radius properties can be defined as follows [24]:

$$\mathrm{Mid}\left(X\right)=\frac{1}{2}\left(x_1+x_2\right)$$

(1)

$$\mathrm{Diameter}\left(X\right)=x_2-x_1$$

(2)

$$\mathrm{Radius}\left(X\right)=\frac{\mathrm{Diameter}\left(X\right)}{2}$$

(3)

Given the interval variables $X=[x_1,x_2]=\{\overline{x}\in \mathbf{R},x_1\le \overline{x}\le x_2\}$ and $Y=[y_1,y_2]=\{\overline{y}\in \mathbf{R},y_1\le \overline{y}\le y_2\}$, arithmetic operations are defined for intervals as follows [24,25]:

$$X+Y=\left[x_1+y_1,x_2+y_2\right]$$

(4)

$$X-Y=\left[x_1-y_1,x_2-y_2\right]$$

(5)

$$X·Y=\left[min(x_1{y}_{1,}{x}_1{y}_2,x_2{y}_1,x_2{y}_2\right),max\left(x_1{y}_{1,}{x}_1{y}_2,x_2{y}_1,x_2{y}_2\right)]$$

(6)

$$\frac{X}{Y}=\left[x_1,x_2\right]\left[\frac{1}{y_1},\frac{1}{y_2}\right]\mathrm{if}0\notin \left[y_1,y_2\right]$$

(7)

The methodology for the evaluation of the ILF requires the solution of an analytical load flow, which is named in this paper as a deterministic load flow. From the input values corresponding to deterministic net injected powers and bus voltages (magnitudes and angles) obtained as solution of the deterministic load flow, the interval values of the net injected powers as well as the magnitudes and angles of the bus voltages are established.

Thus, the convergences of the intervals are analyzed according to the variations of the pre-established diameters. The intervals of the bus voltages are reduced by the Krawczyk method, which in this case consists on the application of the Krawczyk operator [11] at each iteration of the algorithm, as indicated in (8). The Krawczyk operator returns interval values that must be intersected with the current values to reduce overestimation of the interval variables.

$$\mathrm{K}\left(x^h,X^h\right)=x^h-\mathrm{C}\mathrm{f}\left(x^h\right)+\left(\mathrm{I}-\mathrm{C}\mathrm{J}\left(X^h\right)\right)\left(X^h-x^h\right)$$

(8)

In this case, X is a real interval, x is the midpoint of X, f$\left(x\right)$ is the value of the nonlinear function at point x, J$\left(X\right)$ is the Jacobian of the nonlinear function calculated in the range X, C is the matrix of preconditioning which is equal to the Jacobian inverse calculated at the midpoint of X, I is the identity matrix and h is the current iteration. X is the vector of the modules and angles of the interval bus voltages. The main advantage of using the Krawczyk operator is that there is no need to solve a linear system.

3. Interval Load Flow

The most conventional method for solving the load flow problem is the formulation of the power injected at each bus in its polar form [7]. Interval analysis is useful for modeling uncertainties in numerical analyzes [11]; thus, it is used in conjunction with the polar formulation of the Newton-Raphson power flow to create what is called an ILF.

3.1. Algorithm Description

The main steps for obtaining a solution to the power flow considering interval variables are described below.

Step 1: Obtain the bus voltages using a deterministic load flow. In this work, the deterministic load flow corresponds to the solution of the power flow through the Newton-Raphson method with input variables fixed in their deterministic values, i.e., power injections (load and generation) fixed in its predicted maximum values.

Step 2: Define the percentage of variation of demands and calculate their respective interval values according to (9) and (10):

$${Pc}_k^i={Pc}_k^d\left[1-\alpha ,1+\alpha \right]$$

(9)

$${Qc}_k^i={Qc}_k^d\left[1-\alpha ,1+\alpha \right]$$

(10)

where ${Pc}_k^d$ and ${Qc}_k^d$ are the values of the deterministic active and reactive loads, respectively, at bus k, $\alpha$ is the demand percentage variation and ${Pc}_k^i$ and ${Qc}_k^i$ are the interval values of the active and reactive loads, respectively. It is worth mentioning that in this work, the interval representation is only considered for loads; nonetheless, the logic of Equations (9) and (10) can be extended for an interval representation of the generation, if necessary.

Step 3: Calculate the expected interval power through the expected deterministic power and the interval load calculated in (9) and (10), as given by (11) and (12).

$$P_{{exp}_k}^i=P_{{exp}_k}^d-{Pc}_k^d+{Pc}_k^i$$

(11)

$$Q_{{exp}_k}^i=Q_{{exp}_k}^d-{Qc}_k^d+{Qc}_k^i$$

(12)

where $P_{{exp}_k}^d$ and $Q_{{exp}_k}^d$ are, respectively, the active and reactive expected power at bus k used in the deterministic load flow, while $P_{{exp}_k}^i$ and $Q_{{exp}_k}^i$ are, respectively, the interval active and reactive expected power at bus k.

Then, the interval mismatches for the initialization of the bus voltages magnitude and angles can be calculated. The interval mismatches must be obtained with the calculated deterministic powers and the expected interval powers, using expressions (13) and (14) as indicated in [26].

$$\mathsf{\Delta }{P}_k^i=P_{{exp}_k}^i-P_{{calc}_k}^d$$

(13)

$$\mathsf{\Delta }{Q}_k^i=Q_{{exp}_k}^i-Q_{{calc}_k}^d$$

(14)

where $P_{{calc}_k}^d$ and $Q_{{calc}_k}^d$ are the calculated deterministic active and reactive power at bus k, while $\mathsf{\Delta }{P}_k^i$ and $\mathsf{\Delta }{Q}_k^i$ are the interval values of the mismatches at bus k, respectively.

Thus, the interval increments of the bus voltages magnitudes and their respective angles are calculated from the Jacobian matrix of the last iteration of the deterministic power flow and the interval mismatches ($\mathsf{\Delta }{P}^i$ and $\mathsf{\Delta }{Q}^i$) as given by (15).

$$\left[\begin{array}{c}\mathsf{\Delta }{\theta }^i\\ \mathsf{\Delta }{V}^i\end{array}\right]={\left[\mathrm{Jac}\right]}^{-1}\left[\begin{array}{c}\mathsf{\Delta }{P}^i\\ \mathsf{\Delta }{Q}^i\end{array}\right]$$

(15)

where $\mathsf{\Delta }{P}^i$ and $\mathsf{\Delta }{Q}^i$ are the vectors of interval mismatches of active and reactive power, respectively. Also, $\left[\mathrm{Jac}\right]$ is the jacobian matrix of the last iteration of the deterministic load flow and $\mathsf{\Delta }{\theta }^i$ and $\mathsf{\Delta }{V}^i$ are the interval increments of the angles and magnitudes of bus voltages, respectively. Then, the interval voltage magnitudes and phase angles can be calculated using the deterministic solution and the interval increments indicated in (16).

$$\left[\begin{array}{c}{\theta }^i\\ V^i\end{array}\right]=\left[\begin{array}{c}{\theta }^d\\ V^d\end{array}\right]+\left[\begin{array}{c}\mathsf{\Delta }{\theta }^i\\ \mathsf{\Delta }{V}^i\end{array}\right]$$

(16)

where ${\theta }^d$ and $V^d$ are respectively the vectors of the deterministic solution of voltage phase angles and magnitudes, while ${\theta }^i$ and $V^i$ are the vectors of interval values of voltage phase angles and magnitudes, respectively.

Step 4: Obtain the interval active and reactive power calculated using the admittance matrix of the system as well as the current interval magnitudes and angles of the buses voltages, as described in (17) and (18) and, thus recalculate the power mismatches, as proposed by [26]. Note that these equations are similar to those in the deterministic load flow, only now the variables of the problem are replaced by their respective interval variables.

$$P_{{calc}_k}^i={∥\dot{V_k^i}∥}^2{G}_{kk}+∥\dot{V_k^i}∥\sum _{mϵ\mathsf{\Omega }k}^k∥\dot{V_m^i}∥\left(G_{km}cos\left({\theta }_k^i-{\theta }_m^i\right)+B_{km}sin\left({\theta }_k^i-{\theta }_m^i\right)\right)$$

(17)

$$Q_{{calc}_k}^i=-{∥\dot{V_k^i}∥}^2{B}_{kk}+∥\dot{V_k^i}∥\sum _{mϵ\mathsf{\Omega }k}^k∥\dot{V_m^i}∥\left(B_{km}cos\left({\theta }_k^i-{\theta }_m^i\right)+G_{km}sin\left({\theta }_k^i-{\theta }_m^i\right)\right)$$

(18)

where $∥\dot{V_k^i}∥$ and $∥\dot{V_m^i}∥$ are the interval voltage magnitudes at buses k and m, respectively. Also, ${\theta }_k^i$ and ${\theta }_m^i$ are the interval voltage angles at buses k and m, respectively. $G_{ij}$ and $B_{ij}$ are respectively the real and complex components of $Y_{ij}$ of the system admittance matrix and $P_{{calc}_k}^i$ and $Q_{{calc}_k}^i$ are, respectively, the interval active and reactive calculated powers. Thus, the interval mismatches of active and reactive power are updated according to (19) and (20) [26]:

$$\mathsf{\Delta }{P}_k^i=P_{{exp}_k}^i-P_{{calc}_k}^i$$

(19)

$$\mathsf{\Delta }{Q}_k^i=Q_{{exp}_k}^i-Q_{{calc}_k}^i$$

(20)

Step 5: Apply the Krawczyk operator according to (8), where the variables will assume the following values:

$$X^h=\left[\begin{array}{c}{\theta }^i\\ V^i\end{array}\right]$$

(21)

$$x^h=\mathrm{Mid}\left(X^h\right)$$

(22)

$$\mathrm{f}\left(x^h\right)=\mathrm{Mid}\left(\left[\begin{array}{c}\mathsf{\Delta }{P}^i\\ \mathsf{\Delta }{Q}^i\end{array}\right]\right)$$

(23)

Step 6: Update the interval increments for the voltage magnitudes and angles through the intersections between the interval voltages in the previous iteration and the voltage obtained by the Krawczyk operator as indicated in (24).

$${\left[\begin{array}{c}{\theta }^i\\ V^i\end{array}\right]}^{h+1}=\left[\begin{array}{c}{\theta }^i\\ V^i\end{array}\right]\bigcap K\left(x^h,X^h\right)$$

(24)

Step 7: Check if the greater radius variation of the interval magnitude and angle of the bus voltages between the iterations is less than the specified tolerance, if yes, the interval power flow has converged; otherwise, return to step 5.

3.2. Illustrative Example

A detailed explanation of the ILF algorithm using a didactic 3-bus system, shown in Figure 1, is presented. In this figure, the line impedances are in ohms, the loads are in MVA and the bus voltages are in per-unit (p.u). This system has a base power of 1 MVA and a base voltage of 23 kV. In Figure 1, the state variables obtained trough a deterministic load flow in each bus are shown. These calculations correspond to those described in Step 1 of Section 3.1. The Jacobian matrix of the last iteration of the deterministic load flow for this system corresponds to:

$$\left[\begin{array}{c}\mathrm{Jac}\end{array}\right]=\left[\begin{array}{cccccc}{10}^{20}& -1175.96& 0.0& 354.39& -349.20& 0.0\\ -1174.55& 2045.74& -871.19& -353.49& 369.47& -19.22\\ 0.0& -871.14& 871.14& 0.0& -21.15& 19.20\\ -348.76& 348.76& 0.0& {10}^{20}& -1177.47& 0.0\\ 353.49& -372.68& 19.19& -1174.55& 2047.45& -872.67\\ 0.0& 21.12& -21.12& 0.0& -872.26& 871.95\end{array}\right]$$

The Jacobian matrix was considered in its full form. Therefore, high values in the order of ${10}^{20}$ were assumed in the main diagonal corresponding to the slack bus, because in this bus the magnitude and phase angle voltage are fixed in 1.0 p.u and 0 degrees, respectively, and thus iterative updating of these variables is not necessary.

A percentage variation of 5%, i.e., $\alpha$ = 5%, is assumed for the loads of the system to define the interval values, as presented in Step 2, for $P_c^i$ and $Q_c^i$ as:

$$\left[\begin{array}{c}{P}_c^i\end{array}\right]=\left[\begin{array}{c}[0,0]\\ [1.747,1.932]\\ [0.931,1.029]\end{array}\right]$$

$$\left[\begin{array}{c}{Q}_c^i\end{array}\right]=\left[\begin{array}{c}[0,0]\\ [0.437,0.483]\\ [0.323,0.375]\end{array}\right]$$

Next, the interval increments of the bus voltages magnitudes and angles must be calculated, as stated in the Step 3 of the ILF algorithm. For this, the interval mismatches of active and reactive powers at each bus are obtained as:

$$\left[\begin{array}{c}\mathsf{\Delta }{P}^i\\ \mathsf{\Delta }{Q}^i\end{array}\right]=\left[\begin{array}{c}[0,0]\\ [\mathrm{-0.092},0.092]\\ [\mathrm{-0.049},0.049]\\ [0,0]\\ [\mathrm{-0.023},0.023]\\ [\mathrm{-0.017},0.017]\end{array}\right]$$

Using the interval mismatches of the power loads and the Jacobian Matrix, the interval increments of the bus voltages magnitudes and their respective angles are computed as:

$$\left[\begin{array}{c}\mathsf{\Delta }{\theta }^i\\ \mathsf{\Delta }{V}^i\end{array}\right]=\left[\begin{array}{c}[0,0]\\ [\mathrm{-0.0001},0.0001]\\ [\mathrm{-0.0002},0.0002]\\ [0,0]\\ [\mathrm{-6.454},6.454]\\ [\mathrm{-8.541},8.541]\end{array}\right]$$

Using the interval increments of the bus voltages, the interval voltage magnitudes and phase angles are obtained, using Equation (16), as:

$$\left[\begin{array}{c}{\theta }^i\\ V^i\end{array}\right]=\left[\begin{array}{c}[0,0]\\ [\mathrm{-0.0021},\mathrm{-0.0019}]\\ [\mathrm{-0.0033},\mathrm{-0.0029}]\\ [1,1]\\ [0.9986,0.9987]\\ [0.9982,0.9984]\end{array}\right]$$

The interval mismatches of the active and reactive power are updated as presented in Step 4 of the ILF algorithm. For this illustrative example, the updated mismatches correspond to:

$$\left[\begin{array}{c}\mathsf{\Delta }{P}^i\\ \mathsf{\Delta }{Q}^i\end{array}\right]=\left[\begin{array}{c}[0,0]\\ [\mathrm{-0.02116639},0.02116643]\\ [\mathrm{-0.0006336},0.0006334]\\ [0,0]\\ [\mathrm{-0.0065334},0.00065332]\\ [\mathrm{-0.00001280},0.00001275]\end{array}\right]$$

Then, the Krawczyk operator must be applied as shown in Step 5 of the algorithm. The values for $x^h$ and f($x^h$) are calculated using Equations (22)–(24), obtaining:

$$x^h=\left[\begin{array}{c}0\\ -0.00202\\ -0.00313\\ 1\\ 0.9987\\ 0.9983\end{array}\right]$$

$$\mathrm{f}\left(x^h\right)=\left[\begin{array}{c}0\\ 1.91684\\ -6.23102\\ 0\\ -4.98810\\ -2.20536\end{array}\right]\times {10}^{-8}$$

Following the calculations presented in Step 6, the updated values for the voltage magnitudes and angles are obtained as:

$${\left[\begin{array}{c}{\theta }^i\\ V^i\end{array}\right]}^{h+1}=\left[\begin{array}{c}[0,0]\\ [\mathrm{-0.00216},\mathrm{-0.00215}]\\ [\mathrm{-0.0031314},\mathrm{-0.003109}]\\ [1,1]\\ [0.9987137,0.9987141]\\ [0.9982975,0.9982977]\end{array}\right]$$

The convergence of the ILF algorithm is check as stated in Step 7. As the obtained maximum radius variation is 8.52 × ${10}^{-5}$ in the $V_k^i$ of the bus 3, the convergence of the ILF is reached.

4. Probabilistic Load Flow

The PLF was proposed in 1974 [3] and can be solved analytically, as in [2,3,27], numerically as presented in [28] or through a hybrid approach as presented in [2].

The numerical simulation can be performed using the MCS method. In this case, the probability density functions of each of the parameters that are considered as random variables are sampled to proceed with the simulation. The simulation and sampling must be carried out many times so that, with a large sample of input variables, a good result is obtained.

In the MCS performed, uniform probability density functions were adopted for the active and reactive loads varying between $P^{+}$ and $P^{-}$ for active load, and between $Q^{+}$ and $Q^{-}$ for reactive power load as illustrated in Figure 2 and Figure 3. In this case, $P^{+}$, $P^{-}$, $Q^{+}$ and $Q^{-}$ are given by equations (25) to (28).

$$P_k^{+}={Pc}_k^d(1+\alpha )$$

(25)

$$P_k^{-}={Pc}_k^d(1-\alpha )$$

(26)

$$Q_k^{+}={Qc}_k^d(1+\alpha )$$

(27)

$$Q_k^{-}={Qc}_k^d(1-\alpha )$$

(28)

Thus, a large number of samples are carried out in order to ensure that most of the probable values are obtained in the simulation.

5. Test and Results

A program capable of calculating the deterministic power flow using the Newton-Rapshon method [26] with sparse techniques [29] was developed in Python^® and then the IPF was processed as indicated in Section 3.

Medium voltage distribution circuits with radial topologies were simulated with a variation of 5% in the load with a stopping criteria of a variation in the radius of ${10}^{-4}$ between iterations.

Although probabilistic and interval variables are different concepts (interval representation is typically non-probabilistic) [14], the results of the IPF were compared with those of a MCS. For this, the probabilistic load flow was run with a high number of draws to obtain the maximum and minimum values of the simulated quantities, so that results could be compared with the intervals obtained by the interval power flow. Also, a uniform probability distribution was used in the MCS to make the probabilistic approximation closer to the interval, since all the values contained within it are equally representative for the result. In this case, 3000 drawings were used to verify the adherence of the method in the 34-bus and 69-bus test systems and with 30,000 drawings for the 192-bus system. Uniform probability density functions for the net injected power at buses were considered in the MCS, their loads varying by the same percentage as in the ILF.

5.1. IEEE 34-Bus Radial TEST System

The data of this system can be found in [30]. The base values considered are 1000 kVA and 11 kV. The 34-bus system has a total active power demand of 4636.5 kW and total reactive power of 2885.5 kvar. The minimum and maximum values of voltage magnitude at the buses, obtained by a deterministic load flow are 1.0 p.u and 0.942 p.u, respectively.

Table 1. Phase Angle—IEEE 34-bus test system.

Bus	Method	Lower Angle (°)	Upper Angle (°)	Deterministic Angle (°)
5	MCS	0.23615	0.35810	0.31588
	ILF	0.26275	0.36900
9	MCS	0.48285	0.67387	0.61611
	ILF	0.53320	0.69891
23	MCS	0.70353	1.00872	0.90243
	ILF	0.78490	1.01959
27	MCS	0.79289	1.13820	1.00895
	ILF	0.88240	1.13499
33	MCS	0.52957	0.72857	0.66792
	ILF	0.58067	0.75503

and

Table 2. Voltage Magnitude—IEEE 34-bus test system.

Bus	Method	Lower Magnitude (p.u.)	Upper Magnitude (p.u.)	Deterministic Magnitude (p.u.)
5	MCS	0.97556	0.97868	0.97603
	ILF	0.97479	0.97729
9	MCS	0.96113	0.96639	0.96197
	ILF	0.95999	0.96396
23	MCS	0.94417	0.95144	0.94601
	ILF	0.94318	0.94885
27	MCS	0.93963	0.94997	0.94166
	ILF	0.93861	0.94474
33	MCS	0.95896	0.96466	0.95992
	ILF	0.95784	0.96201

show the results of the phase angle and magnitude voltages, respectively, calculated by the ILF and MCS for a variation in the radius of ${10}^{-4}$. Note that the results are shown only for some representative buses of the system. In these tables, the lower and upper values refer to the lower and upper limits of the intervals and, the deterministic values correspond to those obtained trough the execution of the deterministic power flow.

In this test system, the largest variation of the voltage magnitude around the deterministic result with the MCS is 0.88% at bus 27 and the smallest variation is 0.07% at bus 2. For the ILF the largest variation is 0.33% at bus 27 and the smallest one is 0.03% at bus 2. Regarding phase angles the largest variation around the deterministic result for the MCS and ILF are 34.75% and 25.11%, respectively both at bus 2; while the smallest variation for MCS and ILF are 21.32% and 12.54%, respectively at bus 33.

For this system a sensibility analysis was done by changing the percent variation of the loads from 5% to 20% with steps of 5%. Figure 4 depicts the radius of the voltage magnitude and Figure 5 illustrates the radius of the phase angle for each simulated scenario.

As expected, with the increase of the percentage variation on the loads, the radius of the result variables became larger.

5.2. IEEE 69-Bus Distribution Network

The data of this system can be found in [31]. The base values considered are 100 kVA and 12.66 kV. The 69-bus system has a total active power demand of 3802.19 kW and total reactive power of 2694.60 kvar. The minimum and maximum values of voltage magnitude at the buses, obtained by a deterministic load flow are 1.0 p.u and 0.909 p.u, respectively.

Table 3. Phase Angle—IEEE 69-bus test system.

Bus	Method	Lower Angle (°)	Upper Angle (°)	Deterministic Angle (°)
2	MCS	−0.00140	−0.00089	−0.00122
	ILF	−0.00133	−0.00223
29	MCS	−0.00587	−0.00421	−0.00530
	ILF	−0.00578	−0.00482
48	MCS	−0.06559	−0.03440	−0.05247
	ILF	−0.05724	−0.04770
54	MCS	0.00461	0.37110	0.19470
	ILF	0.12409	0.26530
63	MCS	0.33813	187.126	112.526
	ILF	0.90650	134.402

and

Table 4. Voltage Magnitude—IEEE 69-bus test system.

Bus	Method	Lower Magnitude (p.u.)	Upper Magnitude (p.u.)	Deterministic Magnitude (p.u.)
2	MCS	0.99996	0.99997	0.99996
	ILF	0.99996	0.99997
29	MCS	0.99984	0.99987	0.99985
	ILF	0.99984	0.99986
48	MCS	0.99840	0.99883	0.99854
	ILF	0.99847	0.99980
54	MCS	0.96770	0.97738	0.97141
	ILF	0.96985	0.97296
63	MCS	0.89464	0.93411	0.91165
	ILF	0.90679	0.91652

show the results of the state variables calculated by ILF and MCS for a radius of ${10}^{-4}$.

For the IEEE 69-bus test system, the largest variation at the voltage magnitude around the deterministic result with the MCS is 2.61% at bus 65 and the smallest one is 0.002% at bus 29. For the ILF the largest variation is 0.55% at bus 65 and the smallest one is 0.001% at bus 29. Regarding phase angles the largest variation around the deterministic result for MCS and ILF are 307.45% and 116.08%, respectively, both at bus 32; while the smallest variation for MCS and ILF are 21.32% at bus 4 and 21.5% at bus 37, respectively.

5.3. 192-Bus Brazilian Distribution Network

The data of this system can be found in [32]. The base values are 100 kVA and 13.8 kV. The 192-bus system has a total active power demand of 6031 kW and reactive power of 2124 kvar. The minimum and maximum values of voltage magnitude at the buses, obtained by a deterministic load flow are 1.0 p.u and 0.955 p.u, respectively.

Table 5. Phase Angle—192-bus Brazilian distribution network.

Bus	Method	Lower Angle (°)	Upper Angle (°)	Deterministic Angle (°)
8	MCS	−0.24610	−0.21364	−0.24170
	ILF	−0.27692	−0.20648
36	MCS	−0.90111	−0.76996	−0.87943
	ILF	−100.956	−0.75299
64	MCS	−0.44769	−0.38699	−0.43892
	ILF	−0.50307	−0.37477
113	MCS	−0.86144	−0.73693	−0.84702
	ILF	−0.97206	−0.72197
187	MCS	−0.91812	−0.78600	−0.89612
	ILF	−102.858	−0.76365

and

Table 6. Voltage Magnitude—192-bus Brazilian distribution network.

Bus	Method	Lower Magnitude (p.u.)	Upper Magnitude (p.u.)	Deterministic Magnitude (p.u.)
8	MCS	0.98727	0.98888	0.98795
	ILF	0.98732	0.99176
36	MCS	0.95602	0.96027	0.95686
	ILF	0.95461	0.95911
64	MCS	0.97824	0.97994	0.97827
	ILF	0.97713	0.97939
113	MCS	0.95754	0.96171	0.95849
	ILF	0.95632	0.96066
187	MCS	0.95603	0.95959	0.95612
	ILF	0.95612	0.95842

show the results of the state variables calculated by ILF and MCS for a variation in the radius of ${10}^{-4}$.

For this system, the largest variation at the voltage magnitude around the deterministic result with the MCS is 0.37% at bus 51 and the smallest one is 0.01% at bus 2. For the ILF the largest variation is 0.25% at bus 51 and the smallest variation is 0.01% at bus 2. Regarding phase angles, the largest variation around the deterministic result for MCS and ILF are respectively 13.33% at bus 94 and 14.81% at bus 54, while the smallest variation for MCS and ILF are respectively 12.11% at bus 2 and 14.53% at bus 55.

As shown in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6, it can be seen that there is an intersection between the ranges of results obtained through the MCS and the ILF for voltage magnitudes and phase angles in all buses. It is also evident that the MCS and ILF intervals contain the deterministic power flow solution for all cases.

6. Conclusions

Nowadays, power and distribution systems must deal with greater degrees of uncertainty not only related to new generation technologies but also to variations in the demand side. In this context, system operators are in constant search of new tools for power system analyses. This paper contributes in the analysis and discussion of such tools by proposing an ILF approach that takes into account uncertainties in demand. A step-by-step straightforward methodology for calculating an ILF based on the polar form of the Newton-Rapshon method is presented. Furthermore, a comparison of the proposed approach with a MCS is carried out along with a sensibility analysis regarding the input variables of the ILF and its convergence for several benchmark test systems. It was found that with small variations in loads, the ILF presents results as good as those produced by MCS, but with a reduced computational effort. As the operation and planning of electrical energy systems are based on the supply of the load, which in turn is endowed with uncertainties, it can be concluded that the use of the ILF as a methodology for the evaluation of the performance of energy systems is satisfactory and recommended in power flow analyses.