numerically more stable AC optimal power flow

docs/src/manual/acOptimalPowerFlow.md

@@ -32,7 +32,7 @@ system = powerSystem()
 addBus!(system; label = "Bus 1", type = 3, active = 0.1, angle = -0.1)
 addBus!(system; label = "Bus 2", reactive = 0.01, magnitude = 1.1)
 
-@branch(minDiffAngle = -pi, maxDiffAngle = pi, reactance = 0.5, type = 2)
+@branch(minDiffAngle = -pi, maxDiffAngle = pi, reactance = 0.5, type = 1)
 addBranch!(system; label = "Branch 1", from = "Bus 1", to = "Bus 2", maxFromBus = 0.15)
 
 @generator(maxActive = 0.5, minReactive = -0.1, maxReactive = 0.1, status = 0)
@@ -200,14 +200,19 @@ print(system.branch.label, analysis.method.constraint.voltage.angle)
 ---
 
 ##### [Branch Flow Constraints](@id ACBranchFlowConstraintsManual)
-The `flow` field contains references to the inequality constraints associated with the apparent power flow, active power flow, or current flow magnitude limits at the from-bus and to-bus ends of each branch. The type to which one of the constraints will be applied is defined according to the `type` keyword within the [`addBranch!`](@ref addBranch!) function:
-* `type = 1` for the apparent power flow,
-* `type = 2` for the active power flow,
-* `type = 3` for the current flow magnitude.
+The `flow` field refers to inequality constraints that enforce limits on the apparent power flow, active power flow, or current flow magnitude at the from-bus and to-bus ends of each branch. The type of constraint applied is specified using the `type` keyword in the [`addBranch!`](@ref addBranch!) function:
+* `type = 1` active power flow,
+* `type = 2` apparent power flow,
+* `type = 3` apparent power flow with a squared inequality constraint,
+* `type = 4` current flow magnitude,
+* `type = 5` current flow magnitude with a squared inequality constraint.
 
-These limits are specified using the `minFromBus`, `maxFromBus`, `minToBus` and `maxToBus` keywords within the [`addBranch!`](@ref addBranch!) function. By default, these limit keywords are associated with apparent power (`type = 1`).
+!!! tip "Tip"
+    Squared versions of constraints typically make the optimization problem numerically more robust. However, they often result in slower convergence compared to their non-squared counterparts used in the constraints.
+
+These limits are specified using the `minFromBus`, `maxFromBus`, `minToBus` and `maxToBus` keywords within the [`addBranch!`](@ref addBranch!) function. By default, these limit keywords are associated with apparent power (`type = 3`).
 
-However, in the example, we configured it to use active power flow by setting `type = 2`. To access the flow constraints of branches at the from-bus end, we can utilize the following code snippet:
+However, in the example, we configured it to use active power flow by setting `type = 1`. To access the flow constraints of branches at the from-bus end, we can utilize the following code snippet:
 ```@repl ACOptimalPowerFlow
 print(system.branch.label, analysis.method.constraint.flow.from)
 ```

docs/src/manual/powerSystemModel.md

