Adaptive data driven multi period power supply recovery method for distribution networks

In recent years, due to the frequent occurrence of natural disasters and human-induced disruptions, large-scale blackout incidents have posed significant challenges to the safe operation of distribution networks in various countries. To address the crisis of blackouts under different circumstances, countries have begun to establish resilient grids, and enhancing the load recovery capability of power systems is a means to improve grid resilience¹.

With the development of new power systems, the penetration rate of new energy represented by photovoltaic (PV) has gradually increased, and it has become possible to provide short-term support in the power supply recovery of distribution networks². However, compared with traditional fossil fuel power generation, new energy sources are characterized by volatility and randomness, bringing unpredictable disturbance risks to grid operation. Compared with traditional thermal power generation, the output of new energy power generation is highly volatile and difficult to predict, which poses new challenges to the traditional fault recovery methods of power systems^3,4.

At the same time, due to the rapid development of power communication technology in recent years, information and communication technology greatly promotes the efficient operation of power systems, affecting various aspects of power system operations⁵. With the advancement of communication technology, the rapid recovery capability of distribution networks has been greatly improved. The close integration of information systems and power systems is gradually evolving modern power systems into highly intelligent information-physical power systems. This close integration not only enhances the real-time perception and precise control capabilities of distribution networks but also makes the grid more complex and diverse. In the process of power supply recovery in distribution networks, the normal operation of terminal communication equipment corresponding to power nodes plays a critical role^6,7. Therefore, in the above context, considering the factors related to the operation of information networks and utilizing new energy sources and other renewable energy for power system recovery is of great significance.

Power supply recovery of distribution networks can be achieved by restoring the isolated areas resulting from power outages⁸. increased the load recovery of distribution networks by enhancing the utilization rate of distributed generators (DG) after detecting the isolated areas⁹. prioritized power supply recovery for critical loads using dynamic scheduling of DG islands and controllable load alliances, and simplifies the solution process with feeder equivalence and fault equivalence methods¹⁰. proposed constraint simplification methods suitable for the distributed power supply recovery network topology search process after determining the island operation areas. For areas isolated due to faults in the distribution system, soft open point (SOP) technology can effectively provide voltage support for the power outage areas^11,12 established SOP planning and operation models to further enhance power supply recovery levels after analyzing the characteristics and functions of SOP.

In order to cope with the uncertainty of distributed PV output during the power supply recovery process, there are two main types of uncertainty optimization methods for distribution network scheduling: stochastic optimization methods and robust optimization methods^13,14. Robust optimization methods usually use an ensemble form to describe the range of distribution of uncertain parameters, which is getting more and more attention as it does not need to obtain the probability distribution of uncertain parameters compared to stochastic methods and avoids the high-dimensional problem introduced due to a large number of scenarios.

However, different sets have different impacts on the optimization results of distribution network power supply recovery. Therefore, selecting an appropriate set not only reduces the conservativeness of the optimization results but also ensures the robustness of the results. In recent years, some scholars have proposed constructing uncertain sets based on historical data by considering the correlation between uncertain parameters and their random changes over time. In addition to constructing sets using polyhedral sets and ellipsoidal sets, one common method is to construct sets based on extreme scenarios^15,16,17,18. This method first performs envelope processing on the historical data of uncertain parameters, then selects extreme scenarios from the historical data, and constructs sets that include these scenarios. By introducing these extreme scenarios into the optimization model, the model is more robust and includes all possible historical data^19,20 presented an example of constructing such sets. This method addressed the conservativeness of polyhedral sets but still faces the problem of overly conservative results when dealing with extreme scenarios^21,22 proposed a data driven method to improve the conservativeness of polyhedral sets, but constructing sets based on extreme scenarios may still face the problem of computational complexity.

Most of the current research on power supply recovery is limited to topological research and lacks consideration of new uncertainties brought by new energy sources. Therefore, this paper proposes a multi period power supply recovery method of distribution network based on adaptive data drive to address the shortcomings of the above set. This paper first proposes a convex polyhedral set constructed based on historical scenario envelopes, then constructs a convex polyhedral set from the connection vertices of the historical scenario envelopes, and finally addresses the problem of the conservativeness of the convex polyhedral set through contraction to include all historical scenarios. Using the concept of super-planes, this paper constructs a data driven hyperplane polyhedral set. At the same time, considering the establishment of an information-physical model with information nodes, a distribution network power supply recovery model based on the information-physical model is constructed, and the effectiveness of the proposed method is verified through simulations.

1. (1)

   Convex hull polyhedral set.

Traditional uncertainty modeling methods mainly include the probability distribution method and the scenario analysis method. The probability distribution method assumes that the output of new energy sources follows a certain known probability distribution. However, the actual output of new energy sources is affected by a variety of complex factors, and its distribution often does not fully conform to these assumptions. The ellipsoidal uncertainty set and the convex hull polyhedral set are constructed based on actual historical data, which can more accurately reflect the true uncertainty of the new energy output and avoid the modeling bias caused by the wrong assumption of the probability distribution. The scenario analysis method usually simulates the uncertainty by selecting a limited number of representative scenarios. It is very difficult for this method to comprehensively cover all possible situations, and it is easy to miss some extreme but possible scenarios. The convex hull polyhedral set can cover all possible states of new energy output based on historical data, has a wider coverage of scenarios, and can deal with uncertainties more comprehensively.

In order to adapt to the real-time changes in the power output of distributed PV, this paper, based on the power quality detection devices installed at the nodes of the distribution network, captures the real-time signals triggered by sudden weather changes and unexpected power grid faults. By establishing a dedicated data transmission and processing channel, the historical data of the power output of distributed PV under various weather conditions are obtained. Then, according to the scatter plot formed by the historical data of the power output of distributed PV, different envelopes can be used to represent different sets, namely the box set and the ellipsoidal set, as shown in Fig. 1. For different sets, this paper uses the uncertain set \({\mathbf{U}}\) to represent the fluctuation range of distributed PV output.

