1 Introduction

Optimal power flow (OPF) is one of the fundamental static power system calculations. It has wide application in electrical engineering, including scheduling of generators, loss reduction, congestion management, and expansion planning. As distributed generation (DG) and controllable loads (e.g., electric vehicles) proliferate, active network management has been introduced in distribution systems [1]. In this context, OPF is no longer limited to the domain of high voltage transmission networks and has been gradually investigated for application to distribution networks [2].

In general, all the input data of OPF are deterministic. Governed by nonlinear Kirchhoff’s laws, such deterministic optimization problems can be solved by many methods, such as successive linear/quadratic programming [3, 4], trust-region-based methods [5], the Lagrangian Newton method [6] and the interior-point method [7]. However, with increasing internal and external factors of uncertainty, such as the power demand affected by daily economic activities, power generated by renewable energy, and grid parameters obtained by approximate measurements, the input data have increasing uncertainty, which challenges conventional deterministic OPF models. The degree of uncertainty for some factors can be reduced, but for most uncontrollable factors, it is very difficult to decrease the impact of their uncertainty. Hence, OPF should be able to manage uncertainties in power flow performance.

Most conventional methodologies to address uncertainty are based on probabilistic methods that account for the variability and stochastic nature of the input data. Current OPF research on this topic can be divided into two categories, probabilistic OPF (P-OPF) and stochastic OPF (S-OPF). P-OPF [8] is a well-respected approach for characterizing the output of an implicit function whose inputs are random variables, where the cumulant method [9] and the point estimate method [10] are examples of very efficient P-OPF computation. However, the solution of P-OPF is influenced indirectly by the randomness of input variables, and only the probability distributions of control variables can be determined.

In S-OPF problems, the objective function and constraints are usually described by probability equations or inequalities, which means that the randomness of input variables can directly impact the solution [11, 12]. Thus, constraint satisfaction in an uncertain environment can be achieved. Another competitive choice for modeling S-OPF problems often considered is chance-constrained programming (CCP) [13] in which the constraints (some or all) can be violated with a pre-assigned level of probability. These are referred to as “chance constraints”. However, it is generally agreed that closed-form solutions for S-OPF are rarely available due to the different types of parameter distributions and the nonlinear nature of OPF, and thus one has to resort to approximation methods to solve instances of nontrivial size. Researchers have devised solution methods that rely mainly on discretization of the uncertain parameters. Theoretically, these scenario-based approaches to S-OPF can achieve any desired level of accuracy, but the required computational resources for scenario generation and coordinated optimization could be prohibitively expensive.

Although scenario-free approaches to chance-constrained OPF have also been developed, specific uncertainty distributions (e.g., Gaussian distribution) should usually be assumed for analytical reformulation of the chance constraints. Furthermore, to deal with the higher complexity of chance-constrained OPF, the existing approaches either assume a DC-OPF [14], a linearized OPF [15], a convex relaxation-based OPF [16], or they solve iteratively linearized instances of nonlinear OPF [13]. In addition to these approximation or relaxation approaches, accurate probability information in the form of scenarios or pre-defined distributions is always necessary, which is difficult to obtain in practice.

Robust optimization (RO) is another promising approach to model OPF involving uncertainties. It should be noticed that RO is applicable for problems with convex feasible regions. Thus, it is mainly applied to OPF with DC power flow constraints [17] or linearized power flow constraints [18]. To improve its accuracy, second-order cone programming (SOCP) is introduced by relaxing the power flow equality constraints [19] and thus the robust optimal power flow model is converted into a mixed-integer SOCP model through the robust counterpart. However, these approximation-based or relaxation-based methods may still cause a gap to the original exact power flow equations.

Interval arithmetic (IA) offers another approach to model and analyze uncertainties. It can be used to give the variation ranges of output variables and simultaneously verify the satisfaction of constraints. One of its successful applications in electrical engineering is to determine the strict solution bounds of power flow [20]. Ranges of uncertainty expressed by intervals are easily available in practice, so programming methods incorporating IA have received widespread attention. Early attempts in this area were focused on interval linear programming (ILP), [21, 22]. Since the ILP model can be transformed into several deterministic linear programming models according to the signs of the interval coefficients and variables, it can be solved in an efficient way. However, due to the nonlinear and non-convex characteristics of OPF, the best and worst optimums are not necessarily on the bounds of the intervals, resulting in an interval nonlinear programming (INLP) problem that is much more difficult to solve. To the best of our knowledge, the general idea for solving INLP problems in existing literature is to use approximate approaches via order relations of interval numbers or Taylor expansions [23,24,25].

To conquer the drawbacks including the “dependency problem” and the “wrapping effect” in IA, affine arithmetic (AA) is proposed in [26]. By using AA, two self-validated computation approaches for the power flow calculation and OPF with uncertainties are proposed in [27] and [28]. To satisfy constraints fully and to achieve accurate operational bounds of the solution under uncertainties, the authors in [29] further combined the linear approximation method and AA-based power flow calculation to address the uncertain reactive power optimization problem in transmission grid. It is worth noting that OPF incorporating intervals is a typical INLP problem. Although the above solution algorithms have been applied successfully in some applications, an analytical OPF model incorporating intervals like the deterministic OPF model is still unavailable.

Motivated by the requirements of OPF in distribution networks, namely, full satisfaction of constraints and determining the strict operational bounds of the solution under uncertainties, an interval OPF model incorporating robustness is proposed in this paper; the main accomplishments are as follows.

  1. 1)

    An interval OPF (I-OPF) model applied to distribution networks is derived to deal with uncertainties in power flow and constraint satisfaction, by employing the affine arithmetic and interval Taylor expansion.

  2. 2)

    To improve solution accuracy, combining second-order cone relaxation and successive linear approximation, an enhanced method accounting for the equivalent error caused by high-order items is also derived.

  3. 3)

    The proposed I-OPF model and its enhanced method are tested on a modified IEEE 33-bus distribution system and a real 113-bus system, using Monte Carlo simulation and stochastic OPF as comparative analyses.

The remainder of this paper is arranged as follows. Section 2 presents a brief background of the deterministic OPF model using the branch flows, as well as the formulation of OPF incorporating interval uncertainties. The I-OPF model is derived in Section 3. Section 4 derives in detail the enhanced method with improved accuracy. Case studies based on a modified IEEE 33-bus system and a real 113-bus system are presented in Section 5, followed by the conclusion in Section 6.

