Convex Stochastic Approaches for the Optimal Allocation of Distributed Energy Resources in AC Distribution Networks with Measurements Fitted to a Continuous Probability Distribution Function

1. Introduction

The Optimal Power Flow problem (AC-OPF) has been an extensively used tool to find optimal operative set points in the planning and operation of power systems. This analysis has been transforming from a rigid, centralized framework, to a more flexible one by including new economic agents and operational strategies. For instance, with the introduction of Peer-To-Peer trading schemes, a consumer-producer (prosumer) is integrated into energy communities and the overall welfare of the community is optimized by coordinating resources and energy transactions with conventional grids [1], or the aggregation of Distributed Energy Resources (DER) into multi-energy virtual power plants, whose operation is economically optimized to provide services such as upstream reactive power support or frequency control [2]. Similarly, the problem of planning virtual microgrids (optimal partitions of a larger conventional distribution network) with the allocation of Distributed Energy Resources (DER) can be addressed with an integrated approach (instead of separating the problem) with the very same constraints imposed by the original AC-OPF [3]. Although the AC-OPF formulates properly a balanced steady state problem for conventional networks, it has been also modified for the design of microgrids with Photovoltaic (PV) units and storage systems to overcome power intermittence from the utility to supply critical loads[4]. This increase in the number of applications for the AC-OPF has been boosted by more efficient and cheaper technologies in the production of electricity (allowing consumers to meet part of their demand and sell surpluses), the integration of information technologies (IT) advancements in the operation within energy supply chain, and data science analyses (e.g., the implementation of IoT solutions for energy management in microgrids[5], the aggregation of Distributed Energy Resources (DER) via IT[6], the demand charcterization for infrastructure planning[7], the statistical characterization of solar irradiation[8] or solar predictions based on Artificial Neural Networks[9]), requiring grid operators, especially in the distribution, to operate considering uncertainties and multi-directional flows of power, currency and information. Under these conditions, network planners and operators implement the AC-OPF to optimize either individual or in a combination, technical, economical or environmental objectives, such as power losses minimization [10], voltage regulation[11], emissions reduction [12], operational and infrastructure investment cost reductions [13,14], or a combination resulting in techno-economic-environmental multi-objective formulations [15]. This requires the study of newer methodologies involving Active Distribution Network modeling [16], adequate mathematical formulation of optimization problems aiming to increase efficiency and solution quality [17] and proper implementation of optimization techniques to solve the problem [18]. Additionally, recent multidisciplinary research has been directed to increase computational efficiency in the solution of the OPF problem [18], showing results particularly in newer meta-heuristic techniques (e.g., Arithmetic Optimization Algorithm (AOA) [19], hybrid Harris Hawks Optimizer-AOA (hHHO-AOA) [20], Mutation Improved Grey Wolf Optimizer (MIGWO) [21]), and in the convexification of the problem with relaxations based on Second-Order Cones [22,23] or Positive Semi-definite expressions [24]. This is relevant since AC-OPF formulations are more complex than the original power balance problem: it might include non-linear cost functions (in addition to the quadratic expression for the losses found in the original AC-OPF fomulation), Mixed-Integer expressions (as in a reconfiguration [25] or resource allocation problem [26]), and/or uncertainty, which individually or all together aren’t necessarily convex when included to the problem (in addition to the non-convex quadratic constraints in the original AC-Power Flow problem), thus limiting the quality of the solutions and the scope of results. Nevertheless, if the necessary effort is made to make the problem convex, the solution is ensured to be optimal, unique, exact (if feasible), and solvable in polynomial time (if not combinatorial). If convex is not possible, sub-optimal solutions are always obtainable with other tools (meta-heuristics), although the solution can´t be guaranteed to be exact or unique. However, such tools are particularly useful in such cases and in multi-objective Pareto optimization [10].

Uncertainty has become relevant in the analysis of power systems since most common DER technologies and demand have stochastic behavior. Uncertainty is typically included in optimization problems by implementing case analysis [11,27], Montecarlo simulations [28], sensitivity analysis [29,30], or typically by using expected value reductions and treating the stochastic problem deterministically. Depending on the model used to describe uncertainty, the stochastic optimization problem can be formulated as a chance-constrained OPF, in which the knowledge of random variables is expressed by distribution functions, and their associated probabilities are constrained to comply operational constraints, e.g., in [31] demand and PV are modeled to with the normal distribution, and the probabilities constraint the problem to avoid voltage, line power flow, and storage power violations. In [32], the probability of frequency deviation violations constraints the primary frequency control problem. On the other hand, the optimization problem can be formulated in multi-stage approaches, by generating representative cases that follow the distributions of random variables, and solving a deterministic problem for each case. For instance, in [16], a two-stage stochastic program is defined for load restoration in distributed microgrids, whose first stage generates the cases with Monte-Carlo simulations, and reduce them using hierarchical clustering. This information is fed to the second stage that implements the AC-OPF to retrieve the solution set for each case. Similarly, in [33] two optimization stages are proposed, but for the first stage uncertainty is modeled by implementing probability distributions, scenarios are generated with the scenario tree method and reduced with the Kantorovich distance method. A different approach for the two-stage stochastic optimization can be seen in [34], in which for both stages, Monte-Carlo sampling is implemented individually. See also [35] for a review in multiple techniques to handle uncertainty.

In contrast to the two stage stochastic programs, in this paper two single stage approaches are proposed to formulate convex stochastic AC-OPF problems regarding DER sizing and location. The allocation problem aims to install a single DER unit (photovoltaic PV) to minimize power losses in the system, considering demand and irradiance random variables that are modeled from measured data fitted to the logistic continuous distribution. The main contributions of this work are summarized as follows:

• The implementation of single-stage convex stochastic programs, that require no additional software implementations, whose results are expected to be global optimum.
• The uncertainty is modeled using a probabilistic approach, based on data available. Its implementation in the optimization problem is done with continuous logistic distribution functions. With this approach, other distribution functions can be used to calculate probabilities, as long as they can be formulated as concave expressions. With this, the formulation of the optimization problem is more flexible to different probability distributions that might fit better the data.
• The implementation of mathematical expressions for probabilities are tailored to follow the sense of the optimization and the nature of the random variable.
• The results show that the deterministic case (based on expected values) represent a less likely case. On the other hand, both proposed approaches are more likely, thus estimating better the long-term expected behavior.

After this introduction, the rest of the manuscript is organized to present the notation first, followed by the formulation of the proposed optimization problems, an overview of the test systems (one of the systems implements a real-world model from a distribution network in San Andres Island in Colombia), and a summary of the resulting AC-OPF problems (Deterministic and Stochastic) with the computational implementation. Then the results, the discussion, and the conclusions are presented.

2. Methods

In this section, the notation used throughout the article is firstly described, followed by the formulation of non-convex AC-OPF and a convex reformulation. Then, the proposed definition of the objective function regarding power losses and probabilities and, subsequently, constraints regarding the formulated PV-Capacity-Irradiance interface are defined. Finally, an overview of the deterministic and the proposed stochastic problems is shown together with the specifications in the computational implementation.

2.1. Notation