1) Box set.

The specific expression is:

where \({N_{{\text{PV}}}}\) represents the number of PV connections;\({{\mathbf{z}}^{{\text{PV}}}}\)is the uncertain variable of PV output,\({\mathbf{z}}_{{{\text{up}}}}^{{{\text{PV}}}}\)and\({\mathbf{z}}_{{{\text{down}}}}^{{{\text{PV}}}}\)respectively represent the upper and lower bounds of the uncertain variable of PV output.\(\beta\) represents the tuning parameter used to adjust the retention degree of the box set within\((0,1]\).

As shown in Fig. 1 (a), the box set covers all possible values of distributed PV output. However, at different times and locations, the output of distributed PV tends to follow a certain trend and shows spatial correlations. Therefore, PV output data often lie along the y = x and y=-x diagonals. If the uncertain set constructed by the box set fails to reflect the temporal and spatial correlation of distributed PV output, it may result in overly conservative predictions. This is because the box set not only covers the potential fluctuations of PV output but also includes a large number of empty regions where fluctuations are unlikely to occur. Therefore, a more accurate construction method of the uncertain set is required.

Uncertain set.

Full size image

2) Ellipsoid set.

The specific expression is:

where \({\mathbf{c}}\)denotes the centroid of the high-dimensional ellipsoid, and\(\Sigma \in {{\mathbf{R}}^{{N_{{\text{PV}}}} \times {N_{{\text{PV}}}}}}\) is a positive definite matrix that indicates the direction of the axis shift of the high-dimensional ellipsoid.

As shown in Fig. 1 (b), the ellipsoidal set, like the box set, covers all possible values of distributed PV output. At the same time, unlike the box set, the ellipsoidal set reduces the inclusion of empty regions where fluctuations are unlikely to occur, thereby reducing the conservativeness of the decision results. However, because the expression of the ellipsoid set is quadratic, it increases the difficulty of solving the problem in the optimization process.

Based on this¹⁹, proposed a generalized convex hull set that can significantly reduce the conservativeness of the optimization results while maintaining good robustness, and avoids introducing quadratic expressions in the modeling process. Therefore, this paper constructs a data driven uncertain set based on the method in¹⁹, with the modeling process shown in Fig. 2.

Modeling process of convex hull polyhedron uncertainty set.

Full size image

Step (1): First, construct a high-dimensional ellipsoid \({{\mathbf{U}}_{e1}}\) that covers all historical data fluctuations and has the smallest volume. The constructed high-dimensional ellipsoid is shown in Fig. 2 (a), and its specific expression is as follows (1):

Step (2): On the basis of the existing high-dimensional ellipsoid, perform eigenvalue decomposition on the positive definite matrix \({\mathbf{\Sigma }}\):\({\mathbf{\Sigma }}={{\mathbf{P}}^{\mathbf{T}}}{\mathbf{JP}}={{\mathbf{P}}^{ - 1}}{\mathbf{JP}}\). This rotates and translates the ellipsoid so that its center coincides with the origin of the coordinate system, as shown by the green dashed line in Fig. 2 (b). At this point, the uncertain set of the high-dimensional ellipsoid is denoted as \({{\mathbf{U}}_{e2}}\), i.e.,

where \({\mathbf{J}}\) is a diagonal matrix, denoted as \({\mathbf{J}}={\text{diag}}({\lambda _1} \ldots {\lambda _{{N_{{\text{PV}}}}}})\), and P is a transformation matrix representing the rotation angles of the matrix. According to the diagonal matrix \({\mathbf{J}}\), the coordinates of the vertex of the transformed high-dimensional ellipsoid are:

Next, connect each vertex of the high-dimensional ellipsoid to form a high-dimensional polyhedron, as shown by the red dashed line in Fig. 2 (b). At this point, the uncertain set of the high-dimensional polyhedron is denoted as \({{\mathbf{U}}_{p2}}\), which is:

where \({m_i}\) denotes the weight coefficient of the first vertex.

Step (3): Since the high-dimensional polyhedral set obtained in step 2 causes some data points to lie outside the boundary, it is necessary to perform a contraction process on the set. As shown by the red solid line in Fig. 2 (c), the vertices of the contracted high-dimensional polyhedral set are:

At this time, the high-dimensional linear polyhedron after scaling is shown in Fig. 2 (d), and its uncertainty set \({{\mathbf{U}}_{p2}}\) is

where k is the contraction multiple used to adjust the retention degree of the high-dimensional polyhedral set. The calculation method for k can be found in¹⁹. Therefore, there exists a minimum \({k_{\hbox{min} }}\) that allows the contracted polyhedral set to encompass all possible data points. The value range of k is \([0,{k_{\hbox{min} }}]\).The polyhedral sets formed by different k values are shown in Fig. 3.

Range of convex hull polyhedron set under different k.

Full size image

Step (4): Rotate and translate the contracted high-dimensional polyhedral set to encompass the original data points. From (5), it is known that after the rotation and translation, the uncertain set of the high-dimensional polyhedral set is \({{\mathbf{U}}_{p1}}\), which is :

1. (2)

   Hyperplane polyhedral set.