2 Problem formulation

Consider a distribution network with a set of controllable devices denoted by \({\mathcal{G}}\), including DG units, energy storage units and var compensation devices, and loads denoted by \({\mathcal{L}}\). Since the distribution network typically has a radial structure, we model it as a connected tree graph \({(\mathcal{N}},{\mathcal{E})}\), where each node \(i \in {\mathcal{N}}\) represents a bus and each link \(\text{(}i,j\text{)} \in {\mathcal{E}}\) represents a branch. The root of the tree is the feeder with a fixed voltage and flexible power injection.

2.1 Formulation of the deterministic OPF

For the link \(\text{(}i,j\text{)} \in {\mathcal{E}}\), let \(I_{ij} (t)\) and \(S_{ij} \triangleq P_{ij} + {\text{i}}Q_{ij}\) be the complex current and complex power respectively, flowing from buses \(i\) to \(j\), where the corresponding complex impedance is denoted as \(z_{ij} \triangleq r_{ij} + {\text{i}}x_{ij}\). For the bus \(i \in {\mathcal{N}}\), let \(V_{i}\) be its complex voltage. \(s_{i} \triangleq p_{i} + {\text{i}}q_{i}\) is the corresponding net load which is the load minus the generation at bus i, given by:

$$s_{i} = s_{{{\mathcal{L}}_{i} }} - s_{{g_{\text{u}} ,i}} - \sum\limits_{{g_{\text{c}} \in {\mathcal{G}}_{i} }} {s_{{g_{\text{c}} }} }$$
(1)

where \(s_{{{\mathcal{L}}_{i} }}\) is the load at bus i; \(s_{{g_{\text{u}} ,i}}\) is the power injection of the uncontrollable device; \({\mathcal{G}}_{i}\) is the subset of controllable devices connected to bus i; \(s_{{g_{\text{c}} }}\) corresponds to the controllable device \(g_{\text{c}} \in {\mathcal{G}}_{i}\). The power injection \(s_{{g_{\text{c}} }}\) is generally limited to a pre-specified set \(S_{g,i}\), representing for example the capacity and power factor constraints. For brevity, we just denote it as:

$$s_{{g_{\text{c}} }} \in S_{g,i} \,\,\,\,\,\,\,\,\,\, i \in {\mathcal{N}}$$
(2)

The steady-state power flows can be described by the Distflow model for radial networks [30, 31]: \(\forall \text{(}i,j\text{)} \in {\mathcal{E}}\)

$$s_{j} = S_{ij} - z_{ij} \ell_{ij} - \sum\limits_{{k:(j,k) \in {\mathcal{E}}}} {S_{jk} }$$
(3)
$$v_{i} - v_{j} = 2\text{Re} (\bar{z}_{ij} S_{ij} ) - \left| {z_{ij} } \right|^{2} \ell_{ij}$$
(4)
$$\ell_{ij} = \frac{{\left| {S_{ij} } \right|^{2} }}{{v_{i} }}$$
(5)

where \(\ell_{ij} \triangleq \left| {I_{ij} } \right|^{2}\) and \(v_{i} \triangleq \left| {V_{i} } \right|^{2}\).

Regulating the voltages and currents to lie within pre-specified lower and upper bounds is another important goal of control, specified by:

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{v}_{i} \le v_{i} \le \bar{v}_{i} \,\,\,\,\,\,\,\,\,\, i \in {\mathcal{N}}$$
(6)
$$\ell_{ij} \le \bar{\ell }_{ij} \,\,\,\,\,\,\,\,\,\,\forall \text{(}i,j\text{)} \in {\mathcal{E}}$$
(7)

Let \(\varvec{x}\) and \(\varvec{u}\) denote the state variables and control variables in above equality and inequality constraints, where \(\varvec{u}\) is the vector of power injections of controllable devices, i.e., \(\varvec{u}{ = }\left\{ {s_{{g_{\text{c}} }} |\forall g_{\text{c}} \in {\mathcal{G}}} \right\}\). Considering that there may be multiple applications of OPF, such as voltage control to prevent overvoltage, var optimization to minimize loss and expansion planning, the objective function is denoted as \(f(\varvec{x},\varvec{u})\) without loss of generality. Thus, deterministic OPF subject to above constraints can be formulated as:

$$\left\{\begin{array}{lll} {\text{min}}&{f(\varvec{x},\varvec{u})}\\ {\text{s.t. }}&{\varvec{h}(\varvec{x},\varvec{u}) = \varvec{s}_{0}}\\ {}&{\varvec{g}(\varvec{x},\varvec{u}) \le {\mathbf{0}}} \end{array}\right.$$
(8)

where \(\varvec{h}(\varvec{x},\varvec{u})\) and \(\varvec{g}(\varvec{x},\varvec{u})\) represent the equality and inequality constraints, respectively; \(\varvec{s}_{0}\) is the given constant vector constructed by loads and uncontrollable power injections.

2.2 OPF model with interval uncertainty

The most significant uncertainties in active distribution networks are associated with the complex dynamics of active/reactive power supply and demand. Therefore, we take into account the uncertainties of intermittent uncontrollable DG output and load demand in this paper. By IA, \(s_{{{\mathcal{L}}_{i} }}\) and \(s_{{g_{\text{u}} ,i}}\) can be extended to \(s_{{{\mathcal{L}}_{i} }}^{ \pm } = [s_{{{\mathcal{L}}_{i} }}^{ - } ,s_{{{\mathcal{L}}_{i} }}^{ + } ]\) and \(s_{{g_{\text{u}} ,i}}^{ \pm } { = }[s_{{g_{\text{u}} ,i}}^{ - } ,s_{{g_{\text{u}} ,i}}^{ + } ]\) respectively, of which "+" represents the upper bound of variable and "−" represents the lower bound. The net load at each bus \(i \in {\mathcal{N}}\) can be expressed as an interval:

$$s_{i}^{ \pm } = s_{{{\mathcal{L}}_{i} }}^{ \pm } - s_{{g_{\text{u}} ,i}}^{ \pm } - \sum\limits_{{g_{\text{c}} \in {\mathcal{G}}_{i} }} {s_{{g_{\text{c}} }} } { = }[s_{i}^{ - } ,s_{i}^{ + } ] = [p_{i}^{ - } + {\text{i}}q_{i}^{ - } ,p_{i}^{ + } + {\text{i}}q_{i}^{ + } ]$$
(9)

After the net loads in interval form are introduced into the balance equations (3)-(5), the state variables will reflect interval uncertainties, thus producing a new form of OPF model incorporating interval data. It has the form:

$$\left\{\begin{array}{lll} {\text{min}}&{f(\varvec{X}^{ \pm } ,\varvec{u})}\\ {\text{s.t. }}&{\varvec{h}(\varvec{X}^{ \pm } ,\varvec{u}) = [\varvec{s}_{0}^{ - } ,\varvec{s}_{0}^{ + } ]}\\ {}&{\varvec{g}(\varvec{X}^{ \pm } ,\varvec{u}) \le {\mathbf{0}}} \end{array}\right.$$
(10)

where \([\varvec{s}_{0}^{ - } ,\varvec{s}_{0}^{ + } ]\) is an interval vector of the power injection for uncontrollable DGs and loads; \(\varvec{X}^{ \pm }\) is the state variable vector in interval form, which is uniquely determined by the control variable vector \(\varvec{u}\) and net load uncertainties.

Due to the wrapping effect of IA, (10) cannot be fully consistent with the actual physical system described by real values. The existence of interval data also makes (10) hard to solve directly with the classical algorithms used in deterministic OPF problems. For these two problems, we need to transform (10) using some additional measures to obtain a more practical and direct description, which are introduced in the following section.

3 Formulation of I-OPF model

To transform (10) into equivalent deterministic one, we first replace interval arithmetic with AA and further introduce an interval Taylor expansion [25] to reconstruct the nonlinear constraints with intervals.

3.1 Introduction of AA and interval Taylor expansion

AA provides a good solution to keep track of correlations between the input and computed quantities, which allows much tighter bounding of solutions than IA in the computing process. The interval \(x^{ \pm } = [x^{ - } ,x^{ + } ]\) can be rewritten as an affine number, a linear combination of the central value and a series of noise terms:

$$\hat{x} = x_{0} + x_{1} \varepsilon_{1} + x_{2} \varepsilon_{2} + \cdots + x_{n} \varepsilon_{n}$$
(11)

where \(x_{0}\) is the central value; \(\varepsilon_{i}\) represents a noise symbol whose value is unknown but lies in the range [−1,1]; \(x_{i}\) is the known real coefficient of partial deviation.

The linear operation \(w = R(x,y)\) may be transformed to the affine operation \(\hat{w} = R(\hat{x},\hat{y})\) which is a composite procedure, computing the affine form \(\hat{w}\) for \(w\) that is consistent with the affine forms \(\hat{x}\) and \(\hat{y}\). Specially, the addition, subtraction, and scalar multiplication operations are given by:

$$\hat{x} \pm \hat{y} = (x_{0} \pm y_{0} ) + (x_{1} \pm y_{1} )\varepsilon_{1} + \cdots (x_{n} \pm y_{n} )\varepsilon_{n}$$
(12)
$$\alpha \hat{x} = (\alpha x_{0} ) + (\alpha x_{1} )\varepsilon_{1} + (\alpha x_{2} )\varepsilon_{2} + \cdots + (\alpha x_{n} )\varepsilon_{n}$$
(13)

where \(\alpha\) is a constant number.

There are non-affine operations that cannot be expressed exactly as an affine combination of noise symbols. Although Chebyshev approximation [26] and other similar approaches may generate approximate results, the approximation process using compound non-affine operations is quite complicated and the result is hard to express in an intuitive form. For easier engineering practice, we compromise by introducing the first-order interval Taylor expansion to construct the affine approximation. For a function \(R(\varvec{X}^{ \pm } )\) containing an interval variable vector, its first-order interval Taylor expansion about the midpoint \(\varvec{X}^{\text{C}}\) in \(\varvec{X}^{ \pm }\) is given by:

$$\begin{aligned} R(\varvec{X}^{ \pm } ) & \approx R(\varvec{X}^{\text{C}} ) + \sum\limits_{i = 1}^{n} {\frac{{\partial R(\varvec{X}^{\text{C}} )}}{{\partial x_{i} }}\left( {X_{i}^{ \pm } - X_{i}^{\text{C}} } \right)} \\ & = R(\varvec{X}^{\text{C}} ) + \sum\limits_{i = 1}^{n} {\frac{{\partial R(\varvec{X}^{\text{C}} )}}{{\partial x_{i} }}\Delta X_{i}^{ \pm } e} \\ \end{aligned}$$
(14)

where \(\varvec{X}^{ \pm } = [X_{1}^{ \pm } ,X_{2}^{ \pm } , \cdots ,X_{n}^{ \pm } ]^{\text{T}}\); \(\varvec{X}^{\text{C}} = [X_{1}^{\text{C}} ,X_{2}^{\text{C}} , \cdots ,X_{n}^{\text{C}} ]^{\text{T}}\); \(e = [ - 1,1]\); \(\Delta X_{i}^{ \pm } = [X_{i}^{\text{C}} - \Delta X_{i} ,X_{i}^{\text{C}} + \Delta X_{i} ]\).

3.2 AA-based I-OPF model

By AA, the net load in interval form can be expressed as:

$$\hat{s}_{i} \triangleq \hat{p}_{i} + {\text{i}}\hat{q}_{i} = (p_{i}^{\text{c}} + p_{i}^{\text{r}} \varepsilon_{p,i} ) + {\text{i}}(q_{i}^{\text{c}} + q_{i}^{\text{r}} \varepsilon_{q,i} )$$
(15)

where \(p_{i}^{\text{c}} \triangleq {{(p_{i}^{ - } + p_{i}^{ + } )} \mathord{\left/ {\vphantom {{(p_{i}^{ - } + p_{i}^{ + } )} 2}} \right. \kern-0pt} 2}\) and \(q_{i}^{\text{c}} \triangleq {{(q_{i}^{ - } + q_{i}^{ + } )} \mathord{\left/ {\vphantom {{(q_{i}^{ - } + q_{i}^{ + } )} 2}} \right. \kern-0pt} 2}\) are the central values of active and reactive net load at the ith bus; \(p_{i}^{\text{r}} \triangleq {{(p_{i}^{ + } - p_{i}^{ - } )} \mathord{\left/ {\vphantom {{(p_{i}^{ + } - p_{i}^{ - } )} 2}} \right. \kern-0pt} 2} \ge 0\) and \(q_{i}^{\text{r}} \triangleq {{(q_{i}^{ + } - q_{i}^{ - } )} \mathord{\left/ {\vphantom {{(q_{i}^{ + } - q_{i}^{ - } )} 2}} \right. \kern-0pt} 2} \ge 0\) correspond to the partial deviations; \(\varepsilon_{p,i}\) and \(\varepsilon_{q,i}\) are noise symbols of the active and reactive net load at the ith bus.

Because (3) and (4) involve only affine operations, namely, a linear form that can be expressed exactly as an affine combination of noise symbols, after uncertainty modeling for state variables, (3) and (4) can be expressed exactly in affine form as:

$$\hat{p}_{j} = \hat{P}_{ij} - r_{ij} \hat{\ell }_{ij} - \sum\limits_{{k:(j,k) \in {\mathcal{E}}}} {\hat{P}_{jk} }$$
(16)
$$\hat{q}_{j} = \hat{Q}_{ij} - x_{ij} \hat{\ell }_{ij} - \sum\limits_{{k:(j,k) \in {\mathcal{E}}}} {\hat{Q}_{jk} }$$
(17)
$$\hat{v}_{j} = \hat{v}_{i} - 2(r_{ij} \hat{P}_{ij} + x_{ij} \hat{Q}_{ij} ) + (r_{ij}^{2} + x_{ij}^{2} )\hat{\ell }_{ij}$$
(18)

where \(\hat{P}_{ij} \triangleq P_{ij}^{\text{c}} + P_{ij}^{\text{r}} \varepsilon_{P,ij}\), \(\hat{Q}_{ij} \triangleq Q_{ij}^{\text{c}} + Q_{ij}^{\text{r}} \varepsilon_{Q,ij}\), \(\hat{\ell }_{ij} \triangleq \ell_{ij}^{\text{c}} + \ell_{ij}^{\text{r}} \varepsilon_{\ell ,ij}\) and \(\hat{v}_{i} \triangleq v_{i}^{\text{c}} + v_{i}^{\text{r}} \varepsilon_{v,i}\) correspond to the affine forms of \(P_{ij}\), \(Q_{ij}\), \(\ell_{ij}\) and \(v_{i}\) respectively.

Notice also that there still exist non-affine operations in (5) and thus it cannot be expressed as affine combinations of the noise symbols. Using the first-order interval Taylor expansion, (5) can be approximated as:

$$\begin{aligned} \hat{\ell }_{ij} & = \ell_{ij}^{\text{c}} + \ell_{ij}^{\text{r}} \varepsilon_{\ell ,ij} = \frac{{\hat{P}_{ij}^{2} + \hat{Q}_{ij}^{2} }}{{\hat{v}_{i} }} \approx \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{v_{i}^{\text{c}} }} \\ \, & \;\; + \frac{{2P_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}P_{ij}^{\text{r}} \varepsilon_{P,ij} + \frac{{2Q_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}Q_{ij}^{\text{r}} (t)\varepsilon_{Q,ij} - \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{(v_{i}^{\text{c}} )^{2} }}v_{i}^{\text{r}} \varepsilon_{v,i} \\ \end{aligned}$$
(19)

Ignoring the noise symbols, the affine equalities described in (16)-(19) can be further transformed into equivalent decoupled equalities as follows:

$$p_{j}^{\text{c}} = P_{ij}^{\text{c}} - r_{ij} \ell_{ij}^{\text{c}} - \sum\limits_{{k:(j,k) \in {\mathcal{E}}}} {P_{jk}^{\text{c}} }$$
(20)
$$q_{j}^{\text{c}} = Q_{ij}^{\text{c}} - x_{ij} \ell_{ij}^{\text{c}} - \sum\limits_{{k:(j,k) \in {\mathcal{E}}}} {Q_{jk}^{\text{c}} }$$
(21)
$$v_{j}^{\text{c}} = v_{i}^{\text{c}} - 2(r_{ij} P_{ij}^{\text{c}} + x_{ij} Q_{ij}^{\text{c}} ) + (r_{ij}^{2} + x_{ij}^{2} )\ell_{ij}^{\text{c}}$$
(22)
$$\ell_{ij}^{\text{c}} = \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{v_{i}^{\text{c}} }}$$
(23)
$$p_{j}^{\text{r}} = P_{ij}^{\text{r}} - r_{ij} \ell_{ij}^{\text{r}} - \sum\limits_{k:j \to k} {P_{jk}^{\text{r}} }$$
(24)
$$q_{j}^{\text{r}} = Q_{ij}^{\text{r}} - x_{ij} \ell_{ij}^{\text{r}} - \sum\limits_{k:j \to k} {Q_{jk}^{\text{r}} }$$
(25)
$$v_{j}^{\text{r}} = v_{i}^{\text{r}} - 2(r_{ij} P_{ij}^{\text{r}} + x_{ij} Q_{ij}^{\text{r}} ) + (r_{ij}^{2} + x_{ij}^{2} )\ell_{ij}^{\text{r}}$$
(26)
$$v_{i}^{\text{c}} \ell_{ij}^{\text{r}} = 2P_{ij}^{\text{c}} P_{ij}^{\text{r}} + 2Q_{ij}^{\text{c}} Q_{ij}^{\text{r}} - \ell_{ij}^{\text{c}} v_{i}^{\text{r}}$$
(27)

The inequality constraint sets found in (6) and (7) can be easily transformed into the following affine forms:

$$\left\{ {\begin{array}{*{20}c} {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{v}_{i} \le v_{i}^{\text{c}} - v_{i}^{\text{r}} \le \bar{v}_{i} } \\ {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{v}_{i} \le v_{i}^{\text{c}} + v_{i}^{\text{r}} \le \bar{v}_{i} } \\ \end{array} } \right.$$
(28)
$$\left\{ {\begin{array}{*{20}c} {\ell_{ij}^{\text{c}} - \ell_{ij}^{\text{r}} \le \bar{\ell }_{ij} } \\ {\ell_{ij}^{\text{c}} + \ell_{ij}^{\text{r}} \le \bar{\ell }_{ij} } \\ \end{array} } \right.$$
(29)

The comparison of interval objectives generally uses their midpoints and radii, where the midpoints of intervals represent their identified basic locations and their radii represent their uncertainty. In (10), the uncertainty of the objective function mainly depends on uncertainties of load demands and of outputs of uncontrollable devices, and this information is the input data. Thus, we can use the midpoint function as an alternative to the interval objective function. Subject to the above constraints in affine form, the AA-based I-OPF model can be finally obtained:

$$\left\{ {\begin{array}{*{20}l} {\hbox{min}\, f(\varvec{X}^{\text{C}} ,\varvec{u}) \, } \hfill \\ {{\text{s.t. }} { (2), (15), (20)-(29)}} \hfill \\ \end{array} } \right.$$
(30)

It should be pointed out that the proposed I-OPF model in this paper is investigated for a balanced distribution network, being a single-phase formulation. Load unbalance is a typical characteristic of distribution networks, especially in low-voltage distribution networks. When I-OPF is required for an unbalanced system, a three-phase I-OPF model is required, which can be obtained by using the three-phase Distflow model derived in [32] and [33] and applying the same modeling process as described in this paper.

I-OPF depends on the special nature of a radial distribution network, namely that its power flows can be specified by a simple set of linear and quadratic equalities if voltage angles are eliminated, and constraints constructed with central values are exact. In addition to satisfying constraints fully under uncertainty, the interval extensions of state variables in the model also make the operational bounds of solution available, which is one of the advantages of I-OPF or interval programming.

4 Enhanced I-OPF model

The I-OPF model in (30) is deterministic so that it can be solved directly by using conventional nonlinear algorithms, such as the interior-point method. Note that its accuracy mainly relies on the approximation of the interval power flow equalities. When the interval uncertainties increase, the I-OPF model will be less accurate due to neglecting the higher order terms in (27). To solve this problem, an enhanced model is developed through modifications discussed in this section.

4.1 Accurate approximation

A more accurate expression equivalent to (19) but including higher order items is:

$$\begin{aligned} \hat{\ell }_{ij} & = \ell_{ij}^{\text{c}} + \ell_{ij}^{\text{r}} \varepsilon_{\ell ,ij} = \frac{{\hat{P}_{ij}^{2} + \hat{Q}_{ij}^{2} }}{{\hat{v}_{i} }} \\ & = \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{v_{i}^{\text{c}} }} + \frac{{2P_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}P_{ij}^{\text{r}} \varepsilon_{P,ij} + \frac{{2Q_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}Q_{ij}^{\text{r}} \varepsilon_{Q,ij} \\ & \;\; - \frac{{\ell_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}v_{i}^{\text{r}} \varepsilon_{v,i} + \frac{{(P_{ij}^{\text{r}} )^{2} }}{{v_{i}^{\text{c}} }}\sigma_{P,ij} + \frac{{(Q_{ij}^{\text{r}} )^{2} }}{{v_{i}^{\text{c}} }}\sigma_{Q,ij} + \ell_{ij}^{\text{c}} \sum\limits_{d = 1}^{\infty } {\left( {\frac{{v_{i}^{\text{r}} }}{{v_{i}^{\text{c}} }}} \right)^{2d} \sigma_{v,i} } \\ & \;\; - \ell_{ij}^{\text{c}} \sum\limits_{d = 1}^{\infty } {\left( {\frac{{v_{i}^{\text{r}} }}{{v_{i}^{\text{c}} }}} \right)^{2d + 1} } \varepsilon_{v,i} + o(\hat{\ell }_{ij} )\varepsilon_{o} + o(\hat{\ell }_{ij} )\sigma_{o} \\ \end{aligned}$$
(31)

where \(\sigma_{P,ij} ,\sigma_{Q,ij}\) and \(\sigma_{v,i}\) represent new noise terms whose values are unknown but lie in the range [0,1]. Since \(\left| {({{v_{i}^{\text{r}} )^{2} } \mathord{\left/ {\vphantom {{v_{i}^{\text{r}} )^{2} } {(v_{i}^{\text{c}} )^{2} }}} \right. \kern-0pt} {(v_{i}^{\text{c}} )^{2} }}} \right| < 1\), the summation terms in (31) can be simplified to:

$$\begin{aligned} \ell_{ij}^{\text{c}} \sum\limits_{d = 1}^{\infty } {\left( {\frac{{v_{i}^{\text{r}} }}{{v_{i}^{\text{c}} }}} \right)^{2d} \sigma_{v,i} } & = \ell_{ij}^{\text{c}} \mathop {\lim }\limits_{d \to \infty } \frac{{{{(v_{i}^{\text{r}} } \mathord{\left/ {\vphantom {{(v_{i}^{\text{r}} } {v_{i}^{\text{c}} )^{2} }}} \right. \kern-0pt} {v_{i}^{\text{c}} )^{2} }}\left[1 - {{(v_{i}^{\text{r}} } \mathord{\left/ {\vphantom {{(v_{i}^{\text{r}} } {v_{i}^{\text{c}} )}}} \right. \kern-0pt} {v_{i}^{\text{c}} )}}^{2d} \right]}}{{1 - {{(v_{i}^{\text{r}} } \mathord{\left/ {\vphantom {{(v_{i}^{\text{r}} } {v_{i}^{\text{c}} )^{2} }}} \right. \kern-0pt} {v_{i}^{\text{c}} )^{2} }}}}\sigma_{v,i} \\ & = \ell_{ij}^{\text{c}} \frac{{(v_{i}^{\text{r}} )^{2} }}{{(v_{i}^{\text{c}} )^{2} - (v_{i}^{\text{r}} )^{2} }}\sigma_{v,i} \\ \end{aligned}$$
(32)
$$\begin{aligned} \ell_{ij}^{\text{c}} \sum\limits_{d = 1}^{\infty } {\left( {\frac{{v_{i}^{\text{r}} }}{{v_{i}^{\text{c}} }}} \right)^{2d + 1} } \varepsilon_{v,i} & = \ell_{ij}^{\text{c}} \mathop {\lim }\limits_{d \to \infty } \frac{{{{(v_{i}^{\text{r}} } \mathord{\left/ {\vphantom {{(v_{i}^{\text{r}} } {v_{i}^{\text{c}} )^{3} }}} \right. \kern-0pt} {v_{i}^{\text{c}} )^{3} }}\left[1 - {{(v_{i}^{\text{r}} } \mathord{\left/ {\vphantom {{(v_{i}^{\text{r}} } {v_{i}^{\text{c}} )}}} \right. \kern-0pt} {v_{i}^{\text{c}} )}}^{2d} \right]}}{{1 - {{(v_{i}^{\text{r}} } \mathord{\left/ {\vphantom {{(v_{i}^{\text{r}} } {v_{i}^{\text{c}} )^{2} }}} \right. \kern-0pt} {v_{i}^{\text{c}} )^{2} }}}}\varepsilon_{v,i} \\ & = \frac{{\ell_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}\frac{{(v_{i}^{\text{r}} )^{3} }}{{(v_{i}^{\text{c}} )^{2} - (v_{i}^{\text{r}} )^{2} }}\varepsilon_{v,i} \\ \end{aligned}$$
(33)

To eliminate the adverse effect of new noise symbols that couple some terms in (31), we transform (31) into:

$$\begin{aligned} \hat{\ell }_{ij} & = \ell_{ij}^{\text{c}} + \ell_{ij}^{\text{r}} \varepsilon_{\ell ,ij} = \frac{{\hat{P}_{ij}^{2} + \hat{Q}_{ij}^{2} }}{{\hat{v}_{i} }} \\ & = \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{v_{i}^{\text{c}} }} + \frac{{\ell_{ij}^{\text{c}} }}{2}\frac{{(v_{i}^{\text{r}} )^{2} }}{{(v_{i}^{\text{c}} )^{2} - (v_{i}^{\text{r}} )^{2} }} + \frac{{(P_{ij}^{\text{r}} )^{2} }}{{2v_{i}^{\text{c}} }} + \frac{{(Q_{ij}^{\text{r}} )^{2} }}{{2v_{i}^{\text{c}} }} \\ & \;\; + \frac{{2P_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}P_{ij}^{\text{r}} \varepsilon_{P,ij} + \frac{{2Q_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}Q_{ij}^{\text{r}} \varepsilon_{Q,ij} - \frac{{\ell_{ij}^{\text{c}} }}{{v_{i}^{\text{c}} }}v_{i}^{\text{r}} \varepsilon_{v,i} \\ & \;\; + \frac{{(P_{ij}^{\text{r}} )^{2} }}{{2v_{i}^{\text{c}} }}\varepsilon_{P,ij} + \frac{{(Q_{ij}^{\text{r}} )^{2} }}{{2v_{i}^{\text{c}} }}\varepsilon_{Q,ij} + \frac{{\ell_{ij}^{\text{c}} (v_{i}^{\text{c}} - 2v_{i}^{\text{r}} )}}{{2v_{i}^{\text{c}} }}\frac{{(v_{i}^{\text{r}} )^{2} }}{{(v_{i}^{\text{c}} )^{2} - (v_{i}^{\text{r}} )^{2} }}\varepsilon_{v,i} \\ & \;\; + o(\hat{\ell }_{ij} )\varepsilon_{o} + o(\hat{\ell }_{ij} )\sigma_{o} \\ \end{aligned}$$
(34)

If we neglect the remainder \(o(\hat{\ell }_{ij} )\varepsilon_{o}\) and \(o(\hat{\ell }_{ij} )\sigma_{o}\) in (34), then (23) and (27) can be rewritten as:

$$\begin{aligned} \ell_{ij}^{\text{c}} & = \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{v_{i}^{\text{c}} }} + \frac{1}{2}\left[ {\frac{{(P_{ij}^{\text{r}} )^{2} + (Q_{ij}^{\text{r}} )^{2} }}{{v_{i}^{\text{c}} }} + \frac{{\ell_{ij}^{\text{c}} (v_{i}^{\text{r}} )^{2} }}{{(v_{i}^{\text{c}} )^{2} - (v_{i}^{\text{r}} )^{2} }}} \right] \\ & = \frac{{(P_{ij}^{\text{c}} )^{2} + (Q_{ij}^{\text{c}} )^{2} }}{{v_{i}^{\text{c}} }} + \rho^{\text{c}} \\ \end{aligned}$$
(35)
$$\begin{aligned} v_{i}^{\text{c}} \ell_{ij}^{\text{r}} & = 2P_{ij}^{\text{c}} P_{ij}^{\text{r}} + 2Q_{ij}^{\text{c}} Q_{ij}^{\text{r}} - \ell_{ij}^{\text{c}} v_{i}^{\text{r}} \\ & \;\; + \frac{1}{2}\left[ {(P_{ij}^{\text{r}} )^{2} + (Q_{ij}^{\text{r}} )^{2} + \frac{{\ell_{ij}^{\text{c}} (v_{i}^{\text{c}} - 2v_{i}^{\text{r}} )(v_{i}^{\text{r}} )^{2} }}{{(v_{i}^{\text{c}} )^{2} - (v_{i}^{\text{r}} )^{2} }}} \right] \\ & = 2P_{ij}^{\text{c}} P_{ij}^{\text{r}} + 2Q_{ij}^{\text{c}} Q_{ij}^{\text{r}} - \ell_{ij}^{\text{c}} v_{i}^{\text{r}} + \rho^{\text{r}} \\ \end{aligned}$$
(36)

4.2 Successive linear approximation-based modification

Replacing (26) and (27) with (35) and (36), the I-OPF model can be modified to take into account the higher order information of partial deviations for state variables. However, this modification also makes the I-OPF model more non-convex and harder to solve. To solve it efficiently, we introduce successive linear approximation (SLA). By assigning fixed initial values to \(v_{i}^{\text{c}}\), \(P_{ij}^{\text{c}}\), \(Q_{ij}^{\text{r}}\), \(\rho^{\text{r}}\)and \(\rho^{\text{c}}\) in (35) and (36), (36) can be converted to a first-order equality and (35) can be relaxed to the following second-order cone constraint:

$$\left\| {\begin{array}{*{20}c} {2P_{ij}^{\text{c}} } \\ {2Q_{ij}^{\text{c}} } \\ {\ell_{ij}^{\text{c}} - \rho^{\text{c}} - v_{i}^{\text{c}} } \\ \end{array} } \right\|_{2} \le \ell_{ij}^{\text{c}} - \rho^{\text{c}} + v_{i}^{\text{c}}$$
(37)

It has been shown in [31] that, when the objective function is convex, (5) can be exactly relaxed to a second order cone constraint. Roughly speaking, such a relaxation for radial networks is exact if the power injection at each bus is not too large and the voltages are kept around their nominal values [34]. For the relaxation shown in (37), the introduction of a small constant does not affect this feature. Due to space limitations, we do not present the detailed proof process, which is available on request from the authors.

Thus, a second-order cone programming (SOCP) model for I-OPF with stepwise modification, called enhanced I-OPF, can be formulated as:

$$\left\{ {\begin{array}{*{20}l} {\hbox{min}\, f(\varvec{X}^{\text{C}} ,\varvec{u}) \, } \hfill \\ {{\text{s}} . {\text{t}} . { (2), (15), (20)-(26),(28),(29),(36),(37)}} \hfill \\ \end{array} } \right.$$
(38)

Based on SLA, the flowchart of the modification method is shown in Fig. 1. The iterative procedure continues until the ∞-norm of the vector difference for the state and control variables between two successive iterations is small enough:

$$\left\| {(\varvec{X}_{k}^{ \pm } ,\varvec{u}_{k} ) - (\varvec{X}_{{k - \text{1}}}^{ \pm } ,\varvec{u}_{{k - \text{1}}} )} \right\|_{\infty } \le \delta$$
(39)

where \(\delta\) is the convergence threshold. To improve the efficiency of the iterative procedure, the (k−1)th optimization result \((\varvec{X}_{{k\text{ - 1}}}^{ \pm } ,\varvec{u}_{{k\text{ - 1}}} )\) is suggested as the initial search point of kth optimization for enhanced I-OPF.

Fig. 1
figure 1

Flow chart of enhanced I-OPF

5 Case studies

This section illustrates the solutions obtained with the proposed I-OPF model and its enhanced method for two test cases. All the tests were implemented using the OPTI and SDPT3 packages in MATLAB, running on an Intel 4-Core i7-CPU 3.40 GHz personal computer, and they were based on a modified IEEE 33-bus system and a real 10 kV 113-bus distribution system.

5.1 Modified IEEE 33-bus system

The standard IEEE 33-bus system [35] was modified by adding two photovoltaic panels (PV1 and PV2) with the installed capacity of 1 MW at buses 18 and 30, one static var compensator (SVC) with a reactive compensation range of [−100, 200] kvar at bus 8, and one energy storage system (ESS) with a power limit of 240 kW at bus 16, as shown in Fig. 2. The system is connected to root bus 1. Reducing network loss was the control objective in this analysis. The tolerances of voltage magnitude were limited within [0.95, 1.05] p.u. and the active power injected at the root bus was constrained within [0, 2.40] MW. Moreover, the predicted value of PV maximum output was set to 0.8 MW, with the power factor of PV output fixed as 0.95, and a tolerance of ±10% was assumed for the load demand at each bus.

Fig. 2
figure 2

Single-line diagram of modified IEEE 33-bus system

Table 1 shows the solutions of OPF (without considering load uncertainties), I-OPF, and its enhanced version with convergence threshold \(\delta\) set to 1 × 10−5. In addition, stochastic OPF [12] was used to calculate a comparison case, denoted by S-OPF, where 50 scenarios were generated for load uncertainty of ±10%. Note that the CPU time in Table 1 represents the iterative time.

Table 1 Optimization results of I-OPF, enhanced I-OPF and S-OPF

The I-OPF and its enhanced version are more efficient than S-OPF. To illustrate their effectiveness on constraint satisfaction, Monte Carlo simulation (MCS) for loads was conducted using forward-backward power flow calculations, where the load at each bus was uniformly distributed in a ±10% interval for 20000 simulation cases. Figure 3 presents the simulated distribution of active power injected at the root bus according to different solutions. Due to neglecting load uncertainty, the deterministic OPF solution may make the active power injected exceed the upper limit. The scenario-based S-OPF does have the ability to handle the uncertainty caused by load variations, but the simulated maximum active power is still over the upper limit. In comparison, I-OPF model treats the uncertainty of power injection as an interval, thus satisfying the boundary constraint well, especially enhanced I-OPF.

Fig. 3
figure 3

Distribution of active power injected at root bus obtained by MCS using solutions of OPF, S-OPF, I-OPF and enhanced I-OPF

The intervals of node voltages and line currents estimated by I-OPF and enhanced I-OPF are shown in Figs. 4 and 5 respectively, both compared with the MCS results. Note that the node voltages and line currents are valued in the per-unit system, and the line number is equal to its terminal node number minus one. Observe that the I-OPF model performs well in boundary estimation of node voltages and line currents, while its enhanced version can obtain more accurate intervals for the state variables, especially the upper bounds of line currents.

Fig. 4
figure 4

Intervals of voltages and currents estimated by I-OPF and MCS

Fig. 5
figure 5

Intervals of voltages and currents estimated by enhanced I-OPF and MCS

To measure the accuracy of solutions for state variables, the upper bound error E +B and the lower bound error E B are defined as follows:

$$\left\{ {\begin{array}{*{20}c} {E_{\text{B}}^{ + } (x) = \left| {X^{ + } - X_{\text{MCS}}^{ + } } \right|} \\ {E_{\text{B}}^{ - } (x) = \left| {X^{ - } - X_{\text{MCS}}^{ - } } \right|} \\ \end{array} } \right.$$
(40)

where \(X_{\text{MCS}}^{ + }\) and \(X_{\text{MCS}}^{ - }\) are the upper and lower bound of state variable \(x\) obtained by MCS, which is used as a reference solution.

Figure 6 presents the lower bound errors (LBE) and upper bound errors (UBE) of the I-OPF model and its enhanced version with modification. These results further demonstrate that the enhanced method provides more accurate estimated results of state variables, especially for the node voltages.

Fig. 6
figure 6

Lower and upper bound errors of I-OPF and its enhanced version with modification

5.2 Real 113-bus system

The 113-bus system that is part of a distribution system in Beijing was used to test the applicability of I-OPF model to larger distribution networks. It is shown in Fig. 7 including some modification from the real system. The compensation range of the SVC was [−300, 600] kvar, the power limit of the ESS was set as 600 kW, and the predicted value of PV maximum output was set to 1.2 MW. Two control objectives, minimum network loss and maximum utilization ratio of DGs, were investigated as discussed in the following subsections.

Fig. 7
figure 7

Single-line diagram of modified real 113-bus system

1) Case 1: network loss objective with the 113-bus system. To illustrate the effect of load uncertainty on the accuracy of the I-OPF model, a tolerance sequence {±5%, ±10%, ±15%, ±20%} was successively assigned to the load demand. The voltage magnitudes were limited within [0.95, 1.05] p.u. and the active power traded between the main grid and the distribution network was constrained within [0, 5.5] MW.

Figure 8 depicts the intervals of active power transmitted in each line solved by enhanced I-OPF, where the load variation at each bus was ±20%, as well as a profile obtained by deterministic OPF without considering load uncertainties. The distributions of active power injected at root bus corresponding to the load variation of ±20% are summarized in Fig. 9. These were obtained by MCS with 20000 simulation cases using the solutions of I-OPF and its enhanced version.

Fig. 8
figure 8

Intervals of active power transmitted in each line obtained by enhanced I-OPF and a profile obtained by deterministic OPF

Fig. 9
figure 9

Distribution of active power injected at root bus obtained by MCS using solution of I-OPF and enhanced I-OPF

The MCS results show that I-OPF and its enhanced version are both applicable to a large distribution system and they satisfy boundary constraints well, in particular, the constraint on the active power injected at the root bus. Also, we notice that there are still some errors in constraint satisfaction for I-OPF shown in Fig. 9. These errors are caused the approximation that ignores higher order information of partial deviations for state variables, which becomes noticeable for the large load variation of ±20%. However, the enhanced I-OPF modified to include higher order terms still works well.

For different degrees of load variation, Table 2 presents the CPU time and the mean network loss of I-OPF and its enhanced version, as well as the mean bound errors (Emean) of node voltages and line currents, which is defined as:

$$E_{\text{mean}} = {\text{mean}}\left\{ {E_{\text{B}}^{ - } (x) + E_{\text{B}}^{ + } (x)} \right\}$$
(41)
Table 2 Mean network loss, CPU time and Emean of I-OPF and enhanced I-OPF

The I-OPF model and its enhanced version are both computationally efficient for the large system. Although the mean bound errors of node voltages and line currents increase as the load variation increases, the enhanced method can maintain the error level within an acceptable range. These characteristics further demonstrate the superiority of the proposed I-OPF model and its enhanced version. Note that I-OPF without modification is better for application to the large system because it is more efficient. On the other hand, enhanced I-OPF can achieve higher accuracy, especially when the load uncertainty is significant. Therefore, when the system has significant scale or the load variation level is not large, we recommend using the original I-OPF, otherwise the enhanced version is recommended.

2) Case 2: utilization ratio objective with the 113-bus system. The primary target of active network management is to improve the utilization ratio of the renewable energy sources. OPF introduced for dispatch in an active distribution network is generally operating on a 5-minute control cycle or longer. When considering the uncertainties of uncontrolled DG and loads within this control interval, some control instructions may allow overvoltage to occur. The I-OPF model can be employed to prevent overvoltage at the system level and simultaneously to make the most of DGs.

This functionality was demonstrated by setting the loads of each bus at 40% of the peak load, given that the overvoltage generally occurs at the access nodes of DGs at low load levels. The uncertainties of loads including the prediction error within the control interval were all set to ±10%, which corresponds to ±4% of the peak load. The upper limit of node voltage for normal operation was set to 1.042 p.u.. For brevity, a simple objective function is defined as:

$$\hbox{min} \left( {\sum\limits_{{g \in \varvec{\mathcal{G}}}} ( {P_{g}^{\text{Pre}} - P_{g} } }) \right)$$
(42)

where \(P_{g}\) is the output power of DG g and \(P_{g}^{\text{Pre}}\) is the maximum value of predicted DG g output.

Figure 10 shows the range of node voltages solved by enhanced I-OPF, as well as a profile obtained by deterministic OPF. To illustrate the accuracy of solved voltage intervals, we also employ MCS with 20000 simulation cases using the solutions of deterministic OPF and its enhanced version, and the voltage distributions for bus 28 are compared in Fig. 11. The solution of OPF creates a condition that voltages at some nodes reach the upper voltage limit. In this case, load variation will easily make the node voltage exceed the limit, as shown in the upper graph of Fig. 11. I-OPF takes account of the boundary constraint considering load uncertainty, so that the node voltages even under extreme conditions are kept below the upper limit. This confirms the effectiveness of I-OPF in overvoltage prevention at the system level, while also satisfying other boundary constraints.

Fig. 10
figure 10

Range of node voltages solved by enhanced I-OPF with modification and a profile obtained by deterministic OPF

Fig. 11
figure 11

Voltage distribution at bus 28 simulated by using the solutions of deterministic OPF and enhanced I-OPF

6 Conclusion

This paper proposes an I-OPF model that incorporates uncertainties in distribution networks, representing uncertainties as intervals. To improve accuracy, an enhanced method based on AA and successive linear approximation is developed. The proposed technique is tested on a modified IEEE 33-bus system and a real 113-bus distribution network. Some conclusions are drawn as follows.

  1. 1)

    Compared with MCS, the I-OPF model and its enhanced version perform well in bounds estimation and constraint satisfaction of state variables. Enhanced I-OPF can offer a more accurate estimate of state variables and thus meet optimization objectives more effectively than unmodified I-OPF.

  2. 2)

    Compared with scenario-based approaches, the I-OPF model and its enhanced version are more efficient. Unmodified I-OPF is better to use for large distribution systems due to its concise form and computational efficiency, while enhanced I-OPF requires additional computing resources to satisfy constraints accurately when there are large load variations.