A distribution network consisting in the nodes set $N=\{1,2,\dots ,n\}$, and the set of lines $L=\{1,2,\dots ,l\}$ whose elements consist in a pair of nodes $ij|\forall i,j\in N\wedge ij\in L$, can be represented as an undirected graph $G=\{N,L\}$. The reference parameters and are defined for the first node (busbar) of the system, that also hosts the main generator injecting the apparent power $s\in \mathbb{C}$ to supply the demand ${\displaystyle \sum _{i\in N}}{d}_i|d_i\in \mathbb{C}\forall i\in N$ and cover the losses. A node $j\in N$ is a neighbour of node $i\in N$ if the pair $i,j$ are connected through a line $ij\in L$ whose complex admittance, $y_{ij}\in \mathbb{C}$, remains constant. The neighbourhood ${\mathfrak{N}}_i$ of the node $i\in N$, is the set of neighbouring nodes $j|j\in N\wedge ij\in L$. The admittance matrix $Y_G$ (size is $n\times n$) of the network is composed of each line admittance $y_{ij}$ in the neighbourhood ${\mathfrak{N}}_i$. Its diagonal elements are calculated with $Y_{ii}={\displaystyle \sum _{j\in {\mathfrak{N}}_{\mathfrak{i}}}}{y}_{ij}$, while the off-diagonal elements are calculated with $Y_{ij}=\left\{\begin{array}{c}\hfill -y_{ij}|\forall ij\in L\hfill \\ \hfill 0|\forall ij\notin L\hfill \end{array}\right.$. The ∗ operator stands for the complex conjugate. The complex voltage is denoted with v. The squared magnitude of the voltage is represented with the variable u, while for two adjacent nodes $ij\in L$ the product between $v_j^{*}$ and $v_i$ is represented with $w_{ij}$. The Euclidean norm is represented with the operator $\left|\right|\left|\right|$. The symbols $\overline{X}$ and $\underline{X}$ stand for the upper and lower bound respectively. For multiperiod analyses the set $H=\{1,2,\dots ,h\}$ is introduced to extend the size of affected variables, i.e., $X_i\left(h\right)$ for the h period of analysis. In this paper, multiperiod analyses are carried out in hourly steps within a day. Probabilities are assigned to the letter p. In the optimization problems, $f_o$ stands for the objective function. Additional variables and symbols are presented throughout the paper and are explained as they appear. The bolded expression $\mathbf{1}$ represents a row vector of ones, whereas ${\mathbf{1}}^T$ represents the transpose of the row-vector of ones (column vector of ones). The functions $real\left(\right)$ and $imag\left(\right)$ extract the real and imaginary components of the complex input argument.

2.2. AC-OPF formulations

The AC-OPF problem is, as originally formulated by Carpentier, a non-convex Quadratically Constrained Quadratic Program QCQP [36] (omitting objective function definition) since the admittance matrix is not Hermitian. The balance of power can be expressed as in Equation 1.

$$\frac{s_i-d_i}{v_i}=\sum _{j\in N}{Y}_{ij}^{*}{v}_j\ast |\forall i\in N$$

(1)

An alternative formulation can be implemented by taking advantage of the sparsity in the incidence matrix of radial networks. The Bus Injection Model (BIM) considers the current flow through each line (sending and receiving flows) connected to the node as can be observed in Equation 2 [37]. This formulation continues being non-convex.

$$\frac{s_i-d_i}{v_i}=\sum _{j\in {\mathfrak{N}}_i}{(v_i-v_j)}^{*}{y}_{ij}^{*},|\forall i\in N$$

(2)

As can be observed, in Equation 1 the sum takes into account each element within the admittance matrix, whereas in Equation 2 the summation is only done for each line in the neighborhood of the $i^{th}$ node with its respective admittance (it does not consider the admittance matrix, but single element admittance).

To convexify the problem, as proposed in [38] (See [22,39] for more comprehensive information regarding convex relaxations), a linearization is implemented by substituting quadratic expressions with affine ones, and relaxing the equivalent term to form Second Order Conic constraints. The linearization takes place by introducing a new complex variable that replaces quadratic voltage terms, as seen in Equation 3. From these relaxations, convexity can be ensured as long as the convex feasible set lies within the non-convex feasible set.

$$w_{ij}=\left\{\begin{array}{cc}\hfill v_i{v}_j^{*}\hfill & \hfill |\forall i,j\in L\hfill \\ \hfill 0\hfill & \hfill |\forall i,j\notin L\hfill \end{array}\right.$$

(3)

If one builds a matrix W with every $w_{ij}$ element, it is possible to see that W must be hermitian, and the diagonal elements correspond to each node’s squared voltage magnitude. To the diagonal elements of W, a new variable u is assigned. These transformations can be observed in Equations 4,5,6,7 and 8, where the index i is replaced by the index k to leave the i as the imaginary unit.

$$w_{kj}=(v_k^{re}+i{v}_k^{im})(v_j^{re}-i{v}_j^{im})=(v_k^{re}{v}_j^{re}+v_k^{im}{v}_j^{im})+i(v_j^{re}{v}_k^{im}-v_k^{re}{v}_j^{im})$$

(4)

$$w_{jk}=(v_j^{re}+i{v}_j^{im})(v_k^{re}-i{v}_k^{im})=(v_j^{re}{v}_k^{re}+v_j^{im}{v}_k^{im})+i(v_k^{re}{v}_j^{im}-v_j^{re}{v}_k^{im})$$

(5)

$$w_{kj}=w_{jk}^{*}$$

(6)

$$w_{kk}=u_k=(v_k^{re}+i{v}_k^{im})(v_k^{re}-i{v}_k^{im})=(v_k^{re}{v}_k^{re}+v_k^{im}{v}_k^{im})$$

(7)

$$\sqrt{u_k}=\sqrt{{\left(v_k^{re}\right)}^2+{\left(v_k^{im}\right)}^2}=\left|\right|v_k\left|\right|$$

(8)

These transformations are used to create a set of conic constraints (one for each line) taking advantage of the Euclidean norm of complex numbers (Equation 9) and creating a separated expression for the term $u_i{u}_j$ (Equation 10). Combining those equalities, and relaxing them into inequalities, the set of convex conic constraints, namely Second Order Conic Constraints (SOC), is expressed as in Equation 11.

$$\left|\right|w_{ij}{\left|\right|}^2=\left|\right|v_i{\left|\right|}^2\left|\right|v_j{\left|\right|}^2=u_i{u}_j$$

(9)

$$u_i{u}_j=\frac{1}{4}(u_i+u_j)-\frac{1}{4}(u_i-u_j)$$

(10)

$$∥\begin{array}{c}2w_{ij}\\ u_i-u_j\end{array}∥\le u_i+u_j|\forall i,j\in L$$

(11)

The basic convex AC-OPF based on the BIM is then formulated as in Equation 12.

$$\begin{array}{cc}\hfill minimize:& \hfill \hfill \\ \hfill & \hfill f_o\hfill \\ \hfill subjectto:& \hfill \hfill \\ \hfill & \hfill s_i-d_i=\sum _{j\in {\mathfrak{N}}_i}(u_i-w_{ij})y_{ij}^{*},\hfill & \hfill |\forall i\in N\\ \hfill & \hfill ∥\begin{array}{c}2w_{ij}\\ u_i-u_j\end{array}∥\le u_i+u_j\hfill & \hfill |\forall i,j\in L\end{array}$$

(12)

If a multi-period analysis is carried out, the relaxed AC-OPF formulation changes to the expressions in Equation 13.

$$\begin{array}{cc}\hfill minimize:& \hfill \hfill \\ \hfill & \hfill f_o\hfill \\ \hfill subjectto:& \hfill \hfill \\ \hfill & \hfill s_i\left(h\right)-d_i\left(h\right)=\sum _{j\in {\mathfrak{N}}_i}(u_i\left(h\right)-w_{ij}\left(h\right))y_{ij}^{*},\hfill & \hfill |\forall i\in N\wedge \forall h\in H\\ \hfill & \hfill ∥\begin{array}{c}2w_{ij}\left(h\right)\\ u_i\left(h\right)-u_j\left(h\right)\end{array}∥\le u_i\left(h\right)+u_j\left(h\right)\hfill & \hfill |\forall i,j\in L\wedge \forall h\in H\end{array}$$

(13)

Additional operational constraints can be included in the formulation, i.e., voltage limits, generator capacity, and line power flow limits, but in this paper, such restrictions are not considered because uncertainties could make the problem infeasible. Nevertheless, the generation in the reference node is only bounded to be non-negative, as expressed in Equation 14.

$$real\left(s_i\right)\ge 0\wedge imag\left(s_i\right)\ge 0|i=1$$

(14)

2.3. Objective function

Now that the AC-OPF is constrained to convex terms, the objective function must be defined with convex expressions as well. Therefore, a single objective function is defined for two objectives: Minimize power losses and Maximize probabilities.

2.3.1. Energy Losses

The expression for Energy losses $E_{loss}$ is typically non-convex, but by using the aforementioned transformations, it is now convex (affine), as shown in Equation 15.

$$E_{loss}=real\left(\sum _{h\in H}\sum _{i\in N}\sum _{ij\in L}(u_i\left(h\right)-w_{ij}\left(h\right))y_{ij}^{*}\right)$$

(15)

2.3.2. Probabilities

In this paper, the maximization of probabilities for the logistic probability distribution is carried out. The continuous logistic distribution is a two-parameter distribution and has both probability density (pdf) $f(\mu ,\sigma )$, and cumulative probability $F(\mu ,\sigma )$ (cdf) functions, whose definitions are shown in Equations 16 and 17 respectively.

$$f(x;\mu ,\sigma )=\frac{e^{\frac{-x+\mu }{\sigma }}}{\sigma {(1+e^{\frac{-x+\mu }{\sigma }})}^2}$$

(16)

$$F(x;\mu ,\sigma )=\frac{1}{1+e^{\frac{-x+\mu }{\sigma }}}$$

(17)

The first proposed approach to handle uncertainty is called Range Probability Optimization (hereafter named RPO), and is defined to find the value x of the random variable X within an uncertainty range by maximizing its probability of occurrence. Therefore, the maximization of this probability is included in the objective function. To calculate the probability for the random variable following the RPO approach, the cdf is implemented as depicted in Equation 18.

$$maximize:p(a\le X\le b)=F(X\le b)-F(X\le a)$$

(18)

Where $a=(1-ϵ)x$ and $b=(1+ϵ)x$ define the uncertainty range in which the probability is evaluated, and $ϵ$ is the parameter that establishes the half width of the allowed uncertainty. The expression in Equation 19 results from the simplification of Equation 18, and after natural logarithm is applied, it is converted into a sum of concave expressions, as it is shown in equation 20, which are convenient for maximization. The natural logarithm is a monotonically increasing function, thus the maximization of the probability’s logarithm, represents the maximization of the probability itself.

$$p(a\le X\le b)=\frac{1-e^{\frac{a-b-\mu }{\sigma }}}{(1+e^{\frac{-b+\mu }{\sigma }})(1+e^{\frac{-a+\mu }{\sigma }})\left(e^{\frac{a-\mu }{\sigma }}\right)}$$

(19)

$$ln\left(p(a\le X\le b)\right)=ln(1-e^{\frac{a-b-\mu }{\sigma }})-ln(1+e^{\frac{-b+\mu }{\sigma }})-ln(1+e^{\frac{-a+\mu }{\sigma }})-\frac{a-\mu }{\sigma }$$

(20)

Besides the RPO, a different approach named Value-Probability Optimization (hereafter called VPO) is proposed. It includes in the objective function the maximization of the random variable value and its probability of being greater than the value, as described in Equations 21,22 and 23, and the minimization of the value with the maximization of the variable’s probability of being less or equal than its value, as described in Equations 24 and 25.

$$maximize:xp(X>x)$$

(21)

$$p(X>x)=1-\frac{1}{1+e^{\frac{-x+\mu }{\sigma }}}$$

(22)

$$xp(X>x)=x(1-\frac{1}{1+e^{\frac{-x+\mu }{\sigma }}})$$

(23)

$$maximize:p(X\le x)-x,$$

(24)

$$p(X\le x)=\frac{1}{1+e^{\frac{-x+\mu }{\sigma }}}$$

(25)

After simplifications in Equation 23 and by applying the natural logarithm to it and to Equation 25, the resultant concave expressions are shown in Equation 26 and 27.

$$ln\left(xp(X>x)\right)=ln\left(x\right)-ln(1+e^{\frac{x-\mu }{\sigma }})$$

(26)

$$ln\left(p(X\le x)\right)-x=-ln(1+e^{\frac{-x+\mu }{\sigma }})-x$$

(27)

Since the aim of these formulations is to increase probabilities and decrease power losses, the random variables and the calculation of their probabilities must be assigned properly when considering VPO analyses. Generally speaking, if one is determined to reduce power losses, the demand is to be reduced either by managing it directly with demand response schemes or curtailment, or by supplying it with distributed resources. Consequently, and to avoid trivialityin the solution, the demand coefficient is to be minimized and its probability maximized following Equation 25, while the irradiance coefficient and its probability (Equation 22) are both to be maximized. Calculating probabilities otherwise, would end in trivial solutions: if the probability for the demand coefficient is calculated with Equation 22, the maximum probability is obtained when the variable reaches its lower bound while minimizing variable’s value, and if the probability for the irradiance coefficient is calculated with Equation 25, the maximum probability is obtained when the variable reaches its upper bound, while maximizing the variable’s value. Therefore, the demand minimization and its probability maximization are included in the objective function with Equation 27, whereas the maximization of the irradiance coefficient and its probability is included with Equation 26.

2.4. DER and Demand modelling

In this article, the generation through Photovoltaic DER (PV) and the demand are modeled as stochastic processes from real data. The research group Electric Machines and Drives (EM&D) at the National University in Bogotá, Colombia, worked on demand characterization from data collected from Distribution Companies. The collected data for demand was measured hourly from 2018 until 2020. Whereas, the Irradiance (Global Solar Irradiance) hourly measurements were provided by the Colombian Institute of Hydrology, Meteorology, and Environmental Studies (IDEAM) for the period 2015-2018.

To obtain the distribution parameters for both variables, the data were fitted hourly to a Logistic distribution. In Figure 1, the fitted parameters are illustrated (See Table A1 in the appendix for mean $\mu$ and scale $\sigma$ parameter values).

Consequently, irradiance ($I\left(h\right)$) and demand coefficients ($D\left(h\right)$) replace the random variable X seen through Equations 16-27, and modulate the behavior of the PV generation and the demand respectively. Although the demand is typically defined constant for each node in Power Flow Analysis $d_i$, in this article the demand through every node is regulated evenly by the random coefficient$D\left(h\right)$ as it is described in Equation 28.

$$d_i\left(h\right)=d_{i}D\left(h\right)\forall i\in N\wedge \forall h\in H$$

(28)

2.4.1. PV Location modelling

The optimization problem in this article includes the allocation (size and locate) of a PV unit in the system. The location is represented by a boolean variable that indicates if a PV unit is to be installed in a node or not (Equation 29), whereas the size of the PV unit is modeled to deliver active power modulated linearly by the irradiance coefficient ($I\left(h\right)$) with a constrained capacity. The rated capacity of the PV unit can be fully utilized if the irradiance measures $I=1000W/m^2$ (in per unit this would be the base value), although such level was not registered in the measurements. The linearized expression for the power delivered by the PV unit including its location is shown in Equations 30, 31, and 32.

$$\sum _{i\in N}{z}_i=NGmax$$

(29)

$$pp{v}_i\left(h\right)<=z_{i}PGmax$$

(30)

$$pp{v}_i\left(h\right)<=pvh\left(h\right)$$

(31)

$$pp{v}_i\left(h\right)>=pvh\left(h\right)-PGmax(1-z_i)$$

(32)

Where the variables $z_i,pvh\left(h\right)$ and $pp{v}_i\left(h\right)$ represent the location of the PV, the modulated power and the effective power of the located PV unit, respectively, and the parameters $NGmax$ and $PGmax$ stand for the maximum number of PV units to be installed and their maximum capacity, respectively.

2.4.2. PV-Capacity-Irradiance interface

The interface between the PV unit power $pvh\left(h\right)$, its capacity $pvc$ and the irradiance $I\left(h\right)$ is expressed in Equation 33. This interface is named hereafter PCI.

$$pvh\left(h\right)=pvc\times I\left(h\right)$$

(33)

This non-linear expression is linearized using piece-wise linear approximations. The variables $pvc$ and $I\left(h\right)$ are sectioned in k and m equidistant points from their minimum until their maximum values, giving place to $k\times m$ coefficients ${\lambda }_{km}$ interpolating the sectioned values in each variable, as it is shown in Equations 34, 35 and 36.

$$pvc=\sum _k\sum _{m}pv{c}_k{\lambda }_{km}\left(h\right)$$

(34)

$$I\left(h\right)=\sum _k\sum _m{I}_m\left(h\right){\lambda }_{km}\left(h\right)$$

(35)

$$pvh\left(h\right)=\sum _k\sum _{m}pv{h}_{km}\left(h\right){\lambda }_{km}\left(h\right)$$

(36)

Where $pv{c}_k,I_m\left(h\right)$ and $pv{h}_{km}$ are the set of sectioned values for the PV capacity, Irradiance and the resulting modulated power, respectively. The set of ${\lambda }_{km}$ is a Special Ordered Set 2 (SOS2), for which only two adjacent values (in rows and columns) can be non-zero and its whole sum is 1 (Equation 37). To ensure that only two adjacent rows and columns can be non-zero two additional SOS2 sets are defined (one set for rows, the other for columns) such that the sum of their values is also 1, as observed in equations 38 and 39.

$$\sum _k\sum _m{\lambda }_{km}\left(h\right)=1$$

(37)

$${\xi }_k\left(h\right)=\sum _m{\lambda }_{km}\left(h\right)|\sum _k{\xi }_k\left(h\right)=1$$

(38)

$${\omega }_m\left(h\right)=\sum _k{\lambda }_{km}\left(h\right)|\sum _m{\omega }_m\left(h\right)=1$$

(39)

Finally, the two SOS2 (${\xi }_k,{\omega }_m$) include each a set of binary variables (${\delta }_n,{\gamma }_n$) to model the adjacency in the ${\lambda }_{km}$ set. This is shown in Equations 40, 41 42 and 43.

$$\Xi \left(h\right)=\left[\begin{array}{cccccccccc}{\xi }_1\left(h\right)& & & & & -{\delta }_1\left(h\right)& & & & \\ & {\xi }_2\left(h\right)& & & & -{\delta }_1\left(h\right)& -{\delta }_2\left(h\right)& & & \\ & & \ddots & & & & & \ddots & & \\ & & & {\xi }_{n-1}\left(h\right)& & & & & -{\delta }_{n-2}\left(h\right)& -{\delta }_{n-1}\left(h\right)\\ & & & & {\xi }_n\left(h\right)& & & & & -{\delta }_{n-1}\left(h\right)\end{array}\right]\le {\mathbf{0}}^T$$

(40)

$$\sum _{n=1}^{n-1}{\delta }_n\left(h\right)=1$$

(41)

$$\Omega \left(h\right)=\left[\begin{array}{cccccccccc}{\omega }_1\left(h\right)& & & & & -{\gamma }_1\left(h\right)& & & & \\ & {\omega }_2\left(h\right)& & & & -{\gamma }_1\left(h\right)& -{\gamma }_2\left(h\right)& & & \\ & & \ddots & & & & & \ddots & & \\ & & & {\omega }_{n-1}\left(h\right)& & & & & -{\gamma }_{n-2}\left(h\right)& -{\gamma }_{n-1}\left(h\right)\\ & & & & {\omega }_n\left(h\right)& & & & & -{\gamma }_{n-1}\left(h\right)\end{array}\right]\le {\mathbf{0}}^T$$

(42)

$$\sum _{n=1}^{n-1}{\gamma }_n\left(h\right)=1$$

(43)

Mind that the inequality sign ≤ in Equations 39 and 41 is not intended to represent a semi-definite negative constraint, but row-wise inequalities.

2.5. Test systems

For this article, the data of two distribution systems parameterize the AC-OPF formulations: The IEEE33-bus system (hereafter named I33) and the modelled distribution network for the neighborhood "Juan 23" (hereafter named J23) located in San Andres Island (Colombia). J23 system was modelled from data provided by the EM&D research group. The I33 system has a total demand of 3715 $kW$ and 2300 $kVAr$ with 210.9876 $kW$ power losses (see system data in [40]) with a rated voltage of 12.66 $kV$, whereas the J23 system, with 13.6 $kV$ rated voltage, totals 1302 $kW$ and 593 $kVAr$ in demand and 31.30 $kW$ in power losses (see system data in Table A3 in the Appendix). Each respective rated voltage and 1 $MW$ (for both systems) are used as base values. If a multi-period Power Flow Analysis (see Equation 13) is carried out with the demand following the mean averaged hourly profile (${\mu }_D$ in Table 1), then the losses sum up 1.51 $MWh/day$ for the I33 system and 337.25 $kWh/day$ for the J23 system. As can be observed in Figure 2, the I33 and J23 systems have radial and radial/ring configurations respectively.

2.6. Deterministic AC-OPF model

To set reference OPF results, the problem of locating and sizing a PV unit will be also considered deterministically. To do so, the mean profiles for irradiance and demand coeficients ($D\left(h\right)={\mu }_D,I\left(h\right)={\mu }_I$) will be included in the formulation to modulate demand and PV power injections. This model is much simpler than the stochastic one since the AC-OPF would not include the probability expressions in the objective function nor the PV-Irradiance-Capacity interface defined before, but only the PV location modelling. This model is described in Equation 44.

$$\begin{array}{cc}\hfill minimize:& \hfill \hfill \\ \hfill & \hfill E_{loss}\hfill \\ \hfill subjectto:& \hfill \hfill \\ \hfill & \hfill s_i\left(h\right)+pp{v}_i\left(h\right)-d_i\left(h\right)=\sum _{j\in {\mathfrak{N}}_i}(u_i\left(h\right)-w_{ij}\left(h\right))y_{ij}^{*},\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill ∥\begin{array}{c}2w_{ij}\left(h\right)\\ u_i\left(h\right)-u_j\left(h\right)\end{array}∥\le u_i\left(h\right)+u_j\left(h\right)\hfill & \hfill \forall i,j\in L\wedge \forall h\in H\\ \hfill & \hfill d_i\left(h\right)=d_{i}D\left(h\right)\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i\left(h\right)=pp{v}_{i}I\left(h\right)\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i<=z_{i}PGmax\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i<=pvc\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i>=pvc-PGmax(1-z_i)\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill \sum _{i\in N}{z}_i=NGmax\hfill \\ \hfill & \hfill pvc\le PGmax\hfill \\ \hfill & \hfill real\left(s_i\right)\ge 0\hfill & \hfill i=1\\ \hfill & \hfill imag\left(s_i\right)\ge 0\hfill & \hfill i=1\end{array}$$

(44)

2.7. Stochastic AC-OPF model

The formulation including objective function with Equation 20 (RPO) is intended to obtain the maximum probability values within a 10 % range, whereas the AC-OPF with Equations 26 and 27 in the objective (VPO) is intended to find an operating point by minimizing the demand, and maximizing DER injections while maximizing probabilities in both random variables. Therefore, two OPF problems are defined (RPO and VPO), for two test systems (I33 and J23). The complete OPF considering each case’s objective function and some additional operational constraints is detailed in Equation 45.

$$\begin{array}{cc}\hfill minimize:& \hfill \hfill \\ \hfill & \hfill E_{loss}-\sum _{h\in H}(\mathtt{Eq.\; 20}\mathtt{or}\mathtt{Eq.\; 27}.|X=D\left(h\right))-\sum _{h\in H}(\mathtt{Eq.20}.\mathtt{or}\mathtt{Eq.26}|X=I\left(h\right))\hfill \\ \hfill subjectto:& \hfill \hfill \\ \hfill & \hfill s_i\left(h\right)+pp{v}_i\left(h\right)-d_i\left(h\right)=\sum _{j\in {\mathfrak{N}}_i}(u_i\left(h\right)-w_{ij}\left(h\right))y_{ij}^{*},\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill ∥\begin{array}{c}2w_{ij}\left(h\right)\\ u_i\left(h\right)-u_j\left(h\right)\end{array}∥\le u_i\left(h\right)+u_j\left(h\right)\hfill & \hfill \forall i,j\in L\wedge \forall h\in H\\ \hfill & \hfill d_i\left(h\right)=d_{i}D\left(h\right)\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill \underline{u}\le u_i\left(h\right)\le \overline{u}\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i\left(h\right)<=z_{i}PGmax\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i\left(h\right)<=pvh\left(h\right)\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill pp{v}_i\left(h\right)>=pvh\left(h\right)-PGmax(1-z_i)\hfill & \hfill \forall i\in N\wedge \forall h\in H\\ \hfill & \hfill \sum _{i\in N}{z}_i=NGmax\hfill \\ \hfill & \hfill pvc\le PGmax\hfill \\ \hfill & \hfill pvc=\sum _k\sum _{m}pv{c}_k{\lambda }_{km}\left(h\right)\hfill & \hfill \forall h\in H\\ \hfill & \hfill I\left(h\right)=\sum _k\sum _m{I}_m\left(h\right){\lambda }_{km}\left(h\right)\hfill & \hfill \forall h\in H\\ \hfill & \hfill pvh\left(h\right)=\sum _k\sum _{m}pv{h}_{km}\left(h\right){\lambda }_{km}\left(h\right)\hfill & \hfill \forall h\in H\\ \hfill & \hfill \sum _k\sum _m{\lambda }_{km}\left(h\right)=1\hfill & \hfill \forall h\in H\\ \hfill & \hfill {\xi }_k\left(h\right)=\sum _m{\lambda }_{km}\left(h\right)\hfill & \hfill \forall h\in H\\ \hfill & \hfill \sum _k{\xi }_k\left(h\right)=1\hfill & \hfill \forall h\in H\\ \hfill & \hfill {\omega }_m\left(h\right)=\sum _k{\lambda }_{km}\left(h\right)\hfill & \hfill \forall h\in H\\ \hfill & \hfill \sum _m{\omega }_m\left(h\right)=1\hfill & \hfill \forall h\in H\\ \hfill & \hfill \sum _{n=1}^{n-1}{\delta }_n\left(h\right)=1\hfill & \hfill \forall h\in H\\ \hfill & \hfill \sum _{n=1}^{n-1}{\gamma }_n\left(h\right)=1\hfill & \hfill \forall h\in H\\ \hfill & \hfill \Xi \left(h\right)\le 0\hfill & \hfill \forall h\in H\\ \hfill & \hfill \Omega \left(h\right)\le 0\hfill & \hfill \forall h\in H\end{array}$$

(45)

2.8. Computational Implementation

The optimization problems defined in this article are implemented in CVXPY [41,42] through Anaconda’s Python distribution [43] and solved with MOSEK [44]. The following parameters (expressed in per unit) were implemented for the applicable simulation: $ϵ=0.1$, $NGmax=1$, $PGmax=2$, $k=10$ for $0\le pv{c}_k\le PGmax$, $m=10$ for $0\le I_m\left(h\right)\le 1$. Solver parameters were kept in their default values.

3. Results

In this section the main results obtained from the defined optimization problems are shown. For each system, the results for the deterministic AC-OPF problem are first outlined, followed by the results for the stochastic formulations (RPO and VPO) respectively. The RPO and VPO probabilities for the mean coefficients are available in Table A5 in the Appendix.

3.1. Results for the I33 system

3.1.1. Deterministic AC-OPF

An overview on the results of the AC-OPF formulated in Equation 41 for the I33 system is shown in Table 1.

3.1.2. Stochastic OPF

Simulation results for the I33 test system under RPO (Eq. 19 in the objective function) are summarized in Table 2. The demand coefficients profile and their probabilities are depicted in Figure 3 while the irradiance coefficients profile and their probabilities are shown in Figure 4 (See Table A6 in the appendix for detailed hourly modulating coefficients and their associated probabilities).

For simulations under VPO (Equations 24 and 21 for demand and irradiance respectively in the objective function), results are summarized in Table 3 and detailed with coefficients and probabilities in Table A7. Demand coefficients profile and their probabilities are illustrated in Figure 5 while Irradiance coefficients profile and their probabilities are depicted in Figure 6.

3.2. Results for the J23 system

3.2.1. Deterministic AC-OPF

An overview on the results of the AC-OPF formulated in Equation 41 for the I33 system is shown in Table 4.

3.2.2. Stochastic OPF

Simulation results for the J23 test system under RPO (Eq. 19 in the objective function) are summarized in Table 5. The profiles outlined by Demand coefficients and their probabilities are depicted in Figure 7 while the profiles for irradiance and probabilities are shown in Figure 8 (See Table A8 in the appendix for detailed hourly modulating coefficients and their associated probabilities).

For simulations under VPO (Equations 24 and 21 for demand and irradiance respectively in the objective function), results are summarized in Table 6 and detailed with coefficients and probabilities in Table A9. Profiles shaped from Demand coefficients and their probabilities are illustrated in Figure 9 while the profiles constructed from Irradiance coefficients and their probabilities are shown in Figure 10.

4. Discussion

For the I33 system, it could be observed from results that the averaged probabilities during the day are increased by around 1.4 % (RPO) and 35 % (VPO) for demand and by around 1.1 % (RPO) and 20 % (VPO) for irradiance compared to the deterministic case. The location of the PV unit changed from node 8 to nodes 17 and 12 under RPO and VPO respectively. Its capacity was kept near the upper bound in the deterministic and VPO cases, while it suffered an $11\%$ reduction with the RPO, suggesting that, even under uncertainty, DER units should be located in the most congested branch (the active power load is distributed to $11.57\%$ in branch 2-6, $28.93\%$ in branch 7-18, $9.69\%$ in branch 19-22, $25.03\%$ in branch 23-25 and $24.76\%$ in branch 26-33). The minimum power loss was achieved in the deterministic case, followed by a $37\%$ increase for the RPO and a $85\%$ increase for the VPO, approximately, showing how sensitive can be the analysis to uncertain values: Although demand and irradiance profiles suffered small deviations compared to the expected values, the system losses were affected much more negatively compared to the improvements gained in probabilities, even when deploying DER optimally. The results for the VPO can be interpreted as an approximation to the worst-case scenario since, as observed in Figure 5 and Figure 6, the obtained demand and irradiance coefficients are considerably higher and lower than the mean values, respectively, and both probabilities (to have lower demand and greater irradiance coefficients) are very high ($>80\%$ for demand and $>70\%$ for irradiance on average). In other words, the probability of the losses being lesser than $1.91MWh/day$ is on an hourly average greater than $70\%$, compared to the $50\%$ on hourly average resulted from estimating the random variables to the mean values. In contrast to the RPO analysis, a drastic increase in power losses was accompanied with increases in probabilities comparing with the deterministic case. These results allow also the analysis of other operating constraints, such as voltage magnitudes. It could be observed, that the operational voltage limit constraint might be violated (especially the lower bound) if considering scenarios with worsened conditions, as in VPO, i.e., if the demand coefficients from VPO were considered for power flow analysis (without PV), the lower $10\%$ bound would not be complied (see Table 3). However, the installation of the PV unit kept voltage magnitudes within the lower $10\%$ bound. Evidently, tighter bounds in voltages might take the problem to infeasibility regions or, in a practical scenario, make the grid non-compliant of operative constraints under uncertainty, which might represent economic losses due to penalties. For the J23 system, results were somewhat similar to the I33 ones. RPO and VPO probabilities increased by $1.4\%$ and $38\%$ for demand and by $0.9\%$ and $21\%$ for irradiance respectively compared to the deterministic case. This indicates that probabilities are not very affected by the system even though the random variables affect directly the power losses, which at the same time depends on the grid’s configuration and its parameters (size, redundancy, line impedances, load distribution, voltage level). As in the I33, the PV units were located in the most congested branch, in contrast to the deterministic case that located it in the second most congested branch (the active power load is distributed to $20.12\%$ in branch 2-11, $35.71\%$ in branch 12-25, $26.69\%$ in branch 26-43 and $17.46\%$ in branch 28-32). Regarding PV capacity, an important reduction happened when comparing both systems. In the I33 system the PV capacity was kept capped but in the RPO case when it was reduced $11\%$, yet in the J23, PV capacity reached near the $40\%$ of the maximum capacity in the deterministic case, and decreased under uncertainty to the $12\%$ of the maximum with both RPO and VPO approaches. Similar to the results for the I33 system, the J23 system gathered the minimum energy losses from the deterministic case, and were progressively worsening under uncertainty ($26\%$ and $30\%$ increases for the RPO and VPO respectively). Probabilities were increased likewise as in the analyses for the I33. Regarding voltages, the J23 system has significantly lower losses than the I33 system in the power flow analysis (PF). This is also reflected in the voltage envelopes shown in the different OPF analysis, where the minimum voltage remained almost unaltered within the $1\%$ range even in cases with greater demand and lower irradiance (VPO). This explains that the J23 grid has the capability to sustain greater loads under uncertainty, possibly with greater hosting capacities, while operating safely. In Figure A1 and Figure A2 for the I33, and Figure A3 and Figure A4 for the J23 on the Appendix, it is possible to observe how the probability of coefficient profiles under RPO and VPO is compared to the probability under the opposite approach, i.e. the probability of the demand and irradiance profiles obtained from RPO compared to the probability of those profiles under VPO and vice-versa. Those illustrations show that probabilities are greater when calculated with the profiles obtained from their respective analysis than with the opposite. However, for irradiance profiles, the differences between the probabilities are noticeable, especially when comparing VPO probabilities, indicating that RPO irradiance coefficients and mean values are overestimated. Whereas, the VPO profiles show lower RPO probabilities but establish a case with higher probabilities of being more favorable, since the range in which the random variable could fall is much wider. To summarize, mean value estimations have lower probabilities, resulting in important power loss underestimations, whereas, under the RPO approach, a better estimation is obtained even though the irradiance is clearly overestimated. Under the VPO approach, an even more probable case is defined by estimating worsened conditions for the objective, showing that both Deterministic and RPO approaches underestimate the losses with less likely conditions. Finally, the formulation of the proposed convex stochastic programs, allowed to obtain global optimal solutions. Although most of the variables in the proposed stochastic programs are integer, the problem size (directly related to the grid’s size) hadn’t had any considerable impact on the solving time. On the contrary, for the I33 stochastic AC-OPFs, it took longer to find the solution. An explanation for this is that variations in the variables within the PCI (which contains the majority of integer variables and has the same dimensions for each test case), produce greater effects in power losses due to line lengths (impedance) and distribution of greater demands, requiring more intensive Branch-and-Bound steps until the difference between the objective upper and lower bounds, respectively set by the improvements in root relaxation and the incumbent solution, lies within the parameterized relative gap.

5. Conclusions

In this article, two alternatives are proposed to handle uncertainties on convex AC-OPF formulations, whose random variables are based on real measurements of demand and irradiance fitted to the logistic distribution. The effective PV power capability is related to the irradiance and to its capacity, by means of a non-linear expression. To linearize it, a PV-Capacity-Irradiance (PCI) interface was implemented using piece-wise linear approximations. One of the proposed approaches is called Range-Probability Optimization (RPO), and intends to maximize the probability of occurrence for random variables within a range of uncertainty, which is defined as a parameter ($ϵ=\pm 0.1$). The second approach is called Value-Probability Optimization (VPO), in which the values of the random variables are optimized along with the maximization of their probabilities. For the VPO to have sense, and to avoid trivial solutions, random variables and the probability functions are assigned following the sense of the main objective (minimize power losses). Thus, The demand is minimized while maximizing the probability of having lower or equal demand and the irradiance is to be maximized while maximizing the probability of having greater irradiance. The proposed approaches were implemented for two distribution networks: The IEEE 33-bus test system and the modelled distribution network in the "Juan 23" neighborhood in San Andres Island (Colombia), and tested against the deterministic case whose parameters follow the expected values of the distribution. The problem of allocating DER in distribution networks was modified to include uncertainty under RPO and VPO approaches, and as a result, it was possible to observe the problem’s sensitivity to uncertainty. With the RPO, although the variation in the averaged probabilities for demand and irradiance coefficients is very small and the location of the PV unit remained in the most congested branch, the capacity of the PV unit was decreased due to overestimated irradiance and underestimated demand, resulting in increased losses compared to the deterministic case. Meanwhile, with the VPO approach, by increasing the probabilities, the resulting profiles for demand and irradiance coefficients showed negative behaviors towards the minimization of power losses, suggesting that VPO approximates the problem to an scenario with worsened conditions without being the worst-case scenario, which provides robust information for stochastic decision making without underestimating (or overestimating) random variables, or over-sizing DER installations, as the worst-case scenario would suggest. Regarding computational efficiency, the proposed stochastic formulations resulted in the collection of global optimal solutions due to convexity, and the inclusion of mixed-integer expressions hadn’t the expected impact on computation times for combinatorial problems, since systems with bigger size, i.e. J23, were actually solved faster, which is explained by a greater sensitivity of the objectives to the integer variables in systems with greater base losses and deviations in bus voltages, such as in the I33. As future work, it is proposed to perform different RPO analyses varying the range width ($ϵ$), and test the VPO probability of the resulting irradiance profiles to see if the differences between both approaches can be narrowed. Considering that Energy Storage Systems (ESS) are usually implemented in power systems to overcome uncertainties in DER, it is proposed to test both VPO and RPO approaches in the allocation of DER and ESS resources, and test if ESS has effects on probabilities. Additionally, it can be interesting to assign different demand profiles following geographic demand characteristics in the network and to fit the demand and irradiance measured data to other distribution functions to see its impact on the convexity of the formulations and on results themselves. Since the proposed formulations are convex MI-SOCP formulations, it can be interesting to find solutions using meta-heuristic techniques for each system, with convex and non-convex formulations, and see the differences in solution quality and computation times.

Appendix A.1

Table A1. Fitted parameters of logistic Distribution for Demand and Irradiance

h	${\mathit{\mu }}_{\mathit{D}}$	${\mathit{\sigma }}_{\mathit{D}}$	${\mathit{\mu }}_{\mathit{I}}$	${\mathit{\sigma }}_{\mathit{I}}$
1	0.1196	0.0567	0	0
2	0.1165	0.0551	0	0
3	0.1152	0.0541	0	0
4	0.1138	0.0534	0	0
5	0.1151	0.0541	0	0
6	0.1181	0.0573	0	0
7	0.1321	0.0704	0	0
8	0.2111	0.1248	0.0197	0.0082
9	0.4919	0.2152	0.1228	0.0311
10	0.8634	0.1519	0.3214	0.0646
11	1.0000	0.1208	0.5064	0.0963
12	0.9578	0.1064	0.6537	0.1102
13	0.7464	0.1427	0.7540	0.1243
14	0.5029	0.1902	0.7804	0.1315
15	0.6278	0.1950	0.7478	0.1304
16	0.8580	0.1220	0.6398	0.1256
17	0.8989	0.1170	0.4761	0.1060
18	0.8428	0.1346	0.2965	0.0713
19	0.6953	0.1918	0.1180	0.0370
20	0.5249	0.1938	0.0135	0.0068
21	0.2884	0.1268	0	0
22	0.1706	0.0833	0	0
23	0.1397	0.0669	0	0
24	0.1237	0.0592	0	0

Table A1. Fitted parameters of logistic Distribution for Demand and Irradiance

h	${\mathit{\mu }}_{\mathit{D}}$	${\mathit{\sigma }}_{\mathit{D}}$	${\mathit{\mu }}_{\mathit{I}}$	${\mathit{\sigma }}_{\mathit{I}}$
1	0.1196	0.0567	0	0
2	0.1165	0.0551	0	0
3	0.1152	0.0541	0	0
4	0.1138	0.0534	0	0
5	0.1151	0.0541	0	0
6	0.1181	0.0573	0	0
7	0.1321	0.0704	0	0
8	0.2111	0.1248	0.0197	0.0082
9	0.4919	0.2152	0.1228	0.0311
10	0.8634	0.1519	0.3214	0.0646
11	1.0000	0.1208	0.5064	0.0963
12	0.9578	0.1064	0.6537	0.1102
13	0.7464	0.1427	0.7540	0.1243
14	0.5029	0.1902	0.7804	0.1315
15	0.6278	0.1950	0.7478	0.1304
16	0.8580	0.1220	0.6398	0.1256
17	0.8989	0.1170	0.4761	0.1060
18	0.8428	0.1346	0.2965	0.0713
19	0.6953	0.1918	0.1180	0.0370
20	0.5249	0.1938	0.0135	0.0068
21	0.2884	0.1268	0	0
22	0.1706	0.0833	0	0
23	0.1397	0.0669	0	0
24	0.1237	0.0592	0	0

Table A3. J23 distribution network model data.

i	j	${\mathit{R}}_{\mathit{ij}}\left[\mathbf{\Omega }\right]$	${\mathit{X}}_{\mathit{ij}}\left[\mathbf{\Omega }\right]$	${\mathit{P}}_{\mathit{j}}\left[\mathit{MW}\right]$	${\mathit{Q}}_{\mathit{j}}\left[\mathit{MW}\right]$
1	2	0.16715	0.10971	0.03129	0.01425
2	3	0.03768	0.02472	0.03129	0.01425
3	4	0.01704	0.01118	0.03129	0.01425
4	5	0.02891	0.01897	0.02086	0.00950
5	6	0.01808	0.01186	0.03129	0.01425
6	7	0.02551	0.01674	0	0
7	8	0.01723	0.01130	0.04694	0.02138
8	9	0.02391	0.01569	0.03129	0.01425
11	10	0.01067	0.00700	0	0
7	11	0.02905	0.01906	0.03129	0.01425
10	12	0.01181	0.00775	0.03129	0.01425
12	13	0.01027	0.00674	0.04694	0.02138
13	14	0.01506	0.00988	0.04694	0.02138
14	15	0.02613	0.01715	0.03129	0.01425
15	16	0.00610	0.00400	0	0
16	17	0.00528	0.00346	0.03129	0.01425
17	18	0.00278	0.00182	0	0
18	19	0.01193	0.00783	0.04694	0.02138
19	20	0.00741	0.00486	0.04694	0.02138
16	21	0.00827	0.00543	0.04694	0.02138
23	22	0.21027	0.13798	0.04694	0.02138
18	23	0.03503	0.02298	0.04694	0.02138
23	24	0.01012	0.00664	0.03129	0.01425
24	25	0.00371	0.00243	0	0
25	26	0.02282	0.01497	0.02086	0.00950
26	27	0.01808	0.01186	0	0
25	28	0.01047	0.00687	0.04694	0.02138
28	29	0.01921	0.01261	0.04694	0.02138
29	30	0.00544	0.00357	0.01877	0.00855
29	31	0.00751	0.00493	0.06259	0.02850
31	32	0.00699	0.00459	0.04694	0.02138
27	33	0.00699	0.00459	0.04694	0.02138
27	34	0.03341	0.02193	0.03129	0.01425
37	35	0.00931	0.00611	0.03129	0.01425
35	36	0.01037	0.00680	0.03129	0.01425
34	37	0.01782	0.01169	0	0
37	38	0.00618	0.00405	0.03129	0.01425
38	39	0.01395	0.00915	0.02086	0.00950
39	40	0.00917	0.00602	0.04694	0.02138
40	41	0.02503	0.01643	0	0
41	42	0.01629	0.01069	0.04694	0.02138
41	43	0.01463	0.00960	0.03129	0.01425
10	43	0.01384	0.00908	0.03129	0.01425

Table A3. J23 distribution network model data.

i	j	${\mathit{R}}_{\mathit{ij}}\left[\mathbf{\Omega }\right]$	${\mathit{X}}_{\mathit{ij}}\left[\mathbf{\Omega }\right]$	${\mathit{P}}_{\mathit{j}}\left[\mathit{MW}\right]$	${\mathit{Q}}_{\mathit{j}}\left[\mathit{MW}\right]$
1	2	0.16715	0.10971	0.03129	0.01425
2	3	0.03768	0.02472	0.03129	0.01425
3	4	0.01704	0.01118	0.03129	0.01425
4	5	0.02891	0.01897	0.02086	0.00950
5	6	0.01808	0.01186	0.03129	0.01425
6	7	0.02551	0.01674	0	0
7	8	0.01723	0.01130	0.04694	0.02138
8	9	0.02391	0.01569	0.03129	0.01425
11	10	0.01067	0.00700	0	0
7	11	0.02905	0.01906	0.03129	0.01425
10	12	0.01181	0.00775	0.03129	0.01425
12	13	0.01027	0.00674	0.04694	0.02138
13	14	0.01506	0.00988	0.04694	0.02138
14	15	0.02613	0.01715	0.03129	0.01425
15	16	0.00610	0.00400	0	0
16	17	0.00528	0.00346	0.03129	0.01425
17	18	0.00278	0.00182	0	0
18	19	0.01193	0.00783	0.04694	0.02138
19	20	0.00741	0.00486	0.04694	0.02138
16	21	0.00827	0.00543	0.04694	0.02138
23	22	0.21027	0.13798	0.04694	0.02138
18	23	0.03503	0.02298	0.04694	0.02138
23	24	0.01012	0.00664	0.03129	0.01425
24	25	0.00371	0.00243	0	0
25	26	0.02282	0.01497	0.02086	0.00950
26	27	0.01808	0.01186	0	0
25	28	0.01047	0.00687	0.04694	0.02138
28	29	0.01921	0.01261	0.04694	0.02138
29	30	0.00544	0.00357	0.01877	0.00855
29	31	0.00751	0.00493	0.06259	0.02850
31	32	0.00699	0.00459	0.04694	0.02138
27	33	0.00699	0.00459	0.04694	0.02138
27	34	0.03341	0.02193	0.03129	0.01425
37	35	0.00931	0.00611	0.03129	0.01425
35	36	0.01037	0.00680	0.03129	0.01425
34	37	0.01782	0.01169	0	0
37	38	0.00618	0.00405	0.03129	0.01425
38	39	0.01395	0.00915	0.02086	0.00950
39	40	0.00917	0.00602	0.04694	0.02138
40	41	0.02503	0.01643	0	0
41	42	0.01629	0.01069	0.04694	0.02138
41	43	0.01463	0.00960	0.03129	0.01425
10	43	0.01384	0.00908	0.03129	0.01425

Table A5. Operation coefficients and Probabilities for mean profiles (Table A1) according to Equations 19 and 22 for demand and irradiance respectively.

h	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}={\mathit{\mu }}_{\mathit{D}}$	$\mathit{p}(\mathit{X}\le \mathit{x})$ $\mathit{x}={\mathit{\mu }}_{\mathit{D}}$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}={\mathit{\mu }}_{\mathit{I}}$	$\mathit{p}(\mathit{X}>\mathit{x})$ $\mathit{x}={\mathit{\mu }}_{\mathit{I}}$
1	0.1050	0.5000	1.0000	1.0000
2	0.1052	0.5000	1.0000	1.0000
3	0.1060	0.5000	1.0000	1.0000
4	0.1061	0.5000	1.0000	1.0000
5	0.1059	0.5000	1.0000	1.0000
6	0.1026	0.5000	1.0000	1.0000
7	0.0935	0.5000	1.0000	1.0000
8	0.0843	0.5000	0.1195	0.5000
9	0.1137	0.5000	0.1950	0.5000
10	0.2766	0.5000	0.2436	0.5000
11	0.3915	0.5000	0.2570	0.5000
12	0.4219	0.5000	0.2882	0.5000
13	0.2556	0.5000	0.2943	0.5000
14	0.1314	0.5000	0.2888	0.5000
15	0.1596	0.5000	0.2790	0.5000
16	0.3375	0.5000	0.2493	0.5000
17	0.3662	0.5000	0.2209	0.5000
18	0.3031	0.5000	0.2050	0.5000
19	0.1792	0.5000	0.1580	0.5000
20	0.1346	0.5000	0.0990	0.5000
21	0.1131	0.5000	1.0000	1.0000
22	0.1020	0.5000	1.0000	1.0000
23	0.1040	0.5000	1.0000	1.0000
24	0.1040	0.5000	1.0000	1.0000
${\displaystyle \sum _{h\in H}}p$	4.3037	12.0000	2.8975*	6.5*
${\mu }_p$	0.1793	0.5000	0.2228*	0.5*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A5. Operation coefficients and Probabilities for mean profiles (Table A1) according to Equations 19 and 22 for demand and irradiance respectively.

