Patent Publication Number: US-8977524-B2

Title: Interior point method for reformulated optimal power flow model

Description:
TECHNICAL FIELD 
     This disclosure is directed to methods of analyzing and modeling electric power grids, in particular smart grids. 
     DISCUSSION OF THE RELATED ART 
     An Optimal Power Flow (OPF) model of an electric power network is a model whose solution provides the conditions which give the lowest cost per kilowatthour delivered. Electricity flowing through an AC network is subject to Kirchhoff&#39;s Laws, and transmission lines are subject to thermal limits as well as voltage and electrical stability constraints. An OPF model should be capable to quickly and efficiently analyze an electric power network and yield a feasible power flow solution under a wide variety of adverse conditions and determine the most efficient set of actions to take. 
     OPF models present challenges even for cutting-edge mathematical programming techniques due to a sheer size of contingency cases that fully couple all network components for all adverse scenarios. Decoupling presents issues as it is well known that contingency analysis of power networks lacks even the most basic intuitive properties such as monotonicity, in that if tripping a branch results in an infeasible network, then tripping this branch and an extra branch should also be infeasible, but often is not. 
     A basecase OPF model seeks to produce an optimal power flow solution balancing generation/consumption for all buses with power distribution across all connected branches such that the overall system state, as measured by underlying physics parameters including voltages and phase angles, is feasible and safe. 
     A basecase OPF model can control AC and DC energy dispatches, regulate AC bus voltages and phase angles, select AC transformer/phasor tap and shunt switches, and control the AC/DC converter. 
     In a basecase OPF model, an electric power network is made up of the following components and their controls/states:
         AC Bus: represented by a node in the network graph. An AC bus has the following controls &amp; states:
           Load: Utilizes active (MW) and/or reactive (MVar) power;   Generator: Provides active (MW) and reactive (MVar) power;   Voltage and Phase angle: restricted to certain range around nominal;   Shunt capacitor/inductor: are optional, can be fixed or switching, can be positive (capacitor) or negative (inductor) and moderates active and reactive net power produced at the bus;   The main AC bus parameters are: bounds for Load and Generation, and bounds for Voltage and Phase Angle.   
           AC Branch/Line: represented by an arc in the network graph and transmits both active and reactive power. An AC branch has the following controls and states:
           Power flow: active and reactive power injected into the branch by both of its connected buses;   Adjusted impedance: is optional and represents branch impedance adjusted by an attached transformer;   The main AC branch parameters are: Admittance matrix, calculated based on the branch Impedance and the branch fixed shunt admittances, and Flow Capacity Ratings, which are the minimum and maximum total flow as constrained by thermal ratings.   
           AC Transformer: is an extra add-on to the AC branch that may either transform its connected buses&#39; voltages, (referred to herein as Transformer mode) or shift the difference between its connected buses&#39; phase angles (referred to herein below as Phasor mode). The AC Transformer has the following controls and states:
           Tap choice;   The main AC Transformer parameters are Tap Turns Ratio for Transformer mode, and Tap Phase Shift for Phasor mode.   
           DC Bus: represented by a node in the network graph. A DC bus has the following controls &amp; states:
           DC Voltage and Current Injection;   The main DC bus parameters are: bounds for Voltage and Current.   
           DC Branch/Line: represented by an arc in the network graph and transmits a constant current. A DC branch has the following controls and states:
           DC Current;   The main DC bus parameter is: Resistance.   
           Converter: represented by an arc in the network graph that connects AC and DC buses. Converters can be of two types: Rectifier, for AC to DC, and Inverter for DC to AC. A converter has the following controls and states:
           Firing and Overlap angles, AC side voltage, since the converter comes with a Transformer;   Power Factor, Active Power injection;   The main Converter parameters are: bounds for all converter controls and states and Commutating Impedance.   
               

     Other model components include a Voltage Source Converter (VSC), which is a special type of the converter that does not involve firing/overlap angles and directly controls voltage, power factor and active power injection, and FACTS devices, which are not modeled as distinct components but rather define tight bounds on certain controls/states of network components that are regulated by FACTS. Buses may be grouped in areas for defining area power exchange and inter-area power transfer constraints. Area involves several Interfaces that are corresponding groupings of branches that connect one area to another. 
       FIG. 1  is labelled to show each of the elements that may exist in an electric power system network. Their symbolic representation may also be seen. 
     An OPF mathematical formulation according to an embodiment of the invention is focused on making a base-case model formulation suitable for optimization solvers, in particular, utilizing as few as possible non-linear equations and as few as possible integer variables, and is extendable to contingency cases. 
     A power network is modeled as a standard graph N={V(N), E(N)} where V(N) is the set of buses and E(N) is the set of branches. Within V(N), one can distinguish the following subsets: the AC bus set V(AC) and the DC bus set V(DC). Within E(N) one can distinguish the set of AC branches E(AC), the set of DC branches E(DC), and the set of converters V(C). Also, one cane distinguish inter-area interface branch groupings E(I). 
     Every control and state variable mentioned in the previous section has a box constraint that are represented by lower and upper bound constraints, not all of which are explicitly mentioned below. An example of box constraint is the active power generation control variable P ig , which is subject to generation bound box constraints:
 
 P   ig min   ≦P   ig   ≦P   ig max .
 
Note that all choice variables in an OPF model according to an embodiment of the invention, which reflect the choice of each transformer/phasor tap or each shunt switch, denoted w k , not only are subject to be within the [0; 1] box but also should satisfy the Exclusive Choice constraint:
 
                 ∑     k   =   1     K     ⁢           ⁢     w   k       =   1         
The set of AC bus constraints comprises the fundamental set of network constraints, called the Power Flow Conservation set (PFC):
 
                 P   ig     -     P   il     +     P     i   ,   shunt         =         ∑     e   ∈       δ   +     ⁡     (   i   )           ⁢           ⁢     P   e       +       ∑     c   ∈       δ   +     ⁡     (   i   )           ⁢           ⁢     P   c                         Q   ig     -     Q   il     +     Q     i   ,   shunt         =         ∑     e   ∈       δ   +     ⁡     (   i   )           ⁢           ⁢     Q   e       +       ∑     c   ∈       δ   +     ⁡     (   i   )           ⁢           ⁢       Q   c     .               
The l.h.s. reflects the net power injection generation and the r.h.s. reflects the net power flow coming from both V(AC) and V(C). The shunt term in the PFC constraint is given by:
 
               P     i   ,   shunt       =       ∑     k   =   1       K   i       ⁢           ⁢       w   ik     ·     P   ik                       Q     i   ,   shunt       =       ∑     k   =   1       K   i       ⁢           ⁢       w   ik     ·     Q   ik               
The set of AC branch constraints includes:
 
               Voltage   ⁢           ⁢   Law   ⁢           ⁢     (   VL   )       :             P   ij     =         V   ii     ⁢     G   ii       +       V   ij     ⁡     (         G   ij     *     cos   ⁡     (     θ   ij     )         +       B   ij     *     sin   ⁡     (     θ   ij     )           )                       Q   ij     =         -     V   ii       ⁢     B   ii       +       V   ij     ⁡     (         G   ij     *     sin   ⁡     (     θ   ij     )         -       B   ij     *     cos   ⁡     (     θ   ij     )           )                       
Note that the variables V ij  are given by the Voltage Tensor product represented by bilinear constraints:
 
 V   ij   ={circumflex over (V)}   i   ·{circumflex over (V)}   j  
 
In addition, the variables θ ij  are given by the Angle Difference relation represented by:
 
               θ   ij     =       θ   i     -     θ   j     -       ∑     k   =   1       K   p       ⁢           ⁢       w     p   ⁢           ⁢   k       ⁢     Δ     p   ⁢           ⁢   k                   
where
 
               ∑     k   =   1       K   p       ⁢           ⁢       w     p   ⁢           ⁢   k       ⁢     Δ     p   ⁢           ⁢   k               
is the phase shift due to a Phasor (if any).
 
     Note also that for the sake of convenience for every physical undirected AC branch an OPF model according to an embodiment of the invention model has two logical directed AC branches: (i,j) and (j,i). These branches are not coupled but their parameters are calculated based on the same physical branch characteristics, such as impedance. 
     The set of Transformer constraints includes: 
     Joint mixed-integer-linear tap voltage box constraints:
 
 w   tk   · V     tk   ≦{circumflex over (V)}   tk   ≦w   tk   ·  V     tk .
 
     Link with Bus Voltage: 
     
       
         
           
             
               V 
               i 
             
             = 
             
               
                 ∑ 
                 
                   k 
                   = 
                   1 
                 
                 
                   K 
                   t 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   
                     V 
                     ^ 
                   
                   tk 
                 
                 · 
                 
                   
                     t 
                     tk 
                   
                   . 
                 
               
             
           
         
       
     
     Link with Branch (Tapped) Voltage: 
     
       
         
           
             
               
                 V 
                 ^ 
               
               i 
             
             = 
             
               
                 ∑ 
                 
                   k 
                   = 
                   1 
                 
                 
                   K 
                   t 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   
                     V 
                     ^ 
                   
                   tk 
                 
                 . 
               
             
           
         
       
     
     Note that in a model according to an embodiment of the invention there are two logical transformers corresponding to one physical transformer: one for (I,j) line and the other for (j,i) line. These two logical transformers not only have their parameters generated by the same underlying set of taps of the physical transformer but they also are explicitly coupled through an artificial component referred to as Coupling. Coupling holds any pair of component IDs and couples all of their choice variables one-to-one with each other such as:
 
for ( t   1   ,t   2 )εcoupling:  w   t     1     k   =w   t     2     k .
 
Similarly, Phasors are coupled.
 
     Note that the tap choice variables w tk  are restricted to be in [0; 1] box but in addition can be restricted to be binary by a Mixed-Integer integrability constraint: w t     1     k ε{0; 1} These constraints implies using an MIP solver and is expensive computationally while providing more responsive transformer control. 
     Impedance Adjustment constraints produce a choice of the admittance matrix (G, B) that depends on transformer/phasor tap position. For each Impedance Adjustment record a choice variable is introduced that is one-to-one (i.e. do not require extra integrability constraints) linked to transformer/phasor tap choice through the following linear constraints:
 
 w   am   ≧w   tk   ·I (   t     am   ≦t   tk   ≦  t     am )
 
The chosen admittance matrix is then given by the sum:
 
     
       
         
           
             
               ( 
               
                 G 
                 , 
                 B 
               
               ) 
             
             = 
             
               
                 ∑ 
                 
                   k 
                   = 
                   1 
                 
                 
                   K 
                   a 
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   w 
                   ak 
                 
                 · 
                 
                   
                     ( 
                     
                       G 
                       , 
                       B 
                     
                     ) 
                   
                   ak 
                 
               
             
           
         
       
     
     The set of Area constraints involve aggregation of branch power flows into area interface flows: 
     
       
         
           
             
               P 
               f 
             
             = 
             
               
                 ∑ 
                 
                   e 
                   ∈ 
                   
                     I 
                     ⁡ 
                     
                       ( 
                       f 
                       ) 
                     
                   
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 P 
                 e 
               
             
           
         
       
       
         
           
             
               Q 
               f 
             
             = 
             
               
                 ∑ 
                 
                   e 
                   ∈ 
                   
                     I 
                     ⁡ 
                     
                       ( 
                       f 
                       ) 
                     
                   
                 
               
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   Q 
                   e 
                 
                 . 
               
             
           
         
       
     
     The set of DC bus constraints comprises the fundamental set of network constraints called the Current Flow Conservation set (CFC): 
               I   i     =         ∑     e   ∈       δ   +     ⁡     (   i   )           ⁢           ⁢     I   e       -       ∑     e   ∈       δ   -     ⁡     (   i   )           ⁢           ⁢     I   e               
Note that I i  reflects external current injection from AC/DC converters and thus if no converters are attached to the DC node I then the box constraint for I i  is: 0≦I i ≦0
 
     The set of DC branch constraints comprises the DC version of the Voltage Law:
 
 I   ij   ·R   ij   =V   j   −V   i  
 