@@ -263,7 +263,7 @@ The functions [`addBus!`](@ref addBus!), [`addBranch!`](@ref addBranch!), and [`
 ##### Default Keyword Values
 Regarding the [`addBus!`](@ref addBus!) function, the bus type is automatically configured as a demand bus with `type = 1`. The initial bus voltage magnitude is set to `magnitude = 1.0` per-unit, while the base voltage is established as `base = 138e3` volts. Additionally, the minimum and maximum bus voltage magnitudes are set to `minMagnitude = 0.9` per-unit and `maxMagnitude = 1.1` per-unit, respectively.
 
-Transitioning to the [`addBranch!`](@ref addBranch!) function, the default operational status is `status = 1`, indicating that the branch is in-service. The off-nominal turns ratio for the transformer is specified as `turnsRatio = 1.0`, and the phase shift angle is set to `shiftAngle = 0.0`, collectively defining the line configuration with these standard settings. The flow rating is also configured as `type = 1`. Moreover, the minimum and maximum voltage angle differences between the from-bus and to-bus ends are set to `minDiffAngle = -2pi` and `maxDiffAngle = 2pi`, respectively.
+Transitioning to the [`addBranch!`](@ref addBranch!) function, the default operational status is `status = 1`, indicating that the branch is in-service. The off-nominal turns ratio for the transformer is specified as `turnsRatio = 1.0`, and the phase shift angle is set to `shiftAngle = 0.0`, collectively defining the line configuration with these standard settings. The flow rating is also configured as `type = 3`. Moreover, the minimum and maximum voltage angle differences between the from-bus and to-bus ends are set to `minDiffAngle = -2pi` and `maxDiffAngle = 2pi`, respectively.
 
 Similarly, the [`addGenerator!`](@ref addGenerator!) function designates an operational generator by employing `status = 1`, and it sets `magnitude = 1.0` per-unit, denoting the desired voltage magnitude setpoint.
 

docs/src/tutorials/acOptimalPowerFlow.md

@@ -14,7 +14,7 @@ system = powerSystem()
 addBus!(system; label = 1, type = 3, active = 0.1, angle = -0.1)
 addBus!(system; label = 2, reactive = 0.01, magnitude = 1.1)
 
-@branch(minDiffAngle = -pi, maxDiffAngle = pi, reactance = 0.5, type = 2)
+@branch(minDiffAngle = -pi, maxDiffAngle = pi, reactance = 0.5, type = 1)
 addBranch!(system; label = 1, from = 1, to = 2, maxFromBus = 0.15, maxToBus = 0.15)
 
 @generator(maxActive = 0.5, minReactive = -0.1, maxReactive = 0.1)
@@ -286,7 +286,9 @@ print(analysis.method.constraint.voltage.angle)
 ##### Branch Flow Constraints
 The inequality constraints related to the branch flow ratings can be associated with the limits on apparent power flow, active power flow, or current magnitude at the from-bus and to-bus ends of each branch. The type of constraint applied is determined by the `type` keyword within the [`addBranch!`](@ref addBranch!) function.
 
-Specifically, `type = 1` is used for apparent power flow, `type = 2` for active power flow, and `type = 3` for current magnitude. These constraints can be expressed using the equations ``h_{ij}(\mathbf {V}, \bm{\Theta})`` and ``h_{ji}(\mathbf {V}, \bm{\Theta})``, representing the rating constraints at the from-bus and to-bus ends of each branch, respectively:
+The `type` value defines the constraint as follows: `type = 1` applies to active power flow; `type = 2` and `type = 3` apply to apparent power flow; and `type = 4` and `type = 5` apply to current magnitude. When `type = 3` or `type = 5` is selected, squared inequality constraints are used. These constraints typically provide a more numerically robust optimization problem but often result in slower convergence compared to the non-squared versions.
+
+These constraints are mathematically expressed through the equations ``h_{ij}(\mathbf {V}, \bm{\Theta})`` and ``h_{ji}(\mathbf {V}, \bm{\Theta})``, representing the rating constraints at the from-bus and to-bus ends of each branch, respectively:
 ```math
 \begin{aligned}
     F_{ij}^{\text{min}} \leq h_{ij}(\mathbf {V}, \bm{\Theta}) \leq F_{ij}^{\text{max}},\;\;\; \forall (i,j) \in \mathcal{E} \\
@@ -300,7 +302,32 @@ The branch flow limits at the from-bus and to-bus ends, denoted as ``\mathbf{F}_
 𝐅ₜ = [system.branch.flow.minToBus system.branch.flow.maxToBus]
 ```
 
-By default, JuliaGrid employs the rating constraints linked with the apparent power flow (`type = 1`). This constraint at the from-bus is specified as:
+---
+
+The first option is to define the limit keywords for active power flow constraints (`type = 1`) at the from-bus and to-bus ends of each branch:
+```math
+  \begin{aligned}
+    h_{ij}(\mathbf {V}, \bm{\Theta}) &=
+    \cfrac{ g_{ij} + g_{\text{s}ij}}{\tau_{ij}^2} V_{i}^2 -
+    \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) + b_{ij}\sin(\theta_{ij} - \phi_{ij})\right]V_{i}V_{j} \\
+    h_{ji}(\mathbf {V}, \bm{\Theta}) &= (g_{ij} + g_{\text{s}ij}) V_{j}^2 -
+    \cfrac{1}{\tau_{ij}} \left[g_{ij} \cos(\theta_{ij} - \phi_{ij}) - b_{ij} \sin(\theta_{ij}- \phi_{ij})\right] V_{i} V_j.
+  \end{aligned}
+```
+
+In our example, we have chosen to utilize this type of flow constraints. To access the flow constraints of branches at the from-bus end, you can use the following code snippet:
+```@repl ACOptimalPowerFlow
+print(analysis.method.constraint.flow.from)
+```
+
+Similarly, to access the to-bus end flow constraints of branches you can use the following code snippet:
+```@repl ACOptimalPowerFlow
+print(analysis.method.constraint.flow.to)
+```
+
+---
+
+The second option applies constraints to the apparent power flow (`type = 2`). This constraint at the from-bus is specified as:
 ```math
     h_{ij}(\mathbf {V}, \bm{\Theta}) = \sqrt{  A_{ij} V_i^4 + B_{ij} V_i^2 V_j^2 - 2 [C_{ij} \cos(\theta_{ij} - \phi_{ij}) - D_{ij} \sin(\theta_{ij} - \phi_{ij})]V_i^3 V_j},
 ```
@@ -324,28 +351,7 @@ where:
   \end{gathered}
 ```
 