h	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}={\mathit{\mu }}_{\mathit{D}}$	$\mathit{p}(\mathit{X}\le \mathit{x})$ $\mathit{x}={\mathit{\mu }}_{\mathit{D}}$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}={\mathit{\mu }}_{\mathit{I}}$	$\mathit{p}(\mathit{X}>\mathit{x})$ $\mathit{x}={\mathit{\mu }}_{\mathit{I}}$
1	0.1050	0.5000	1.0000	1.0000
2	0.1052	0.5000	1.0000	1.0000
3	0.1060	0.5000	1.0000	1.0000
4	0.1061	0.5000	1.0000	1.0000
5	0.1059	0.5000	1.0000	1.0000
6	0.1026	0.5000	1.0000	1.0000
7	0.0935	0.5000	1.0000	1.0000
8	0.0843	0.5000	0.1195	0.5000
9	0.1137	0.5000	0.1950	0.5000
10	0.2766	0.5000	0.2436	0.5000
11	0.3915	0.5000	0.2570	0.5000
12	0.4219	0.5000	0.2882	0.5000
13	0.2556	0.5000	0.2943	0.5000
14	0.1314	0.5000	0.2888	0.5000
15	0.1596	0.5000	0.2790	0.5000
16	0.3375	0.5000	0.2493	0.5000
17	0.3662	0.5000	0.2209	0.5000
18	0.3031	0.5000	0.2050	0.5000
19	0.1792	0.5000	0.1580	0.5000
20	0.1346	0.5000	0.0990	0.5000
21	0.1131	0.5000	1.0000	1.0000
22	0.1020	0.5000	1.0000	1.0000
23	0.1040	0.5000	1.0000	1.0000
24	0.1040	0.5000	1.0000	1.0000
${\displaystyle \sum _{h\in H}}p$	4.3037	12.0000	2.8975*	6.5*
${\mu }_p$	0.1793	0.5000	0.2228*	0.5*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A6. Operation coefficients and Probabilities according Equation 19 obtained from the RPO-OPF on the I33.

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{D}\left(\mathit{h}\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{I}\left(\mathit{h}\right)$
1	0.1613	0.1240	0	1
2	0.1570	0.1242	0	1
3	0.1547	0.1248	0	1
4	0.1528	0.1249	0	1
5	0.1546	0.1248	0	1
6	0.1608	0.1219	0	1
7	0.1875	0.1141	0	1
8	0.310	0.1063	0.0253	0.1368
9	0.6238	0.1314	0.1375	0.2064
10	0.9001	0.2839	0.3473	0.2528
11	1.0184	0.3960	0.5436	0.2656
12	0.9745	0.4262	0.6917	0.2958
13	0.7914	0.2642	0.7946	0.3018
14	0.6166	0.1474	0.8229	0.2959
15	0.7236	0.1731	0.7921	0.2869
16	0.8848	0.3436	0.6886	0.2583
17	0.9207	0.3714	0.5224	0.2310
18	0.8734	0.3097	0.3333	0.2157
19	0.7701	0.1910	0.1388	0.1718
20	0.6282	0.1500	0.0186	0.1188
21	0.3736	0.1311	0	1
22	0.2323	0.1214	0	1
23	0.1890	0.1232	0	1
24	0.1674	0.1232	0	1
	${\displaystyle \sum _{h\in H}}p$	4.6529		3.0383*
	${\mu }_p$	0.1938		0.2337*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A6. Operation coefficients and Probabilities according Equation 19 obtained from the RPO-OPF on the I33.

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{D}\left(\mathit{h}\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{I}\left(\mathit{h}\right)$
1	0.1613	0.1240	0	1
2	0.1570	0.1242	0	1
3	0.1547	0.1248	0	1
4	0.1528	0.1249	0	1
5	0.1546	0.1248	0	1
6	0.1608	0.1219	0	1
7	0.1875	0.1141	0	1
8	0.310	0.1063	0.0253	0.1368
9	0.6238	0.1314	0.1375	0.2064
10	0.9001	0.2839	0.3473	0.2528
11	1.0184	0.3960	0.5436	0.2656
12	0.9745	0.4262	0.6917	0.2958
13	0.7914	0.2642	0.7946	0.3018
14	0.6166	0.1474	0.8229	0.2959
15	0.7236	0.1731	0.7921	0.2869
16	0.8848	0.3436	0.6886	0.2583
17	0.9207	0.3714	0.5224	0.2310
18	0.8734	0.3097	0.3333	0.2157
19	0.7701	0.1910	0.1388	0.1718
20	0.6282	0.1500	0.0186	0.1188
21	0.3736	0.1311	0	1
22	0.2323	0.1214	0	1
23	0.1890	0.1232	0	1
24	0.1674	0.1232	0	1
	${\displaystyle \sum _{h\in H}}p$	4.6529		3.0383*
	${\mu }_p$	0.1938		0.2337*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A7. Operation coefficients and Probabilities according Equations 25 and 22 (for demand and irradiance respectively) from the VPO-OPF on the I33

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}\le \mathit{D}(\mathit{h}\left)\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}>\mathit{I}(\mathit{h}\left)\right)$
1	0.2731	0.9373	0	1
2	0.2675	0.9391	0	1
3	0.2645	0.9403	0	1
4	0.2620	0.9411	0	1
5	0.2644	0.9404	0	1
6	0.2726	0.9366	0	1
7	0.3054	0.9212	0	1
8	0.4319	0.8542	0.0181	0.5472
9	0.7025	0.7268	0.0990	0.6821
10	1.0613	0.7862	0.2553	0.7355
11	1.1890	0.8269	0.4025	0.7462
12	1.1455	0.8537	0.5164	0.7765
13	0.9625	0.8196	0.5910	0.7876
14	0.7426	0.7790	0.6072	0.7885
15	0.8557	0.7628	0.5844	0.7776
16	1.0581	0.8373	0.5064	0.7430
17	1.0912	0.8380	0.3819	0.7086
18	1.0395	0.8115	0.2395	0.6899
19	0.8941	0.7381	0.0996	0.6214
20	0.7339	0.7461	0.0135	0.4977
21	0.5057	0.8471	0	1
22	0.3584	0.9049	0	1
23	0.3079	0.9251	0	1
24	0.2810	0.9343	0	1
	${\displaystyle \sum _{h\in H}}p$	20.5486		9.1023*
	${\mu }_p$	0.8561		0.7001*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A7. Operation coefficients and Probabilities according Equations 25 and 22 (for demand and irradiance respectively) from the VPO-OPF on the I33

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}\le \mathit{D}(\mathit{h}\left)\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}>\mathit{I}(\mathit{h}\left)\right)$
1	0.2731	0.9373	0	1
2	0.2675	0.9391	0	1
3	0.2645	0.9403	0	1
4	0.2620	0.9411	0	1
5	0.2644	0.9404	0	1
6	0.2726	0.9366	0	1
7	0.3054	0.9212	0	1
8	0.4319	0.8542	0.0181	0.5472
9	0.7025	0.7268	0.0990	0.6821
10	1.0613	0.7862	0.2553	0.7355
11	1.1890	0.8269	0.4025	0.7462
12	1.1455	0.8537	0.5164	0.7765
13	0.9625	0.8196	0.5910	0.7876
14	0.7426	0.7790	0.6072	0.7885
15	0.8557	0.7628	0.5844	0.7776
16	1.0581	0.8373	0.5064	0.7430
17	1.0912	0.8380	0.3819	0.7086
18	1.0395	0.8115	0.2395	0.6899
19	0.8941	0.7381	0.0996	0.6214
20	0.7339	0.7461	0.0135	0.4977
21	0.5057	0.8471	0	1
22	0.3584	0.9049	0	1
23	0.3079	0.9251	0	1
24	0.2810	0.9343	0	1
	${\displaystyle \sum _{h\in H}}p$	20.5486		9.1023*
	${\mu }_p$	0.8561		0.7001*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A8. Operation coefficients and Probabilities according Equation 19 obtained from the RPO-OPF on the J23

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{D}\left(\mathit{h}\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{D}\left(\mathit{h}\right)$
1	0.1615	0.1240	0	1
2	0.1572	0.1242	0	1
3	0.1549	0.1248	0	1
4	0.1530	0.1249	0	1
5	0.1548	0.1248	0	1
6	0.1611	0.1219	0	1
7	0.1879	0.1141	0	1
8	0.3157	0.1063	0.0253	0.1368
9	0.6415	0.1317	0.1374	0.2064
10	0.9149	0.2846	0.3333	0.2503
11	1.0295	0.3968	0.5419	0.2657
12	0.9820	0.4267	0.6904	0.2958
13	0.7985	0.2644	0.7946	0.3018
14	0.7209	0.1383	0.8888	0.2804
15	0.7338	0.1733	0.7777	0.2861
16	0.8925	0.3439	0.6666	0.2566
17	0.9294	0.3720	0.5208	0.2310
18	0.8849	0.3103	0.3286	0.2159
19	0.7902	0.1916	0.1385	0.1718
20	0.6448	0.1503	0.0186	0.1188
21	0.3773	0.1312	0	1
22	0.2332	0.1214	0	1
23	0.1894	0.1232	0	1
24	0.1677	0.1232	0	1
	${\displaystyle \sum _{h\in H}}p$	4.6490		3.0180*
	${\mu }_p$	0.1937		0.2321*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A8. Operation coefficients and Probabilities according Equation 19 obtained from the RPO-OPF on the J23

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{D}\left(\mathit{h}\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{a}\le \mathit{X}\le \mathit{b})$ $\mathit{x}=\mathit{D}\left(\mathit{h}\right)$
1	0.1615	0.1240	0	1
2	0.1572	0.1242	0	1
3	0.1549	0.1248	0	1
4	0.1530	0.1249	0	1
5	0.1548	0.1248	0	1
6	0.1611	0.1219	0	1
7	0.1879	0.1141	0	1
8	0.3157	0.1063	0.0253	0.1368
9	0.6415	0.1317	0.1374	0.2064
10	0.9149	0.2846	0.3333	0.2503
11	1.0295	0.3968	0.5419	0.2657
12	0.9820	0.4267	0.6904	0.2958
13	0.7985	0.2644	0.7946	0.3018
14	0.7209	0.1383	0.8888	0.2804
15	0.7338	0.1733	0.7777	0.2861
16	0.8925	0.3439	0.6666	0.2566
17	0.9294	0.3720	0.5208	0.2310
18	0.8849	0.3103	0.3286	0.2159
19	0.7902	0.1916	0.1385	0.1718
20	0.6448	0.1503	0.0186	0.1188
21	0.3773	0.1312	0	1
22	0.2332	0.1214	0	1
23	0.1894	0.1232	0	1
24	0.1677	0.1232	0	1
	${\displaystyle \sum _{h\in H}}p$	4.6490		3.0180*
	${\mu }_p$	0.1937		0.2321*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A9. Operation coefficients and Probabilities according Equations 25 and 22 (for demand and irradiance respectively) from the VPO-OPF on the J23

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}\le \mathit{D}(\mathit{h}\left)\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}\le \mathit{I}(\mathit{h}\left)\right)$
1	0.2790	0.9431	0	1
2	0.2732	0.9447	0	1
3	0.2700	0.9457	0	1
4	0.2673	0.9464	0	1
5	0.2698	0.9458	0	1
6	0.2786	0.9425	0	1
7	0.3137	0.9293	0	1
8	0.4536	0.8745	0.0181	0.5477
9	0.7675	0.7825	0.0986	0.6850
10	1.1220	0.8457	0.2221	0.8227
11	1.2378	0.8773	0.3968	0.7572
12	1.1827	0.8922	0.4444	0.8698
13	1.0004	0.8556	0.5554	0.8316
14	0.7765	0.8081	0.6667	0.7035
15	0.9014	0.8026	0.5850	0.7769
16	1.0970	0.8762	0.5019	0.7498
17	1.1335	0.8813	0.3767	0.7186
18	1.0912	0.8634	0.2221	0.7395
19	0.9682	0.8057	0.0989	0.6259
20	0.7985	0.8040	0.0135	0.4983
21	0.5322	0.8722	0	1
22	0.3701	0.9163	0	1
23	0.3159	0.9329	0	1
24	0.2874	0.9406	0	1
	${\displaystyle \sum _{h\in H}}p$	21.2298		9.3271*
	${\mu }_p$	0.8845		0.7174*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Table A9. Operation coefficients and Probabilities according Equations 25 and 22 (for demand and irradiance respectively) from the VPO-OPF on the J23

