Bus-Privacy-Preserving Distributed Economic Dispatch Method for Power Systems Using Iteration Functions

The sufficiency of electricity market competition and the economy of power system operation require that the economic dispatch should not only preserve the bus privacy but consider network losses. The existing economic dispatch methods, however, can’t preserve the privacy of each bus from leakage and even complicate the solution when they reasonably consider network losses, due to the coupling complexity of network losses. To address these issues, this paper presents a novel distributed economic dispatch (DED) method based on KKT optimality conditions and iteration functions. It focuses on an economic dispatch model that contains bus active power balance equations and generator power limit inequalities. First, the model’s KKT optimality conditions with inequalities are simplified to equivalent pure equality-type optimality conditions. Then a set of iteration functions is created using the Taylor series, which overcomes the difficulty caused by the fact that unknown bus voltage phases are completely contained in the sine and cosine functions of optimality conditions. Finally, a fast fixed-point iteration procedure is proposed for bus agents to solve the economic dispatch model. The proposed method reasonably considers network losses by using the set of bus active power balance equations. It not only preserves the privacy of each bus from leakage but is straightforward in solution (does not complicate the solution). In addition, it is fully bus-level distributed. Numerical simulation results of systems with different sizes verified the effectiveness and advantages of the proposed method.

Subscript letting its related variable take the (k + 1)-th step estimated solution if available.