----
-
-The second option is to define the limit keywords for active power flow constraints (`type = 2`) at the from-bus and to-bus ends of each branch:
-```math
-  \begin{aligned}
-    h_{ij}(\mathbf {V}, \bm{\Theta}) &=
-    \cfrac{ g_{ij} + g_{\text{s}ij}}{\tau_{ij}^2} V_{i}^2 -
-    \cfrac{1}{\tau_{ij}} \left[g_{ij}\cos(\theta_{ij} - \phi_{ij}) + b_{ij}\sin(\theta_{ij} - \phi_{ij})\right]V_{i}V_{j} \\
-    h_{ji}(\mathbf {V}, \bm{\Theta}) &= (g_{ij} + g_{\text{s}ij}) V_{j}^2 -
-    \cfrac{1}{\tau_{ij}} \left[g_{ij} \cos(\theta_{ij} - \phi_{ij}) - b_{ij} \sin(\theta_{ij}- \phi_{ij})\right] V_{i} V_j.
-  \end{aligned}
-```
-
-In our example, we have chosen to utilize this type of flow constraints. To access the flow constraints of branches at the from-bus end, you can use the following code snippet:
-```@repl ACOptimalPowerFlow
-print(analysis.method.constraint.flow.from)
-```
-
-Similarly, to access the to-bus end flow constraints of branches you can use the following code snippet:
-```@repl ACOptimalPowerFlow
-print(analysis.method.constraint.flow.to)
-```
+If users choose `type = 3`, it means that the equations are squared (i.e., the square root is omitted), and the limit values will also be squared accordingly.
 
 ---
 