Finally, the set of Converter constraints comprises:
 
     
       
         
           
             AC 
             ⁢ 
             
               / 
             
             ⁢ 
             DC 
             ⁢ 
             
                 
             
             ⁢ 
             Voltage 
             ⁢ 
             
                 
             
             ⁢ 
             relation 
             ⁢ 
             
               : 
             
             ⁢ 
             
                 
             
             ⁢ 
             
               
                 
                   
                     
                       V 
                       DC 
                     
                     = 
                     
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         
                           V 
                           AC 
                         
                         · 
                         
                           ( 
                           
                             
                               cos 
                               ⁡ 
                               
                                 ( 
                                 α 
                                 ) 
                               
                             
                             + 
                             
                               cos 
                               ⁡ 
                               
                                 ( 
                                 μ 
                                 ) 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   
                     
                       V 
                       DC 
                     
                     = 
                     
                       
                         
                           V 
                           AC 
                         
                         · 
                         
                           cos 
                           ⁡ 
                           
                             ( 
                             
                               α 
                               + 
                               μ 
                             
                             ) 
                           
                         
                       
                       - 
                       
                         
                           I 
                           DC 
                         
                         · 
                         
                           R 
                           com 
                         
                       
                     
                   
                 
               
             
           
         
       
       
         
           
             Power 
             ⁢ 
             
                 
             
             ⁢ 
             factor 
             ⁢ 
             
                 
             
             ⁢ 
             relation 
             ⁢ 
             
               : 
             
             ⁢ 
             
                 
             
             ⁢ 
             
               
                 
                   
                     
                       cos 
                       ⁡ 
                       
                         ( 
                         R 
                         ) 
                       
                     
                     = 
                     
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         ( 
                         
                           
                             cos 
                             ⁡ 
                             
                               ( 
                               α 
                               ) 
                             
                           
                           + 
                           
                             cos 
                             ⁡ 
                             
                               ( 
                               μ 
                               ) 
                             
                           
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   
                     Q 
                     = 
                     
                       P 
                       · 
                       
                         tan 
                         ⁡ 
                         
                           ( 
                           R 
                           ) 
                         
                       
                     
                   
                 
               
             
           
         
       
         
         
           
             Active power injection from AC side: P=V DC ·I DC . 
           
         
       
    
     Note that the sign of active power injection is opposite for Rectifier and Inverter. 
     Finally, an OPF model can handle an objective function of arbitrary form. A basecase OPF objective function according to an embodiment of the invention is a generic convex quadratic cost of all buses&#39; generation and load: 
               min   ⁢       ∑     i   ∈     V   ⁡     (   AC   )           ⁢           ⁢     [           C   ig   a     ⁡     (       P   ^     ig     )       2     +         C   il   a     ⁡     (       P   ^     il     )       2     +         C   ig   r     ⁡     (       Q   ^     ig     )       2     +         C   il   r     ⁡     (       Q   ^     il     )       2     +       B   ig   a     ⁢       P   ^     ig       +       B   il   a     ⁢       P   ^     il       +       B   ig   r     ⁢       Q   ^     ig       +       B   il   r     ⁢       Q   ^     il         ]         ,         
where each active/reactive load or generation term reflects a linear transformation, that is, a centering and scaling with respect to a given reference value, of the corresponding model variables. Alternatively, another more complex objective function is one that minimizes the active power losses, defined as:
 
             min   ⁢       ∑     i   =   1     n     ⁢           ⁢       ∑     j   =   1     n     ⁢           ⁢       G   ij     ⁡     (       V   i   2     +       (       V   j       ri   j       )     2     -     2   ⁢         V   i     ⁢     V   j         r   ij       ⁢     cos   ⁡     (       θ   i     -     θ   j     +     φ   ij       )           )                 
where V i , V j  and θ i , θ j  represent the voltages and angles of the i-th and j-th buses, respectively, and G ij  represents the conductances of the branch connecting buses i and j. The ratio of the branch (i,j) is denoted by r ij , and the additional angle difference between buses i and j is denoted by φ ij . Note that φ ij =0 in the case where there is no phase shifter transformer in the branch (i, j). According an embodiment of the invention, a quadratic approximation of a non-linear objection function will be used.
 
     A basecase OPF model according to an embodiment of the invention as defined above is very complex for numerical optimization due to: 
     a) Non-convex Voltage Law constraints; 
     b) Fully coupled Power Flow Conservation network constraints; 
     c) Integer choice variables for transformer/phasor taps &amp; switching shunts; 
     d) Large network size in realistic examples; 
     From the above presentation it is evident that solving OPF models is a very challenging task. Its solution demands the development of sophisticated optimization algorithms which can handle both the large number of decision variables and the complex (non-convex) type of constraints. However, embodiments of the invention provide algorithms that have been successfully applied to the solution of very large scale electrical networks and the smart grid. 
     SUMMARY 
     Exemplary embodiments of the invention as described herein generally include systems and methods for approximating an optimal power flow of a smart electric power grid. 
     According to an aspect of the invention, there is provided a method for approximating an optimal power flow of a smart electric power grid, including providing a cost function that models a smart electric power grid having buses iε[1,N] connected by branches ijε[1,N]×[1,N] with connection matrix 
               R   ij     =     {           1   ,             ∃     (     i   ,   j     )       ,               0   ,             o   /   w     ,                   
in terms of bus control vector u i  that include voltage magnitudes V i , phase angles θ i , and shunt switch positions Q i• , bus dependent variables vector x i  that include net active x i   P  and reactive x i   Q  power injection, branch control vector u ij , that includes transformer tap selections t ij  and phase shifts φ ij , and dependent branch variable vector x ij , that includes net active x ij   P  and reactive x ij   Q  power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements, deriving a set of linear equations that minimize the cost function subject to the constraints, from an expression of an extremum of the cost function with respect to all arguments, reducing a dimension of the linear equations by solving for a subset of the linear equations, re-organizing the reduced dimension linear equations into primal and dual parts, and decomposing the re-organized reduced dimensional linear equations into two systems of block matrix equations, where a solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid.
 
     According to a further aspect of the invention, the cost function is 
                 f   ⁡     (     x   ,   u     )       =         ∑     i   =   1     N     ⁢           ⁢     [         1   2     ⁢     x   i   T     ⁢     V   i   X     ⁢     x   i       +       1   2     ⁢     u   i   T     ⁢     V   i   U     ⁢     u   i       +       S   i   X     ⁢     x   i       +       S   i   U     ⁢     u   i         ]       +       ∑     i   ,     j   =   1       N     ⁢           ⁢       R   ij     ·     [         1   2     ⁢     x   ij   T     ⁢     W   ij   X     ⁢     x   ij       +       1   2     ⁢     u   ij   T     ⁢     W   ij   U     ⁢     u   ij       +       T   ij   X     ⁢     x   ij       +       T   ij   U     ⁢     u   ij         ]             ,         
where bus control vector u i  includes voltage magnitudes V i , phase angles θ i , and shunt switch positions Q i• , bus dependent variables vector x i  includes net active x i   P  and reactive x i   Q  power injection, branch control vector u ij  includes transformer tap selections t ij  and phase shifts φ ij , and dependent branch variable vector x ij  includes net active x ij   P  and reactive x ij   Q  power injections, and connection matrix
 
               R   ij     =     {                 1   ,             ∃     (     i   ,   j     )       ,               0   ,             o   /   w     ,           ⁢           ⁢   where   ⁢           ⁢     x   i       =       ∑     j   =   1     N     ⁢           ⁢       R   ij     ·     x   ij           ,             
where V X  is a Hessian of the objective function with respect to x i , V U  is a Hessian of the objective function with respect to u i , S X  is a gradient of the objective function with respect to x i , S U  is a gradient of the objective function with respect to u i , W X  is a Hessian of the objective function with respect to x ij , W U  is a Hessian of the objective function with respect to u ij , T X  is a gradient of the objective function with respect to x ij , and T U  is a gradient of the objective function with respect to u ij , subject to the constraints
 
 x   i   LB   ≦x   i   ≦x   i   UB   , x   ij   LB   ≦x   ij   ≦x   ij   UB , and
 
 u   i   LB   ≦u   i   ≦u   i   UB   , u   ij   LB   ≦u   ij   ≦u   ij   UB ,
 
satisfying a voltage law x ij =R ij ·(A ij   0 u i +B ij   0 u j +C ij   0 u ij +D ij   0 u ji +E ij   0 ), where the matrices A ij   0 , B ij   0 , C ij   0 , D ij   0 , E ij   0  can be explicitly pre-calculated for every link (i, j) for a given anchor point (u i   0 , u j   0 , u ij   0 ).
 
     According to a further aspect of the invention, the linear equations that minimize the cost function are: 
                   [           H   11         0         -   I           H   41   T           H   51   T           -   I         0           0       0         -   I         0       0       0         -   I               -   I           -   I         0       0       0       0       0             H   41         0       0       0       0       0       0             H   51         0       0       0       0       0       0             diag   ⁡     (   Zx   )           0       0       0       0         diag   ⁡     (   X   )           0           0         diag   ⁡     (   Zs   )           0       0       0       0         diag   ⁡     (   S   )             ]     ·     [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   S               Δ   ⁢           ⁢   Y               Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2                 Δ   ⁢           ⁢     Z   x                 Δ   ⁢           ⁢     Z   s             ]       =     -     [             ∇     f   ⁡     (   X   )         -     Z   x     -   Y   -       y   1   T     ⁢     ∇     o   ⁡     (       x   ij     ,     x   i       )           -       y   2   T     ⁢     ∇     v   ⁡     (   X   )                         -   Y     -     Z   s                 -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -     v   ⁡     (   X   )                       Z   x     ·   X     -     μ   ⁢           ⁢   e                     Z   e     ·   S     -     μ   ⁢           ⁢   e             ]         ,         
where [ΔX, ΔS, ΔY, Δy 1 , Δy 2 , ΔZ x , ΔZ s ] are unknowns,
 
                     ⁢         H   11     =     diag   ⁡     (       V   x     ,     W   x     ,     V   u     ,     W   U       )         ,     
     ⁢           ⁢       H   41     =     [       -   I     ,   R   ,   0   ,   0     ]       ,     
     ⁢     R   =     (           [       R   11     ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R     1   ⁢           ⁢     N   1         ⁢   I     ]         0       …       0           0         [       R   21     ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R     2   ⁢           ⁢     N   2         ⁢   I     ]         …       0           0       ⋮       ⋱       ⋮           0       …       0         [       R       N   1     ⁢   1       ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R       N   1     ⁢     N   1         ⁢   I     ]           )       ,           
a row of branch (i, j) of H 51  is
 
                   (     H   51     )     ij     =     [     0   ,     -   I     ,     …   ⁢           ⁢     R   ij     ⁢     A   ij   0       ,     …   ⁢           ⁢     R   ij     ⁢     B   ij   0       ,     …   ⁢           ⁢     R   ij     ⁢     C   ij   0       ,     …   ⁢           ⁢     R   ij     ⁢     D   ij   0       ,     …   ⁢           ⁢   0       ]       ,     
     ⁢     X   =       (       {     x   il     }     ,     {     x   ijl     }     ,     {     u   il     }     ,     {     u   ijl     }       )     T       ,     
     ⁢       slack   ⁢           ⁢   variables   ⁢           ⁢   S     =       (       {     s   xi     }     ,     {     s   xij     }     ,     {     s   ui     }     ,     {     s   uij     }       )     T       ,     
     ⁢       o   ⁡     (       x   ij     ,     x   i       )       =       x     i   ⁢           ⁢   1       -       ∑     j   =   1     N     ⁢           ⁢       R   ij     ·     x     ij   ⁢           ⁢   1           -     b   x         ,     
     ⁢       v   ⁡     (   X   )       =       x     ij   ⁢           ⁢   1       -       R   ij     ·     (         A   ij   0     ·     u     i   ⁢           ⁢   1         +       B   ij   0     ·     u     j   ⁢           ⁢   1         +       C   ij   0     ·     u     ij   ⁢           ⁢   1         +       D   ij   0     ·     u     ji   ⁢           ⁢   1           )       -     b   v         ,         
Y, y 1 , and y 2  are dual variables, and Z X =μ/X, Z S =μ/S where μ is a barrier parameter.
 
     According to a further aspect of the invention, reducing a dimension of the linear equations further comprises solving for
 
Δ Z   X =diag( X ) −1 (− r   Z     X   −diag( Z   X )Δ X ),
 
where
 
Δ Z   S =diag( S ) −1 (− r   Z     S   −diag( Z   S )Δ S ),
 
 r   Z     X     =Z   X   ·X−μe, r   Z     S     =Z   S   ·S−μe.  
 
     According to a further aspect of the invention, the reduced dimension linear equations re-organized into primal and dual parts are 
     
       
         
           
             
               
                 [ 
                 
                   
                     
                       
                         H 
                         41 
                       
                     
                     
                       0 
                     
                     
                       | 
                     
                     
                       0 
                     
                     
                       0 
                     
                     
                       0 
                     
                   
                   
                     
                       
                         H 
                         51 
                       
                     
                     
                       0 
                     
                     
                       | 
                     
                     
                       0 
                     
                     
                       0 
                     
                     
                       0 
                     
                   
                   
                     
                       
                         - 
                         I 
                       
                     
                     
                       
                         - 
                         I 
                       
                     
                     
                       | 
                     
                     
                       0 
                     
                     
                       0 
                     
                     
                       0 
                     
                   
                   
                     
                       - 
                     
                     
                       - 
                     
                     
                       
                           
                       
                     
                     
                       - 
                     
                     
                       - 
                     
                     
                       - 
                     
                   
                   
                     
                       
                         
                           H 
                           11 
                         
                         + 
                         
                           
                             
                               diag 
                               ⁡ 
                               
                                 ( 
                                 X 
                                 ) 
                               
                             
                             
                               - 
                               1 
                             
                           
                           ⁢ 
                           
                             diag 
                             ⁡ 
                             
                               ( 
                               
                                 Z 
                                 x 
                               
                               ) 
                             
                           
                         
                       
                     
                     
                       0 
                     
                     
                       | 
                     
                     
                       
                         - 
                         I 
                       
                     
                     
                       
                         H 
                         41 
                         T 
                       
                     
                     
                       
                         H 
                         51 
                         T 
                       
                     
                   
                   
                     
                       0 
                     
                     
                       
                         
                           
                             diag 
                             ⁡ 
                             
                               ( 
                               S 
                               ) 
                             
                           
                           
                             - 
                             1 
                           
                         
                         ⁢ 
                         
                           diag 
                           ⁡ 
                           
                             ( 
                             
                               Z 
                               s 
                             
                             ) 
                           
                         
                       
                     
                     
                       | 
                     
                     
                       
                         - 
                         I 
                       
                     
                     
                       0 
                     
                     
                       0 
                     
                   
                 
                 ] 
               
               ⁡ 
               
                 [ 
                 
                   
                     
                       
                         Δ 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         X 
                       
                     
                   
                   
                     
                       
                         Δ 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         S 
                       
                     
                   
                   
                     
                       - 
                     
                   
                   
                     
                       
                         Δ 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         Y 
                       
                     
                   
                   
                     
                       
                         Δ 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           y 
                           1 
                         
                       
                     
                   
                   
                     
                       
                         Δ 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           y 
                           2 
                         
                       
                     
                   
                 
                 ] 
               
             
             = 
             
                 
               
                 
                   [ 
                   
                     
                       
                         
                           - 
                           
                             r 
                             
                               y 
                               1 
                             
                           
                         
                       
                     
                     
                       
                         
                           - 
                           
                             r 
                             
                               y 
                               2 
                             
                           
                         
                       
                     
                     
                       
                         
                           - 
                           
                             r 
                             Y 
                           
                         
                       
                     
                     
                       
                         - 
                       
                     
                     
                       
                         
                           
                             - 
                             
                               r 
                               X 
                             
                           
                           - 
                           
                             
                               
                                 diag 
                                 ⁡ 
                                 
                                   ( 
                                   X 
                                   ) 
                                 
                               
                               
                                 - 
                                 1 
                               
                             
                             ⁢ 
                             
                               
                                 r 
                                 Z 
                               
                               . 
                             
                           
                         
                       
                     
                     
                       
                         
                           
                             - 
                             
                               r 
                               S 
                             
                           
                           - 
                           
                             
                               
                                 diag 
                                 ⁡ 
                                 
                                   ( 
                                   S 
                                   ) 
                                 
                               
                               
                                 - 
                                 1 
                               
                             
                             ⁢ 
                             
                               
                                 r 
                                 Z 
                               
                               . 
                             
                           
                         
                       
                     
                   
                   ] 
                 
                 , 
                 
                   
 
                 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   
                     where 
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     [ 
                     
                       
                         
                           
                             r 
                             X 
                           
                         
                       
                       
                         
                           
                             r 
                             S 
                           
                         
                       
                       
                         
                           
                             r 
                             Y 
                           
                         
                       
                       
                         
                           
                             r 
                             
                               y 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               1 
                             
                           
                         
                       
                       
                         
                           
                             r 
                             
                               y 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               2 
                             
                           
                         
                       
                     
                     ] 
                   
                   = 
                   
                     
                       [ 
                       
                         
                           
                             
                               
                                 ∇ 
                                 
                                   f 
                                   ⁡ 
                                   
                                     ( 
                                     X 
                                     ) 
                                   
                                 
                               
                               - 
                               
                                 Z 
                                 x 
                               
                               - 
                               Y 
                               - 
                               
                                 
                                   y 
                                   1 
                                   T 
                                 
                                 ⁢ 
                                 
                                   ∇ 
                                   
                                     o 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         
                                           x 
                                           ij 
                                         
                                         , 
                                         
                                           x 
                                           i 
                                         
                                       
                                       ) 
                                     
                                   
                                 
                               
                               - 
                               
                                 
                                   y 
                                   2 
                                   T 
                                 
                                 ⁢ 
                                 
                                   ∇ 
                                   
                                     v 
                                     ⁡ 
                                     
                                       ( 
                                       X 
                                       ) 
                                     
                                   
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 - 
                                 Y 
                               
                               - 
                               
                                 Z 
                                 s 
                               
                             
                           
                         
                         
                           
                             
                               - 
                               
                                 ( 
                                 
                                   X 
                                   + 
                                   S 
                                   - 
                                   
                                     X 
                                     UB 
                                   
                                 
                                 ) 
                               
                             
                           
                         
                         
                           
                             
                               - 
                               
                                 o 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       x 
                                       ij 
                                     
                                     , 
                                     
                                       x 
                                       i 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                         
                           
                             
                               - 
                               
                                 v 
                                 ⁡ 
                                 
                                   ( 
                                   X 
                                   ) 
                                 
                               
                             
                           
                         
                       
                       ] 
                     
                     . 
                   
                 
               
             
           
         
       
     
     According to a further aspect of the invention, the two systems of block matrix equations are 
     
       
         
           
             
               
                 
                   [ 
                   
                     
                       
                         
                           H 
                           41 
                         
                       
                       
                         0 
                       
                     
                     
                       
                         
                           H 
                           51 
                         
                       
                       
                         0 
                       
                     
                     
                       
                         
                           - 
                           I 
                         
                       
                       
                         
                           - 
                           I 
                         
                       
                     
                   
                   ] 
                 
                 ⁡ 
                 
                   [ 
                   
                     
                       
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           X 
                         
                       
                     
                     
                       
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           S 
                         
                       
                     
                   
                   ] 
                 
               
               = 
               
                 [ 
                 
                   
                     
                       
                         - 
                         
                           r 
                           
                             y 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             1 
                           
                         
                       
                     
                   
                   
                     
                       
                         - 
                         
                           r 
                           
                             y 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             2 
                           
                         
                       
                     
                   
                   
                     
                       
                         - 
                         
                           r 
                           Y 
                         
                       
                     
                   
                 
                 ] 
               
             
             , 
             
               
                 
                   
                     and 
                     ⁢ 
                     
                       
 
                     
                     [ 
                     
                       
                         
                           
                             
                               H 
                               ~ 
                             
                             11 
                           
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           
                             
                               H 
                               ~ 
                             
                             22 
                           
                         
                       
                     
                     ] 
                   
                   ⁡ 
                   
                     [ 
                     
                       
                         
                           
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             X 
                           
                         
                       
                       
                         
                           
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             S 
                           
                         
                       
                     
                     ] 
                   
                 
                 + 
                 
                   
                     [ 
                     
                       
                         
                           
                             - 
                             I 
                           
                         
                         
                           
                             H 
                             41 
                             T 
                           
                         
                         
                           
                             H 
                             51 
                             T 
                           
                         
                       
                       
                         
                           
                             - 
                             I 
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                     
                     ] 
                   
                   ⁡ 
                   
                     [ 
                     
                       
                         
                           
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             Y 
                           
                         
                       
                       
                         
                           
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               y 
                               1 
                             
                           
                         
                       
                       
                         
                           
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               y 
                               2 
                             
                           
                         
                       
                     
                     ] 
                   
                 
               
               = 
               
                 [ 
                 
                   
                     
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           X 
                         
                       
                     
                   
                   
                     
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           S 
                         
                       
                     
                   
                 
                 ] 
               
             
             , 
             
               
 
             
             ⁢ 
             where 
           
         
       
       
         
           
             
               
                 
                   H 
                   ~ 
                 
                 11 
               
               = 
               
                 
                   H 
                   11 
                 
                 + 
                 
                   
                     
                       diag 
                       ⁡ 
                       
                         ( 
                         X 
                         ) 
                       
                     
                     
                       - 
                       1 
                     
                   
                   ⁢ 
                   
                     diag 
                     ⁡ 
                     
                       ( 
                       
                         Z 
                         X 
                       
                       ) 
                     
                   
                 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 
                   H 
                   ~ 
                 
                 22 
               
               = 
               
                 
                   
                     diag 
                     ⁡ 
                     
                       ( 
                       S 
                       ) 
                     
                   
                   
                     - 
                     1 
                   
                 
                 ⁢ 
                 
                   diag 
                   ⁡ 
                   
                     ( 
                     
                       Z 
                       S 
                     
                     ) 
                   
                 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 
                   r 
                   ~ 
                 
                 X 
               
               = 
               
                 
                   r 
                   X 
                 
                 + 
                 
                   
                     
                       diag 
                       ⁡ 
                       
                         ( 
                         X 
                         ) 
                       
                     
                     
                       - 
                       1 
                     
                   
                   ⁢ 
                   
                     r 
                     
                       Z 
                       X 
                     
                   
                 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 
                   r 
                   ~ 
                 
                 S 
               
               = 
               
                 
                   r 
                   S 
                 
                 + 
                 
                   
                     
                       diag 
                       ⁡ 
                       
                         ( 
                         S 
                         ) 
                       
                     
                     
                       - 
                       1 
                     
                   
                   ⁢ 
                   
                     r 
                     
                       Z 
                       S 
                     
                   
                 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 Δ 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 X 
               
               = 
               
                 [ 
                 
                   
                     
                       
                         { 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             x 
                             i 
                           
                         
                         } 
                       
                     
                   
                   
                     
                       
                         { 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             x 
                             ij 
                           
                         
                         } 
                       
                     
                   
                   
                     
                       
                         { 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             u 
                             i 
                           
                         
                         } 
                       
                     
                   
                   
                     
                       
                         { 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             u 
                             ij 
                           
                         
                         } 
                       
                     
                   
                 
                 ] 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 Δ 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   y 
                   1 
                 
               
               = 
               
                 { 
                 
                   Δ 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     y 
                     1 
                     
                       x 
                       i 
                     
                   
                 
                 } 
               
             
             , 
             
               
 
             
             ⁢ 
             
               
                 Δ 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   y 
                   2 
                 
               
               = 
               
                 
                   { 
                   
                     Δ 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       y 
                       2 
                       
                         x 
                         ij 
                       
                     
                   
                   } 
                 
                 . 
               
             
           
         
       
     
     According to a further aspect of the invention, the first block matrix equation is expanded for each bus i and branch(i, j) as follows: 
                 Δ   ⁢           ⁢     x   i       -       ∑   j     ⁢           ⁢       R   ij     ⁢   Δ   ⁢           ⁢     x   ij           =     r     y   ⁢           ⁢   1       x   i                       Δ   ⁢           ⁢     x   ij       -       R   ij     ⁡     (         A   ij     ⁢   Δ   ⁢           ⁢     u   i       +       B   ij     ⁢   Δ   ⁢           ⁢     u   j       +       C   ij     ⁢   Δ   ⁢           ⁢     u   ij       +       D   ij     ⁢   Δ   ⁢           ⁢     u   ji         )         =     r     y   ⁢           ⁢   2       x   ij                       Δ   ⁢           ⁢   X     =       Δ   ⁢           ⁢   S     =   rY       ,         
and the second block matrix equation is expanded for each bus i and branch(i, j) as follows:
 
                     V   ~     i   X     ⁢   Δ   ⁢           ⁢     x   i       -     Δ   ⁢           ⁢     Y     x   i         -     Δ   ⁢           ⁢     y   1     x   i           =     -       r   ~       x   i                             W   ~     ij   X     ⁢   Δ   ⁢           ⁢     x   ij       -     Δ   ⁢           ⁢     Y     x   ij         +       R   ij     ⁢   Δ   ⁢           ⁢     y   1     x   i         -     Δ   ⁢           ⁢     y   2     x   ij           =     -       r   ~       x   ij                             V   ~     i   U     ⁢   Δ   ⁢           ⁢     u   i       -     Δ   ⁢           ⁢     Y     u   i         +       ∑   j     ⁢       R   ij     ⁢     A   ij   T     ⁢   Δ   ⁢           ⁢     y   2     x   ij           +       ∑   j     ⁢       R   ji     ⁢     B   ji   T     ⁢   Δ   ⁢           ⁢     y   2     x   ji             =     -       r   ~       u   i                             W   ~     ij   U     ⁢   Δ   ⁢           ⁢     u   ij       -     Δ   ⁢           ⁢     Y     u   ij         +       R   ij     ⁢     C   ij   T     ⁢   Δ   ⁢           ⁢     y   2     x   ij         +       R   ji     ⁢     D   ji   T     ⁢   Δ   ⁢           ⁢     y   2     x   ji           =     -       r   ~       u   ij                               H   ~     22     ⁢   Δ   ⁢           ⁢   S     -     Δ   ⁢           ⁢   Y       =     -       r   ~     S         ,         
where the method further comprises solving for the other variables in the block matrix equations in the following order:
 
{Δ u   i   },{Δu   ij   }→{Δx   ij   }→{Δx   i   }→S→ΔY→Δy   1   →Δy   2 .
 
     According to a further aspect of the invention, the method includes initializing unknowns [ΔX k , ΔS k , ΔY k , Δy 1   k , Δy 2   k ] where k is an iteration counter, updating the barrier parameter as μ k+1 ←ε·μ k , εε(0.1, . . . , 0.9), evaluating (ΔX k , ΔS k , ΔY k , Δy i   k , Δy 2   k ) using the two block matrix equations, calculating (ΔZ X   k , ΔZ S   k ), determining a step length Δβ k  from a line search, and updating
 
 X   k+1   ←X   k +Δβ k   ΔX   k ,
 
 S   k+1   ←S   k +Δβ k   ΔS   k ,
 
 Y   k+1   ←Y   k +Δβ k   ΔY   k ,
 
 y   1   k+1   ←y   1   k +Δβ k   Δy   1   k ,
 
 Z   X   k+1   ←Z   X   k +Δβ k   ΔZ   X   k ,
 
 Z   S   k+1   ←Z   S   k +Δβ k   ΔZ   S   k ,
 
where the steps of updating the barrier parameter, evaluating (ΔX k , ΔS 4 , ΔY k , Δy i   k , Δy 2   k ), calculating (ΔZ X   k , ΔZ S   k ), determining a step length Δβ k , and updating X k+1 , S k+1 , Y k+1 , y 1   k+1 , y 2   k+1 , Z X   k+1  and Z S   k+1  are repeated until either a dual gap or a maximum of (r X , r S ) is less than a pre-specified minimum.
 
     According to a further aspect of the invention, the cost function is, 
                 f   ⁡     (     x   ,   u     )       =         ∑     i   =   1     N     ⁢     [         1   2     ⁢     x   i   T     ⁢     V   i   X     ⁢     x   i       +       1   2     ⁢       ∑   Ki     ⁢       u   i     k   ⁢           ⁢   T       ⁢     V   i   U     ⁢     u   i   k           +       S   i   X     ⁢     x   i       +       S   i   U     ⁢       ∑   Ki     ⁢     u   i   k           ]       +       ∑     i   ,     j   =   1       N     ⁢       R   ij     ·     [         1   2     ⁢     x   ij   T     ⁢     W   ij   X     ⁢     x   ij       +       1   2     ⁢       ∑   Kij     ⁢       u   ij     k   ⁢           ⁢   T       ⁢     W   ij   U     ⁢     u   ij   k           +       T   ij   X     ⁢     x   ij       +       T   ij   U     ⁢       ∑   Kik     ⁢     u   ij   k           ]             ,         
where bus control vector u i  includes voltage magnitudes V i , phase angles θ i , and shunt switch positions Q i• , bus dependent variables vector x i  includes net active x i   P  and reactive x i   Q  power injection, branch control vector u ij  includes transformer tap selections t ij  and phase shifts φ ij , and dependent branch variable vector x ij  includes net active x ij   P  and reactive x ij   Q  power injections, and connection matrix
 
               R   ij     =     {                 1   ,             ∃     (     i   ,   j     )       ,               0   ,             o   /   w     ,           ⁢     
     ⁢   where   ⁢           ⁢     x   i       =       ∑     j   =   1     N     ⁢       R   ij     ·     x   ij           ,             
where V X  is a Hessian of the objective function with respect to x i , V U  is a Hessian of the objective function with respect to u i , S X  is a gradient of the objective function with respect to x i , S U  is a gradient of the objective function with respect to u i  W X  is a Hessian of the objective function with respect to x ij , W U  is a Hessian of the objective function with respect to u ij , T X  is a gradient of the objective function with respect to x ij , and T U  is a gradient of the objective function with respect to u ij , subject to constraints
 
                 x   i   LB     ≤     x   i     ≤     x   i   UB       ,       u   i   LB     ≤       ∑   Ki     ⁢     u   i   k       ≤     u   i   UB       ,     
     ⁢       x   ij   LB     ≤     x   ij     ≤     x   ij   UB       ,       u   ij   LB     ≤       ∑   Kij     ⁢     u   ij   k       ≤     u   ij   UB       ,         
where a control grid is discretized by binary variables zε{0,1} as
 
 z   i   k   u   ik   LB   ≦u   i   k   ≦z   i   k   u   ik   UB ,
 
 z   ij   k   u   ijk   LB   ≦u   ij   k   ≦z   ij   k   u   ijk   UB ,
 
where the binary variables are subject to constraints
 
               z   i   LB     ≤       ∑   Ki     ⁢     z   i   k       ≤     z   i   UB                   z   i   LB     ≤       ∑   Ki     ⁢     z   i   k       ≤     z   i   UB           
and satisfying a voltage law x ij =R ij ·(Σ Ki A ij   k u i   k +Σ Kj B ij   k u j   k +Σ Kij C ij   k u ij   k +Σ Kji D ij   k u ji   k ) where the matrices A ij   0 , B ij   0 , C ij   0 , D ij   0 , E ij   0  can be explicitly pre-calculated for every link (i, j) for a given anchor point (u i   0 , u j   0 , u ij   0 ).
 
     According to a further aspect of the invention, the cost function is
 
ƒ( {tilde over (x)},u,y,v )=½ [{tilde over (x)}   T   V   X   {tilde over (x)}+u   T   V   U   u+y   T   W   Y   y+v   T   W   V   v]+S   X   {tilde over (x)}+S   U   u+T   Y   y+T   V   v,  
 
where bus control vector u includes voltage magnitudes V i , phase angles θ i , and shunt switch positions Q i• , bus dependent variables vector x includes net active X i   P  and reactive X i   Q  power injection, branch control vector v includes transformer tap selections t ij  and phase shifts φ ij , and dependent branch variable vector y includes net active X ij   P  and reactive X ij   Q  power injections, where bus dependent vector is defined as
 
               x   =           A   ~     ⁢     x   ~       ⇔     x   i       =         ∑     k   =   1       G   i       ⁢     x     i   ,   k     g       -       ∑     k   =   1       L   i       ⁢     x     i   ,   k     l       +       ∑     k   =   1       SH   i       ⁢     x     i   ,   k     sh             ,         
where G i , L, and SH i  are the number of generators, loads and shunts placed on bus i and x i,k   g , x i,k   l  and x i,k   sh  are the k th  active and reactive generation, load and shunt variables at bus i, respectively, and Ãε{−1,0,1} N×2(G+L+SH) , where V X  is a Hessian of the objective function with respect to x, V U  is a Hessian of the objective function with respect to u, S X  is a gradient of the objective function with respect to x, S U  is a gradient of the objective function with respect to u, W Y  is a Hessian of the objective function with respect to y, W V  is a Hessian of the objective function with respect to v, T Y  is a gradient of the objective function with respect to y, and T V  is a gradient of the objective function with respect to v, the cost function having connection matrix
 
               R   ij     =     {               1   ,             ∃     (     i   ,   j     )       ,               0   ,             o   /   w     ,           ⁢     
     ⁢   where   ⁢           ⁢     a   i     ⁢     x   ~       =       ∑     j   =   1     N     ⁢       R   ij     ⁢     y   ij                   
where a i  is i th  is row of matrix Ã, subject to the constraints
 
 {tilde over (x)}   LB   ≦{tilde over (x)}≦{tilde over (x)}   UB   , y   LB   ≦y≦y   UB ,
 
 u   LB   ≦u≦u   UB   , v   LB   ≦v≦v   UB ,
 
and satisfying a voltage law y ij =R ij ·(A ij   0 ·u i +B ij   0 ·u j +C ij   0 ·v ij +D ij   0 ·v ji +E ij   0 ), where the matrices A ij   0 , B ij   0 , C ij   0 , D ij   0 , E ij   0  can be explicitly pre-calculated for every link (i, j) for a given anchor point (u i   0 , u j   0 , u ij   0 ).
 
     According to another aspect of the invention, there is provided a method for approximating an optimal power flow of a smart electric power grid, including providing a cost function that models a smart electric power grid having buses iε[1,N] connected by branches ijε[1, N]×[1, N] with connection matrix 
               R   ij     =     {           1   ,             ∃     (     i   ,   j     )       ,               0   ,             o   /   w     ,                   
in terms of bus control vector u i  that include voltage magnitudes V i , phase angles θ i , and shunt switch positions Q i• , bus dependent variables vector x i  that include net active X i   P  and reactive X i   Q  power injection, branch control vector u ij , that includes transformer tap selections t ij  and phase shifts φ ij , and dependent branch variable vector x ij , that includes net active x ij   P  and reactive x ij   Q  power injections, constraining the bus and branch control vectors, and the bus and branch dependent variables vectors based on physical and operational requirements and actual measurements, deriving a set of non-linear vector equations
 
               L   ~     =       [             ∇     f   ⁡     (   X   )         -     μ   /   X     -   Y   -       y   1   T     ⁢     ∇     o   ⁡     (       x   ij     ,     x   i       )           -       y   2   T     ⁢     ∇     v   ⁡     (   X   )                           -   μ     /   S     -   Y               -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -     v   ⁡     (   X   )               ]     =   0           
that minimize the cost function subject to the constraints, from an expression of an extremum of the cost function with respect to all arguments, where X=({x i }, {x ij }, {u i }, {u ij }) T , bus control vector u i , includes voltage magnitudes V i , phase angles θ i , and shunt switch positions Q i• , bus dependent variables vector x, includes net active x i   P  and reactive x i   Q  power injection, branch control vector u ij  includes transformer tap selections t ij  and phase shifts φ ij , and dependent branch variable vector x ij  includes net active x ij   P  and reactive x ij   Q  power injections, S=({s xi }, {s xij }, {s ui }, {s uij }) T , where s xi , s ui , s xij , s uij  are slack variables,
 
                 o   ⁡     (       x   ij     ,     x   i       )       =       x     i   ⁢           ⁢   1       -       ∑     j   =   1     N     ⁢       R   ij     ·     x     ij   ⁢           ⁢   1           -     b   x         ,     
     ⁢       v   ⁡     (   X   )       =       x     ij   ⁢           ⁢   1       -       R   ij     ·     (         A   ij   0     ·     u     i   ⁢           ⁢   1         +       B   ij   0     ·     u     j   ⁢           ⁢   1         +       C   ij   0     ·     u     ij   ⁢           ⁢   1         +       D   ij   0     ·     u     ji   ⁢           ⁢   1           )       -     b   v               
b x  and b v  are predetermined constants,
 
               R   ij     =     {           1   ,             ∃     (     i   ,   j     )       ,               0   ,             o   /   w     ,                   
is a connection matrix where
 
                 x     i   ⁢               =       ∑     j   =   1     N     ⁢       R   ij     ·     x     ij   ⁢                     ,     
     ⁢       x     ij   ⁢               =       R   ij     ·     (         A   ij   0     ⁢     u     i   ⁢                 +       B   ij   0     ⁢     u     j   ⁢                 +       C   ij   0     ⁢     u     ij   ⁢                 +       D   ij   0     ⁢     u     ji   ⁢                 +     E   ij   0       )               
is a voltage law where the matrices A ij   0 , B ij   0 , C ij   0 , D ij   0 , E ij   0  can be explicitly pre-calculated for every link (i, j) for a given anchor point (u i   0 , u j   0 , u ij   0 ), Y, y 1 , and y 2  are dual variables, and μ is a barrier parameter, re-organizing the vector {tilde over (L)} as
 
                 L   ~     =     [           {     L   i     }               {     L   ij     }           ]       ,         
where each L i  represent the constraints of bus i, and each L ij  represents the constraints of branch (i, j), applying Newton&#39;s method to every bus L i  and branch L ij , to obtain
 
                     H   i     ·     [           Δ   ⁢           ⁢     X   i                 Δ   ⁢           ⁢     S   i                 Δ   ⁢           ⁢     Y   i                 Δ   ⁢           ⁢     y   1   i                 Δ   ⁢           ⁢     y   2   i             ]       =     -     L   i         ;         H   ij     ·     [           Δ   ⁢           ⁢     X   ij                 Δ   ⁢           ⁢     S   ij                 Δ   ⁢           ⁢     Y   ij                 Δ   ⁢           ⁢     y   1   ij                 Δ   ⁢           ⁢     y   2   ij             ]       =     -     L   ij           ,         
where H i  and H ij  are Jacobian matrices, ΔX i =({Δx i }, {Δu i }) T , ΔX ij =({Δx ij }, {Δu ij }), ΔS i =({Δs xi }, {Δs ui }) T , ΔS ij =({Δs xij }, {Δs uij }), and [ΔX, ΔS, ΔY, Δy 1 , Δy 2 ] are unknowns, and using cutting planes to solve these block-wise Newton&#39;s equations, where a solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid.
 
     According to another aspect of the invention, there is provided a non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executed by the computer to perform the method steps for approximating an optimal power flow of a smart electric power grid. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a diagram of an exemplary, non-limiting electric power system network, according to an embodiment of the invention. 
         FIG. 2  is a flow chart of an interior point method for solving an optimal power flow model of a smart power grid, according to an embodiment of the invention. 
         FIG. 3  is a flow chart of an iterative method for solving an optimal power flow model of a smart power grid, according to an embodiment of the invention. 
         FIG. 4  is a block diagram of an exemplary computer system for implementing an interior point method for solving an optimal power flow model of a smart power grid, according to an embodiment of the invention. 
     
    
    
     DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS 
     Exemplary embodiments of the invention as described herein generally include systems and methods for estimating a solution of an optimal power flow model of a smart power grid. Accordingly, while the invention is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that there is no intent to limit the invention to the particular forms disclosed, but on the contrary, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention. 
     An optimal power flow (OPF) model seeks to produce an optimal power flow solution that balances generation/consumption for all buses with power distribution across all connected branches such that the overall system state, as measured by underlying physical parameters including voltages and phase angles, is feasible and safe. A basecase OPF model can be used to control AC and DC energy dispatch (i.e. bus generation/load and branch power flow), regulate AC bus voltages and phase angles, select AC transformer/phasor taps and shunt switches, and control an AC/DC converter. 
     Solver Model Specification 
     An electric grid Optimal Power Flow (OPF) model according to an embodiment of the invention includes buses iε[1,N] connected by branches ijε[1,N]×[1,N]. The AC-AC branch (i.e., a branch obeying alternating current physics laws) connection matrix R is given as: 
               R   ij     =     {           1   ,             ∃     (     i   ,   j     )       ,               0   ,           o   /     w   .                     
Minimal AC OPF Model Specification
 
     A minimal OPF model according to an embodiment of the invention assumes an AC-only grid with continuous-only state-dependent variables x and controls u. In addition to that, there are no thermal capacity limits enforced on branch power flows. Finally, discrete logic z defining a linearization region of the originally non-convex AC voltage law constraints is assumed to be a fixed binary vector z 0 . 
     Box-Constraints: 
     According to an embodiment of the invention, the following two box-constraint families are defined that independently limit bus control vector u i  and dependent bus variable vector x i :
 
 x   i   LB   ≦x   i   ≦x   i   UB ,
 
 u   i   LB   ≦u   i   ≦u   i   UB ,  (1)
 
where the superscripts LB and UB refer to lower bound and upper bound, respectively. Bus control vector u i  typically includes voltage magnitude V i , log-scaled in this model, phase angle θ i , and shunt switch position Q i• . Bus dependent variables vector x i  typically includes net active x i   P  and reactive x i   Q  power injection. Similarly, according to an embodiment of the invention, a box-constraint family is defined for branch variables that independently limit branch control vector u ij , which typically includes transformer tap selection t ij  and phase shift φ ij , and dependent branch variable vector x ij , which typically includes net active x ij   P  and reactive x ij   Q  power injections:
 
 x   ij   LB   ≦x   ij   ≦x   ij   UB ,
 
 u   ij   LB   ≦u   ij   ≦u   ij   UB ,  (2)
 
     Flow Balance: 
     At every bus a net power injection balance should be maintained between a bus and connected branches. Branch dependent variables vector x ij , in which typically x ij =(x ij   P , x ij   Q ), can be defined to specify this balance as follows: 
     
       
         
           
             
               
                 
                   
                     x 
                     i 
                   
                   = 
                   
                     
                       ∑ 
                       
                         j 
                         = 
                         1 
                       
                       N 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         R 
                         ij 
                       
                       · 
                       
                         
                           x 
                           ij 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     The matrix R is always highly sparse since the connection density of real-world power grids is typically 1-5 branches per bus. 
     Discrete Control Grid: 
     Some of the components in bus control vector u i  may be discrete, while others may be continuous. According to an embodiment of the invention, there is a minimal discretization of the overall Cartesian product space of the vector u i  into intervals kε[1,K i ] each defined by a box [u ik   LB , u ik   LB ], a fixed anchor vector u ik   0  and a binary choice variable z ik . The latter is assumed fixed at z ik   0  in this minimal AC OPF specification. Then, a discrete control grid is given by the following discrete box-constraint family: 
                       ∑     k   =   1       K   i       ⁢           ⁢       Z   ik   0     ⁢     u   ik   LB         ≤     u   i     ≤       ∑     k   =   1       K   i       ⁢           ⁢       Z   ik   0     ⁢       u   ik   UB     .                 (   4   )               
Similarly, according to an embodiment of the invention, a discrete control grid is defined for the branch control vector u ij  defined by kε└1,K ij ┘ intervals each defined by a box [u ijk   LB , u ijk   UB ], a fixed anchor u ijk   0  and a binary choice variable z ijk , assumed fixed at z ijk   0  in this minimal AC OPF specification:
 
     
       
         
           
             
               
                 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       
                         K 
                         ji 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         Z 
                         ijk 
                         0 
                       
                       ⁢ 
                       
                         u 
                         ijk 
                         LB 
                       
                     
                   
                   ≤ 
                   
                     u 
                     ij 
                   
                   ≤ 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         1 
                       
                       
                         K 
                         ij 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         Z 
                         ijk 
                         0 
                       
                       ⁢ 
                       
                         
                           u 
                           ijk 
                           UB 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     Discretized Voltage Law: 
     Each interval k implies a new feasibility set X ij   k  for the branch power flow vector x ij  represented as a rectangle generated by intervals [u ik   LB , u ik   UB ], [u jk   LB , u jk   UB ], and [u ijk   LB , u ijk   UB ], following the voltage law gradient lines drawn at the combined anchor point (u ik   0 , u jk   0 , u ijk   0 ). According to an embodiment of the invention, the feasibility set X ij   k  is not explicitly included in the model specification. Instead, a discretized version of the AC voltage law is provided that defines the relation between branch power flow vectors x ij  and control vectors u i , u j , u ij  that are linearized with respect to a given anchor point (u i   0 , u j   0 , u ij   0 ) that is typically given by the grid anchor point (u ik   0 , u jk   0 , u ijk   0 ) corresponding to z ik   0  and z ijk   0 . Then, an embodiment of the invention assumes the following form of the discretized AC voltage law:
 
 x   ij   =R   ij ·( A   ij   0   u   i   +B   ij   0   u   j   +C   ij   0   u   ij   +D   ij   0   u   ji   +E   ij   0 ).  (6)
 
Note that this constraint exhibits two types of sparsity: one inherited from R ij  and the other due to dependence on only two control vectors (u i , u j ) and (u ij , u ji ) instead of the whole-grid control vector U=└(u i ) i=1   N , (u ij ) i,j=1   N ┘. The matrices A ij   0 , B ij   0 , C ij   0 , D ij   0 , E ij   0  can be explicitly pre-calculated for every link (i, j) for the given anchor point (u i   0 , u j   0 , u ij   0 ).
 
     Objective Function: 
     According to an embodiment of the invention, it is assumed that there exists a generic convex quadratic cost function to be minimized. This function includes terms related to managing the control vector u i  that minimizes control adjustment costs, bus dependent variables vector x i  that minimize generation costs, and branch dependent variables vector x ij  that minimize transmission losses: 
                       f   ⁡     (     x   ,   u     )       =         ∑     i   =   1     N     ⁢           ⁢     [         1   2     ⁢     x   i   T     ⁢     V   i   X     ⁢     x   i       +       1   2     ⁢     u   i   T     ⁢     V   i   U     ⁢     u   i       +       S   i   X     ⁢     x   i       +       S   i   U     ⁢     u   i         ]       +       ∑     i   ,     j   =   1       N     ⁢           ⁢       R   ij     ·     [         1   2     ⁢     x   ij   T     ⁢     W   ij   X     ⁢     x   ij       +       1   2     ⁢     u   ij   T     ⁢     W   ij   U     ⁢     u   ij       +       T   ij   X     ⁢     x   ij       +       T   ij   U     ⁢     u   ij         ]             ,           (   7   )               
where:
 
     V X  is a Hessian of the objective function with respect to x i , 
     V U  is a Hessian of the objective function with respect to u i , 
     S X  is a gradient of the objective function with respect to x i , 
     S U  is a gradient of the objective function with respect to u i , 
     W X  is a Hessian of the objective function with respect to x ij , 
     W U  is a Hessian of the objective function with respect to u ij , 
     T X  is a gradient of the objective function with respect to x ij , and 
     T U  is a gradient of the objective function with respect to u ij . 
     Discrete AC OPF Model Specification 
     A discrete OPF model according to an embodiment of the invention assumes that a discrete logic z defines a linearization region of the originally non-convex AC voltage law constraints is also a model variable that is constrained to be a binary zε{0,1}. A discrete OPF model according to an embodiment of the invention has a form similar to the minimal OPF except for the following changes: 
     Box-Constraints: 
     The dependent variables box-bounds remain unchanged:
 
 x   i   LB   ≦x   i   ≦x   i   UB ,
 
 x   ij   LB   ≦x   ij   ≦x   ij   UB .  (8)
 
However, the box-bound constraints for control variables u i  and u ij  are changed to have the following simplex form:
 
                       u   i   LB     ≤       ∑   Ki     ⁢           ⁢     u   i   k       ≤     u   i   UB       ⁢     
     ⁢       u   ij   LB     ≤       ∑   Kij     ⁢           ⁢     u   ij   k       ≤       u   ij   UB     .               (   9   )               
Also, a simplex form is introduced for the discrete variables:
 
                       z   i   LB     ≤       ∑   Ki     ⁢           ⁢     z   i   k       ≤     z   i   UB       ⁢     
     ⁢       z   i   LB     ≤       ∑   Ki     ⁢           ⁢     z   i   k       ≤     z   i   UB               (   10   )               
Discrete Control Grid:
 
 z   i   k   u   ik   LB   ≦u   i   k   ≦z   i   k   u   ik   UB ,
 
 z   ij   k   u   ijk   LB   ≦u   ij   k   ≦z   ij   k   u   ijk   UB .  (11)
 
     Note that the vector u ij   k  of branch control variables, constrained within a discrete control grid interval k, now also includes power flow anchor points (P ij   k , Q ij   k ) for every branch. 
     Discretized Voltage Law.
 
 x   ij   =R   ij ·(Σ Ki   A   ij   k   u   i   k +Σ Ki   B   ij   k   u   j   k +Σ Kij   C   ij   k   u   ij   k +Σ Kji   D   ij   k   u   ji   k )  (12)
 
     Objective Function: 
                     f   ⁡     (     x   ,   u     )       =         ∑     i   =   1     N     ⁢           ⁢     [         1   2     ⁢     x   i   T     ⁢     V   i   X     ⁢     x   i       +       1   2     ⁢       ∑   Ki     ⁢           ⁢       u   i   kT     ⁢     V   i   U     ⁢     u   i   k           +       S   i   X     ⁢     x   i       +       S   i   U     ⁢       ∑   Ki     ⁢           ⁢     u   i   k           ]       +       ∑     i   ,     j   =   1       N     ⁢           ⁢       R   ij     ·     [         1   2     ⁢     x   ij   T     ⁢     W   ij   X     ⁢     x   ij       +       1   2     ⁢       ∑   Kij     ⁢           ⁢       u   ij   kT     ⁢     W   ij   U     ⁢     u   ij   k           +       T   ij   X     ⁢     x   ij       +       T   ij   U     ⁢       ∑   Kik     ⁢           ⁢     u   ij   k           ]                   (   13   )               
Interior Point Method Solver Math for a Minimal AC OPF Model
 
     First, a model according to an embodiment of the invention is standardized to the following format: 
                     min   ⁢           ⁢     f   ⁡     (     x   ,   u     )         =         ∑     i   =   1     N     ⁢           ⁢     [         1   2     ⁢       x     i   ⁢           ⁢   1     T     ·     V   i   X     ·     x     i   ⁢           ⁢   1     T         +       1   2     ⁢       u     i   ⁢           ⁢   1     T     ·     V   i   U       ⁢     u     i   ⁢           ⁢   1     T       +       S     i   ⁢           ⁢   1     X     ·     x     i   ⁢           ⁢   1         +       S     i   ⁢           ⁢   1     U     ·     u     i   ⁢           ⁢   1           ]       +       ∑     i   ,     j   =   1       N     ⁢           ⁢       R   ij     ·     [               1   2     ⁢       x     ij   ⁢           ⁢   1     T     ·     W     ij   ⁢           ⁢   1     X     ·     x     ij   ⁢           ⁢   1     T         +                   1   2     ⁢       u   ij   T     ·     W   ij   U     ·     u     ij   ⁢           ⁢   1     T         +       T     ij   ⁢           ⁢   1     X     ·     x     ij   ⁢           ⁢   1         +       T     ij   ⁢           ⁢   1     U     ·     u     ij   ⁢           ⁢   1                 ]                   (   14.1   )                       ⁢       X   +   S     =       X   UB     ⁢     {               x     i   ⁢           ⁢   1       +     s   xi       =     x     i   ⁢           ⁢   1     UB                     u     i   ⁢           ⁢   1       +     s   ui       =     u     i   ⁢           ⁢   1     UB                     x     ij   ⁢           ⁢   1       +     s   xij       =     x     ij   ⁢           ⁢   1     UB                     u     ij   ⁢           ⁢   1       +     s   uij       =     u     ij   ⁢           ⁢   1     UB                           (   14.2   )                       ⁢         x     i   ⁢           ⁢   1       -       ∑     j   =   1     N     ⁢           ⁢       R   ij     ·     x     ij   ⁢           ⁢   1             =     b   x               (   14.3   )                       ⁢         x     ij   ⁢           ⁢   1       -       R   ij     ·     (         A   ij   0     ·     u     i   ⁢           ⁢   1         +       B   ij   0     ·     u     j   ⁢           ⁢   1         +       C   ij   0     ·     u     ij   ⁢           ⁢   1         +       D   ij   0     ·     u     ji   ⁢           ⁢   1           )         =     b   v               (   14.4   )                       ⁢       x     i   ⁢           ⁢   1       ,     x     ij   ⁢           ⁢   1       ,     u     i   ⁢           ⁢   1       ,     u     ij   ⁢           ⁢   1       ,     s   xi     ,     s   ui     ,     s   xij     ,       s   uij     ≥   0     ,             (   14.5   )               
where X=({x i1 }, {x ij1 }, {u u1 }, {u ij1 }) T , s xi , s ui , s xij , s uij  are slack variables, S=({s xi }, {s xij }, {s ui }, {s uij }) T , and
 
                           ⁢         x     i   ⁢           ⁢   1       =       x   i     -     x   i   LB         ,     
     ⁢           ⁢       x     ij   ⁢           ⁢   1       =       x   ij     -     x   ij   LB         ,     
     ⁢           ⁢       x     i   ⁢           ⁢   1     UB     =       x   i   UB     -     x   i   LB         ,     
     ⁢           ⁢       x     ij   ⁢           ⁢   1     UB     =       x   ij   UB     -     x   ij   LB         ,             (   15.1   )                       ⁢         u     i   ⁢           ⁢   1       =       u   i     -     max   ⁢     {       u   i   LB     ,       ∑     k   =   1       K   i       ⁢           ⁢       z   ik   0     ·     u   ik   LB           }           ,     
     ⁢           ⁢       u     ij   ⁢           ⁢   1       =       u   ij     -     max   ⁢     {       u   ij   LB     ,       ∑     k   =   1       K   ij       ⁢           ⁢       z   ijk   0     ·     u   ijk   LB           }           ,             (   15.2   )                       ⁢         u     i   ⁢           ⁢   1     UB     =       min   ⁢     {       u   i   UB     ,       ∑     k   =   1       K   i       ⁢           ⁢       z   ik   0     ·     u   ik   UB           }       -     max   ⁢     {       u   i   LB     ,       ∑     k   =   1       K   i       ⁢           ⁢       z   ik   0     ·     u   ik   LB           }           ,     
     ⁢           ⁢       u     ij   ⁢           ⁢   1     UB     =       min   ⁢     {       u   ij   UB     ,       ∑     k   =   1       K   ij       ⁢           ⁢       z   ijk   0     ·     u   ijk   UB           }       -     max   ⁢     {       u   ij   LB     ,       ∑     k   =   1       K   ij       ⁢           ⁢       z   ijk   0     ·     u   ijk   LB           }           ,             (   15.3   )                       ⁢         b   x     =         ∑     j   =   1     N     ⁢           ⁢       R   ij     ·     x   ij   LB         -     x   i   LB         ,             (   15.4   )                 b   v     =         R   ij     ·     (           A   ij   0     ·   max     ⁢     {       u   i   LB     ,       ∑     k   =   1       K   i       ⁢           ⁢       z   ik   0     ·     u   ik   LB           }       +         B   ij   0     ·   max     ⁢     {       u   j   LB     ,       ∑     k   =   1       K   i       ⁢           ⁢       z   ik   0     ·     u   jk   LB           }       +         C   ij   0     ·   max     ⁢     {       u   ij   LB     ,       ∑     k   =   1       K   ij       ⁢           ⁢       z   ijk   0     ·     u   ijk   LB           }       +         D   ij   0     ·   max     ⁢     {       u   ij   LB     ,       ∑     k   =   1       K   ij       ⁢           ⁢       z   ijk   0     ·     u   jik   LB           }       +     E     ij   ⁢           ⁢   1     0       )       -       x   ij   LB     .               (   15.5   )               
The coefficients of the objective function can also be changed:
 
     
       
         
           
             
               
                 
                   
                     
                       S 
                       
                         i 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         1 
                       
                       X 
                     
                     = 
                     
                       
                         S 
                         i 
                         X 
                       
                       + 
                       
                         
                           x 
                           i 
                           LB 
                         
                         · 
                         
                           V 
                           i 
                           X 
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       S 
                       
                         i 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         1 
                       
                       U 
                     
                     = 
                     
                       
                         S 
                         i 
                         U 
                       
                       + 
                       
                         max 
                         ⁢ 
                         
                           
                             { 
                             
                               
                                 u 
                                 j 
                                 LB 
                               
                               , 
                               
                                 
                                   ∑ 
                                   
                                     k 
                                     = 
                                     1 
                                   
                                   
                                     K 
                                     i 
                                   
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   
                                     z 
                                     ik 
                                     0 
                                   
                                   · 
                                   
                                     u 
                                     jk 
                                     LB 
                                   
                                 
                               
                             
                             } 
                           
                           · 
                           
                             V 
                             i 
                             U 
                           
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       T 
                       
                         ij 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         1 
                       
                       X 
                     
                     = 
                     
                       
                         T 
                         ij 
                         X 
                       
                       + 
                       
                         
                           x 
                           i 
                           LB 
                         
                         · 
                         
                           W 
                           ij 
                           X 
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       T 
                       
                         ij 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         1 
                       
                       U 
                     
                     = 
                     
                       
                         T 
                         ij 
                         U 
                       
                       + 
                       
                         max 
                         ⁢ 
                         
                           
                             { 
                             
                               
                                 u 
                                 j 
                                 LB 
                               
                               , 
                               
                                 
                                   ∑ 
                                   
                                     k 
                                     = 
                                     1 
                                   
                                   
                                     K 
                                     i 
                                   
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   
                                     z 
                                     ik 
                                     0 
                                   
                                   · 
                                   
                                     u 
                                     jk 
                                     LB 
                                   
                                 
                               
                             
                             } 
                           
                           · 
                           
                             W 
                             ij 
                             U 
                           
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     Let 
     
       
         
           
             
               
                 
                   
                     v 
                     ⁡ 
                     
                       ( 
                       X 
                       ) 
                     
                   
                   = 
                   
                     
                       x 
                       
                         ij 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         1 
                       
                     
                     - 
                     
                       
                         R 
                         ij 
                       
                       · 
                       
                         ( 
                         
                           
                             
                               A 
                               ij 
                               0 
                             
                             · 
                             
                               u 
                               
                                 i 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           + 
                           
                             
                               B 
                               ij 
                               0 
                             
                             · 
                             
                               u 
                               
                                 j 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           + 
                           
                             
                               C 
                               ij 
                               0 
                             
                             · 
                             
                               u 
                               
                                 ij 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                           + 
                           
                             
                               D 
                               ij 
                               0 
                             
                             · 
                             
                               u 
                               
                                 ji 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 1 
                               
                             
                           
                         
                         ) 
                       
                     
                     - 
                     
                       
                         b 
                         v 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
             
               
                 
                   
                     o 
                     ⁡ 
                     
                       ( 
                       
                         
                           x 
                           ij 
                         
                         , 
                         
                           x 
                           i 
                         
                       
                       ) 
                     
                   
                   = 
                   
                     
                       x 
                       
                         i 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         1 
                       
                     
                     - 
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         N 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           R 
                           ij 
                         
                         · 
                         
                           x 
                           
                             ij 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             1 
                           
                         
                       
                     
                     - 
                     
                       
                         b 
                         x 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
           
         
       
     
     Then, by introducing logarithmic barrier terms, one can obtain the following Lagrangian function:
 
 L   μ ( X,S )=ƒ( X )−μ T (ln  X +ln  S )− Y   T ( X+S−X   UB )− y   1   T   o ( x   ij   ,x   i )− y   2   T   v ( X ))  (19)
 
where X=({x i1 }, {x ij1 }, {u i1 }, {u ij1 }) T , S=({s xi }, {s xij }, {s ui }, {s uij }) T . The Lagrange multipliers Y and y are called dual variables. According to an embodiment of the invention, the Karush-Kuhn-Tucker (KKT) optimality conditions of the barrier system may be obtained from the gradient of the Lagrangian with respect to all variables, and solving for the extremum:
 
                     [           ∇     L   X                 ∇     L   S                 ∇     L   Y                 ∇     L     y   1                   ∇     L     y   2               ]     =       ⁢       [             ∇     f   ⁡     (   X   )         -     μ   /   X     -   Y   -       y   1   T     ⁢     ∇     o   ⁡     (       x   ij     ,     x   i       )           -       y   2   T     ⁢     ∇     v   ⁡     (   X   )                           -   μ     /   S     -   Y               -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -     v   ⁡     (   X   )               ]     ⁢           =   0.             (   20   )               
Introducing variables Z x =μ/X, Z s =μ/S transforms EQ. (20) into:
 
                       [             ∇     f   ⁡     (     X   ,   Y     )         -     Z   x     -   Y   -       y   1   T     ⁢     ∇     o   ⁡     (       x   ij     ,     x   i       )           -       y   2   T     ⁢     ∇     v   ⁡     (   X   )                         -   Y     -     Z   x                 -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -     v   ⁡     (   X   )                       Z   x     ·   X     -     μ   ⁢           ⁢   e                     Z   s     ·   S     -     μ   ⁢           ⁢   e             ]     =   0     ,           (   21   )               
where e is a unit vector. Newton&#39;s method can be used to solve the above non-linear equation, to yield the following linear equation:
 
                     H   ·     [         ΔX           ΔS           ΔY             Δy   1               Δy   2               ΔZ   x               ΔZ   s           ]       =     -     [             ∇     f   ⁡     (     X   ,   Y     )         -     Z   x     -   Y   -       y   1     ·     ∇     o   ⁡     (       x   ij     ,     x   i       )           -       y   2     ·     ∇     v   ⁡     (   X   )                         -   Y     -     Z   s                 -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -     v   ⁡     (   X   )                       Z   S     ·   X     -     μ   ⁢           ⁢   e                     Z   S     ·   S     -     μ   ⁢           ⁢   e             ]               (   22   )               
in which H is the following Jacobian matrix:
 
                     H   =     [           H   11         0         -   I           H   41   T           H   51   T           -   I         0           0       0         -   I         0       0       0         -   I               -   I           -   I         0       0       0       0       0             H   41         0       0       0       0       0       0             H   51         0       0       0       0       0       0             diag   ⁡     (     Z   X     )           0       0       0       0         diag   ⁡     (   X   )           0           0         diag   ⁡     (     Z   S     )           0       0       0       0         diag   ⁡     (   S   )             ]       ,           (   23   )               
where
 
 H   11 =diag( V   X   ,W   X   v   U   ,W   U ),  (24.1)
 
where diag( ) is a vector formed from the diagonal components of the matrix argument,
 
                           ⁢         H   41     =     [       -   I     ,   R   ,   0   ,   0     ]       ,             (   24.2   )               R   =     (           [       R   11     ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R       1   ⁢     N   1       ⁢               ⁢   I     ]         0       …       0           0         [       R   21     ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R     2   ⁢     N   2         ⁢   I     ]         …       0           0       ⋮       ⋱       ⋮           0       …       0         [       R       N   1     ⁢   1       ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R       N   1     ⁢     N   1         ⁢   I     ]           )             (   24.3   )               
and the row of branch (i, j) of H 51  is
 
( H   51 ) ij =[0 ,−I, . . . R   ij   A   ij   0   , . . . R   ij   B   ij   0   , . . . R   ij   C   ij   0   , . . . R   ij   D   ij   0 , . . . 0].  (24.4)
 
Note that, in EQ. (24.1), the arguments V X , W X , V U , W U  themselves form a block diagonal
 
                           ⁢       matrix   ⁢         :     ⁢           [           V   X         0       0       0           0         W   X         0       0           0       0         V   U         0           0       0       0         W   U           ]     .     
     ⁢           ⁢     Define   ⁢     
     ⁢           [           r   X               r   S               r   Y               r     y   1                 r     y   2                 r     Z   x                 r     Z   s             ]         =       [             ∇     f   ⁡     (   X   )         -     Z   x     -   Y   -       y   1   T     ⁢     ∇           ⁢     o   ⁡     (       x   ij     ,     x   i       )           -       y   2   T     ⁢     ∇           ⁢     v   ⁡     (   X   )                         -   Y     -     Z   s                 -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -           ⁢     v   ⁡     (   X   )                       Z   x     ·   X     -     μ   ⁢           ⁢   e                     Z   s     ·   S     -     μ   ⁢           ⁢   e             ]     .               (   25   )                       ⁢     Then   ,                                 [           H   11         0         -   I           H   41   T           H   51   T           -   I         0           0       0         -   I         0       0       0         -   I               -   I           -   I         0       0       0       0       0             H   41         0       0       0       0       0       0             H   51         0       0       0       0       0       0             diag   ⁡     (     Z   X     )           0       0       0       0         diag   ⁡     (   X   )           0           0         diag   ⁡     (     Z   S     )           0       0       0       0         diag   ⁡     (   S   )             ]     ⁡     [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   S               Δ   ⁢           ⁢   Y               Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2                 Δ   ⁢           ⁢     Z   x                 Δ   ⁢           ⁢     Z   s             ]       =     -           ⁢       [           r   X               r   S               r   Y               r     y   1                 r     y   2                 r     Z   x                 r     Z   s             ]     .               (   26   )               
Note that the dimension of the system of EQS. (26) can be very large. As a result the solution time needed to solve it may be prohibitive for large scale electrical networks. However, a procedure according to an embodiment of the invention that uses a series of back substitutions, described below, can reduce the dimensionality of the system so that the equations require much less computational time to solve.
 
     According to an embodiment of the invention, to solve these Newton equations, the dimension can be reduced as follows. From the last two equations, one can obtain:
 
Δ Z   x =diag( X ) −1 (− r   Z     x   −diag( Z   X )Δ X )
 
Δ Z   s =diag( S ) −1 (− r   Z     s   −diag( Z   S )Δ S )  (27)
 
The solutions of EQS. (27) can be substituted into the other equations of EQS. (26), which can be reorganized into primal and dual parts:
 
                     [           H   41         0       ❘       0       0       0             H   51         0       ❘       0       0       0             -   I           -   I         ❘       0       0       0           -       -                   -       -       -               H   11     +         diag   ⁡     (   X   )         -   1       ⁢     diag   ⁡     (     Z   x     )               0       ❘         -   I           H   41   T           H   51   T             0             diag   ⁡     (   S   )         -   1       ⁢     diag   ⁡     (     Z   s     )             ❘         -   I         0       0         ]     ⁢             [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   S             -             Δ   ⁢           ⁢   Y               Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2             ]     =           [           -     r     y   ⁢           ⁢   1                   -     r     y   ⁢           ⁢   2                   -     r   Y               -               -     r   x       -         diag   ⁡     (   X   )         -   1       ⁢     r     Z   x                       -     r   s       -         diag   ⁡     (   S   )         -   1       ⁢     r     Z   s                 ]                   (   28   )               
By defining
 
                         H   ~     11     =       H   11     +         diag   ⁡     (   X   )         -   1       ⁢     diag   ⁡     (     Z   x     )             ⁢     
     ⁢         H   ~     22     =         diag   ⁡     (   S   )         -   1       ⁢     diag   ⁡     (     Z   S     )           ⁢     
     ⁢         r   ~     x     =       r   x     +         diag   ⁡     (   X   )         -   1       ⁢     r     Z   x             ⁢     
     ⁢         r   ~     s     =       r   s     +         diag   ⁡     (   S   )         -   1       ⁢     r     Z   s             ⁢     
     ⁢       Δ   ⁢           ⁢   X     =     [           {     Δ   ⁢           ⁢     x   i       }               {     Δ   ⁢           ⁢     x   ij       }               {     Δ   ⁢           ⁢     u   i       }               {     Δ   ⁢           ⁢     u   ij       }           ]       ⁢           ⁢       Δ   ⁢           ⁢     y   1       =     {     Δ   ⁢           ⁢     y   1     x   i         }       ⁢     
     ⁢       Δ   ⁢           ⁢     y   2       =     {     Δ   ⁢           ⁢     y   2     x   ij         }               (   29   )               
one can obtain, according to an embodiment of the invention:
 
                         [           H   41         0             H   51         0             -   I           -   I           ]     ⁡     [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   S           ]       =     [           -     r     y   ⁢           ⁢   1                   -     r     y   ⁢           ⁢   2                   -     r   Y             ]       ,           (   30   )             and                               [             H   ~     11         0           0           H   ~     22           ]     ⁡     [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   S           ]       +       [           -   I           H   41   T           H   51   T               -   I         0       0         ]     ⁡     [           Δ   ⁢           ⁢   Y               Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2             ]         =       [           -       r   ~     x                 -       r   ~     s             ]     .             (   31   )               
More specifically, EQS. (30) and (31) can be expanded using the definitions of EQS. (24.1)-(24.4) and the model definitions of EQS. (14.1)-(14.5) and (15.1)-(15.5). Thus, EQ. (30) has the following expanded form for each bus i and branch(i, j):
 
                         Δ   ⁢           ⁢     x   i       -       ∑   j             ⁢           ⁢       R   ij     ⁢   Δ   ⁢           ⁢     x   ij           =     r     y   1       x   i         ⁢     
     ⁢         Δ   ⁢           ⁢     x   ij       -       R   ij     ⁡     (         A   ij     ⁢   Δ   ⁢           ⁢     u   i       +       B   ij     ⁢   Δ   ⁢           ⁢     u   j       +       C   ij     ⁢   Δ   ⁢           ⁢     u   ij       +       D   ij     ⁢   Δ   ⁢           ⁢     u   ji         )         =     r     y   2       x   ij         ⁢     
     ⁢         Δ   ⁢           ⁢   X     +     Δ   ⁢           ⁢   S       =     r   Y               (   32   )               
Similarly, EQ. (31) has the following expanded form:
 
     
       
         
           
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         
                           
                             
                               V 
                               ~ 
                             
                             i 
                             X 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             x 
                             i 
                           
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             Y 
                             
                               x 
                               i 
                             
                           
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             y 
                             1 
                             
                               x 
                               i 
                             
                           
                         
                       
                       = 
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           
                             x 
                             i 
                           
                         
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         
                           
                             
                               W 
                               ~ 
                             
                             ij 
                             X 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             x 
                             ij 
                           
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             Y 
                             
                               x 
                               ij 
                             
                           
                         
                         + 
                         
                           
                             R 
                             ij 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             y 
                             1 
                             
                               x 
                               i 
                             
                           
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             y 
                             2 
                             
                               x 
                               ij 
                             
                           
                         
                       
                       = 
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           
                             x 
                             ij 
                           
                         
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                       
                         
                           
                             
                               V 
                               ~ 
                             
                             i 
                             U 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             u 
                             i 
                           
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             Y 
                             
                               u 
                               i 
                             
                           
                         
                         + 
                         
                           
                             ∑ 
                             j 
                             
                                 
                             
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             
                               R 
                               ij 
                             
                             ⁢ 
                             
                               A 
                               ij 
                               T 
                             
                             ⁢ 
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               y 
                               2 
                               
                                 x 
                                 ij 
                               
                             
                           
                         
                         + 
                         
                           
                             ∑ 
                             j 
                             
                                 
                             
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             
                               R 
                               ji 
                             
                             ⁢ 
                             
                               B 
                               ji 
                               T 
                             
                             ⁢ 
                             Δ 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               y 
                               2 
                               
                                 x 
                                 ji 
                               
                             
                           
                         
                       
                       = 
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           
                             u 
                             i 
                           
                         
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         
                           
                             
                               W 
                               ~ 
                             
                             ij 
                             U 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             u 
                             ij 
                           
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             Y 
                             
                               u 
                               ij 
                             
                           
                         
                         + 
                         
                           
                             R 
                             ij 
                           
                           ⁢ 
                           
                             C 
                             ij 
                             T 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             y 
                             2 
                             
                               x 
                               ij 
                             
                           
                         
                         + 
                         
                           
                             R 
                             ji 
                           
                           ⁢ 
                           
                             D 
                             ji 
                             T 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             y 
                             2 
                             
                               x 
                               ji 
                             
                           
                         
                       
                       = 
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           
                             u 
                             ij 
                           
                         
                       
                     
                     ⁢ 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         
                           
                             
                               H 
                               ~ 
                             
                             22 
                           
                           ⁢ 
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           S 
                         
                         - 
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           Y 
                         
                       
                       = 
                       
                         - 
                         
                           
                             r 
                             ~ 
                           
                           S 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   33 
                   ) 
                 
               
             
           
         
       
     
     According to an embodiment of the invention, by first solving for control variables ΔU=[{Δu i }, {Δu ij }] T  in EQS. (33), the other variables in EQS. (32) and (33) can be solved using the following sequence of substitutions:
 
{Δ u   i },{Δu ij   }→{Δx   ij   }→{Δx   i   }→ΔS→ΔY→Δy   1   →Δy   2 .  (34)
 
       FIG. 2  is a flowchart of a method according to an embodiment of the invention for estimating the optimal power flow of a smart electric power grid. Referring now to the figure, a method according to an embodiment of the invention begins at step  21  by providing a cost function that models a smart electric power grid, such as the cost functions of EQS. (7) and (13), above, and EQ. (64), below. Constraints for the bus and branch control vectors, and the bus and branch dependent variables vectors, are provided at step  22 , based on physical and operational requirements and actual measurements. At step  23 , a set of linear equations is derived that minimize the cost function subject to the constraints. The dimension of the linear equations can be reduced at step  24  by solving for a subset of the linear equations. The reduced dimension linear equations are re-organized into primal and dual parts at step  25 , and the re-organized reduced dimension linear equations are de-composed at step  26  into two systems of block matrix equations that can solved. A solution of the two systems of block matrix equations yields conditions for a lowest cost per kilowatthour delivered through the smart electric power grid. 
     Coding Details 
     According to an embodiment of the invention, a preconditioned conjugate gradient method can be used to solve the above equations, which uses a Jacobi preconditioning matrix composed of the diagonal elements of a matrix M of the linear system Mx=b. For every Δu i , the diagonal parameter can be computed as: 
               M     u   i       =         V   ~     i   U     +   I   +       ∑   j             ⁢           ⁢       R   ij     ·     A   ij   T     ·     [         (         W   ~     ij   X     +       H   ~     22     x   ij         )     ⁢     A   ij       +       (         V   ~     i   X     +       H   ~     22     x   i         )     ⁢     (     ∑       R   ij     ⁢     A   ij         )         ]         +     ∑         R   ji     ·       B   ji   T     ⁡     (         W   ~     ji   X     +       H   ~     22     x   ji       +       V   ~     j   X     +       H   ~     22     x   j         )         ⁢     B   ji                 
Similarly, for every Δu ij , the diagonal parameter can be computed as:
 
     
       
         
           
             
               M 
               
                 u 
                 ij 
               
             
             = 
             
               
                 
                   W 
                   ~ 
                 
                 ij 
                 U 
               
               + 
               
                 
                   H 
                   ~ 
                 
                 22 
                 
                   u 
                   ij 
                 
               
               + 
               
                 
                   
                     R 
                     ij 
                   
                   · 
                   
                     
                       C 
                       ij 
                       T 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             W 
                             ~ 
                           
                           ij 
                           X 
                         
                         + 
                         
                           
                             H 
                             ~ 
                           
                           22 
                           
                             x 
                             ij 
                           
                         
                         + 
                         
                           
                             V 
                             ~ 
                           
                           i 
                           X 
                         
                         + 
                         
                           
                             H 
                             ~ 
                           
                           22 
                           
                             x 
                             i 
                           
                         
                       
                       ) 
                     
                   
                 
                 ⁢ 
                 
                   C 
                   ij 
                 
               
               + 
               
                 
                   
                     R 
                     ji 
                   
                   · 
                   
                     
                       D 
                       ji 
                       T 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           
                             W 
                             ~ 
                           
                           ji 
                           X 
                         
                         + 
                         
                           
                             H 
                             ~ 
                           
                           22 
                           
                             x 
                             ji 
                           
                         
                         + 
                         
                           
                             V 
                             ~ 
                           
                           j 
                           X 
                         
                         + 
                         
                           
                             H 
                             ~ 
                           
                           22 
                           
                             x 
                             j 
                           
                         
                       
                       ) 
                     
                   
                 
                 ⁢ 
                 
                   D 
                   ji 
                   T 
                 
               
             
           
         
       
     
     According to an embodiment of the invention, the quadratic objective function ƒ(X) is written as ƒ(X)=½X T QX+p T X, which has the following associated dual program and duality gap: 
               max   ⁢           ⁢     P   T     ⁢     Q     -   1       ⁢     P   T       +       y   1   T     ⁢     b   x       +       y   2   T     ⁢     b   v       +       Y   T     ⁢     X   UB                   s   .   t   .           ⁢     {               Z   S     +   Y     ≤   0                   Q     -   1       ⁢   P     ≥   0                 Z   S     ≥   0                 Z   X     ≥   0                   
where
 
 P=p−Z   x   −Y−y   1   T   ∇o ( X )− y   2   T   ∇v ( X ).
 
According to an embodiment of the invention, the current duality gap is:
 
ƒ( X )− g ( Y,u   1   ,y   2   ,Z   S   ,Z   X )= Y ·(− r   Y )+ y   1 (− r   y     1   )+ y   2 ·(− r   y     2   )+½ r   X   T   Q   −1   r   X   T +( Z   X   ·X+Z   S   ·S )− S ·( Z   S   +Y ).
 
In this expression of the duality gap, the first four terms come from the primal infeasibility for the current primal variable (X, S), and the last two terms come from the complementary slackness. When an optimal primal and dual solution is reached, the duality gap is 2nμ, where n is the dimension.
 
Interior Point Method Solver Math for a Discrete AC OPF Model
 
     First, a model according to an embodiment of the invention is standardized to the following format: 
                           ⁢       X   +     S   X     -     X   UB       =     0   ⁢           ⁢     {               x   i     +     s     x   i         =     x   i   UB                     x   ij     +     s     x   ij         =     x   ij   UB                           (   35.1   )                       ⁢       MU   +   NZ   +     S   U       =     b   ⁢           ⁢     {                 ∑   Ki     ⁢     u   i   k       -     s   ui     1   ,   k         =     u   i   LB                       ∑   Ki     ⁢     u   i   k       +     s   ui     2   ,   k         =     u   i   UB                       ∑   Ki     ⁢     u   ij   k       -     s   uij     1   ,   k         =     u   ij   LB                       ∑   Ki     ⁢     u   ij   k       +     s   uij     1   ,   k         =     u   ij   UB                       ∑   Ki     ⁢     z   i   k       -     s   zi     1   ,   k         =     z   i   LB                       ∑   Ki     ⁢     z   i   k       +     s   zi     2   ,   k         =     z   i   UB                       ∑   Ki     ⁢     z   ij   k       -     s   zij     1   ,   k         =     z   ij   LB                       ∑   Ki     ⁢     z   ij   k       +     s   zij     2   ,   k         =     z   ij   UB                     -     u   i   k       +       z   i   k     ·     u   ik   LB       +     s   uzi     1   ,   k         =   0                   u   i   k     -       z   i   k     ·     u   ik   UB       +     s   uzi     2   ,   k         =   0                   -     u   ij   k       +       z   ij   k     ·     u   ijk   LB       +     s   uzij     1   ,   k         =   0                   u   i   k     -       z   ij   k     ·     u   ijk   UB       +     s   uzij     2   ,   k         =   0                         (   35.2   )                       ⁢           α   ⁡     (   X   )       ⁢     :     ⁢           ⁢     x   i       -       ∑   j             ⁢       R   ij     ⁢     x   ij         -     b   v       =   0             (   35.3   )                     β   ⁡     (     X   ,   U     )       ⁢     :     ⁢           ⁢     x   ij       -       R   ij     ·     (         ∑     K   i               ⁢       A   ij   k     ·     u   i   k         +       ∑     K   i               ⁢       B   ij   k     ·     u   j   k         +       ∑     K   ij               ⁢       C   ij   k     ·     u   ij   k         +       ∑     K   ji               ⁢       D   ij   k     ·     u   ji   k           )         =   0           (   35.4   )                       ⁢           γ   ⁡     (     X   ,   U   ,   Z     )       ⁢     :     ⁢           ⁢   FX     +   GU   +   HZ   +   J     =   0             (   35.5   )                       ⁢     (     Cutting   ⁢           ⁢   planes     )                             
where the variables X, U, Z, M, N, and b are defined as follows:
 
                         X   =     [           {     x   i     }               {     x   ij     }           ]                             U   =     [           {     u   i   k     }               {     x   ij   k     }           ]                             Z   =     [           {     z   i   k     }               {     z   ij   k     }           ]                               (   36.1   )               M   =     [           -     E   1           0             E   1         0           0         -     E   1               0         E   1             0       0           0       0           0       0           0       0             -   I         0           I       0           0         -   I             0       I         ]             (   36.2   )               N   =     [         0       0           0       0           0       0           0       0             -     E   1           0             E   1         0           0         -     E   1               0         E   1               I   1         0             -     I   2           0           0         I   3             0         I   4           ]                             b   =     [           {     -     u   i   LB       }               {     u   i   UB     }               {     -     u   ij   LB       }               {     u   ij   UB     }               {     -     z   i   LB       }               {     z   i   UB     }               {     -     z   ij   LB       }               {     z   ij   UB     }             0           0           0           0         ]                             
and where
 
     
       
         
           
             
               
                 
                   
                     E 
                     1 
                   
                   = 
                   
                     [ 
                     
                       
                         
                           
                             ⅇ 
                             T 
                           
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           
                             ⅇ 
                             T 
                           
                         
                         
                           … 
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           ⋱ 
                         
                         
                           ⋮ 
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           … 
                         
                         
                           
                             ⅇ 
                             T 
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   36.3 
                   ) 
                 
               
             
             
               
                 
                   
                     I 
                     1 
                   
                   = 
                   
                     diag 
                     ⁢ 
                     
                       { 
                       
                         u 
                         ik 
                         LB 
                       
                       } 
                     
                   
                 
               
               
                 
                   ( 
                   36.4 
                   ) 
                 
               
             
             
               
                 
                   
                     I 
                     2 
                   
                   = 
                   
                     diag 
                     ⁢ 
                     
                       { 
                       
                         u 
                         ik 
                         UB 
                       
                       } 
                     
                   
                 
               
               
                 
                     
                 
               
             
             
               
                 
                   
                     I 
                     3 
                   
                   = 
                   
                     diag 
                     ⁢ 
                     
                       { 
                       
                         u 
                         ijk 
                         LB 
                       
                       } 
                     
                   
                 
               
               
                 
                     
                 
               
             
             
               
                 
                   
                     I 
                     4 
                   
                   = 
                   
                     diag 
                     ⁢ 
                     
                       { 
                       
                         u 
                         ijk 
                         uB 
                       
                       } 
                     
                   
                 
               
               
                 
                     
                 
               
             
           
         
       
     
     Following the above notation, a Lagrangian function according to an embodiment of the invention can be formulated as follows:
 
 L =ƒ( X,U )+μ T (ln  X +ln  U +ln  Z +ln  S   X +ln  S   U )+λ 1   T ( X+S   X   −X   UB )+λ 2   T ( MU+NZ+S   U   −b )+ y   1   T α( X )+ y   2   T β( X,U )+ y   3   T γ( X,U,Z )  (37)
 
The KKT first order optimality conditions of the resulting program are:
 
                     [           ∇     L   X                 ∇     L   U                 ∇     L   Z                 ∇     L     S   X                   ∇     L     S   U                   ∇     L     λ   1                   ∇     L     λ   2                   ∇     L     y   1                   ∇     L     y   2                   ∇     L     y   3               ]     =       ⁢                   [               f   X     ⁡     (     X   ,   U     )       +     μ   /   X     +     λ   1   T     +       y   1   T     ⁢       α   X     ⁡     (   X   )         +       y   2   T     ⁢       β   X     ⁡     (     X   ,   U     )         +                 y   3   T     ⁢       γ   X     ⁡     (     X   ,   U   ,   Z     )                       f   U     ⁡     (     X   ,   U     )       +     μ   /   U     +       y   2   T     ⁢       β   U     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   U     ⁡     (     X   ,   U   ,   Z     )                       μ   /   Z     +       λ   2   T     ⁢   N     +       y   3   T     ⁢       γ   Z     ⁡     (     X   ,   U   ,   Z     )                       μ   /     S   X       +     λ   1   T                   μ   /     S   U       +     λ   2   T                 X   +     S   X     -     X   UB                 MU   +   NZ   +     S   U     -   b               α   ⁡     (   X   )                 β   ⁡     (     X   ,   U     )                 γ   ⁡     (     X   ,   U   ,   Z     )             ]     ⁢           =   0                     (   38   )               
By introducing
 
                   T   =       [           T   X               T   U               T   Z               T     S   X                 T     S   U             ]     =     [           μ   /   X               μ   /   U               μ   /   Z               μ   /     S   X                 μ   /     S   U             ]               (   39   )               
The optimality conditions can be rewritten as
 
                     [                     f   X     ⁡     (     X   ,   U     )       +     T   X     +     λ   1   T     +       y   1   T     ⁢     α   X     ⁢     (   X   )       +                   y   2   T     ⁢       β   X     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   X     ⁡     (     X   ,   U   ,   Z     )                               f   U     ⁡     (     X   ,   U     )       +     T   U     +       y   2   T     ⁢       β   U     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   U     ⁡     (     X   ,   U   ,   Z     )                       T   Z     +       λ   2   T     ⁢   N     +       y   3   T     ⁢       γ   Z     ⁡     (     X   ,   U   ,   Z     )                       T     S   X       +     λ   1   T                   T     S   U       +     λ   2   T                 X   +     S   X     -     X   UB                 MU   +   NZ   +     S   U     -   b               α   ⁡     (   X   )                 β   ⁡     (     X   ,   U     )                 γ   ⁡     (     X   ,   U   ,   Z     )                   T   ·     [         X           U           Z             S   X               S   U           ]       -     μ   ⁢           ⁢   e             ]     =   0.           (   40   )               
Newton&#39;s method can be used to solve the above non-linear equation, to yield the following linear equation:
 
                     H   ·     [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   U               Δ   ⁢           ⁢   Z               Δ   ⁢           ⁢     S   X                 Δ   ⁢           ⁢     S   U                 Δ   ⁢           ⁢     λ   1                 Δ   ⁢           ⁢     λ   2                 Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2                 Δ   ⁢           ⁢     y   3                 Δ   ⁢           ⁢   T           ]       =     -     [                     f   X     ⁡     (     X   ,   U     )       +     T   X     +     λ   1   T     +       y   1   T     ⁢     α   X     ⁢     (   X   )       +                   y   2   T     ⁢       β   X     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   X     ⁡     (     X   ,   U   ,   Z     )                               f   U     ⁡     (     X   ,   U     )       +     T   U     +       y   2   T     ⁢       β   U     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   U     ⁡     (     X   ,   U   ,   Z     )                       T   Z     +       λ   2   T     ⁢   N     +       y   3   T     ⁢       γ   Z     ⁡     (     X   ,   U   ,   Z     )                       T     S   X       +     λ   1   T                   T     S   U       +     λ   2   T                 X   +     S   X     -     X   UB                 MU   +   NZ   +     S   U     -   b               α   ⁡     (   X   )                 β   ⁡     (     X   ,   U     )                 γ   ⁡     (     X   ,   U   ,   Z     )                   T   ·     [         X           U           Z             S   X               S   U           ]       -     μ   ⁢           ⁢   e             ]               (   41   )               
where H is the following Jacobian matrix:
 
                   H   =     [           H     1   ,   1           0       0       0       0       I       0         H     1   ,   8             H     1   ,   9             H     1   ,   10             E   1             0         H     2   ,   2           0       0       0       0         H     2   ,   7           0         H     2   ,   9             H     2   ,   10             E   2             0       0       0       0       0       0         H     3   ,   7           0       0         H     3   ,   10             E   3             0       0       0       0       0       I       0       0       0       0         E   4             0       0       0       0       0       0       I       0       0       0         E   5             I       0       0       I       0       0       0       0       0       0       0           0         H     2   ,   7             H     3   ,   7           0       I       0       0       0       0       0       0             H     1   ,   8           0       0       0       0       0       0       0       0       0       0             H     1   ,   9             H     2   ,   9           0       0       0       0       0       0       0       0       0             H     1   ,   10             H     2   ,   10             H     3   ,   10           0       0       0       0       0       0       0       0             T   1           T   2           T   3           T   4           T   5         0       0       0       0       0         T   6           ]             (   42   )               
where
 
E i  is the i th  block identity matrix where
 
                           ⁢       I   =       ∑     i   =   1     5     ⁢     E   i         ;             (   43.1   )                   T   i     =       E   i     ⁡     [           diag   ⁢     {     T   X     }                                                                             diag   ⁢     {     T   U     }                                                                             diag   ⁢     {     T   Z     }                                                                             diag   ⁢     {     T     S   X       }                                                                             diag   ⁢     {     T     S   U       }             ]         ⁢     
     ⁢           ⁢         for   ⁢           ⁢   i     =   1     ,   2   ,   3   ,   4   ,     5   ;               (   43.2   )                       ⁢         T   6     =     [           diag   ⁢     {   X   }                                                                             diag   ⁢     {   U   }                                                                             diag   ⁢     {   Z   }                                                                             diag   ⁢     {     S   X     }                                                                             diag   ⁢     {     S   U     }             ]       ;             (   43.3   )                       ⁢         H     1   ,   1       =     [       V   x     |     W   x       ]       ;             (   43.4   )                       ⁢         H     1   ,   8       =       [     I   |     -   R       ]     T       ;             (   43.5   )                 R   =     (           [       R   11     ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R     1   ⁢     N   1         ⁢   I     ]         0       …       0           0         [       R   21     ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R     2   ⁢     N   2         ⁢   I     ]         …       0           0       ⋮       ⋱       ⋮           0       …       0         [       R       N   1     ⁢   1       ⁢   I   ⁢           ⁢   …   ⁢           ⁢     R       N   1     ⁢     N   1         ⁢   I     ]           ⁢           )       ⁢           ;           (   43.6   )                       ⁢         H   1.9     =       [     0   |   I     ]     T       ;             (   43.7   )                       ⁢         H   1.10     =     F   T       ;             (   43.8   )                       ⁢         H   2.2     =     [       V   U     |     W   U       ]       ;             (   43.9   )                       ⁢         H   2.7     =     M   T       ;             (   43.10   )                       ⁢         H   2.9     =       [     A   |   B     ]     T       ,             (   43.11   )               
where A ij , B ij  that correspond to the branch (i,j) are defined as
 
 A   ij =[0 , . . . A   ij   1   , . . . A   ij   K     i   , 0 , . . . B   ij   1   , . . . B   ij   K     j   ,0, . . . ],  (43.12)
 
 B   ij =[0 , . . . C   ij   1   , . . . C   ij   K     ij   ,0 , . . . D   ij   1   , . . . D   ij   K     ji   ,0, . . . ],  (43.13)
 
 H   2.10   =N   T ;  (43.14)
 
 H   3.10   =H   T .  (43.15)
 
By defining
 
                     [           r   X               r   U               r   Z               r     S   X                 r     S   U                 r     λ   1                 r     λ   2                 r     y   1                 r     y   2                 r     y   3                 r   T           ]     =           ⁢           [                     f   X     ⁡     (     X   ,   U     )       +     T   X     +     λ   1   T     +       y   1   T     ⁢     α   X     ⁢     (   X   )       +                   y   2   T     ⁢       β   X     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   X     ⁡     (     X   ,   U   ,   Z     )                               f   U     ⁡     (     X   ,   U     )       +     T   U     +       y   2   T     ⁢       β   U     ⁡     (     X   ,   U     )         +       y   3   T     ⁢       γ   U     ⁡     (     X   ,   U   ,   Z     )                       T   Z     +       λ   2   T     ⁢   N     +       y   3   T     ⁢       γ   Z     ⁡     (     X   ,   U   ,   Z     )                       T     S   X       +     λ   1   T                   T     S   U       +     λ   2   T                 X   +     S   X     -     X   UB                 MU   +   NZ   +     S   U     -   b               α   ⁡     (   X   )                 β   ⁡     (     X   ,   U     )                 γ   ⁡     (     X   ,   U   ,   Z     )                   T   ·     [         X           U           Z             S   X               S   U           ]       -     μ   ⁢           ⁢   e             ⁢           ]               (   44   )               
one has
 
     
       
         
           
             
               
                 
                   
                     [ 
                     
                       
                         
                           
                             H 
                             
                               1 
                               , 
                               1 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           I 
                         
                         
                           0 
                         
                         
                           
                             H 
                             
                               1 
                               , 
                               8 
                             
                           
                         
                         
                           
                             H 
                             
                               1 
                               , 
                               9 
                             
                           
                         
                         
                           
                             H 
                             
                               1 
                               , 
                               10 
                             
                           
                         
                         
                           
                             E 
                             1 
                           
                         
                       
                       
                         
                           0 
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               2 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               7 
                             
                           
                         
                         
                           0 
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               9 
                             
                           
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               10 
                             
                           
                         
                         
                           
                             E 
                             2 
                           
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           
                             H 
                             
                               3 
                               , 
                               7 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           
                             H 
                             
                               3 
                               , 
                               10 
                             
                           
                         
                         
                           
                             E 
                             3 
                           
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           I 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           
                             E 
                             4 
                           
                         
                       
                       
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           I 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           
                             E 
                             5 
                           
                         
                       
                       
                         
                           I 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           I 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                       
                         
                           0 
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               7 
                             
                           
                         
                         
                           
                             H 
                             
                               3 
                               , 
                               7 
                             
                           
                         
                         
                           0 
                         
                         
                           I 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                       
                         
                           
                             H 
                             
                               1 
                               , 
                               8 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                       
                         
                           
                             H 
                             
                               1 
                               , 
                               9 
                             
                           
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               9 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                       
                         
                           
                             H 
                             
                               1 
                               , 
                               10 
                             
                           
                         
                         
                           
                             H 
                             
                               2 
                               , 
                               10 
                             
                           
                         
                         
                           
                             H 
                             
                               3 
                               , 
                               10 
                             
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                       
                       
                         
                           
                             T 
                             1 
                           
                         
                         
                           
                             T 
                             2 
                           
                         
                         
                           
                             T 
                             3 
                           
                         
                         
                           
                             T 
                             4 
                           
                         
                         
                           
                             T 
                             5 
                           
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           0 
                         
                         
                           
                             T 
                             6 
                           
                         
                       
                     
                     ] 
                   
                   ⁢ 
                   
                       
                     
                       
                         [ 
                         
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 X 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 U 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 Z 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   S 
                                   X 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   S 
                                   U 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   λ 
                                   1 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   λ 
                                   2 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   y 
                                   1 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   y 
                                   2 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   y 
                                   3 
                                 
                               
                             
                           
                           
                             
                               
                                 Δ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 T 
                               
                             
                           
                         
                         ] 
                       
                       = 
                       
                         - 
                         
                           
                             [ 
                             
                               
                                 
                                   
                                     r 
                                     X 
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     U 
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     Z 
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       S 
                                       X 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       S 
                                       U 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       λ 
                                       1 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       λ 
                                       2 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       y 
                                       1 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       y 
                                       2 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     
                                       y 
                                       3 
                                     
                                   
                                 
                               
                               
                                 
                                   
                                     r 
                                     T 
                                   
                                 
                               
                             
                             ] 
                           
                           . 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   45 
                   ) 
                 
               
             
           
         
       
     
     Similar to EQS. (26), the dimension of the system of EQS. (45) can be very large. By decomposing r T  as follows, 
                       r   T     =     [           r   Tx               r     T   U                 r   Tz               r     Ts   X                 r     Ts   U             ]       ,           (   46   )               
the dimension of EQS. (45) can be reduced as follows. From the last block constraints there is:
 
                     [           Δ   ⁢           ⁢     T   X                 ΔT   U               Δ   ⁢           ⁢     T   Z                 Δ   ⁢           ⁢     T     S   X                   Δ   ⁢           ⁢     T     S   U               ]     =       [               diag     -   1       ⁡     (   X   )       ⁢     (       r     T   X       -       diag   ⁡     (     T   X     )       ⁢   Δ   ⁢           ⁢   X       )                     diag     -   1       ⁡     (   U   )       ⁢     (       r     T   U       -       diag   ⁡     (     T   U     )       ⁢   Δ   ⁢           ⁢   U       )                     diag     -   1       ⁡     (   Z   )       ⁢     (       r     T   Z       -       diag   ⁡     (     T   Z     )       ⁢   Δ   ⁢           ⁢   Z       )                     diag     -   1       ⁡     (     S   X     )       ⁢     (       r     T     S   X         -       diag   ⁡     (     T     S   X       )       ⁢   Δ   ⁢           ⁢     S   X         )                     diag     -   1       ⁡     (     S   U     )       ⁢     (       r     T     S   U         -       diag   ⁡     (     T     S   U       )       ⁢   Δ   ⁢           ⁢     S   U         )             ]     .             (   47   )               
This definition for the ΔTs can be substituted into the other blocks of the matrix on the left hand side of EQ. (45), which can be reorganized into primal and dual parts:
 
                     [         I       0       0       I       0       ❘       0       0       0       0       0       0           0         H     2   ,   7             H     3   ,   7           0       I       ❘       0       0       0       0       0       0             H     1   ,   8           0       0       0       0       ❘       0       0       0       0       0       0             H     1   ,   9             H     2   ,   9           0       0       0       ❘       0       0       0       0       0       0             H     1   ,   10             H     2   ,   10             H     3   ,   10           0       0       ❘       0       0       0       0       0       0           -       -       -       -       -       -       -       -       -       -       -                           H   ~       1   ,   1           0       0       0       0       ❘       I       0         H     1   ,   8             H     1   ,   9             H     1   ,   10             E   1             0           H   ~       2   ,   2           0       0       0       ❘       0         H     2   ,   7           0         H     2   ,   9             H     2   ,   10             E   2             0       0           H   ~       3   ,   3           0       0       ❘       0         H     3   ,   7           0       0         H     3   ,   10             E   3             0       0       0           H   ~       4   ,   4           0       ❘       I       0       0       0       0         E   4             0       0       0       0           H   ~       5   ,   5           ❘       0       I       0       0       0         E   5           ]     ⁢             [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   U               Δ   ⁢           ⁢   Z               Δ   ⁢           ⁢     S   X                 Δ   ⁢           ⁢     S   U                 -   ❘               Δ   ⁢           ⁢     λ   1                 Δ   ⁢           ⁢     λ   2                 Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2                 Δ   ⁢           ⁢     y   3             ]     =       -     [           r     λ   1                 r     λ   2                 r     y   1                 r     y   2                 r     y   3               -               r   ~     X                 r   ~     U                 r   ~     Z                 r   ~       S   X                   r   ~       S   U             ]       ⁢     
     ⁢           ⁢   where                 (   48   )                       ⁢         [             H   ~       1   ,   1                   H   ~       2   ,   2                   H   ~       3   ,   3                   H   ~       4   ,   4                   H   ~       5   ,   5             ]     =     [             H     1   ,   1       -         diag     -   1       ⁡     (   X   )       ⁢     diag   ⁡     (     T   X     )                       H     2   ,   2       -         diag     -   1       ⁡     (   U   )       ⁢     diag   ⁡     (     T   U     )                       -       diag     -   1       ⁡     (   Z   )         ⁢     diag   ⁡     (     T   Z     )                     -       diag     -   1       ⁡     (     S   X     )         ⁢     diag   ⁡     (     T     S   X       )                     -       diag     -   1       ⁡     (     S   U     )         ⁢     diag   ⁡     (     T     S   U       )               ]       ,     
     ⁢           ⁢   and             (   49   )                       ⁢       [             r   ~     X                 r   ~     U                 r   ~     Z                 r   ~       S   X                   r   ~       S   U             ]     =       [             r   X     +         diag     -   1       ⁡     (   X   )       ⁢     r     T   X                       r   U     +         diag     -   1       ⁡     (   U   )       ⁢     r     T   U                       r   Z     +         diag     -   1       ⁡     (   Z   )       ⁢     r     T   Z                       r     S   X       +         diag     -   1       ⁡     (     S   X     )       ⁢     r     S     T   X                         r     S   U       +         diag     -   1       ⁡     (     S   U     )       ⁢     r     S     T   U                   ]     .               (   50   )               
According to an embodiment of the invention, EQ. (45) can be reduced to
 
                           ⁢           [         I       0       0       I       0           0         H     2   ,   7             H     3   ,   7           0       I             H     1   ,   8           0       0       0       0             H     1   ,   9             H     2   ,   9           0       0       0             H     1   ,   10             H     2   ,   10             H     3   ,   10           0       0         ]     ⁡     [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   U               Δ   ⁢           ⁢   Z               Δ   ⁢           ⁢     S   X                 Δ   ⁢           ⁢     S   U             ]       =     -     [           r     λ   1                 r     λ   2                 r     y   1                 r     y   2                 r     y   3             ]         ⁢     
     ⁢           ⁢   and             (   51   )                 [             H   ~       1   ,   1           0       0       0       0           0           H   ~       2   ,   2           0       0       0           0       0           H   ~       3   ,   3           0       0           0       0       0           H   ~       4   ,   4           0           0       0       0       0           H   ~       5   ,   5             ]     ⁢               [           Δ   ⁢           ⁢   X               Δ   ⁢           ⁢   U               Δ   ⁢           ⁢   Z               Δ   ⁢           ⁢     S   X                 Δ   ⁢           ⁢     S   U             ]     +       [         I       0         H     1   ,   8             H     1   ,   9             H     1   ,   10             E   1             0         H     2   ,   7           0         H     2   ,   9             H     2   ,   10             E   2             0         H     3   ,   7           0       0         H     3   ,   10             E   3             I       0       0       0       0         E   4             0       I       0       0       0         E   5           ]     ⁡     [           Δ   ⁢           ⁢     λ   1                 Δ   ⁢           ⁢     λ   2                 Δ   ⁢           ⁢     y   1                 Δ   ⁢           ⁢     y   2                 Δ   ⁢           ⁢     y   3             ]         =     [             r   ~     X                 r   ~     U                 r   ~     Z                 r   ~       S   X                   r   ~       S   U             ]                 (   52   )               
EQS. (51) and (52) can be solved similarly to EQS. (30) and (31) above, by expanding the block components in the matrices of EQS. (51) and (52) using their definitions above, and carrying out a series of substitutions similar to that of EQ. (34).
 
Heuristic Block-Wise Newton Method
 
     It can be seen that a system according to an embodiment of the invention is highly structural in the way that all constraints can be grouped by buses and branches. Without loss of generality, a minimal AC OPF model according to an embodiment of the invention will be considered as an example, and the same heuristic can be applied to a more general discrete AC OPF model according to an embodiment of the invention. 
     For a minimal AC OPF model according to an embodiment of the invention, the following non-linear system will be solved: 
                     L   ~     =       [             ∇     f   ⁡     (   X   )         -     μ   /   X     -   Y   -       y   1   T     ⁢     ∇     o   ⁡     (       x   ij     ,     x   i       )           -       y   2   T     ⁢     ∇     v   ⁡     (   X   )                           -   μ     /   S     -   Y               -     (     X   +   S   -     X   UB       )                 -     o   ⁡     (       x   ij     ,     x   i       )                   -     v   ⁡     (   X   )               ]     =   0.             (   53   )               
According to an embodiment of the invention, the vector {tilde over (L)} can be re-organized and rewritten as follows:
 
                       L   ~     =     [           {     L   i     }               {     L   ij     }           ]       ,           (   54   )               
where each L i  denotes the constraints of bus i, thus is a function of x i  and {x ik }. Similarly each L ij  denotes the constraints of branch (i, j), which is a function of x i , x j , u i , u j , u ij , and u ji . Therefore, according to an embodiment of the invention, the system can be decomposed by different buses and branches. More specifically, instead of directly applying Newton method to {tilde over (L)}, it can be used for every bus L i  and branch L ij :
 
                           H   i     ·     [           Δ   ⁢           ⁢     X   i                 Δ   ⁢           ⁢     S   i                 Δ   ⁢           ⁢     Y   i                 Δ   ⁢           ⁢     y   1   i                 Δ   ⁢           ⁢     y   2   i             ]       =     -     L   i         ;         H   ij     ·     [           Δ   ⁢           ⁢     X   ij                 Δ   ⁢           ⁢     S   ij                 Δ   ⁢           ⁢     Y   ij                 Δ   ⁢           ⁢     y   1   ij                 Δ   ⁢           ⁢     y   2   ij             ]       =     -     L   ij           ,           (   55   )               
where the Jacobian matrices H i  and H ij  can be obtained by the above techniques. Moreover, cutting plane methods can be used to solve these block-wise Newton&#39;s equations since each cutting plane is a function of either L i  or L ij .
 
     One issue with respect to this block-wise Newton method according to an embodiment of the invention is its convergence because some variables will be updated more than once in each iteration. For instance, there will be different Δu i  when solving L ij  and L ik . A straightforward technique would be to add Au, whenever there is an update for it. However, no theoretic proof will guarantee the convergence in this case. 
     Interior Point Method Algorithm 
     Based on the above linear system, given a tolerance ε&gt;0, barrier parameter μ=μ 0 , and iteration counter k=1, an interior point algorithm according to an embodiment of the invention can be designed as follows, with reference to the flowchart of  FIG. 3 . 
     Begin Initialize variables X k =X 0 , S k =S 0 , . . . (Step  31 ) 
     Repeat
         Update μ k+1 ←ε·μu k , εε(0.1, . . . , 0.9); (Step  32 )   Evaluate (ΔX k , ΔS k , ΔY k , Δy i   k , Δy 2   k ) using EQS. (30) and (31); (Step  33 )   Calculate (ΔZ X   k , ΔZ S   k ) using EQS. (27); (Step  34 )   Find a step length Δβ k  by a line search; (Step  35 )   Update: (Step  36 )
 
 X   k+1   ←X   k +Δβ k   ΔX   k ,
 
 S   k+1   ←S   k +Δβ k   ΔS   k ,
 
 Y   k+1   ←Y   k +Δβ k   ΔY   k ,
 
 y   1   k+1   ←y   1   k +Δβ k   Δy   1   k ,
 
 y   2   k+1   ←y   2   k +Δβ k   Δy   2   k ,
 
 Z   X   k+1   ←Z   X   k +Δβ k   ΔZ   X   k ,
 
 Z   S   k+1   ←Z   S   k +Δβ k   ΔZ   S   k ,
       

     Until Stop conditions are satisfied (Step  37 ). 
     The initial values of the variables depend upon the situation. To measure optimality, a dual gap according to an embodiment of the invention is defined as:
 
gap= X   UB     T     Y+b   x   T   y   1   +b   v   T   y   2 −½ X   T   VX −ƒ( x,u ).  (56)
 
A stop condition according to an embodiment of the invention can be one of the following:
 
gap≦ε, or  (57)
 
max{ r   X   ,r   S }≦ε.  (58)
 
     The line search of Step  15  tries to find the suitable step size. Two methods which could be considered are the Trust-region method and the Filter line-search. 
     A trust region method has a region around the current search point, where the quadratic model q k (β k )=ƒ(X k )+∇ƒ(X k ) T β k +½β k   T B k β k  for “local minimization” is “trusted” to be correct and steps are chosen to stay within this region. The size of the region is modified during the search, based on how well the model agrees with actual function evaluations. Once a step β k  is chosen, the function is evaluated at the new point, and the actual function value is checked against the value predicted by the quadratic model. What is actually computed is the ratio of actual to predicted reduction: 
               β   k     =           f   ⁡     (     X   k     )       -     f   ⁡     (       X   k     +     β   k       )               q   k     ⁡     (   0   )       -       q   k     ⁡     (     β   k     )           =         actual   ⁢           ⁢   reduction   ⁢           ⁢   of   ⁢           ⁢   f       predicted   ⁢           ⁢   model   ⁢           ⁢   reduction   ⁢           ⁢   of   ⁢           ⁢   f       .             
If β k  is close to 1, then the quadratic model is quite a good predictor and the region can be increased in size. On the other hand, if β k  is too small, the region is decreased in size. When k is below a threshold η, the step is rejected and recomputed.
 
     A filter line-search interprets EQS. (19) and (14.1)-(14.5) as a bi-objective optimization program with twin goals of minimizing the objective function L μ (X, S) and the constraint violation θ=∥c(X)∥. Following this paradigm, a trial point X k (β k )=X k +β k ΔX k  may be acceptable, if it leads to sufficient progress toward either objective compared to the current iterate, i.e., if
 
θ( X   k (β k ))≦(1−γ θ )θ( X   k ), or
 
 L   μ ( X   k (β k ), S   k (β k ))≦ L   μ ( X   k   ,S   k )−γ L θ( X   k )
 
holds for fixed constant γ θ , γ L ε(0,1).
 
Extended OPF Model
 
     An extended OPF model according to other embodiments of the invention assumes an AC-only grid with continuous-only dependent (non-control) variables x and control variables u, and that a particular bus can be connected to more than one generator or load, or both at the same time. In addition to that, no thermal capacity limits are enforced on the branch power flows. An extended OPF model according to other embodiments of the invention can be solved by the above disclosed methods for solving a discrete AC OPF model and a minimal AC OPF model according to other embodiments of the invention. 
     Bus dependent variables vector x i  typically includes net active x i   P  and reactive x i   Q , and it is computed for each bus i according to EQ. (1): 
                       x   i     =         ∑     k   =   1       G   i       ⁢     x     i   ,   k     g       -       ∑     k   =   1       L   i       ⁢     x     i   ,   k     l       +       ∑     k   =   1       SH   i       ⁢     x     i   ,   k     sh           ,           (   60   )               
where G i , L i  and SH i  are the number of generators, loads and shunts placed on bus i and x i,k   g , x i,k   l  and x i,k   sh  are the k th  active and reactive generation, load and shunt variables at bus i, respectively. Power system models may include loads and shunts as a function of the associated node&#39;s voltage magnitudes. However, in a model according to an embodiment of the invention, load and shunt devices are assumed to operate as constant power devices. EQ. (60) may be rewritten in a vector notation as
 
 x=Ã{tilde over (x)},   (61)
 
where Ãε{−1,0,1} N×2(G+L+SH)  indicates which generators, loads or shunts are connected to the same bus, and {tilde over (x)} is a vector of generators, loads and shunts on the network.
 
     Box Constraints: 
     Two box-constraint families may be defined that independently limit bus control vector u and dependent bus variable column vector {tilde over (x)}:
 
 {tilde over (x)}   LB   ≦{tilde over (x)}≦{tilde over (x)}   UB ,
 
 u   LB   ≦u≦u   UB ,
 
Bus control vector u i  typically includes voltage magnitude V i  (log-scaled according to an embodiment of the invention), phase angle θ i ; and shunt switch position Q i .
 
     Similarly, box-constraints may be defined for branch variables that independently limit branch control vector v ij , which typically includes transformer tap selection t ij  and phase shift φ ij , and dependent branch variable vector y ij , which typically includes net active y ij   P  and reactive y ij   Q  power injections, for each branch ij:
 
 y   LB   ≦y≦y   UB ,
 
 v   LB   ≦v≦v   UB ,
 
     Flow Balance Set: 
     At every bus, a net power injection balance should be maintained between a bus and the connected branches. A branch dependent variables vector y ij  may be defined, in which typically y ij =y ij   P , y ij   Q ), and the balance for each bus i may be specified as follows: 
                       a   i     ⁢     x   ~       =       ∑     j   =   1     N     ⁢       R   ij     ⁢     y   ij                 (   62   )               
where a i  is i th  row of matrix Ã. Note that matrix R is always highly sparse since the connection density in real-world power grids is typically 1 to 5 branches per bus.
 
     Linearized Voltage Law: 
     According to an embodiment of the invention, a first order approximation of the voltage law constraints is used. A linearized version of the AC voltage law defining the relation between branch power flow vectors y ij  and control vectors u i , u j , v ij  which are linearized with respect to a given fixed point (u i   0 , u j   0 , v ij   0 ) is written as
 
 y   ij   =R   ij ·( A   ij   0   ·u   i   +B   ij   0   ·u   i   +C   ij   0   ·v   ij   +D   ij   0   ·v   ji   +E   ij   0 )  (63)
 
     This constraint exhibits two types of sparsity: one inherited from R ij , and the other due to dependence on only two control vectors (u i , u j ) and (v ij , v ji ) instead of the whole-grid control vector U=└(u i ) i=1   N ; (v ij ) i,j=1   N ┘. The matrices A ij   0 , B ij   0 , C ij   0 , D ij   0 , E ij   0  can be explicitly pre-calculated for every link (i,j) for the given point (u i   0 , u j   0 , v ij   0 ). 
     Objective Function: 
     An OPF model according to an embodiment of the invention assumes a generic convex quadratic cost function to be minimized. This function includes terms related to managing the vector u i  that minimizes control adjustment costs, bus dependent variables vector {tilde over (x)} i  that minimize generation costs, and branch dependent variables vector y ij  that minimize transmission losses, and a penalty coefficient for the auxiliary variable β:
 
ƒ( {tilde over (x)},u,y,v )=½ [{tilde over (x)}   T   V   X   {tilde over (x)}+u   T   V   U   u+y   T   W   Y   y+v   T   W   V   v]+S   X   {tilde over (x)}+S   U   u+T   Y   y+T   V   v.   (64)
 
Interior Point Method Solver Math for the Minimal OPF Model
 
     A model according to an embodiment of the invention can be standardized to the following format:
 
½ [x   1   T   V   X   {tilde over (x)}   1   +u   1   T   V   U   u   1   +y   i   T   W   Y   y   1   +v   1   T   W   V   v   1   ]+S   X   {tilde over (x)}   1   +S   U   u   1   +T   Y   y   1   +T   V   v   1 .  (65)
 
Slack variables are introduced for all bound constraints:
 
                           x   ~       i   ⁢           ⁢   1       +     s     x   ⁢           ⁢   i         =         x   ~       i   ⁢           ⁢   1     UB     ⁢     ∀   i         ⁢     
     ⁢         u     i   ⁢           ⁢   1       +     s   ui       =       u     i   ⁢           ⁢   1     UB     ⁢     ∀   i         ⁢     
     ⁢         y     ij   ⁢           ⁢   1       +     s   yij       =       y     ij   ⁢           ⁢   1     UB     ⁢     ∀   ij         ⁢     
     ⁢         v     ij   ⁢           ⁢   1       +     s   vij       =       v     ij   ⁢           ⁢   1     UB     ⁢     ∀   ij                 (   66.1   )                       α   ~     i     ⁢       x   ~     1       -       ∑     j   =   1     N     ⁢       R   ij     ·     y     ij   ⁢           ⁢   1             =     b   x             (   66.2   )                   y     ij   ⁢           ⁢   1       -       R   ij     ·     (         A   ij   0     ·     u     i   ⁢           ⁢   1         +       B   ij   0     ·     u     j   ⁢           ⁢   1         +       C   ij   0     ·     v     ij   ⁢           ⁢   1         +       D   ij   0     ·     v     ij   ⁢           ⁢   1           )         =     b   v             (   66.3   )                   x   ~     1     ,     y   1     ,     u   1     ,     v   1     ,     s     x   ~       ,     s   u     ,     s   y     ,       s   v     ≥   0             (   66.4   )               
where S {tilde over (x)} , S u , S y , S v  are slack variables, and
 
     
       
         
           
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         
                           x 
                           ~ 
                         
                         
                           i 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                       
                       = 
                       
                         
                           
                             x 
                             ~ 
                           
                           i 
                         
                         - 
                         
                           
                             x 
                             ~ 
                           
                           i 
                           LB 
                         
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         y 
                         
                           ij 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                       
                       = 
                       
                         
                           y 
                           ij 
                         
                         - 
                         
                           y 
                           ij 
                           LB 
                         
                       
                     
                     , 
                   
                 
               
               
                 
                   ( 
                   67.1 
                   ) 
                 
               
             
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         
                           x 
                           ~ 
                         
                         
                           i 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                         UB 
                       
                       = 
                       
                         
                           
                             x 
                             ~ 
                           
                           i 
                           UB 
                         
                         - 
                         
                           
                             x 
                             ~ 
                           
                           i 
                           LB 
                         
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         y 
                         
                           i 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           j 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                         UB 
                       
                       = 
                       
                         
                           y 
                           ij 
                           UB 
                         
                         - 
                         
                           y 
                           ij 
                           LB 
                         
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         u 
                         
                           i 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                       
                       = 
                       
                         
                           u 
                           i 
                         
                         - 
                         
                           u 
                           i 
                           LB 
                         
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         v 
                         
                           ij 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                       
                       = 
                       
                         
                           v 
                           ij 
                         
                         - 
                         
                           v 
                           ij 
                           LB 
                         
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         u 
                         
                           i 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                         UB 
                       
                       = 
                       
                         
                           u 
                           i 
                           UB 
                         
                         - 
                         
                           u 
                           i 
                           LB 
                         
                       
                     
                     , 
                     
                       
 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       
                         v 
                         
                           ij 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           1 
                         
                         UB 
                       
                       = 
                       
                         
                           v 
                           ij 
                           UB 
                         
                         - 
                         
                           v 
                           ij 
                           LB 
                         
                       
                     
                     , 
                   
                 
               
               
                 
                   ( 
                   67.2 
                   ) 
                 
               
             
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     
                       
                         b 
                         x 
                       
                       = 
                       
                         
                           
                             ∑ 
                             
                               j 
                               = 
                               1 
                             
                             N 
                           
                           ⁢ 
                           
                             
                               R 
                               ij 
                             
                             · 
                             
                               y 
                               ij 
                               LB 
                             
                           
                         
                         - 
                         
                           
                             
                               a 
                               ~ 
                             
                             i 
                           
                           ⁢ 
                           
                             
                               x 
                               ~ 
                             
                             i 
                             LB 
                           
                         
                       
                     
                     , 
                   
                 
               
               
                 
                   ( 
                   67.3 
                   ) 
                 
               
             
             
               
                 
                   
                     b 
                     v 
                   
                   = 
                   
                     
                       
                         R 
                         ij 
                       
                       · 
                       
                         ( 
                         
                           
                             E 
                             
                               ij 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               1 
                             
                             0 
                           
                           + 
                           
                             
                               A 
                               ij 
                               0 
                             
                             · 
                             
                               u 
                               i 
                               LB 
                             
                           
                           + 
                           
                             
                               B 
                               ij 
                               0 
                             
                             · 
                             
                               u 
                               j 
                               LB 
                             
                           
                           + 
                           
                             
                               C 
                               ij 
                               0 
                             
                             · 
                             
                               v 
                               ij 
                               LB 
                             
                           
                           + 
                           
                             
                               D 
                               ij 
                               0 
                             
                             · 
                             
                               v 
                               ji 
                               LB 
                             
                           
                         
                         ) 
                       
                     
                     - 
                     
                       y 
                       ij 
                       LB 
                     
                   
                 
               
               
                 
                   ( 
                   67.4 
                   ) 
                 
               
             
           
         
       
     
     According to an embodiment of the invention, the coefficients of the objective function ƒ can be changed:
 
 S   i1   X   =S   i   X   +V   i   X   {tilde over (x)}   i   LB  
 
 S   i1   U   =S   i   U   +V   i   U   u   j   LB  
 
 T   ij1   Y   =T   ij   Y   +W   ij   X   {tilde over (x)}   i   LB  
 
 T   ij1   V   =T   ij   V   +W   ij   V   u   j   LB   (68)
 
Let
 
 v ( y   ij   ,u   i   ,v   ij ) ij   =y   ij1   −R   ij ·( A   ij   0   ·u   i1   +B   ij   0   ·u   j1   +C   ij   0   ·v   ij1   +D   ij   0   ·v   ji1 )− b   v .  (69)
 
and
 
     
       
         
           
             
               
                 
                   
                     
                       o 
                       ⁡ 
                       
                         ( 
                         
                           
                             y 
                             ij 
                           
                           , 
                           
                             x 
                             ~ 
                           
                         
                         ) 
                       
                     
                     i 
                   
                   = 
                   
                     
                       
                         
                           a 
                           ~ 
                         
                         i 
                       
                       ⁢ 
                       
                         
                           x 
                           ~ 
                         
                         1 
                       
                     
                     - 
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           1 
                         
                         N 
                       
                       ⁢ 
                       
                         
                           R 
                           ij 
                         
                         · 
                         
                           y 
                           
                             ij 
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             1 
                           
                         
                       
                     
                     - 
                     
                       
                         b 
                         x 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   70 
                   ) 
                 
               
             
           
         
       
     
     According to an embodiment of the invention, notation can be simplified and readability can be improved by redefining the variable vector x and slack variables s so that 
     
       
         
           
             
               
                 
                   x 
                   = 
                   
                     
                       
                         [ 
                         
                           
                             
                               
                                 x 
                                 ~ 
                               
                             
                           
                           
                             
                               u 
                             
                           
                           
                             
                               y 
                             
                           
                           
                             
                               v 
                             
                           
                         
                         ] 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       and 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       s 
                     
                     = 
                     
                       [ 
                       
                         
                           
                             
                               s 
                               
                                 x 
                                 ~ 
                               
                             
                           
                         
                         
                           
                             
                               s 
                               u 
                             
                           
                         
                         
                           
                             
                               s 
                               y 
                             
                           
                         
                         
                           
                             
                               s 
                               v 
                             
                           
                         
                       
                       ] 
                     
                   
                 
               
               
                 
                   ( 
                   71 
                   ) 
                 
               
             
           
         
       
     
     By introducing logarithmic barrier terms, the following Lagrangian function can be derived:
 
 L   μ ( x,s )=ƒ( x )−μ T (ln  x +ln  s )−η T ( x+s−x   UB )−λ 1   T   o ( y   ij   ,{tilde over (x)} )−λ 2   T   v ( y   ij   ,u   i   ,v   ij ).  (72)
 
The Lagrange multiplier vectors η, λ 1 , λ 2 , are referred to as dual variables.
 
     The KKT optimality conditions of the barrier system are: 
                     [           Δ   ⁢           ⁢     L   X                 Δ   ⁢           ⁢     L   S                 Δ   ⁢           ⁢     L   η                 Δ   ⁢           ⁢     L     λ   1                   Δ   ⁢           ⁢     L     λ   2               ]     =       [                   Δ   ⁢           ⁢     f   ⁡     (   x   )         -     μ   /   x     -   η   -         ∇   0     ⁢       (       y     ij   ⁢           ⁢   1       ,       x   ~     1       )     T       ⁢     λ   1       -               v   ⁢       (       y   ij     ,     u   i     ,     v   ij       )     T     ⁢     λ   2                           -   μ     /   s     -   η               -     (     x   +   s   -     x   UB       )                 -     o   ⁡     (       y   ij     ,       x   ~     1       )                   -     v   ⁡     (       y   ij     ,     u   i     ,     v   ij       )               ]     =   0.             (   73   )               
By defining z x =μ/x and z s =μ/s, EQ. (73) becomes equivalent to:
 
                     [             ∇     f   ⁡     (   x   )         -     z   x     -   η   -         ∇   0     ⁢       (       y     ij   ⁢           ⁢   1       ,       x   ~     1       )     T       ⁢     λ   1       -       ∇       v   ⁡     (       y   ij     ,     u   i     ,     v   ij       )       T       ⁢     λ   2                     -   η     -     z   s                 -     (     x   +   s   -     x   UB       )                 -     o   ⁡     (       y   ij     ,       x   ~     1       )                   -     v   ⁡     (       y   ij     ,     u   i     ,     v   ij       )                       z   x   T     ⁢   x     -     μ   ⁢           ⁢   e                     z   s   T     ⁢   s     -     μ   ⁢           ⁢   e             ]     =   0           (   74   )               
Newton&#39;s method can be used to solve the above non-linear equations, which yields the following linear system of equations, where His a Jacobian matrix:
 
     
       
         
           
             
               
                 
                   
                     H 
                     · 
                     
                       [ 
                       
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               x 
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               s 
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               η 
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 λ 
                                 1 
                               
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 λ 
                                 2 
                               
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 z 
                                 x 
                               
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 z 
                                 s 
                               
                             
                           
                         
                       
                       ] 
                     
                   
                   = 
                   
                     - 
                     
                       [ 
                       
                         
                           
                             
                               
                                 
                                   
                                     
                                       ∇ 
                                       
                                         f 
                                         ⁡ 
                                         
                                           ( 
                                           x 
                                           ) 
                                         
                                       
                                     
                                     - 
                                     
                                       z 
                                       x 
                                     
                                     - 
                                     η 
                                     - 
                                     
                                       
                                         
                                           ∇ 
                                           0 
                                         
                                         ⁢ 
                                         
                                           
                                             ( 
                                             
                                               
                                                 y 
                                                 
                                                   ij 
                                                   ⁢ 
                                                   
                                                       
                                                   
                                                   ⁢ 
                                                   1 
                                                 
                                               
                                               , 
                                               
                                                 
                                                   x 
                                                   ~ 
                                                 
                                                 1 
                                               
                                             
                                             ) 
                                           
                                           T 
                                         
                                       
                                       ⁢ 
                                       
                                         λ 
                                         1 
                                       
                                     
                                     - 
                                   
                                 
                               
                               
                                 
                                   
                                     
                                       ∇ 
                                       v 
                                     
                                     ⁢ 
                                     
                                       
                                         ( 
                                         
                                           
                                             y 
                                             ij 
                                           
                                           , 
                                           
                                             u 
                                             i 
                                           
                                           , 
                                           
                                             v 
                                             ij 
                                           
                                         
                                         ) 
                                       
                                       T 
                                     
                                     ⁢ 
                                     
                                       λ 
                                       2 
                                     
                                   
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 - 
                                 η 
                               
                               - 
                               
                                 z 
                                 s 
                               
                             
                           
                         
                         
                           
                             
                               - 
                               
                                 ( 
                                 
                                   x 
                                   + 
                                   s 
                                   - 
                                   
                                     x 
                                     UB 
                                   
                                 
                                 ) 
                               
                             
                           
                         
                         
                           
                             
                               - 
                               
                                 o 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       y 
                                       ij 
                                     
                                     , 
                                     
                                       
                                         x 
                                         ~ 
                                       
                                       1 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                         
                           
                             
                               - 
                               
                                 v 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       y 
                                       ij 
                                     
                                     , 
                                     
                                       u 
                                       i 
                                     
                                     , 
                                     
                                       v 
                                       ij 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                         
                           
                             
                               
                                 
                                   z 
                                   x 
                                   T 
                                 
                                 ⁢ 
                                 x 
                               
                               - 
                               
                                 μ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 e 
                               
                             
                           
                         
                         
                           
                             
                               
                                 
                                   z 
                                   s 
                                   T 
                                 
                                 ⁢ 
                                 s 
                               
                               - 
                               
                                 μ 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 e 
                               
                             
                           
                         
                       
                       ] 
                     
                   
                 
               
               
                 
                   ( 
                   75 
                   ) 
                 
               
             
           
         
       
     
     Let Z x , Z s , X and S be the diagonal matrices of z x , z s , x and s, respectively. Then, the Jacobian matrix H takes the following form: 
                           ⁢       H   =     [           H   11         0         -   I           H   41   T           H   51   T           -   I         0           0       0         -   I         0       0       0         -   I               -   I           -   I         0       0       0       0       0             H   41         0       0       0       0       0       0             H   51         0       0       0       0       0       0             Z   x         0       0       0       0       X       0           0         Z   s         0       0       0       0       S         ]       ⁢     
     ⁢           ⁢   where             (   76   )                       ⁢         H   11     =     diag   ⁡     (       V   X     ,     V   U     ,     W   Y     ,     W   V       )         ,             (   77.1   )                       ⁢         H   41     =     [       -   I     ,   0   ,   R   ,   0     ]       ,             (   77.2   )                       ⁢         H   51     =     [     0   ,     H   51   u     ,     -   I     ,     H   51   v       ]       ,             (   77.3   )                       ⁢       R   =     (           [       R   11     ⁢           ⁢   …   ⁢           ⁢     R     1   ⁢   N         ]         0       …       0           0         [       R   21     ⁢           ⁢   …   ⁢           ⁢     R     2   ⁢   N         ]         …       0           ⋮       ⋮       ⋱       ⋮           0       …       0         [       R     N   ⁢           ⁢   1       ⁢           ⁢   …   ⁢           ⁢     R   NN       ]           )       ,             (   77.4   )                   (     H   51   u     )     ij     =     {           (         0       …           R   ij     ⁢     A   ij   0           0       …           R   ij     ⁢     B   ij   0           0       …         )             if   ⁢           ⁢   i     &lt;   j               (         0       …           R   ij     ⁢     B   ij   0           0       …           R   ij     ⁢     A   ij   0           0       …         )             if   ⁢           ⁢   i     &gt;   j               (         0       …             R   ij     ⁢     A   ij   0       +       R   ij     ⁢     B   ij   0             0       …       0       0       …         )             if   ⁢           ⁢   i     =   j                     (   77.5   )                   (     H   51   v     )     ij     =     {           (         0       …           R   ij     ⁢     C   ij   0           0       …           R   ij     ⁢     D   ij   0           0       …         )             if   ⁢           ⁢   i     &lt;   j               (         0       …           R   ij     ⁢     D   ij   0           0       …           R   ij     ⁢     C   ij   0           0       …         )             if   ⁢           ⁢   i     &gt;   j               (         0       …             R   ij     ⁢     C   ij   0       +       R   ij     ⁢     D   ij   0             0       …       0       0       …         )             if   ⁢           ⁢   i     =   j                                     
where (H 51   u ) ij  or (H 51   v ) ij  is the row corresponding to branch (i,j). Note that all four diagonal blocks in the matrix H 11  are themselves diagonal by design.
 
     According to an embodiment of the invention, one can define 
                       [             ∇     f   ⁡     (   x   )         -     z   x     -   η   -         ∇   0     ⁢       (       y     ij   ⁢           ⁢   1       ,       x   ~     1       )     T       ⁢     λ   1       -       ∇       v   ⁡     (       y   ij     ,     u   i     ,     v   ij       )       T       ⁢     λ   2                     -   η     -     z   s                 -     (     x   +   s   -     x   UB       )                 -     o   ⁡     (       y   ij     ,       x   ~     1       )                   -     v   ⁡     (       y   ij     ,     u   i     ,     v   ij       )                       z   x   T     ⁢   x     -     μ   ⁢           ⁢   e                     z   s   T     ⁢   s     -     μ   ⁢           ⁢   e             ]     =     [           r   x               r   s               r   η               r     λ   1                 r     λ   2                 r     z   x                 r     z   s             ]       ⁢     
     ⁢           ⁢     Then   ,             (   78   )                       ⁢         [           H   11         0         -   I           H   41   T           H   51   T           -   I         0           0       0         -   I         0       0       0         -   I               -   I           -   I         0       0       0       0       0             H   41         0       0       0       0       0       0             H   51         0       0       0       0       0       0             Z   x         0       0       0       0       X       0           0         Z   s         0       0       0       0       S         ]     ⁡     [           Δ   ⁢           ⁢   x               Δ   ⁢           ⁢   s               Δ   ⁢           ⁢   η               Δ   ⁢           ⁢     λ   1                 Δ   ⁢           ⁢     λ   2                 Δ   ⁢           ⁢     z   x                 Δ   ⁢           ⁢     z   s             ]       =     -       [           r   x               r   s               r   η               r     λ   1                 r     λ   2                 r     z   x                 r     z   s             ]     .                 (   79   )               
EQS. (79) are analogous to EQS. (26), above, and like EQS. (26), the dimension of the system of EQS. (79) can be very large. However, a procedure according to an embodiment of the invention that uses a series of back substitutions, described below, similar to that used for solving EQS. (26), can reduce the dimensionality of the system so that the equations require much less computational time to solve.
 
     According to an embodiment of the invention, the above system of Newton equations can be solved by reducing the dimension as follows. From the last two equation there is:
 
Δ z   X   =X   −1 (− r   Z     x     −Z   X   Δx )  (80)
 
Δ z   S   =S   −1 (− r   Z     s     −Z   S   Δs )  (81)
 
EQS. (80) and (81) can be substituted back into the other equations, which can be re-organized into primal and dual parts as follows:
 
     
       
         
           
             
               
                 
                   
                     
                       [ 
                       
                         
                           
                             
                               H 
                               41 
                             
                           
                           
                             0 
                           
                           
                             | 
                           
                           
                             0 
                           
                           
                             0 
                           
                           
                             0 
                           
                         
                         
                           
                             
                               H 
                               51 
                             
                           
                           
                             0 
                           
                           
                             | 
                           
                           
                             0 
                           
                           
                             0 
                           
                           
                             0 
                           
                         
                         
                           
                             
                               - 
                               I 
                             
                           
                           
                             
                               - 
                               I 
                             
                           
                           
                             | 
                           
                           
                             0 
                           
                           
                             0 
                           
                           
                             0 
                           
                         
                         
                           
                             - 
                           
                           
                             - 
                           
                           
                             
                                 
                             
                           
                           
                             - 
                           
                           
                             - 
                           
                           
                             - 
                           
                         
                         
                           
                             
                               
                                 H 
                                 11 
                               
                               + 
                               
                                 
                                   X 
                                   
                                     - 
                                     1 
                                   
                                 
                                 ⁢ 
                                 
                                   Z 
                                   x 
                                 
                               
                             
                           
                           
                             0 
                           
                           
                             | 
                           
                           
                             
                               - 
                               I 
                             
                           
                           
                             
                               H 
                               41 
                               T 
                             
                           
                           
                             
                               H 
                               51 
                               T 
                             
                           
                         
                         
                           
                             0 
                           
                           
                             
                               
                                 S 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 Z 
                                 s 
                               
                             
                           
                           
                             | 
                           
                           
                             
                               - 
                               I 
                             
                           
                           
                             0 
                           
                           
                             0 
                           
                         
                       
                       ] 
                     
                     ⁡ 
                     
                       [ 
                       
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               x 
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               s 
                             
                           
                         
                         
                           
                             - 
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               η 
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 λ 
                                 1 
                               
                             
                           
                         
                         
                           
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 λ 
                                 2 
                               
                             
                           
                         
                       
                       ] 
                     
                   
                   = 
                   
                     [ 
                     
                       
                         
                           
                             - 
                             
                               r 
                               
                                 λ 
                                 1 
                               
                             
                           
                         
                       
                       
                         
                           
                             - 
                             
                               r 
                               
                                 λ 
                                 2 
                               
                             
                           
                         
                       
                       
                         
                           
                             - 
                             
                               r 
                               η 
                             
                           
                         
                       
                       
                         
                           - 
                         
                       
                       
                         
                           
                             
                               - 
                               
                                 r 
                                 x 
                               
                             
                             - 
                             
                               
                                 X 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 r 
                                 
                                   z 
                                   x 
                                 
                               
                             
                           
                         
                       
                       
                         
                           
                             
                               - 
                               
                                 r 
                                 s 
                               
                             
                             - 
                             
                               
                                 S 
                                 
                                   - 
                                   1 
                                 
                               
                               ⁢ 
                               
                                 r 
                                 
                                   z 
                                   s 
                                 
                               
                             
                           
                         
                       
                     
                     ] 
                   
                 
               
               
                 
                   ( 
                   82 
                   ) 
                 
               
             
           
         
       
     
     According to an embodiment of the invention, define: 
                         H   ~     11     =       H   11     +       X     -   1       ⁢     Z   x           ⁢     
     ⁢         H   ~     22     =       S     -   1       ⁢     Z   s                 (   83.1   )                     r   ~     x     =       r   x     +       X     -   1       ⁢     r     z   x             ⁢     
     ⁢         r   ~     s     =       r   s     +       S     -   1       ⁢     r     Z   s                     (   83.2   )                 Δ   ⁢           ⁢   x     =     [           {     Δ   ⁢           ⁢     x   ~       }               {     Δ   ⁢           ⁢   y     }               {     Δ   ⁢           ⁢   u     }               {     Δ   ⁢           ⁢   v     }           ]             (   83.3   )                   Δ   ⁢           ⁢     λ   1       =     {     Δ   ⁢           ⁢     λ   1     x   ~         }       ⁢     
     ⁢       Δ   ⁢           ⁢     λ   2       =     {     Δ   ⁢           ⁢     λ   2   y       }               (   83.4   )               
Therefore, according to an embodiment of the invention, one has:
 
                         [           H   41         0             H   51         0             -   I           -   I           ]     ⁡     [           Δ   ⁢           ⁢   x               Δ   ⁢           ⁢   s           ]       =     [           -     r     λ   1                   -     r     λ   2                   -     r   η             ]       ,     
     ⁢   and           (   84   )                     [             H   ~     11         0           0           H   ~     22           ]     ⁡     [           Δ   ⁢           ⁢   x               Δ   ⁢           ⁢   s           ]       +       [           -   I           H   41   T           H   51   T               -   I         0       0         ]     ⁡     [           Δ   ⁢           ⁢   η               Δ   ⁢           ⁢     λ   1                 Δ   ⁢           ⁢     λ   2             ]         =       [           -       r   ~     x                 -       r   ~     s             ]     .             (   85   )               
Back Substitution Solution Series for Reduced Newton System
 
     More specifically, similar to the solution of EQS. (30) and (31), above, EQS. (84) and (85) can be expanded using the definitions of EQS. (83.1)-(83.4) and the model definitions of EQS. (67.1)-(67.4), (68)-(71), and (77.1)-(77.5). Thus, EQ. (84) has the following expanded form for each bus i and branch (i, j): 
                           ⁢               a   ~     i     ⁢   Δ   ⁢           ⁢       x   ~     1       -       ∑   j     ⁢       R   ij     ⁢   Δ   ⁢           ⁢     y     ij   ⁢           ⁢   1             =       (     r     λ   1       )     i       ⁢     
     ⁢         Δ   ⁢           ⁢     y     ij   ⁢           ⁢   1         -       R   ij     ⁡     (         A   ij   0     ⁢   Δ   ⁢           ⁢     u     i   ⁢           ⁢   1         +       B   ij   0     ⁢   Δ   ⁢           ⁢     u     j   ⁢           ⁢   1         +       C   ij   0     ⁢   Δ   ⁢           ⁢     v     ij   ⁢           ⁢   1         +       D   ij   0     ⁢   Δ   ⁢           ⁢     v     ji   ⁢           ⁢   1           )         =       (     r     λ     2   ⁢                 )     ij       ⁢     
     ⁢           ⁢         Δ   ⁢           ⁢   x     +     Δ   ⁢           ⁢   s       =     r   η                 (   86   )               
Similarly, EQ. (85) has the following expanded form:
 
                           ⁢           Δ   ⁢           ⁢   η     -     Δ   ⁢           ⁢     λ   1         =     -       r   ~       β   ⁢           ⁢   1           ⁢     
     ⁢           ⁢             V   ~     X     ⁢   Δ   ⁢           ⁢       x   ~     1       -     Δ   ⁢           ⁢   η     -         A   ~     T     ⁢   Δ   ⁢           ⁢     λ   1         =     -       r   ~         x   ~     1                     (   87.1   )                       ⁢               W   ~     Y     ⁢   Δ   ⁢           ⁢     y   1       -     Δ   ⁢           ⁢   η     +     R   ⁢           ⁢   Δ   ⁢           ⁢     λ   1       -     Δ   ⁢           ⁢     λ   2         =     -       r   ~       y   1           ⁢     
     ⁢           ⁢             H   ~     22     ⁢   Δ   ⁢           ⁢   s     -     Δ   ⁢           ⁢   η       =     -       r   ~     s         ⁢     
     ⁢           ⁢   and             (   87.2   )                         V   ~     i   U     ⁢   Δ   ⁢           ⁢     u     i   ⁢           ⁢   1         -     Δ   ⁢           ⁢     η   i       +       ∑   j     ⁢       R   ij     ⁢     A   ij   T     ⁢   Δ   ⁢           ⁢     λ   2     y     ij   ⁢           ⁢   1             +       ∑   j     ⁢       R   ji     ⁢     B   ji   T     ⁢   Δ   ⁢           ⁢     λ   2     y     ji   ⁢           ⁢   1               =     -       (       r   ~       u   1       )     i         ⁢     
     ⁢           ⁢             W   ~     ij   V     ⁢   Δ   ⁢           ⁢     v     ij   ⁢           ⁢   1         -     Δ   ⁢           ⁢     η   ij       ⁢           +       R   ij     ⁢     C   ij   T     ⁢   Δ   ⁢           ⁢     λ   2     y     ij   ⁢           ⁢   1           +       R   ji     ⁢     D   ji   T     ⁢   Δ   ⁢           ⁢     λ   2     v     ji   ⁢           ⁢   1             =     -       (       r   ~       v   1       )     ij                 (   87.3   )               
Suppose there is a candidate Newton step solution for the control variables ΔU=[{Δu i1 }, {Δv ij1 }] T . Then, one can solve for all the other variables directly through EQS. (86) and (87.1)-(87.2) forming a back substitution series of direct linear elimination:
 
{Δ u   1   },{Δv   1   }→{Δy   1 }→{Δ{tilde over (x)} 1   }→Δs→Δη→Δλ   1 →Δλ 2   (88)
 
Then, this back substitution series rule of EQ. (88) can be used to define a linear operator L:ΔU→{Δy 1 , Δ{tilde over (x)} 1 , ΔS, Δη, Δλ 1 , Δλ 2 } defined over the span of ΔU=[{Δu i1 }, {Δv ij1 }] T . This operator can be substituted back into the linear system of EQS. (87.3) which then becomes a system of joint linear equations for the ΔU only. This system is typically full-rank, well-defined and its size is a few orders of magnitude smaller than the reduced Newton system of EQS. (85). In the simplest case of {Δv ij1 }=∅, this system is N by N and relatively sparse (density factor ˜0.1%) since the non-zero entries for any row i are given by the set of “double-neighbor” nodes S 2 (i)={j: ∃k s.t. R jk =1, R k1 =1}. This system is best solved via indirect or iterative methods, such as the Preconditioned Conjugate Method, since computing the full matrix form of the linear operator L is too costly for massive-scale problems.
 
System Implementation
 
     It is to be understood that embodiments of the present invention can be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof. In one embodiment, the present invention can be implemented in software as an application program tangible embodied on a computer readable program storage device. The application program can be uploaded to, and executed by, a machine comprising any suitable architecture. 
       FIG. 4  is a block diagram of an exemplary computer system for implementing an interior point method solution for an optimal power flow model of a smart power grid according to an embodiment of the invention. Referring now to  FIG. 4 , a computer system  41  for implementing the present invention can comprise, inter alia, a central processing unit (CPU)  42 , a memory  43  and an input/output (I/O) interface  44 . The computer system  41  is generally coupled through the I/O interface  44  to a display  45  and various input devices  46  such as a mouse and a keyboard. The support circuits can include circuits such as cache, power supplies, clock circuits, and a communication bus. The memory  43  can include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combinations thereof. The present invention can be implemented as a routine  47  that is stored in memory  43  and executed by the CPU  42  to process the signal from the signal source  48 . As such, the computer system  41  is a general purpose computer system that becomes a specific purpose computer system when executing the routine  47  of the present invention. 
     The computer system  41  also includes an operating system and micro instruction code. The various processes and functions described herein can either be part of the micro instruction code or part of the application program (or combination thereof) which is executed via the operating system. In addition, various other peripheral devices can be connected to the computer platform such as an additional data storage device and a printing device. 
     It is to be further understood that, because some of the constituent system components and method steps depicted in the accompanying figures can be implemented in software, the actual connections between the systems components (or the process steps) may differ depending upon the manner in which the present invention is programmed. Given the teachings of the present disclosure provided herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations or configurations of the present invention. 
     While embodiments of the present invention has been described in detail with reference to exemplary embodiments, those skilled in the art will appreciate that various modifications and substitutions can be made thereto without departing from the spirit and scope of the disclosure as set forth in the appended claims.