I. INTRODUCTION
Economic dispatch is one of the most important routines in the power system operation. It aims at optimizing the base power of individual generators by minimizing the total generation cost under some system constraints, to establish the base point for load-frequency control. Conventional economic dispatch is implemented by a power system control center. The control center collects power and cost information (bus privacy) about the generators/loads on all buses of the system at first, then establishes and solves the optimization model for economic dispatch. Thus, the conventional economic dispatch methods are in a centralized manner. See for example the equal incremental cost (λ) method and piecewise linear programming method [1], genetic algorithm approach [2], neural network method [3], convex relaxation methods [4], [5], and so forth. The conventional economic dispatch methods need a control center to collect global information and execute centralized computation. Thus, they not only have the defect of disclosing bus privacy but require a high-performance computer and a reliable control-center-oriented communication network. The recently proposed distributed economic dispatch (DED) methods [6], [7], [8], [9] not only have the potential ability to preserve the bus privacy due to their distributed calculation but can lower the requirements for the computer performance and communication by distributing the computational burden among distributed computers and exchanging information between neighbors. These features are required not only by the sufficiency of electricity market competition but by the dispersity of renewable power sources. As a result, the DED methods are the inevitable trend of power system dispatch in the future.
The earliest distributed economic dispatch method is the incremental cost consensus (λ consensus) method [6], here called the basic λ consensus DED method. It updates each generator's λ using the average consensus dynamics at first, then determines each generator output power following its λ, and lets a specially deployed leader generator collect power information about all generators/loads and correct its λ by the system mismatch power neglecting network losses at last. This process is repeated and does not stop until the λ's of all buses' generators reach the same. The basic λ consensus DED method needs a leader generator for collecting global power information and requires the network-loss-neglected system mismatch power to vanish. Thus, it not only discloses bus privacy but neglects network losses. To overcome the defect of disclosing bus privacy, several improved λ consensus DED methods have been proposed. By creating a local (bus) estimated system mismatch power dynamics (instead of the leader generator's action), an improved λ consensus DED method is proposed in [10], which achieves each bus privacy preservation. By creating a generator power dynamics (to avoid calculating system mismatch power) and accordingly improving the λ dynamics, another improved λ consensus DED method is proposed in [11], which preserves bus privacy. By creating another local (bus) estimated system mismatch power dynamics and accordingly improving the λ dynamics, one more improved λ consensus DED method is proposed in [12], which protects bus privacy from leakage. All these improved λ consensus DED methods preserve bus privacy by removing the leader generator. Developed from the basic λ consensus DED method, they still require the network-loss-neglected system mismatch power to vanish. As a result, although these improved λ consensus DED methods can preserve the privacy of each bus from leakage, they fail to account for network losses.
The economy of power system operation requires that the economic dispatch should consider network losses. So several improved λ consensus DED methods have been proposed by introducing network losses. The simplest practice was taking a certain percent of the value of total loads as the value of total network losses [13], [14], it neglected the impact of locations and outputs of generators, so it is low in rationality. By assuming that each bus knows, from protective relay devices, the currents through the lines connected to it, the bus's local portion of line losses was separately calculated in [15] using these currents. This method improves the economy of dispatch indeed but is irrational (impractical) since the required line-current information is outside of the normal dispatch conditions. By letting the value of network losses induced by a generator be equal to a coefficient times the square of the generator's power, the generator's local portion of losses was separately calculated in [16] and [17]. This method improves the economy compared to the basic λ consensus DED method. However, the value of total network losses in a power system is, based on loss coefficients, a quadratic function of all generators' power [18]. That means it depends upon not only the square of each generator power but the cross-products of two arbitrarily different generators' power. So the value of total network losses in a power system is not locally separable for each bus/generator. Following this idea, loss coefficients were reasonably used in [19] and [20] to calculate the value of total network losses in a power system. Loss coefficients, however, can only be obtained by centralized calculation [21]. As a result, the improved λ consensus DED methods can't preserve the privacy of each bus from leakage when reasonably considering network losses.
The value of total network losses in a power system is equal to the sum of all buses' injection power. So the network losses are reasonably considered by replacing the total system power balance equation constraint (used in the basic λ consensus DED method) with the set of all buses' power balance equations in economic dispatch. In this case, the economic dispatch problem becomes an optimal power flow (OPF) problem. The existing distributed algorithms for OPF have been comprehensively analyzed in [22], [23], and [24]. There are mainly six distributed OPF algorithms applicable for economic dispatch: analytical target cascading (ATC), alternating direction method of multipliers (ADMM), proximal message passing (PMP), auxiliary problem principle (APP), optimality condition decomposition (OCD), and consensus + innovations (C+I). The former four are implemented by relaxing border constraints into the objective function in terms of augmented Lagrangian and then decomposing the global model according to areas. The reformulation of the global model, especially the nonlinearity introduced by the augmented Lagrange terms, complicates the solution. The OCD method is implemented by assigning each coupling constraint to a specific area, formulating each area subproblem model, and letting the objective function of each subproblem contain (by Lagrange term) the coupling constraints assigned to the other areas. The reformulation of the primal model, especially the lack of the theory for how to assign coupling constraints, complicates the solution. The C+I method is implemented via designing an iteration formula for each kind of unknown variable by combining, using artificial tuning parameters, the mismatches of KKT conditions. However, there is at least one artificial tuning parameter that heavily impacts the convergence of one kind of unknown variable. And the lack of a theory for how to select these tuning parameters complicates the solution. In addition, only the area-level distributed methods (excluding C+I) have been simulated and analyzed before, they have the defects of intra-area disclosure of bus privacy and degrading the sufficiency of electricity market competition. The existing OPF-based DED methods not only disclose bus privacy in an intra-area way but complicate the solution (reformulate the model and introduce the tuning parameters).
In summary, economic dispatch should not only preserve the bus privacy but consider network losses. The existing DED methods, however, can't preserve the privacy of each bus from leakage and even complicate the solution when reasonably considering network losses. To settle these problems, this paper presents a novel DED method based on KKT optimality conditions and iteration functions. It views each bus as an entity managed by an agent ''a bus agent'', and assumes that the exchange of information between neighbor buses is implemented by power line carrier communication. The main contributions of this paper are as follows. (1) It creates a novel DED method using iteration functions, the creation overcomes the difficulty caused by the fact that unknown bus voltage phases are completely contained in the sine and cosine functions of KKT conditions. (2) It settles the problems that the existing DED methods can't preserve the privacy of each bus from leakage and even complicate the solution when reasonably considering network losses. (3) It raises the privacypreserving level from the intra-area level of the existing DED methods to the bus level. (4) It needs no leader generator for collecting and processing global information, so it is fully bus-level distributed.