The polyhedral set described in the previous section is constructed by first establishing a high-dimensional polyhedral set based on the vertices of an ellipsoid, and then performing successive contractions to obtain an uncertain high-dimensional polyhedral set. This set encompasses all historical data points. However, if the contraction is uniformly applied, it might lead to some new data points being located outside the boundary due to excessive enlargement of the bounding region. Consequently, the corresponding bounding regions will include more empty space. Therefore^21,22, proposed constructing the uncertain set incrementally based on the positions of the boundary points, forming an irregular polyhedral set through iterative extreme value envelopes, as shown in Fig. 4.

Uncertainty set based on extreme scenarios.

Full size image

The uncertain set based on extreme scenarios is expressed as:

where \({N_{ex}}\) represents the number of extreme scenarios. \({\mathbf{z}}_{i}^{{{\text{PV}}}}\) represents the ith extreme scenario. \({\sigma _i}\) represents the weight coefficient of the extreme scenario.

Compared to the uncertain set of the convex hull, the uncertain set based on extreme scenarios can significantly reduce the inclusion of empty regions with low probability distributions. Therefore, this method maintains the best fidelity. However, as shown in (11), the construction of the uncertain set based on extreme scenarios depends on the number of extreme scenarios, i.e., the number of vertex points. If the number of extreme scenarios is large, it will increase the computational complexity of the optimization process. Therefore, this paper proposes a hyperplane polyhedral set, which cuts and flattens the boundary from the vertices based on the original box set range, significantly reducing empty regions, as shown in Fig. 5.

Hyperplane polyhedron set.

Full size image

In general, for a \({\rm K}\)-dimensional space, the hyperplane is generally expressed as:

where \({{\mathbf{\alpha }}_m}\) is a non-zero \({\rm K}\)-dimensional nonzero vector. \({\beta _m}\) is a scalar. m is the vertex index of the \({\rm K}\)-dimensional box set. Let \({\mathbf{z}}_{m}^{{\text{B}}}\) be the vertex of the \({\rm K}\)-dimensional box set. The vertex produced by this hyperplane cut is \({\mathbf{z}}_{n}^{{\text{H}}}\).Then the relationship between the vertex indices \(m,n\) is:

At this point, the vertices produced by the hyperplane cut can be obtained by solving the following model:

In (14), \({\xi _{ms}}\) represents the Euclidean distance between \({\mathbf{z}}_{n}^{{\text{H}}}\) and the corresponding \({\mathbf{z}}_{m}^{{\text{B}}}\). Since (14) is a nonlinear equation, to reduce computational pressure without affecting the decision variables to be solved, it is transformed into the form of (14-a):

(15) divides the \({\rm K}\)-dimensional space into inside and outside regions, representing the vector direction of the \({\mathbf{z}}_{m}^{{\text{B}}}\) of the box set, with the data vector \({\mathbf{\tau }}\) inside, and \({\rm T}\) representing the data set. (16) indicates that, except for the first \(o=s\)-dimension, other dimensions have coordinates with a standard deviation of 0. (17) shows the calculation formula for the Euclidean distance \({\xi _{ms}}\), where \({\theta _{ms}}\) takes values of 1 or −1, depending on whether \({\mathbf{z}}_{{ms}}^{{\text{B}}}\) in the s-dimension is greater or less than \({\mathbf{z}}_{{ns}}^{{\text{H}}}\). (18) indicates that any hyperplane does not intersect with the \({\rm K}\)-dimensional box set. h and \(H(s)\) represent the order and collection of intervals of the s-dimensional uncertain set of the box, while \({m_1}\) and \({m_2}\) represent the boundaries of the regions.

The above model is a nonlinear model, so the interior point method is used for solving. After determining the coordinates of the vertices of the hyperplane and combining with (11), the uncertain set of the hyperplane is expressed as:

where \({\varepsilon _i}\) denotes the weighting coefficient of the ith vertex.

1. 1)

   Objective function.

This paper constructs a distribution network power supply recovery model that maximizes load recovery and minimizes network loss. The model is formulated as:

where \({\lambda _{i,t}}\) represents the power supply recovery status of node i in time period t, where 0 indicates no recovery and 1 indicates full recovery. \({r_i}\) represents the importance weight of node i. This paper divides the load into three levels²³, with corresponding values of 10, 5, and 1 for each level. \(P_{{i,t}}^{{{\text{INF}}}}\) denotes the power of load recovery of node i in time period t. \({P_{ij,t}}\) and \({P_{ji,t}}\) denote the active power of branch \(ij\) in time period t flow. \({P_{ji,t}}\) denotes the power flow from node i to node j at time t. \({P_{ji,t}}\) denotes the power flow through node i to node j. C represents the total power supply recovery cost of the distribution network. \({c_{{\text{loss}}}}\) denotes the power loss coefficient of the distribution network. D represents the set of distribution network nodes. L denotes the set of all paths in the distribution network.

1. 2)

   Constraint condition.

1) Second-order cone relaxation power flow constraints.

The second-order cone relaxation power flow model under the polar coordinate system can be expressed as:

where \(P_{{i,t}}^{{{\text{sum}}}}\) and \(Q_{{i,t}}^{{{\text{sum}}}}\) denote the active and reactive power injected by the node i at the t time, respectively. \(Q_{{i,t}}^{{\text{D}}}\) denotes the reactive power load of node i at time t. \({Q_{ij,t}}\) denotes the reactive power flowing from node i to node j at time t of branch \(ij\) in the distribution network. \({P_{i,t}}\) and \({Q_{i,t}}\) denote the gateway power injected by node i at time t, respectively. \(P_{{i,t}}^{{{\text{SOP}}}}\) and \(Q_{{i,t}}^{{{\text{SOP}}}}\) represent the active and reactive power of the SOP at node i at time t. \(N(i)\) denotes the set of all nodes connected to node i. \({g_l}\) and \({b_l}\) denote the conductance of branch l and the conductance to ground. \({I_l}\) indicates the current amplitude of branch l. Also introduce \({u_{i,t}}={{V_{{i,t}}^{2}} \mathord{\left/ {\vphantom {{V_{{i,t}}^{2}} {\sqrt 2 }}} \right. \kern-0pt} {\sqrt 2 }}\), \({R_{l,t}}={V_{i,t}}{V_{j,t}}{\text{cos}}{\theta _{ij,t}}\), \({T_{l,t}}={V_{i,t}}{V_{j,t}}{\text{sin}}{\theta _{ij,t}}\), where \({V_{i,t}}\) and \({V_{j,t}}\) denote the voltage amplitude of node i and j at time t, respectively. \({\theta _{ij,t}}\) represents the voltage phase angle difference between both ends of the line \(ij\) at time t. Additionally, since \({u_{i,t}}\) does not include the branch switch variable \({\alpha _{l,t}}\), a new variable \(u_{{i,t}}^{l}\) that includes the branch switch variable is introduced to jointly determine the branch status and the voltage at the branch end nodes.

When the branch l is open at time t both \(u_{{i,t}}^{l}\) and \(u_{{j,t}}^{l}\) are 0. When the branch l is closed at time t, \(u_{{i,t}}^{l}\) and \(u_{{j,t}}^{l}\) are equal to the corresponding node voltages \({u_{i,t}}\) and \({u_{j,t}}\). To express the relationship between \(u_{{i,t}}^{l}\), \(u_{{j,t}}^{l}\), and \({u_{i,t}}\), \({u_{j,t}}\), additional constraints are added, i.e.,

where \(V_{i}^{{{\text{max}}}}\) and \(V_{j}^{{{\text{max}}}}\) represent the maximum value of voltage amplitude at nodes i and j, respectively. Additionally, the constraints for node voltages and branch currents are respectively given by:

where \(I_{l}^{{\hbox{max} }}\) represents the maximum value of current amplitude in branch l.

2) Distributed PV constraints

(30)-(32) represent the operational constraints of distributed PV, where \(P_{{i,t}}^{{{\text{PV,f}}}}\) represents the maximum available power at node x connected to distributed PV before the fluctuation at time t. \(\Delta {P^{{\text{PV,max}}}}\) denotes the maximum fluctuation of the distributed PV. \(Q_{{i,t}}^{{{\text{PV}}}}\) and \(S_{{i,t}}^{{{\text{PV}}}}\) denote the reactive power and capacity of the connected distributed PV at node i at time t.

3) Battery energy storage constraints

(33)-(35) represent the operational constraints of the battery energy storage. Where \(S_{{i,t}}^{{{\text{SOC}}}}\) represents the charging state of the connected battery energy storage at node i at time t. \(D_{{i,t}}^{{{\text{ch}}}}\) and \(D_{{i,t}}^{{{\text{dis}}}}\) are 0–1 variables, which represent the charging and discharging state of the connected battery energy storage at node i at time t, with 1 indicating charging and 0 indicating discharging. \({\eta _{i,{\text{ch}}}}\) and \({\eta _{i,{\text{dis}}}}\) represent the charging and discharging efficiencies of the connected battery energy storage at node i, respectively. \({E_{i,{\text{max}}}}\) represents the maximum amount of battery storage energy at node i. \(S_{{i,\hbox{min} }}^{{\text{S}\text{O}\text{C}}}\) and \(S_{{i,\hbox{max} }}^{{\text{S}\text{O}\text{C}}}\) represent the minimum and maximum states of charge of the battery at node i, respectively. \(\Delta t\) represents the time interval for the battery state adjustment. H indicates the set of charging and discharging time of the battery energy storage.

4) Load model constraints.

This paper classifies loads into two types: controllable loads and uncontrollable loads. For controllable loads, the distribution network power supply recovery can restore part of the load. For uncontrollable loads, the distribution network power supply recovery must fully restore or not restore them at all. The specific model is as follows:

where \({k_{i,t}}\) denotes the load recovery factor of node i at time t, when the load is controllable, \({k_{i,t}} \in [0,1]\). \({k_{i,t}}\) is equal to 1 when the load is equivalent to an uncontrollable load. \(P_{{i,t}}^{{\text{D}}}\) denotes the amount of load of node i at time t.

5) Radial topology structure.

This paper uses the ST constraints to describe the radial topology constraint of the distribution network, as shown in (37).

(37) determines whether the first and last nodes i and j of a branch are child nodes or parent nodes by introducing two 0–1 variables \({\beta _{ij,t}}\) and \({\beta _{ji,t}}\), ensuring that the generated radial structure connects to the root nodes to form a tree. \({\beta _{ij,t}}=1\) indicates that node i is the parent node of node j, and \({\beta _{ij,t}}=0\) indicates that it is not. N represents the set of all nodes in the distribution network, and \({N_s}\) represents the set of root nodes.

6) SOP constraints

where \(P_{{i,t}}^{{{\text{SOP,L}}}}\) and \(P_{{j,t}}^{{{\text{SOP,L}}}}\) represent the losses of the SOP. \({\eta ^{{\text{SOP}}}}\) represents the loss coefficient of the SOP. \(S_{i}^{{{\text{SOP}}}}\) represents the capacity of the SOP.

1. 3)

   Information-physical integration framework.

In the physical-information integrated framework constructed in this paper, the physical network is the distribution network, and the information network refers to the network formed by the terminal devices that detect and control the distribution network. There is a one-to-one correspondence between the nodes of the information network and the nodes of the distribution network, meaning each load node corresponds to an information node²⁴.

1) Energy supply constraints for information nodes.

During the power supply recovery of the distribution network, the information node will only be recovered if the recovered load at the load node is greater than the energy demand of the information node. The specific model is as follows:

where \({\varphi _{i^{\prime},t}}\) is a 0–1 variable indicating the effective situation of energy supply of the information node \(i^{\prime}\), equal to 0 indicates that the energy supply is invalid, and equal to 1 indicates that the energy supply is effective. x indicates the coupling parameter of the information node and the load node, \(x \in [0,1]\). The larger x is, the larger the coupling strength is. \({N_{\text{D}}}\) denotes the set of load nodes connected to the information node.

(43) is the information node validity constraint, which indicates that when the load recovery of a load node is greater than the supply demand of the information node, the node is flagged as a valid state and vice versa. (44) is the definition of the load recovery of the load node.

Since (43) is non-linear, it needs to be linearized as shown in the following equation:

where \(\sigma\) is a very small positive number.

2) Information network traffic constraints.

Whether an information node works or not depends not only on the energy supply status of that node, but also on whether the node has communication traffic with the control node, as modelled below:

where \({\psi _{i^{\prime},t}}\) indicates the working state of information node \(i^{\prime}\) in time period t. \({T_{i^{\prime}j^{\prime},t}}\) and \({T_{j^{\prime}i^{\prime},t}}\) indicate that the communication traffic flows from node \(i^{\prime}\) to node \(j^{\prime}\) and from node \(j^{\prime}\) to node \(i^{\prime}\) in time period t, respectively. \({n_{{\text{C,goal}}}}\) denotes the number of terminals of the information node. \({N_{\text{C}}}\) denotes the set of information nodes. \({N_{{\text{C,0}}}}\) denotes the set of information source nodes. \({B_{i^{\prime}}}\) denotes the bandwidth of the node \(i^{\prime}\).

3) Information node control constraints.

If the information node fails, its corresponding load node will be uncontrollable and its corresponding model is as follows.

(51) indicates that if the information node fails, its corresponding load node will not be recoverable. (52) indicates that the external power supply connected to the load node can be put into operation only when the information node is valid. \(P_{{i,\hbox{max} }}^{{{\text{sum}}}}\) denotes the maximum value of the output power of the external power supply.

1. 1)

   Robust power supply recovery optimization model for distribution network.

Let the power flow constraints variable be vector \({\mathbf{P}}=\left\{ {{P_{ij,t}},{P_{ji,t}},{Q_{ij,t}},{Q_{ji,t}},{u_{i,t}},{R_{l,t}},{T_{l,t}},{I_{l,t}},} \right\}\). Let the PV constraints variable be vector \({{\mathbf{P}}^{{\text{PV}}}}=\left\{ {P_{{i,t}}^{{{\text{PV}}}},Q_{{i,t}}^{{{\text{PV}}}}} \right\}\). Let the battery energy storage operating variable be vector \({{\mathbf{P}}^{{\text{ESS}}}}=\left\{ {S_{{i,t}}^{{{\text{SOC}}}},P_{{i,t}}^{{{\text{ch}}}},P_{{i,t}}^{{{\text{dis}}}}} \right\}\). Let the load model constraints be vector \({{\mathbf{P}}^{\text{D}}}=\left\{ {P_{{i,t}}^{{{\text{INF}}}}} \right\}\). Let the radial topology constraints be vector \({{\mathbf{P}}^{{\text{switch}}}}=\left\{ {{\beta _{ij,t}},{\beta _{ji,t}},{\alpha _{l,t}}} \right\}\). Let the information network constraints be vector \({{\mathbf{P}}^{{\text{INF}}}}=\left\{ {{\varphi _{i,t}},{\lambda _{i,t}},{T_{ij,t}},{\psi _{i,t}}} \right\}\). Let the SOP constraints be \({{\mathbf{P}}^{{\text{SOP}}}}=\left\{ {P_{{i,t}}^{{{\text{SOP}}}},P_{{i,t}}^{{{\text{SOP,L}}}},Q_{{i,t}}^{{{\text{SOP}}}}} \right\}\).

On the basis of establishing the data driven polyhedral set for distributed PV output, this paper establishes a two-stage robust power supply recovery optimization model for distribution networks, whose matrix form is:

where \({\mathbf{x}},{\mathbf{y}}\) are decision variables of the model and \({\mathbf{u}}\) is an uncertainty variable. \({\mathbf{x}}=\left\{ {D_{{i,t}}^{{{\text{ch}}}},D_{{i,t}}^{{{\text{dis}}}},{\varphi _{i,t}},{\psi _{i,t}},{{\mathbf{P}}^{{\text{switch}}}}} \right\}\) is the first stage decision variable. \({\mathbf{y}}=\left\{ {{\mathbf{P}},{{\mathbf{P}}^{{\text{PV}}}},{{\mathbf{P}}^{{\text{ESS}}}},{{\mathbf{P}}^{\text{D}}},{{\mathbf{P}}^{{\text{INF}}}},{{\mathbf{P}}^{{\text{SOP}}}}} \right\}\) is the second stage decision variable. \({\mathbf{u}}=\left\{ {\tilde {P}_{{i,t}}^{{{\text{PV}}}}} \right\}\) is the second stage uncertainty variable. Constant matrix \({\mathbf{A}}\) denotes the coefficient matrix associated with decision variable \({\mathbf{x}}\), and column vector \({\mathbf{d}}\) is a constant denoting the coefficient vector associated with decision variable \({\mathbf{x}}\). Constant matrices \({\mathbf{G}}\),\({\mathbf{E}}\),\({\mathbf{R}}\) denote the coefficient matrix associated with decision variable \({\mathbf{y}}\). Column vectors \({\mathbf{h}}\) and \({\mathbf{r}}\) are constant vectors denoting the coefficient vector associated with decision variable \({\mathbf{y}}\). The constant matrix \({\mathbf{M}}\) represents the matrix of coefficients associated with the uncertain variable \({\mathbf{u}}\). \(\Omega \left( {{\mathbf{y}},{\mathbf{u}}} \right)\) is the feasible ___domain of continuous variable \({\mathbf{y}}\) given \(\left( {{\mathbf{x}},{\mathbf{u}}} \right)\). \({{\mathbf{c}}^{\text{T}}}{\mathbf{y}}\) denotes the objective function of the second stage, corresponding to (16). (53-a) corresponds to the constraints associated with variable \({\mathbf{x}}\) in the first stage. (53-b) corresponds to the constraints associated with variable \({\mathbf{y}}\) in the second stage. (53-c) corresponds to the second-order cone constraints associated with variable \({\mathbf{y}}\) in the second stage.

