Current voltage form (lanl-ansi#622)

* add current-voltage formulation in rectangular coordinates
* add tests for OPF and PF problem types
* add initial docs
* constraint formulation tweak
* remove unneeded constraints
* adding doc strings
* Debugged and unit tested
* IVRPowerModel subtype of ACRPowerModel
* dcline power limits through variable definitions
* separate series and total branch current variables
* Update CHANGELOG.md
* add check variable bounds test

Author: frederikgeth

CHANGELOG.md

@@ -9,6 +9,8 @@ PowerModels.jl Change Log
 - Removed `Inf` bounds from variables (#630)
 - Removal of unused functions in `solution.jl`: `get_solution`, `add_generator_power_setpoint`, `add_storage_setpoint`, `add_branch_flow_setpoint`, `add_dcline_flow_setpoint` (breaking) (#637)
 - Fixed a solution reporting bug when the optimizer result count is zero
+- Added IVRForm (rectangular current-voltage formulation) (#622)
+
 
 ### v0.13.1
 - Added DCMPPowerModel for replication of Matpower's DC model (#612)
@@ -17,6 +19,7 @@ PowerModels.jl Change Log
 - Fixed quadratic terms in nonlinear objectives (#621)
 - Fixed incorrect flow bounds based on current flow limits (#631)
 
+
 ### v0.13.0
 - Added SOC formulation of the storage model
 - Added support for optional line flow constraints (#576)

CONTRIBUTING.md

@@ -1,6 +1,6 @@
 ## Comments, Questions, Feature Requests, Bug Reports
 
-At this time, any kind of comment, question, feature request or bug report is welcome in the issue tracker.  If you are asking a question, please search the issues for an answer before making a post. 
+At this time, any kind of comment, question, feature request or bug report is welcome in the issue tracker.  If you are asking a question, please search the issues for an answer before making a post.
 
 ## Code Contributions
 
@@ -9,7 +9,7 @@ Community-driven development and enhancement of PowerModels is welcomed and enco
 That said it is important to keep in mind the code-quality requirements and scope of PowerModels in mind, before preparing a contribution.
 
 ### PowerModels Scope
-The primary motivation for PowerModels is to help researchers quickly and easily replicate _established_ results in power system optimization. 
+The primary motivation for PowerModels is to help researchers quickly and easily replicate _established_ results in power system optimization.
 More preliminary research and development efforts are best left to other packages, which extend PowerModels.
 That said, the PowerModels package should have a software design that makes building such extensions as easy as possible.
 
@@ -19,7 +19,7 @@ To that end, the following code quality controls are necessary for all contribut
 - all features must have unit tests
 - naming conventions must be consistent
 - core features must have documentation
-- code style must be consistent (e.g. 4 spaces for indetnation)
+- code style must be consistent (e.g. 4 spaces for indentation)
 
 Due to these quality control requirements, pull request to PowerModels may require lengthy discussion and code review (e.g. see [#111](https://github.com/lanl-ansi/PowerModels.jl/pull/111)).
 

README.md

@@ -28,6 +28,8 @@ This enables the definition of a wide variety of power network formulations and
 * SDP Relaxation (W-space)
 * SOC Relaxation (W-space)
 * QC Relaxation (W+L-space)
+* IV (rectangular coordinates)
+
 
 **Network Data Formats**
 * Matpower ".m" files
@@ -58,7 +60,7 @@ The primary developer is Carleton Coffrin(@ccoffrin) with support from the follo
 - Hakan Ergun (@hakanergun) KU Leuven, HVDC lines
 - David Fobes (@pseudocubic) LANL, PSS(R)E v33 data support
 - Rory Finnegan (@rofinn) Invenia, Memento Logging
-- Frederik Geth (@frederikgeth) CSIRO, storage modeling advise, Branch Flow formulation
+- Frederik Geth (@frederikgeth) CSIRO, storage modeling advise, Branch Flow and current-voltage formulation
 - Jonas Kersulis (@kersulis) University of Michigan, Sparse SDP formulation
 - Miles Lubin (@mlubin) MIT, Julia/JuMP advise
 - Yeesian Ng (@yeesian) MIT, Documenter.jl setup
@@ -69,13 +71,13 @@ The primary developer is Carleton Coffrin(@ccoffrin) with support from the follo
 
 If you find PowerModels useful in your work, we kindly request that you cite the following [publication](https://ieeexplore.ieee.org/document/8442948/):
 ```
-@inproceedings{8442948, 
-  author = {Carleton Coffrin and Russell Bent and Kaarthik Sundar and Yeesian Ng and Miles Lubin}, 
-  title = {PowerModels.jl: An Open-Source Framework for Exploring Power Flow Formulations}, 
-  booktitle = {2018 Power Systems Computation Conference (PSCC)}, 
+@inproceedings{8442948,
+  author = {Carleton Coffrin and Russell Bent and Kaarthik Sundar and Yeesian Ng and Miles Lubin},
+  title = {PowerModels.jl: An Open-Source Framework for Exploring Power Flow Formulations},
+  booktitle = {2018 Power Systems Computation Conference (PSCC)},
   year = {2018},
   month = {June},
-  pages = {1-8}, 
+  pages = {1-8},
   doi = {10.23919/PSCC.2018.8442948}
 }
 ```

docs/src/formulations.md

@@ -7,6 +7,7 @@ AbstractACPModel <: AbstractPowerModel
 AbstractDCPModel <: AbstractPowerModel
 AbstractWRModel <: AbstractPowerModel
 AbstractWModel <: AbstractPowerModel
+AbstractIVRModel <: AbstractPowerModel
 ```
 
 ## Power Models
@@ -16,6 +17,8 @@ ACPPowerModel <: AbstractACPForm
 
 DCPPowerModel <: AbstractDCPForm
 
+IVRPowerModel <: AbstractIVRModel
+
 SOCWRPowerModel <: AbstractWRForm
 QCRMPowerModel <: AbstractWRForm
 

docs/src/math-model.md

@@ -124,3 +124,42 @@ Note that constraints $\eqref{eq_line_losses} - \eqref{eq_ohms_bfm}$ replace $\e
 - Eq. $\eqref{eq_series_power_flow}$ - implicit, substituted out before implementation
 - Eq. $\eqref{eq_complex_power_definition}$ - `constraint_model_voltage` in `constraint_template.jl`
 - Eq. $\eqref{eq_ohms_bfm}$ - `constraint_voltage_magnitude_difference` in `constraint_template.jl`
+
+
+## AC Optimal Power Flow in Current-Voltage Variables
+A variable $I^{s}_{ij}$, representing the current in the direction $i$ to $j$, through the series part of the pi-section, is used.
+The mathematical structure for a current-voltage formulation is conceived as:
+
+```math
+\begin{align}
+%
+\mbox{variables: } & \nonumber \\
+& I^g_k \;\; \forall k\in G \nonumber \\
+& V_i \;\; \forall i\in N \nonumber \\
+& I^{s}_{ij} \;\; \forall (i,j) \in E \cup E^R  \mbox{ - branch complex (series) current}\\
+& I_{ij} \;\; \forall (i,j) \in E \cup E^R  \mbox{ - branch complex (total) current} \label{var_total_current}\\
+%
+\mbox{minimize: } & \sum_{k \in G} c_{2k} (\Re(S^g_k))^2 + c_{1k}\Re(S^g_k) + c_{0k} \nonumber\\
+%
+\mbox{subject to: } & \nonumber \\
+& \angle V_{r} = 0  \;\; \forall r \in R \nonumber \\
+& S^{gl}_k \leq \Re(V_i (I^g_k)^*) + j \Im(V_i (I^g_k)^*) \leq S^{gu}_k \;\; \forall k \in G   \label{eq_complex_power_definition_gen}\\
+& v^l_i \leq |V_i| \leq v^u_i \;\; \forall i \in N \nonumber\\
+& \sum_{\substack{k \in G_i}} I^g_k - \sum_{\substack{k \in L_i}} (S^d_k/V_i)^{*} - \sum_{\substack{k \in S_i}} Y^s_k V_i = \sum_{\substack{(i,j)\in E_i \cup E_i^R}} I_{ij} \;\; \forall i\in N  \label{eq_kcl_current} \\
+& I_{ij} =  \frac{I^{s}_{ij}}{T_{ij}^*} + Y^c_{ij} \frac{V_i}{|T_{ij}|^2}  \;\; \forall (i,j)\in E \label{eq_current_from} \\
+& I_{ji} = -I^{s}_{ij} + Y^c_{ji} V_j  \;\; \forall (i,j)\in E \label{eq_current_to} \\
+& \frac{V_i}{{T}_{ij}} = V_j + z_{ij} I^{s}_{ij}  \;\; \forall (i,j) \in E \label{eq_ohms_iv} \\
+& |S_{ij}| = |V_{i}| |I_{ij}| \leq s^u_{ij} \;\; \forall (i,j) \in E \cup E^R \nonumber\\
+& |I_{ij}| \leq i^u_{ij} \;\; \forall (i,j) \in E \cup E^R \nonumber\\
+& \theta^{\Delta l}_{ij} \leq \angle (V_i V^*_j) \leq \theta^{\Delta u}_{ij} \;\; \forall (i,j) \in E \nonumber
+%
+\end{align}
+```
+
+### Mapping to function names
+- Eq. $\eqref{var_total_current}$ - total current flow into a branch on either end `variable_branch_current`
+- Eq. $\eqref{eq_complex_power_definition_gen}$  - models active and reactive power range of a generator `constraint_gen`
+- Eq. $\eqref{eq_kcl_current}$  - Kirchhoff's current law in current variables  `constraint_current_balance`
+- Eq. $\eqref{eq_current_from}$  - branch from-side current constraint in `constraint_current_from`
+- Eq. $\eqref{eq_current_to}$  - branch to-side current constraint in `constraint_current_to`
+- Eq. $\eqref{eq_ohms_iv}$  - Ohm's law `constraint_voltage_difference`

src/PowerModels.jl

@@ -53,6 +53,8 @@ include("core/admittance_matrix.jl")
 
 include("io/json.jl")
 
+include("form/iv.jl")
+
 include("form/acp.jl")
 include("form/acr.jl")
 include("form/act.jl")
@@ -67,8 +69,10 @@ include("form/shared.jl")
 include("prob/opb.jl")
 include("prob/pf.jl")
 include("prob/pf_bf.jl")
+include("prob/pf_iv.jl")
 include("prob/opf.jl")
 include("prob/opf_bf.jl")
+include("prob/opf_iv.jl")
 include("prob/ots.jl")
 include("prob/tnep.jl")
 include("prob/test.jl")

src/core/constraint_template.jl

@@ -52,6 +52,16 @@ function constraint_model_voltage_ne(pm::AbstractPowerModel; nw::Int=pm.cnw, cnd
     constraint_model_voltage_ne(pm, nw, cnd)
 end
 
+"""
+This constraint captures problem agnostic constraints that define limits for
+voltage magnitudes (where variable bounds cannot be used)
+
+Notable examples include IVRPowerModel and ACRPowerModel
+"""
+function constraint_voltage_magnitude_bounds(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    bus = ref(pm, nw, :bus, i)
+    constraint_voltage_magnitude_bounds(pm, nw, cnd, i, bus["vmin"][cnd], bus["vmax"][cnd])
+end
 
 ### Current Constraints ###
 
@@ -82,12 +92,26 @@ function constraint_reactive_gen_setpoint(pm::AbstractPowerModel, i::Int; nw::In
     constraint_reactive_gen_setpoint(pm, nw, cnd, gen["index"], gen["qg"][cnd])
 end
 
+""
 function constraint_generation_on_off(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
     gen = ref(pm, nw, :gen, i)
 
     constraint_generation_on_off(pm, nw, cnd, i, gen["pmin"][cnd], gen["pmax"][cnd], gen["qmin"][cnd], gen["qmax"][cnd])
 end
 
+"defines limits on active power output of a generator where bounds can't be used"
+function constraint_gen_active_power_limits(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    gen = ref(pm, nw, :gen, i)
+    bus = gen["gen_bus"]
+    constraint_gen_active_power_limits(pm, nw, cnd, i, bus, gen["pmax"][cnd], gen["pmin"][cnd])
+end
+
+"defines limits on reactive power output of a generator where bounds can't be used"
+function constraint_gen_reactive_power_limits(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    gen = ref(pm, nw, :gen, i)
+    bus = gen["gen_bus"]
+    constraint_gen_reactive_power_limits(pm, nw, cnd, i, bus, gen["qmax"][cnd], gen["qmin"][cnd])
+end
 
 ### Bus - Setpoint Constraints ###
 
@@ -202,6 +226,31 @@ function constraint_power_balance_ne(pm::AbstractPowerModel, i::Int; nw::Int=pm.
     constraint_power_balance_ne(pm, nw, cnd, i, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_arcs_ne, bus_gens, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
 end
 
+""
+function constraint_current_balance(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    if !haskey(con(pm, nw, cnd), :kcl_cr)
+        con(pm, nw, cnd)[:kcl_cr] = Dict{Int,JuMP.ConstraintRef}()
+    end
+    if !haskey(con(pm, nw, cnd), :kcl_ci)
+        con(pm, nw, cnd)[:kcl_ci] = Dict{Int,JuMP.ConstraintRef}()
+    end
+
+    bus = ref(pm, nw, :bus, i)
+    bus_arcs = ref(pm, nw, :bus_arcs, i)
+    bus_arcs_dc = ref(pm, nw, :bus_arcs_dc, i)
+    bus_gens = ref(pm, nw, :bus_gens, i)
+    bus_loads = ref(pm, nw, :bus_loads, i)
+    bus_shunts = ref(pm, nw, :bus_shunts, i)
+
+
+    bus_pd = Dict(k => ref(pm, nw, :load, k, "pd", cnd) for k in bus_loads)
+    bus_qd = Dict(k => ref(pm, nw, :load, k, "qd", cnd) for k in bus_loads)
+
+    bus_gs = Dict(k => ref(pm, nw, :shunt, k, "gs", cnd) for k in bus_shunts)
+    bus_bs = Dict(k => ref(pm, nw, :shunt, k, "bs", cnd) for k in bus_shunts)
+
+    constraint_current_balance(pm, nw, cnd, i, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs)
+end
 
 ### Branch - Ohm's Law Constraints ###
 
@@ -277,6 +326,39 @@ function constraint_ohms_y_to(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cn
 end
 
 
+""
+function constraint_current_from(pm::AbstractIVRModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    branch = ref(pm, nw, :branch, i)
+    f_bus = branch["f_bus"]
+    t_bus = branch["t_bus"]
+    f_idx = (i, f_bus, t_bus)
+
+    tr, ti = calc_branch_t(branch)
+    g_fr = branch["g_fr"][cnd]
+    b_fr = branch["b_fr"][cnd]
+    tm = branch["tap"][cnd]
+
+    constraint_current_from(pm, nw, cnd, f_bus, f_idx, g_fr, b_fr, tr[cnd], ti[cnd], tm)
+end
+
+""
+function constraint_current_to(pm::AbstractIVRModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    branch = ref(pm, nw, :branch, i)
+    f_bus = branch["f_bus"]
+    t_bus = branch["t_bus"]
+    f_idx = (i, f_bus, t_bus)
+    t_idx = (i, t_bus, f_bus)
+
+    tr, ti = calc_branch_t(branch)
+    g_to = branch["g_to"][cnd]
+    b_to = branch["b_to"][cnd]
+    tm = branch["tap"][cnd]
+
+    constraint_current_to(pm, nw, cnd, t_bus, f_idx, t_idx, g_to, b_to)
+end
+
+
+
 ### Branch - On/Off Ohm's Law Constraints ###
 
 ""
@@ -363,6 +445,21 @@ function constraint_ohms_yt_to_ne(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw
 end
 
 
+""
+function constraint_voltage_drop(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    branch = ref(pm, nw, :branch, i)
+    f_bus = branch["f_bus"]
+    t_bus = branch["t_bus"]
+    f_idx = (i, f_bus, t_bus)
+
+    tr, ti = calc_branch_t(branch)
+    r = branch["br_r"][cnd,cnd]
+    x = branch["br_x"][cnd,cnd]
+    tm = branch["tap"][cnd]
+
+    constraint_voltage_drop(pm, nw, cnd, i, f_bus, t_bus, f_idx, r, x, tr[cnd], ti[cnd], tm)
+end
+
 
 ### Branch - Current ###
 
@@ -621,8 +718,8 @@ function constraint_flow_losses(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw,
     f_idx = (i, f_bus, t_bus)
     t_idx = (i, t_bus, f_bus)
 
-    r = branch["br_r"][cnd]
-    x = branch["br_x"][cnd]
+    r = branch["br_r"][cnd,cnd]
+    x = branch["br_x"][cnd,cnd]
     tm = branch["tap"][cnd]
     g_sh_fr = branch["g_fr"][cnd]
     g_sh_to = branch["g_to"][cnd]
@@ -640,8 +737,8 @@ function constraint_voltage_magnitude_difference(pm::AbstractPowerModel, i::Int;
     f_idx = (i, f_bus, t_bus)
     t_idx = (i, t_bus, f_bus)
 
-    r = branch["br_r"][cnd]
-    x = branch["br_x"][cnd]
+    r = branch["br_r"][cnd,cnd]
+    x = branch["br_x"][cnd,cnd]
     g_sh_fr = branch["g_fr"][cnd]
     b_sh_fr = branch["b_fr"][cnd]
     tm = branch["tap"][cnd]
@@ -760,7 +857,7 @@ function constraint_storage_on_off(pm::AbstractPowerModel, i::Int; nw::Int=pm.cn
     storage = ref(pm, nw, :storage, i)
     charge_ub = storage["charge_rating"]
     discharge_ub = storage["discharge_rating"]
-    
+
     inj_lb, inj_ub = ref_calc_storage_injection_bounds(ref(pm, nw, :storage), ref(pm, nw, :bus), cnd)
     pmin = inj_lb[i]
     pmax = inj_ub[i]
@@ -797,3 +894,33 @@ function constraint_active_dcline_setpoint(pm::AbstractPowerModel, i::Int; nw::I
 
     constraint_active_dcline_setpoint(pm, nw, cnd, f_idx, t_idx, pf, pt)
 end
+
+
+""
+function constraint_dcline_power_limits_from(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    dcline = ref(pm, nw, :dcline, i)
+    f_bus = dcline["f_bus"]
+    t_bus = dcline["t_bus"]
+    f_idx = (i, f_bus, t_bus)
+
+    pmax = dcline["pmaxf"][cnd]
+    pmin = dcline["pminf"][cnd]
+
+    qmax = dcline["qmaxf"][cnd]
+    qmin = dcline["qminf"][cnd]
+    constraint_dcline_power_limits_from(pm, nw, cnd, i, f_bus, f_idx, pmax, pmin, qmax, qmin)
+end
+
+""
+function constraint_dcline_power_limits_to(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    dcline = ref(pm, nw, :dcline, i)
+    f_bus = dcline["f_bus"]
+    t_bus = dcline["t_bus"]
+    t_idx = (i, t_bus, f_bus)
+
+    pmax = dcline["pmaxt"][cnd]
+    pmin = dcline["pmint"][cnd]
+    qmax = dcline["qmaxt"][cnd]
+    qmin = dcline["qmint"][cnd]
+    constraint_dcline_power_limits_to(pm, nw, cnd, i, t_bus, t_idx, pmax, pmin, qmax, qmin)
+end