II. ECONOMIC DISPATCH MODEL CONSIDERING LOSSES AND GENERATOR POWER LIMITS
The economic dispatch model formulated in [1] is most accepted and familiar to power engineers, which is aiming at solving the problem of (i) minimizing the total generation cost, (ii) letting the total system power balance equation constraint be satisfied (network losses may be considered using B coefficients), (iii) letting each generator operate within its lower and upper limits.
In this paper, the set of bus active power balance equations is used to replace the total system power balance equation constraint in [1] to reasonably consider network losses. Considering that all bus voltage magnitudes are close to 1.0 p.u. in actual power systems, the economic dispatch model to be studied in this paper is formulated as follows.
A reference bus voltage phase must be selected and set to zero, θ r = 0, in solving model (1)-(3) since the bus voltage phases are not physically independent.
In model (1)-(3), (1) is the quadratic objective function ''total generation cost'' to be minimized, (2) is bus i power balance equation constraint which reasonably and implicitly contains network losses that vary with system load level, (3) is bus i generator power limit inequality constraint. Compared to the model in [1], model (1)-(3) contains the set of nonlinear bus power balance equations rather than the total system power balance equation, and all buses' voltage phases θ i s are newly added unknowns. Thus, the number of equation constraints increases by n (counting θ r = 0), so does the number of unknowns. Consequently, model (1)-(3) is a large-scale nonlinear programming just simpler than AC OPF models.

III. KKT OPTIMALITY CONDITIONS AND SIMPLIFICATIONS
First, the Lagrangian of the model (1)-(3) is given in this section. Then the basic KKT optimality conditions including inequalities are presented. Finally, a set of equivalent pure equality-type optimality conditions is derived by simplification to avoid the inherent defect of inequality optimality conditions failing to generate iteration functions. VOLUME 10, 2022

A. BASIC KKT OPTIMALITY CONDITIONS
According to the optimization theory, the Lagrangian L of the model (1)-(3) is as follows.
The dual variable λ i is associated with bus i power balance equation constraint (2) and is also called a Lagrangian multiplier. It is the marginal cost of bus i active power. The dual variables γ Li and γ Ui are associated with the lower and upper power limit inequalities in (3) respectively, they are the penalty costs of bus i generator active power limits.
According to the optimization theory, the optimal solution to an optimization problem must satisfy the KKT optimality necessary conditions. Based on (4), the KKT optimality conditions for model (1)-(3) are as follows.

B. PURE EQUALITY-TYPE KKT OPTIMALITY CONDITIONS
An iteration function f i (x i ) is the one that satisfies the equation f i (x i ) = x i . Naturally, it can only be derived from equations rather than inequalities. The KKT conditions (8), (9), and (12) are, however, inequalities. They inherently fail to generate iteration functions. So (7)- (12) are going to be simplified by removing variables γ Li and γ Ui to produce pure equality-type conditions.
First, look at (7). If the active power P i of the generator on bus i violates/equal a limit, it can either be the lower limit P Li or upper limit P Ui (namely, P i ≤ P Li or P i ≥ P Ui ), and not both simultaneously. Thus, either inequality condition (8) or (9) is active at a time, that is, either γ Li or γ Ui exists (means nonzero), but never both. So equation (7) can be rewritten as follows by replacing γ Li and γ Ui with γ i .
It is evident that (14) at any feasible solution, with λ i from (6), is the negative gradient −∂F ∂P i of the primal objective function F in (1) with respect to P i . It is γ i that guarantees the correct direction for minimizing F in solving for unknowns.
At last, merging (14) and (15), or finding their intersection parts, yields From (16), we know that γ i = γ Li = γ Ui = 0 and P i = (λ i − β i ) / (2α i ) if P Li < P i < P Ui . So we may take (λ i − β i ) / (2α i ) as a guess solution of P i , denoting it by f Pi produces f Pi = (λ i − β i ) / (2α i ). In addition, −2α i P i − β i + λ i ≤ 0 indeed if P i = P Li and −2α i P i − β i + λ i ≥ 0 indeed if P i = P Ui since the cost function α i P 2 i + β i P i + µ i of the generator on bus i is convex. Consequently, (16) can be rewritten as follows after removing γ i .
According to (17), P i is determined by the following equation based on the local (bus i) variable λ i .
The operator P P i projects the value (λ i − β i ) / (2α i ) onto the feasible interval [P Li , P Ui ]. Consequently, combined constraints (7)- (12) can be equivalently simplified into the equality-type condition (18) that removes variables γ Li and γ Ui . The variables γ Li and γ Ui can be computed by (17) at last. Based on the above analysis, the KKT optimality conditions (5)-(12) for solving model (1)-(3) are equivalent to (5), (6), and (18), which remove variables γ Li and γ Ui . So θ i , λ i and P i can be solved without consideration of γ Li and γ Ui , i = 1, 2, · · · , n. This practice will be adopted throughout the following parts of this paper.