@@ -357,6 +363,8 @@ The last option involves defining the limit keywords for current magnitude const
   \end{aligned}
 ```
 
+If users choose `type = 5`, it means that the equations are squared (i.e., the square root is omitted), and the limit values will also be squared accordingly.
+
 ---
 
 ##### [Generator Power Capability Constraints](@id ACPowerCapabilityConstraintsTutorials)

src/backend/expressions.jl

@@ -221,6 +221,19 @@ function Iij(
     sqrt(p.A * Vi^2 + p.B * Vj^2 - 2 * Vi * Vj * (p.C * cosθij - p.D * sinθij))
 end
 
+function Iij2(
+    system::PowerSystem,
+    Vi::VariableRef,
+    Vj::VariableRef,
+    sinθij::NonlinearExpr,
+    cosθij::NonlinearExpr,
+    idx::Int64
+)
+    p = IijCoefficient(system.branch, system.model.ac, idx)
+
+    p.A * Vi^2 + p.B * Vj^2 - 2 * Vi * Vj * (p.C * cosθij - p.D * sinθij)
+end
+
 ##### To-Bus End Current Magnitude #####
 function Iji(system::PowerSystem, method::LAV, idx::Int64)
     s = method.state
@@ -250,6 +263,19 @@ function Iji(
     sqrt(p.A * Vi^2 + p.B * Vj^2 - 2 * Vi * Vj * (p.C * cosθ + p.D * sinθ))
 end
 
+function Iji2(
+    system::PowerSystem,
+    Vi::VariableRef,
+    Vj::VariableRef,
+    sinθ::NonlinearExpr,
+    cosθ::NonlinearExpr,
+    idx::Int64
+)
+    p = IjiCoefficient(system.branch, system.model.ac, idx)
+
+    p.A * Vi^2 + p.B * Vj^2 - 2 * Vi * Vj * (p.C * cosθ + p.D * sinθ)
+end
+
 ##### From-Bus End Current Angle #####
 function ψij(system::PowerSystem, method::LAV, idx::Int64)
     IijRe, IijIm = ReImIij(system, method, idx)
@@ -278,6 +304,19 @@ function Sij(
     sqrt(p.A * Vi^4 + p.B * Vi^2 * Vj^2 - 2 * Vi^3 * Vj * (p.C * cosθ - p.D * sinθ))
 end
 
+function Sij2(
+    system::PowerSystem,
+    Vi::VariableRef,
+    Vj::VariableRef,
+    sinθ::NonlinearExpr,
+    cosθ::NonlinearExpr,
+    idx::Int64
+)
+    p = IijCoefficient(system.branch, system.model.ac, idx)
+
+    p.A * Vi^4 + p.B * Vi^2 * Vj^2 - 2 * Vi^3 * Vj * (p.C * cosθ - p.D * sinθ)
+end
+
 ##### To-Bus End Apparent Power #####
 function Sji(
     system::PowerSystem,
@@ -292,6 +331,19 @@ function Sji(
     sqrt(p.A * Vi^2 * Vj^2 + p.B * Vj^4 - 2 * Vi * Vj^3 * (p.C * cosθ + p.D * sinθ))
 end
 
+function Sji2(
+    system::PowerSystem,
+    Vi::VariableRef,
+    Vj::VariableRef,
+    sinθ::NonlinearExpr,
+    cosθ::NonlinearExpr,
+    idx::Int64
+)
+    p = IjiCoefficient(system.branch, system.model.ac, idx)
+
+    p.A * Vi^2 * Vj^2 + p.B * Vj^4 - 2 * Vi * Vj^3 * (p.C * cosθ + p.D * sinθ)
+end
+
 ##### Real and Imaginary Components of Bus Voltage Phasor #####
 function ReImVi(method::LAV, idx::Int64)
     s = method.state

src/backend/internal.jl

@@ -403,7 +403,7 @@ macro default(mode::Symbol)
         template.branch.label = "?"
         template.branch.turnsRatio = 1.0
         template.branch.status = Int8(1)
-        template.branch.type = Int8(1)
+        template.branch.type = Int8(3)
     end
 
     if mode == :template || mode == :generator

src/definition/internal.jl

@@ -116,7 +116,7 @@ Base.@kwdef mutable struct BranchTemplate
     label::String = "?"
     turnsRatio::Float64 = 1.0
     status::Int8 = Int8(1)
-    type::Int8 = Int8(1)
+    type::Int8 = Int8(3)
 end
 
 Base.@kwdef mutable struct GeneratorTemplate

src/optimalPowerFlow/acOptimalPowerFlow.jl

@@ -642,23 +642,32 @@ function addFlow(
 
         if branch.flow.type[idx] == 1
             if from
-                expr = Sij(system, Vi, Vj, sinθij, cosθij, idx)
+                expr = Pij(system, Vi, Vj, sinθij, cosθij, idx)
                 refFrom[idx] = @constraint(jump, minFrom <= expr <= maxFrom)
             end
             if to
-                expr = Sji(system, Vi, Vj, sinθij, cosθij, idx)
+                expr = Pji(system, Vi, Vj, sinθij, cosθij, idx)
                 refTo[idx] = @constraint(jump, minTo <= expr <= maxTo)
             end
         elseif branch.flow.type[idx] == 2
             if from
-                expr = Pij(system, Vi, Vj, sinθij, cosθij, idx)
+                expr = Sij(system, Vi, Vj, sinθij, cosθij, idx)
                 refFrom[idx] = @constraint(jump, minFrom <= expr <= maxFrom)
             end
             if to
-                expr = Pji(system, Vi, Vj, sinθij, cosθij, idx)
+                expr = Sji(system, Vi, Vj, sinθij, cosθij, idx)
                 refTo[idx] = @constraint(jump, minTo <= expr <= maxTo)
             end
         elseif branch.flow.type[idx] == 3
+            if from
+                expr = Sij2(system, Vi, Vj, sinθij, cosθij, idx)
+                refFrom[idx] = @constraint(jump, minFrom^2 <= expr <= maxFrom^2)
+            end
+            if to
+                expr = Sji2(system, Vi, Vj, sinθij, cosθij, idx)
+                refTo[idx] = @constraint(jump, minTo^2 <= expr <= maxTo^2)
+            end
+        elseif branch.flow.type[idx] == 4
             if from
                 expr = Iij(system, Vi, Vj, sinθij, cosθij, idx)
                 refFrom[idx] = @constraint(jump, minFrom <= expr <= maxFrom)
@@ -667,6 +676,15 @@ function addFlow(
                 expr = Iji(system, Vi, Vj, sinθij, cosθij, idx)
                 refTo[idx] = @constraint(jump, minTo <= expr <= maxTo)
             end
+        elseif branch.flow.type[idx] == 5
+            if from
+                expr = Iij2(system, Vi, Vj, sinθij, cosθij, idx)
+                refFrom[idx] = @constraint(jump, minFrom^2 <= expr <= maxFrom^2)
+            end
+            if to
+                expr = Iji2(system, Vi, Vj, sinθij, cosθij, idx)
+                refTo[idx] = @constraint(jump, minTo^2 <= expr <= maxTo^2)
+            end
         end
     end
 end
@@ -772,7 +790,7 @@ function updateBalance(
 end
 
 function checkLimit(system::PowerSystem, minFlow::Float64, maxFlow::Float64, i::Int64)
-    if system.branch.flow.type[i] == 1 || system.branch.flow.type[i] == 3
+    if system.branch.flow.type[i] != 1
         if minFlow < 0.0
             minFlow = 0.0
         end

src/powerSystem/branch.jl

@@ -33,9 +33,11 @@ The branch is defined with the following keywords:
 * `minToBus` (pu, VA, W, or A): Minimum branch flow rating at the to-bus end.
 * `maxToBus` (pu, VA, W, or A): Maximum branch flow rating at the to-bus end.
 * `type`: Types of `minFromBus`, `maxFromBus`, `minToBus`, and `maxToBus` branch flow ratings:
-  * `type = 1`: apparent power flow (pu or VA),
-  * `type = 2`: active power flow (pu or W),
-  * `type = 3`: current magnitude flow (pu or A).
+  * `type = 1`: active power flow (pu or W),
+  * `type = 2`: apparent power flow (pu or VA),
+  * `type = 3`: apparent power flow (pu or VA), internally constraint involves its square,
+  * `type = 4`: current magnitude flow (pu or A),
+  * `type = 5`: current magnitude flow (pu or A), internally constraint involves its square.
 
 Note that when powers are given in SI units, they correspond to three-phase power.
 
@@ -47,7 +49,7 @@ are transmitted to the  `Analysis` type.
 
 # Default Settings
 By default, certain keywords are assigned default values: `status = 1`, `turnsRatio = 1.0`,
-`type = 1`, `minDiffAngle = -2pi`, and `maxDiffAngle = 2pi`. The rest of the keywords are
+`type = 3`, `minDiffAngle = -2pi`, and `maxDiffAngle = 2pi`. The rest of the keywords are
 initialized with a value of zero. However, the user can modify these default settings by
 utilizing the [`@branch`](@ref @branch) macro.
 
@@ -566,10 +568,12 @@ macro branch(kwargs...)
                 elseif parameter in [:shiftAngle; :minDiffAngle; :maxDiffAngle]
                     pfxLive = pfx.voltageAngle
                 elseif parameter in [:minFromBus; :maxFromBus; :minToBus; :maxToBus]
-                    if template.branch.type in [1, 3]
-                        pfxLive = pfx.apparentPower
-                    elseif template.branch.type == 2
+                    if template.branch.type == 1
                         pfxLive = pfx.activePower
+                    elseif template.branch.type in [2, 3]
+                        pfxLive = pfx.apparentPower
+                    elseif template.branch.type in [4, 5]
+                        pfxLive = pfx.currentMagnitude
                     end
                 end
                 if pfxLive != 0.0
@@ -605,14 +609,14 @@ function flowType(system::PowerSystem, pfx::PrefixLive, basePowerInv::Float64, i
     baseVoltg = system.base.voltage
 
     if branch.flow.type[i] == 1
-        pfxLive = pfx.apparentPower
+        pfxLive = pfx.activePower
         baseInvFrom = basePowerInv
         baseInvTo = basePowerInv
-    elseif branch.flow.type[i] == 2
-        pfxLive = pfx.activePower
+    elseif branch.flow.type[i] == 2 || branch.flow.type[i] == 3
+        pfxLive = pfx.apparentPower
         baseInvFrom = basePowerInv
         baseInvTo = basePowerInv
-    elseif branch.flow.type[i] == 3
+    elseif branch.flow.type[i] == 4 || branch.flow.type[i] == 5
         pfxLive = pfx.currentMagnitude
         baseInvFrom = baseCurrentInv(
             basePowerInv, baseVoltg.value[branch.layout.from[i]] * baseVoltg.prefix

src/powerSystem/load.jl

@@ -441,7 +441,7 @@ function loadBranch(system::PowerSystem, branchLine::Vector{String})
     branch.flow.maxFromBus = similar(branch.parameter.conductance)
     branch.flow.minToBus = similar(branch.parameter.conductance)
     branch.flow.maxToBus = similar(branch.parameter.conductance)
-    branch.flow.type = fill(Int8(1), branch.number)
+    branch.flow.type = fill(Int8(3), branch.number)
 
     branch.layout.from = fill(0, branch.number)
     branch.layout.to = similar( branch.layout.from)