For the two-stage robust optimization model in (53), it includes continuous variables and integer variables. Since the second stage contains uncertain variables \({\mathbf{u}}\), it cannot be solved directly. Therefore, this paper adopts the \({\text{C\& CG}}\) method^25,26 to transform it into a master-slave problem for solution. The master problem seeks the worst-case scenario that maximizes the load recovery amount. The problem is then decomposed into solving the integer variables of the master problem first (including the charging and discharging state of the storage device, switch status, and the number of information nodes, etc.), optimizing the objective function for the second-stage variables, and iteratively updating the feasible region to obtain the total load recovery amount under the worst-case scenario.

1. 2)

   C&CG solution methods.

2. 1)

   Main and subproblem model.

The main and subproblem model corresponding to (53) are:

First, solve the main problem \({\text{MP1}}\) corresponding to (54). \({\text{MP1}}\) is a mixed-integer second-stage robust optimization problem. The first stage variable solution \({{\mathbf{x}}^ * }\) corresponding to \({\text{MP1}}\) and the auxiliary variable \(\eta\) introduced in the \(k+1\) iteration are \({\text{C\& CG}}\) cuts. Then, the first-stage decision variable \({{\mathbf{x}}^ * }\) is brought into the second-stage subproblem \({\text{SP1}}\) to seek the worst-case scenario \({{\mathbf{u}}^l}\), where l is the historical iteration number and k is the current iteration number. Finally, the worst-case scenario \({{\mathbf{u}}^l}\) of the second-stage subproblem \({\text{SP1}}\) is brought back into the first-stage main problem \({\text{MP1}}\) for solving. The last three constraints in (54) are the set of extreme cut and non-intersecting cut formed by the worst-case scenario variables \({\mathbf{\pi }}\), \({{\mathbf{\tau }}^a}\), \({{\mathbf{\tau }}^b}\).

2) Subproblem solving methods.

(55) is a \(\hbox{max} - \hbox{min}\) optimization problem, we use the duality theorem to convert the inner \(\hbox{min}\) problem of (55) into its dual form to merge it into a maximization problem, which is shown in the form of (56):

where there exists a bilinear term \({\left( {{\mathbf{Mu}}} \right)^{\text{T}}}{\mathbf{\pi }}\), which is solved using the outer approximation method for bilinear terms²⁶. The main problem \({\text{MP2}}\) and subproblem \({\text{SP2}}\) are obtained, as shown in (57) and (58):

where \({\text{MP2}}\) and \({\text{SP2}}\) are used to solve the upper and lower bounds of (56), respectively. m is the number of historical iterations, and n is the number of current iterations. An auxiliary variable \({\mathbf{\beta }}\) is introduced to replace the bilinear term in the original equation, and a bilinear term exists in (58) \({G^m}\left( {{\mathbf{u}},{\mathbf{\pi }}} \right)={\left( {{\mathbf{Mu}}} \right)^{\text{T}}}{\mathbf{\pi }}\). Therefore, it is necessary to use the outer approximation method for linearization, and the linearization formula is shown in (59):

1. 3)

   Solution steps and processes.

The specific solving steps of the \({\text{C\& CG}}\) algorithm are as follows:

1) Set the initial values of the upper and lower bounds of the master-slave problem to be \(U1=+\infty\) and \(L1= - \infty\), respectively, and the initial number of iterations to be \(k=1\), given that the convergence value is \(\varepsilon {1^{\hbox{max} }}\).

2) Solve the worst-case main problem, where the constraints of the main problem do not contain \({\text{C\& CG}}\) cut, and obtain the integer solution \({{\mathbf{x}}^ * }\)

3) Based on the integer solution \({{\mathbf{x}}^ * }\) obtained from solving the main problem, solve the subproblem to obtain the worst case scenario \({({{\mathbf{u}}^{l+1}})^ * }\), and the objective function value of the subproblem, and update the upper limit \(U1=\hbox{min} \left( {U1,{{\mathbf{c}}^{\text{T}}}{\mathbf{y}}} \right)\).

4) Determine whether \(U1 - L1\) is smaller than the convergence value \(\varepsilon {1^{\hbox{max} }}\), if it is valid, then end. If not, make \(k=k+1\), and add a new set of scenario variables \({{\mathbf{u}}^k}\) and \({\text{C\& CG}}\) cuts for the main problem. Solve the main problem to get \({\eta ^ * }\) to update the lower bound \(L1=\hbox{max} \left( {L1,{\eta ^ * }} \right)\) and go to step 3.

The specific solution steps are shown in Fig. 6.

Solution flow chart of C&CG.

Full size image

To verify the feasibility of the proposed adaptive data driven multi period robust power supply recovery method for distribution networks, this section uses the IEEE-33 bus system and a distribution network in the northeastern city of Hubei Province for case analysis. For convenience, it is assumed that the availability factors of the two PV systems before the fault are equal. According to the calculation method provided in¹⁹, the value of \({k_{\hbox{min} }}\) is 1.41.