A fixed point for an iteration function
We are going to use iteration functions to simultaneously compute the solutions of the equivalent KKT equations. So a local bus KKT equation should be used to generate an iteration function for computing one of this bus's unknowns. Let x i = [θ i , λ i , P i ] T be bus i vector-valued variable, the iteration function for computing the θ i is going to be derived from one of (5), (6), and (18) in this section. So are the iteration functions for computing the λ i and P i .
First, look at (18). It is obvious that (18) can directly be written as follows.
Formula (19) is an explicit function that maps bus i vectorvalued variable [θ i , λ i , P i ] T to its one component P i . Thus, it is the expected iteration function for computing variable P i .
Second, look at (5). Equation (5) contains only two local (bus i) unknown variables, θ i and P i . The iteration function for computing variable P i has been obtained (see (19)). So the iteration function for computing variable θ i should be generated from (5). The variable θ i in (5) is, however, completely contained in the sine and cosine functions rather than in polynomial functions. So we can't write out the iteration function for θ i directly from (5).
Considering that the coefficients of sin θ ij and cos θ ij are b ij and g ij respectively, and b ij is generally bigger than g ij in actual power systems, we expand sin θ ij in (5) at θ ij = θ i − θ j = 0 using the Taylor series as follows. where By (20), even if θ ij = 90 • , the relative error of sin θ ij is −6.6278×10 −10 , very small. So we can view (20) as an exact expression of sin θ ij . As a result, the ''='' is used instead of ''≈'' in (20). Substituting (20) into (5) and manipulating yield Formula (22) is an explicit function that maps bus i vectorvalued variable [θ i , λ i , P i ] T to its one component θ i . Thus, it is the expected iteration function for computing variable θ i . The derivation of (22) overcomes the difficulty caused by the fact that bus voltage phases are completely contained in the sine and cosine functions.
At last, since the iteration functions (19) and (22) for computing θ i and P i have been derived from (18) and (5) respectively, the iteration function for computing λ i has to be generated from the rest KKT condition (6). Separating λ i in (6) and manipulating yield.
Formula (23) is an explicit function that maps bus i vectorvalued variable [θ i , λ i , P i ] T to its one component λ i . Thus, it is the expected iteration function for computing variable λ i . The right-hand side functions of (22), (23), and (19) constitute the vector-valued iteration function f i (x i ), which maps bus i vector-valued variable

V. ITERATION PROCEDURES AND ALGORITHM
Based on the vector-valued iteration function f i (x i ) that consists of the right-hand side functions of (22), (23), and (19), a basic fixed-point iteration procedure is established for bus agents to solve the economic dispatch model (1)-(3) in this section. Then a fast version of the basic procedure is put forward. At last, the algorithm for the fast fixed-point iteration procedure and the way of variable initialization are presented. VOLUME 10, 2022

A. BASIC FIXED-POINT ITERATION PROCEDURE
Let the variables on the right-hand side of (22), (23), and (19) be the (k)-th step estimated solution, and let the variables on the left-hand side of them be the (k + 1)-th step estimated solution. Then (22), (23), and (19) become Updating θ i , λ i , and P i , using (24a)-(26a), independently by each of all bus agents and repeating the updating will produce the estimated solution sequence {θ 1(k) , λ 1(k) , P 1(k) , . . . , θ n(k) , λ n(k) , P n(k) } ∞ k=0 , which approaches the exact solution of (1)-(3). The set of formulas (24a)-(26a) originates from iteration functions (22), (23), and (19). Thus, it is a fixed-point iteration procedure. The fixed-point iteration procedure (24a)-(26a) is directly derived from the KKT conditions of (1)-(3), so we call it the basic fixed-point iteration procedure for bus agents to solve the economic dispatch model (1)- (3). Compared to the existing algorithms for economic dispatch, the basic fixed-point iteration procedure is different in both derivation method and result expressions, so it is novel.