h	$\mathit{D}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}\le \mathit{D}(\mathit{h}\left)\right)$	$\mathit{I}\left(\mathit{h}\right)$	$\mathit{p}(\mathit{X}\le \mathit{I}(\mathit{h}\left)\right)$
1	0.2790	0.9431	0	1
2	0.2732	0.9447	0	1
3	0.2700	0.9457	0	1
4	0.2673	0.9464	0	1
5	0.2698	0.9458	0	1
6	0.2786	0.9425	0	1
7	0.3137	0.9293	0	1
8	0.4536	0.8745	0.0181	0.5477
9	0.7675	0.7825	0.0986	0.6850
10	1.1220	0.8457	0.2221	0.8227
11	1.2378	0.8773	0.3968	0.7572
12	1.1827	0.8922	0.4444	0.8698
13	1.0004	0.8556	0.5554	0.8316
14	0.7765	0.8081	0.6667	0.7035
15	0.9014	0.8026	0.5850	0.7769
16	1.0970	0.8762	0.5019	0.7498
17	1.1335	0.8813	0.3767	0.7186
18	1.0912	0.8634	0.2221	0.7395
19	0.9682	0.8057	0.0989	0.6259
20	0.7985	0.8040	0.0135	0.4983
21	0.5322	0.8722	0	1
22	0.3701	0.9163	0	1
23	0.3159	0.9329	0	1
24	0.2874	0.9406	0	1
	${\displaystyle \sum _{h\in H}}p$	21.2298		9.3271*
	${\mu }_p$	0.8845		0.7174*

^* When irradiance is zero, its probability is not accounted in the sum or the average.

Figure A1. RPO probabilities for profiles obtained from RPO and VPO analyses for the I33. Left: Demand. Right: Irradiance. system

Figure A1. RPO probabilities for profiles obtained from RPO and VPO analyses for the I33. Left: Demand. Right: Irradiance. system

Figure A2. VPO probabilities for profiles obtained from RPO and VPO analyses for the I33 system. Left: Demand. Right: Irradiance.

Figure A2. VPO probabilities for profiles obtained from RPO and VPO analyses for the I33 system. Left: Demand. Right: Irradiance.

Figure A3. RPO probabilities for profiles obtained from RPO and VPO analyses for the J23. Left:Demand. Right: Irradiance. system

Figure A3. RPO probabilities for profiles obtained from RPO and VPO analyses for the J23. Left:Demand. Right: Irradiance. system

Figure A4. VPO probabilities for profiles obtained from RPO and VPO analyses for the J23 system. Left: Demand. Right: Irradiance.

Figure A4. VPO probabilities for profiles obtained from RPO and VPO analyses for the J23 system. Left: Demand. Right: Irradiance.