1. 1)

   33-bus system example analysis.

The connection diagram of the modified IEEE-33 bus system is shown in Fig. 7. In Fig. 7, with the substation disconnected, the distribution network is in a state of total power outage. The system load is powered by diesel generators, PV systems and BESS. Table 1 provides the parameter settings for the integrated system’s PV, BESS, and SOP. The system’s base voltage is 12.66 kV, and the base capacity is 10MVA. The active power range for the connection is 0–5000 kW, the reactive power range is 0-5000kvar, the upper limit for branch current magnitude is 0.5 p.u., and the node voltage magnitude range is 0.95–1.05 p.u.

Modified IEEE-33 bus test system.

Full size image

1. 1)

   Power supply recovery results when using hyperplane set.

Figs.8,9,10 illustrate the load recovery situation and the energy supply of the system during the power outage phase when using the hyperplane uncertain set. According to the diagrams, the energy supply of the 33-bus case study system during the power outage phase relies mainly on distributed PV systems, BESS, and diesel generators. When there is no PV output (1–7 h), the energy supply mainly relies on diesel generators and BESS. During the PV output period (13–14 h), the large-scale access to PV energy leads to a significant reduction in the load recovery ratio, and the system’s load is more adequately supplied. For areas with adequate PV resources, the large-scale PV output can significantly reduce the load recovery ratio. When the PV output is insufficient, the system compensates for the insufficient energy supply with BESS, which meets the load recovery needs, while diesel generators operate to meet the system’s basic operating requirements. During the grid fault phase, the intelligent switch S2 at node 4 is triggered and opened, causing node 3 to lose power. When the PV output increases, the system begins to supply power from node 4 to node 9, reversing the power flow from node 4 to node 9. This change is due to the large-scale PV energy injection.

Power supply restoration situation.

Full size image

Energy output situation.

Full size image

SOP transmission status

Full size image

1. 2)

   Effect of scaling factor k on power supply recovery

Table 2 shows the load recovery situation of the 33-bus system under different sets. According to the table data, when the convex hull set is used, as the number of contractions increases, the load recovery ratio of the system gradually decreases, the number of recovered loads also decreases, and network losses are also reduced. However, the rate of reduction slows down as the number of contractions increases. This is because the convex hull set expands the historical maximum data boundary, making the worst-case scenario more likely to occur. To achieve a balance between reducing network loss and the uncertainty introduced by PV injection, the system needs to filter out part of the power injected by distributed PV, reducing the load recovery ratio and the number of recovered loads. For network loss, the power loss and load loss both decrease due to the reduction in new energy injection, resulting in a certain degree of stability in the system. Compared with the convex hull set and the box set, the recovered weighted load of the hyperplane set can be increased by 0.55% and 0.56% respectively. This verifies that the hyperplane set can improve the robustness of the optimization results when considering all the power outputs of photovoltaics.

When using the hyperplane set, it applies the mathematical model of the hyperplane to the historical maximum data boundary, significantly reducing the coverage range. This method can effectively balance the system’s robustness and reduce system loss. Therefore, using the hyperplane set, the load recovery ratio and the number of recovered loads are the smallest. Compared with the convex hull set, it can better respond to the worst-case scenario of PV output. Due to the balance of robustness and efficiency, the load recovery ratio and the number of recovered loads of the system are higher.

1. 3)

   Impact of three uncertainty sets on power supply recovery.

The load recovery situation at each time period is shown in Figs. 11 and 12. From the figures, it can be seen that the 33-bus system load recovery reaches its peak at 13:00, which is consistent with the PV output situation. When using the three types of uncertain sets for system fault recovery, it is observed that the hyperplane set results in a higher load recovery at each time period compared to the convex hull set and the box set. This is consistent with the previous analysis.

Load recovery percentage under three uncertain sets.

Full size image

:Restoring weighted load cases using three uncertain sets. (From left to right are the convex hull set, the hyperplane set, and the box set).

Full size image

1. 4)

   Effect of coupling factor \(\gamma\) on power supply recovery results

The Figs. 13, 14 and 15 give the impact of taking each of the three uncertainty sets on power supply recovery as the coupling degree changes from 0.2 to 1. From the Fig. 13, it can be seen that the weighted load of 33-bus system recovery decreases as the coupling degree increases. This is due to the fact that when the coupling degree increases, the information node has a higher recovery demand on the power node, resulting in the power node not being able to meet the energy supply demand of the information node’s end devices during power supply recovery, which leads to a decrease in the amount of load recovered. In the Fig. 14, when the coupling degree increases, in order to increase the energy supply to the information nodes as much as possible, the power nodes increase the transmission flow of power to each node, which leads to an increase in the network loss of the system. In the Fig. 15, when the coupling coefficient is small, the information nodes are less dependent on the power nodes, and the information network is more resistant to interference, so more information nodes can be recovered. However, when the coupling coefficient increases, since the premise of information node restoration is that the power node load must be completely restored, and it can be seen from the attached Fig. 13 that at this time the amount of load restored by the power system is less, which leads to a reduction in the number of information nodes restored as well, and makes the information network relatively fragile.

Restoration of weighted load under different. (From left to right are the convex hull set, the hyperplane set, and the box set).

Full size image

Network loss situation under different.

Full size image

Information node recovery under different.

Full size image

1. 2)

   Analysis of the 26-bus example in an urban area in Hubei.