B. FAST FIXED-POINT ITERATION PROCEDURE
Examining the iteration procedure (24a)-(26a), the variables on the right-hand sides of (24a)-(26a) all have the subscript k. So the (k + 1)-th step updating for bus i unknown variables is completely based on the (k)-th step estimated solutions. When the bus i agent updates bus i unknown variables in the (k + 1)-th step, the (k + 1)-th step estimated solutions θ j(k+1) and λ j(k+1) of neighbor buses j ∈ N i may have been obtained, such as the case when j < i in the laboratory simulation. Since θ j(k+1) and λ j(k+1) are closer to their exact solutions than θ j(k) and λ j(k) , replacing θ j(k) and λ j(k) with θ j(k+1) and λ j(k+1) in the (k + 1)-th step updating for bus i unknown variables will certainly speed up the convergence of the iteration. So we have the following modified iteration procedure.
The subscriptk implies its related variable takes the (k +1)-th step estimated solution if available. This modification speeds up the convergence of the iteration. The (k + 1)-th step estimated solutions θ i(k+1) and λ i(k+1) (instead of θ i(k) and λ i(k) ) of bus i are also used on the right-hand sides of (24b)-(26b). It further speeds up the convergence of the iteration. So the new iteration procedure (24b)-(26b) is called the fast fixed-point iteration procedure for bus agents to solve the economic dispatch model (1)-(3).
The updating of bus i unknown variables by the fast fixedpoint iteration procedure (24b)-(26b) needs estimated solutions θ i(k) , λ i(k) , and P i(k) of bus i itself, and θ j(k) and λ j(k) of its neighbor buses ∈ N i (θ j(k) and λ j(k) are delivered via power line carrier communication). Thus, it is bus-level distributed. Moreover, it needs no leader generator for collecting and processing global information, so it is fully buslevel distributed (or bus-level decentralized). The updating of bus i unknown variables by (24b)-(26b) does not need private information about neighbor buses j ∈ N i (such as the load power, generator power, cost coefficients, and generator limits), so the fast fixed-point iteration procedure preserves the privacy of each bus from leakage (bus-privacypreserving), and consequently raises the privacy-preserving level from the intra-area level of the existing DED methods to the bus level. Unlike the existing DED methods, it considers network losses by bus active power balance equations and does not need model reformulation and tuning parameters (does not complicate the solution), so it is a straightforward distributed optimization method reasonably considering network losses. Since the fast fixed-point iteration procedure is fully bus-level distributed and the exchange of information between neighbor buses is assumed to be implemented by power line carrier communication, it is easily scalable. The above conclusions or advantages are made for the structures of the model/expressions, so they are always true. Thus, the fast fixed-point iteration procedure (24b)-(26b), as the final algorithm of this paper, has significant advantages over the existing DED methods.

C. INITIALIZATION AND ALGORITHM
To start the calculation using the fast fixed-point iteration procedure (24b)-(26b), all variables must be assigned with initial values. For the convenience of comparison and analysis, the way of variable initialization is given below.
(1) The initial value of each bus voltage phase θ i(0) is assigned with zero. (2) The initial value of each generator's active power P i (0) is assigned proportionally to the mean value of its lower and upper power limits. is assigned with the mean value of initial marginal prices 2α i P i(0) + β i of all generators. The pseudo-code type of algorithm for solving (1)-(3) by the fast fixed-point iteration procedure (24b)-(26b) is designed and given in the following. The stopping criterion is ., n is the max relative correction of variables, r(x (k+1) ) ∞ is the max residual of KKT equations (5)- (7). To further improve the convergence, the updating by (24b)-(26b) for each bus unknown variable in each iteration step is repeated up to 4 times and broken if the max relative correction of these bus variables is less than ε.   Tables 1 and 2 show the economy and convergence related results of the 9-, 30-and 57-bus systems obtained by the fast FPI procedure and the basic FPI procedure respectively.