The connection diagram of the distribution network system in a city in Hubei Province is shown in Fig. 16. Similar to the 33-bus system, in Fig. 16, the substation is disconnected, and the distribution network is in a state of total power outage. The system load is supplied by diesel generators, distributed PV systems, and BESS. Table 3 provides the parameter settings for the integrated system’s PV, BESS, and SOP. The system’s base voltage is 12.66 kV, and the base capacity is 10 MVA. The active power range for the connection is 0–5000 kW, the reactive power range is 0–5000 kvar, the upper limit for branch current magnitude is 0.5 p.u., and the node voltage magnitude range is 0.95–1.05 p.u.

Testing system for a certain urban area in Hubei.

Full size image

1. 1)

   Fault recovery results when using hyperplane aggregation.

The load recovery situation of the 26-bus power supply system in a city in Hubei is shown in Figs. 17, 18 and 19. It is worth mentioning that due to the almost negligible total load of the third-level loads in this city’s system, the third-level loads are not recovered during fault recovery. The recovery results are similar to the 33-bus system, with the load recovery reaching its peak during the PV output period. During other periods, diesel generators operate at full power to compensate for the system’s load deficiency.

Restoration of power supply in a certain urban area system.

Full size image

Energy output of a certain urban area system.

Full size image

SOP transmission status.

Full size image

1. 2)

   Effect of scaling factor k on optimization results

The effect of deflation multiplier k on the fault recovery in this urban area is shown in Table 4, from which it can be seen that when convex packet aggregation is used, the deflation multiplier decreases as it gradually increases, both in terms of the total amount of weighted loads recovered, the number of loads and the network loss, which is consistent with the results of the 33-node system fault recovery, while the use of hyperplane aggregation is still one of the uncertainty aggregations with the highest amount of weighted loads recovered. Compared with the convex hull set and the box set, the recovered weighted load of the hyperplane set can be increased by 0.33% and 0.35% respectively. This verifies that the hyperplane set can improve the robustness of the optimization results when considering all the power outputs of photovoltaics.

1. 3)

   Impact of three uncertainty sets on power supply recovery.

The results of the three uncertainty sets for the recovery of a 24-node fault in an urban area are shown in the Figs. 20 and 21. The only difference with the 33-node system is that the 24-node system, due to its particular load distribution, therefore, the total amount of recovered tertiary loads is almost zero, while the rest of the variations are similar to that of the 33-node system.

Load recovery percentage under three uncertain sets.

Full size image

Restoring weighted load cases using three uncertain sets. (From left to right are the convex hull set, the hyperplane set, and the box set).

Full size image

1. 4)

   Effect of coupling factor \(\gamma\) on optimization results

The power supply recovery results of coupling coefficient \(\gamma\) for 24 nodes in an urban area are shown in Figs. 22, 23 and 24, and it can be seen from Fig. 22 that when the coupling coefficient is changed from 0.2 to 0.8, the change of the recovered weighted load is not so obvious compared with that of the 33-node system, which is because of the difference of the load distribution of the actual example and the load distribution of the standard example, and the load recovery is not so obvious due to the fact that the load distribution of many electric power nodes is zero for the actual example. node loads are 0, which leads to the fact that when the coupling coefficient changes, the restored node loads may not change, so the change in load restoration is not so obvious, while when \(\gamma =1\) is in place, the information nodes have a high demand for restoration of the power nodes due to the complete coupling of information and power networks, so the impact on the load restoration of the power system is also the largest. The variation of network loss and node recovery under different coupling coefficients is similar to that of 33-node.

Restoration of weighted load under different. (From left to right are the convex hull set, the hyperplane set, and the box set).

Full size image

Network loss situation under different.

Full size image

Information node recovery under different.

Full size image

This paper constructs a distribution network fault recovery model based on adaptive data driven methods and employs the \({\text{C\& CG}}\) method for solving. Finally, through the case studies of the IEEE-33 bus system and a 24-bus system in a city in Hubei Province, the power supply recovery models under three types of uncertain sets are calculated, and the results are as follows:

(1) Under the same coupling coefficient, in the IEEE-33 bus system, the recovered weighted load of the hyperplane set is increased by 0.55% and 0.56% respectively compared with that of the convex hull set and the box set. In the actual 26-bus example system, the recovered weighted load of the hyperplane set is increased by 0.33% and 0.35% respectively compared with that of the convex hull set and the box set. This verifies that the hyperplane set can improve the robustness of the optimization results when considering all the power outputs of photovoltaics. Furthermore, with the increase in the coupling factor, the load recovery ratio and the number of recovered loads of the system using the convex hull set gradually decrease, increasing the system’s robustness.

(2) As the system load increases, the dependence of the information network and the power network on each other continues to increase. Therefore, the power supply recovery ratio of the system shows a downward trend. However, the hyperplane set performs better, especially in accurately depicting the parameter range of the uncertain set, reducing the empty regions with low probability and reducing the system loss. Hence, the hyperplane set based on adaptive data driven methods provides a more robust and reliable fault recovery method than the convex hull set, showing stronger robustness.

The data driven multi period robust fault recovery method for distribution networks mentioned in this paper is aimed at the hourly fault recovery scenarios of an actual distribution network in a certain place. In time-sensitive real-time fault recovery scenarios, due to the iterative characteristics of the C&CG algorithm, there may be the problem of relatively large computational complexity. Future research can take advantage of the strengths of modern multi-core processors and distributed computing architectures, conduct in-depth research on the parallelizable parts of the C&CG algorithm. By designing a reasonable parallel algorithm framework, multiple computing nodes can simultaneously handle different subtasks, thereby significantly reducing the overall computing time. At the same time, in order to demonstrate the applicability of this method to large-scale networks in the real world, in our subsequent research, we will actively cooperate with relevant power enterprises to obtain representative large-scale actual distribution network data, providing a reliable basis for testing our method.