Algorithm of the Fast Fixed-Point Iteration Procedure
Comparing the values of the total cost in Tables 1 and 2, it is obvious that the values of the total cost of the same system by these two procedures are the same (their relative error is close to zero). So are the values of the total losses and the values of the loss rate. This is because the expressions of the two iteration procedures are identical in structure except for the number of iterations. Consequently, the generators' power by the fast FPI procedure and the basic FPI procedure are the same. So are the generator bus lambdas. Tables 3 and 4 show the generators' power and generator bus lambdas by the fast FPI procedure. Fig. 1 and Fig. 2 show the generators' power versus the number of iterations (P − k) curves and bus lambdas versus the number of iterations (λ − k) curves of the 9-bus system. Fig. 3 and Fig. 4 show the same curves of the 30-bus system. Fig. 5 and Fig. 6 show the same curves of the 57-bus system. All these curves are drawn by the fast FPI procedure based on the estimated solutions of the first 1/2 iterations to make oscillation parts of these curves clearer. In the simulation, the buses selected as a voltage phase reference for the 9-, 30-and 57-bus systems are bus 5, bus 21, and bus 38, respectively. They are selected by trial to let the number of iterations needed to reach convergence be the smallest.
The significant advantages or conclusions of the fast FPI procedure (24b)-(26b), given in sub-section V.B, over the existing DED methods are made for the structures of the model/expressions, so they are always true and no simulation verification required. So other characteristics of the fast FPI procedure are further revealed based on simulations in the following.    Table 1 shows that the maximum residual of KKT equations of the 9-bus system is 9.95724e-06<1e-5, small enough. Namely, on the one hand, the bus power balance equation constraints ∂L ∂λ i = 0, i = 1, 2, · · · , n, are all satisfied, and the obtained solution of model (1)-(3) is feasible. On the other hand, the optimality conditions ∂L ∂θ i = 0 and ∂L ∂P i = 0 and ∂L ∂λ i = 0 are all satisfied, and the obtained solution is optimal or suboptimal, owing to the practical meaning of the OPF problems [26], [27]. This feasibility and optimality are true for the solutions of the other systems by the fast FPI procedure. Consequently, the solution of the

2) ON THE RAPIDITY OF THE FAST FPI PROCEDURE
Comparing the results in Table 1 and Table 2, the numbers of iterations needed for solving the 9-bus system by the fast and basic FPI procedures are 656 and 1232 respectively, the latter is 1.878 times as big as the former, so the former is much faster than the latter. This is true for the other systems. Consequently, the fast FPI procedure is much faster than the basic one. This is because the former updates variables using the most up-to-date estimated solutions rather than the last step estimated solutions used by the latter. Table 1 shows that the transmission loss rates of the 9-, 30-and 57-bus systems are 1.24918%, 1.33801%, and 1.42167% respectively, that of the 30-bus system is in the middle. However, Table 4 shows that the minimum and maximum lambdas of generator buses in the 30-bus system are 3.74454 and 3.89495 respectively, their relative deviation reaches 4.01692%. Since the lambdas of generator buses are the locational marginal prices of these buses, their relative deviation of 4.01692% is economically notable. This is true   for the other systems. Consequently, the impact of network losses on bus lambdas is notable.

4) ON THE DISCREPANCY BETWEEN P − k CURVES AND λ − k CURVES
According to (26), when a generator output power P is not binding on its lower/upper limit, P is proportional to λ, that means the generator's P − k curve and λ − k curve should have the same shape. Fig. 1 shows that the P 2 − k curve of the generator on bus 2 in the 9-bus system is the top one which has the shape of rise-fall-rise-fall. And the first rise finishes in a very short period. Fig. 2 shows, however, the λ 2 −k curve of bus 2 in the 9-bus system is the bottom one which has the shape of fall-rise-fall. So bus-2-related P 2 −k curve and λ 2 −k curve do not have the same shape. This is true for all buses of all systems. Consequently, the approaching processes of the same bus's generator power and lambda to their exact solutions are not in step. This is because λ depends on θ (see (25)), θ depends on P (see (24)), and the impacts of the latter on the former exist time delay due to distributed calculation.

VII. CONCLUSION
The newly proposed method ''fast fixed-point iteration procedure method'' for economic dispatch is derived based on KKT optimality conditions and iteration functions. It not only reasonably considers network losses by the set of bus power balance equations but preserves the privacy of each bus from leakage by bus-level distributed optimization (without the need for exchange of neighbor bus privacy). It is a fast and straightforward distributed economic dispatch method since it neither needs to reformulate the model nor needs to introduce the tuning parameters. It is fully bus-level distributed since it needs no leader generator for collecting and processing global information. It raises the privacy-preserving level from the intra-area level of the existing distributed economic dispatch methods to the bus level. Due to the bus-privacy preservation and distributed calculation, it provides a more efficient tool for the economic operation of future power systems with increasing distributed renewable power sources.
The proposed method does not consider line capacity limit constraints, which limits its application. The disposal of these constraints is challenging since the complementary slackness conditions associated with them include inequalities that fail to directly generate iteration functions. So extending this method to consider line capacity limits and so forth is our future work.