Adding var to GenericPowerModel (#149)

* Added var to GenericPowerModel and updated all models to use it
* moved ACR formulation to separate file
* fixed name of dcline voltage setpoint
* updating how dc models ignore reactive power, closes #147
* removing need for model.ext from dc models
* moving WRM pm.model.ext to pm.ext
* adding extra qc ots case for extra code coverage
* minor tweak to dcline flow variable constructor
* adding changelog notes

CHANGELOG.md

@@ -2,21 +2,23 @@ PowerModels.jl Change Log
 =================
 
 ### Staged
-- Made solution default units, per unit (breaking)
+- Added pm.var and made all JuMP variables anonymous (breaking)
+- Eliminated usage of pm.model.ext, for details see [#149](https://github.com/lanl-ansi/PowerModels.jl/pull/149)
+- Made solution default units per-unit (breaking)
 - Removed deprecated bus-less constraint_theta_ref function (breaking)
 - Moved check_cost_models into the objective building function
 - Fixed out of range bug in calc_theta_delta_bounds
 
 ### v0.3.4
 - Added support for Matpower data with dclines (thanks to @frederikgeth, @hakanergun)
 - Added support for QC-OTS
-- Added support for multiple refrence buses
+- Added support for multiple reference buses
 - Added rectangular voltage formulation of AC-OPF
 - Added w-theta formulation of AC-OPF
 - Added data units checks to update_data
 - Made branch flow parameter names consistent with Matpower
 - Fixed bug in constants for w-space phase angle difference constraints
-- Fixed bug when no refrence bus was specified
+- Fixed bug when no reference bus was specified
 - Fixed dcline parsing bug
 
 ### v0.3.3

docs/src/specifications.md

@@ -124,8 +124,8 @@ for (i,branch) in pm.ref[:branch]
     constraint_ohms_yt_to(pm, branch)
 end
 for (i,dcline) in pm.ref[:dcline]
-    constraint_active_dc_line_setpoint(pm, dcline)
-    constraint_dcline_voltage(pm, dcline; epsilon = 0.00001)
+    constraint_active_dcline_setpoint(pm, dcline)
+    constraint_voltage_dcline_setpoint(pm, dcline; epsilon = 0.00001)
 end
 ```
 

src/PowerModels.jl

@@ -21,6 +21,7 @@ include("core/objective.jl")
 include("core/solution.jl")
 
 include("form/acp.jl")
+include("form/acr.jl")
 include("form/act.jl")
 include("form/dcp.jl")
 include("form/wr.jl")

src/core/base.jl

@@ -18,7 +18,9 @@ type GenericPowerModel{T<:AbstractPowerFormulation}
     data::Dict{String,Any}
     setting::Dict{String,Any}
     solution::Dict{String,Any}
+    var::Dict{Symbol,Any} # model variable lookup
     ref::Dict{Symbol,Any} # reference data
+    ext::Dict{Symbol,Any} # user extentions
 end
 ```
 where
@@ -41,7 +43,9 @@ type GenericPowerModel{T<:AbstractPowerFormulation}
     setting::Dict{String,Any}
     solution::Dict{String,Any}
 
-    ref::Dict{Symbol,Any}
+    ref::Dict{Symbol,Any} # data reference data
+
+    var::Dict{Symbol,Any} # JuMP variables
 
     # Extension dictionary
     # Extensions should define a type to hold information particular to
@@ -62,26 +66,14 @@ function GenericPowerModel(data::Dict{String,Any}, T::DataType; setting = Dict{S
         setting, # setting
         Dict{String,Any}(), # solution
         build_ref(data), # refrence data
+        Dict{Symbol,Any}(), # vars
         Dict{Symbol,Any}() # ext
     )
 
     return pm
 end
 
 
-#
-# Just seems too hard to maintain with the default constructor
-#
-#function setdata{T}(pm::GenericPowerModel{T}, data::Dict{String,Any})
-#    data, sets = process_raw_data(data)
-
-#    pm.model = Model()
-#    pm.set = sets
-#    pm.solution = Dict{String,Any}()
-#    pm.data = data
-
-#end
-
 # TODO Ask Miles, why do we need to put JuMP. here?  using at top level should bring it in
 function JuMP.setsolver(pm::GenericPowerModel, solver::MathProgBase.AbstractMathProgSolver)
     setsolver(pm.model, solver)

src/core/constraint.jl

@@ -7,32 +7,32 @@
 # Generic thermal limit constraint
 "`p[f_idx]^2 + q[f_idx]^2 <= rate_a^2`"
 function constraint_thermal_limit_from(pm::GenericPowerModel, f_idx, rate_a)
-    p_fr = getindex(pm.model, :p)[f_idx]
-    q_fr = getindex(pm.model, :q)[f_idx]
+    p_fr = pm.var[:p][f_idx]
+    q_fr = pm.var[:q][f_idx]
     c = @constraint(pm.model, p_fr^2 + q_fr^2 <= rate_a^2)
     return Set([c])
 end
 
 "`p[t_idx]^2 + q[t_idx]^2 <= rate_a^2`"
 function constraint_thermal_limit_to(pm::GenericPowerModel, t_idx, rate_a)
-    p_to = getindex(pm.model, :p)[t_idx]
-    q_to = getindex(pm.model, :q)[t_idx]
+    p_to = pm.var[:p][t_idx]
+    q_to = pm.var[:q][t_idx]
     c = @constraint(pm.model, p_to^2 + q_to^2 <= rate_a^2)
     return Set([c])
 end
 
 "`norm([p[f_idx]; q[f_idx]]) <= rate_a`"
 function constraint_thermal_limit_from{T <: AbstractConicPowerFormulation}(pm::GenericPowerModel{T}, f_idx, rate_a)
-    p_fr = getindex(pm.model, :p)[f_idx]
-    q_fr = getindex(pm.model, :q)[f_idx]
+    p_fr = pm.var[:p][f_idx]
+    q_fr = pm.var[:q][f_idx]
     c = @constraint(pm.model, norm([p_fr; q_fr]) <= rate_a)
     return Set([c])
 end
 
 "`norm([p[t_idx]; q[t_idx]]) <= rate_a`"
 function constraint_thermal_limit_to{T <: AbstractConicPowerFormulation}(pm::GenericPowerModel{T}, t_idx, rate_a)
-    p_to = getindex(pm.model, :p)[t_idx]
-    q_to = getindex(pm.model, :q)[t_idx]
+    p_to = pm.var[:p][t_idx]
+    q_to = pm.var[:q][t_idx]
     c = @constraint(pm.model, norm([p_to; q_to]) <= rate_a)
     return Set([c])
 end
@@ -41,58 +41,73 @@ end
 
 "`p[f_idx]^2 + q[f_idx]^2 <= (rate_a * line_z[i])^2`"
 function constraint_thermal_limit_from_on_off(pm::GenericPowerModel, i, f_idx, rate_a)
-    p_fr = getindex(pm.model, :p)[f_idx]
-    q_fr = getindex(pm.model, :q)[f_idx]
-    z = getindex(pm.model, :line_z)[i]
+    p_fr = pm.var[:p][f_idx]
+    q_fr = pm.var[:q][f_idx]
+    z = pm.var[:line_z][i]
     c = @constraint(pm.model, p_fr^2 + q_fr^2 <= rate_a^2*z^2)
     return Set([c])
 end
 
 "`p[t_idx]^2 + q[t_idx]^2 <= (rate_a * line_z[i])^2`"
 function constraint_thermal_limit_to_on_off(pm::GenericPowerModel, i, t_idx, rate_a)
-    p_to = getindex(pm.model, :p)[t_idx]
-    q_to = getindex(pm.model, :q)[t_idx]
-    z = getindex(pm.model, :line_z)[i]
+    p_to = pm.var[:p][t_idx]
+    q_to = pm.var[:q][t_idx]
+    z = pm.var[:line_z][i]
     c = @constraint(pm.model, p_to^2 + q_to^2 <= rate_a^2*z^2)
     return Set([c])
 end
 
 "`p_ne[f_idx]^2 + q_ne[f_idx]^2 <= (rate_a * line_ne[i])^2`"
 function constraint_thermal_limit_from_ne(pm::GenericPowerModel, i, f_idx, rate_a)
-    p_fr = getindex(pm.model, :p_ne)[f_idx]
-    q_fr = getindex(pm.model, :q_ne)[f_idx]
-    z = getindex(pm.model, :line_ne)[i]
+    p_fr = pm.var[:p_ne][f_idx]
+    q_fr = pm.var[:q_ne][f_idx]
+    z = pm.var[:line_ne][i]
     c = @constraint(pm.model, p_fr^2 + q_fr^2 <= rate_a^2*z^2)
     return Set([c])
 end
 
 "`p_ne[t_idx]^2 + q_ne[t_idx]^2 <= (rate_a * line_ne[i])^2`"
 function constraint_thermal_limit_to_ne(pm::GenericPowerModel, i, t_idx, rate_a)
-    p_to = getindex(pm.model, :p_ne)[t_idx]
-    q_to = getindex(pm.model, :q_ne)[t_idx]
-    z = getindex(pm.model, :line_ne)[i]
+    p_to = pm.var[:p_ne][t_idx]
+    q_to = pm.var[:q_ne][t_idx]
+    z = pm.var[:line_ne][i]
     c = @constraint(pm.model, p_to^2 + q_to^2 <= rate_a^2*z^2)
     return Set([c])
 end
 
 "`pg[i] == pg`"
 function constraint_active_gen_setpoint(pm::GenericPowerModel, i, pg)
-    pg_var = getindex(pm.model, :pg)[i]
+    pg_var = pm.var[:pg][i]
     c = @constraint(pm.model, pg_var == pg)
     return Set([c])
 end
 
 "`qq[i] == qq`"
 function constraint_reactive_gen_setpoint(pm::GenericPowerModel, i, qg)
-    qg_var = getindex(pm.model, :qg)[i]
+    qg_var = pm.var[:qg][i]
     c = @constraint(pm.model, qg_var == qg)
     return Set([c])
 end
 
+"""
+Creates Line Flow constraint for DC Lines (Matpower Formulation)
+
+```
+p_fr + p_to == loss0 + p_fr * loss1
+```
+"""
+function constraint_dcline{T}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, loss0, loss1)
+    p_fr = pm.var[:p_dc][f_idx]
+    p_to = pm.var[:p_dc][t_idx]
+
+    c1 = @constraint(pm.model, (1-loss1) * p_fr + (p_to - loss0) == 0)
+    return Set([c1])
+end
+
 "`pf[i] == pf, pt[i] == pt`"
 function constraint_active_dcline_setpoint(pm::GenericPowerModel, i, f_idx, t_idx, pf, pt, epsilon)
-    p_fr = getindex(pm.model, :p_dc)[f_idx]
-    p_to = getindex(pm.model, :p_dc)[t_idx]
+    p_fr = pm.var[:p_dc][f_idx]
+    p_to = pm.var[:p_dc][t_idx]
 
     if epsilon == 0.0
         c1 = @constraint(pm.model, p_fr == pf)

src/core/constraint_template.jl

@@ -91,31 +91,17 @@ function constraint_dcline(pm::GenericPowerModel, dcline)
     return constraint_dcline(pm, f_bus, t_bus, f_idx, t_idx, loss0, loss1)
 end
 
-"""
-Creates Line Flow constraint for DC Lines (Matpower Formulation)
-
-```
-p_fr + p_to == loss0 + p_fr * loss1
-```
-"""
-function constraint_dcline{T}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, loss0, loss1)
-    p_fr = getindex(pm.model, :p_dc)[f_idx]
-    p_to = getindex(pm.model, :p_dc)[t_idx]
-
-    c1 = @constraint(pm.model, (1-loss1) * p_fr + (p_to - loss0) == 0)
-    return Set([c1])
-end
 
 ""
-function constraint_dcline_voltage(pm::GenericPowerModel, dcline; epsilon = 0.0)
+function constraint_voltage_dcline_setpoint(pm::GenericPowerModel, dcline; epsilon = 0.0)
     @assert epsilon >= 0.0
     i = dcline["index"]
     f_bus = dcline["f_bus"]
     t_bus = dcline["t_bus"]
     vf = dcline["vf"]
     vt = dcline["vt"]
 
-    return constraint_dcline_voltage(pm, f_bus, t_bus, vf, vt, epsilon)
+    return constraint_voltage_dcline_setpoint(pm, f_bus, t_bus, vf, vt, epsilon)
 end
 
 function constraint_active_dcline_setpoint(pm::GenericPowerModel, dcline; epsilon = 0.0)
@@ -126,6 +112,7 @@ function constraint_active_dcline_setpoint(pm::GenericPowerModel, dcline; epsilo
     t_idx = (i, t_bus, f_bus)
     pf = dcline["pf"]
     pt = dcline["pt"]
+
     return constraint_active_dcline_setpoint(pm, i, f_idx, t_idx, pf, pt, epsilon)
 end
 

src/core/objective.jl

@@ -30,8 +30,9 @@ end
 function objective_min_fuel_cost(pm::GenericPowerModel)
     check_cost_models(pm)
 
-    pg = getindex(pm.model, :pg)
-    dc_p = getindex(pm.model, :p_dc)
+    pg = pm.var[:pg]
+    dc_p = pm.var[:p_dc]
+
     from_idx = Dict(arc[1] => arc for arc in pm.ref[:arcs_from_dc])
 
     return @objective(pm.model, Min, 
@@ -44,16 +45,24 @@ end
 function objective_min_fuel_cost{T <: AbstractConicPowerFormulation}(pm::GenericPowerModel{T})
     check_cost_models(pm)
 
-    pg = getindex(pm.model, :pg)
-    dc_p = getindex(pm.model, :p_dc)
+    pg = pm.var[:pg]
+    dc_p = pm.var[:p_dc]
     from_idx = Dict(arc[1] => arc for arc in pm.ref[:arcs_from_dc])
 
-    @variable(pm.model, pm.ref[:gen][i]["pmin"]^2 <= pg_sqr[i in keys(pm.ref[:gen])] <= pm.ref[:gen][i]["pmax"]^2)
+    pg_sqr = pm.var[:pg_sqr] = @variable(pm.model, 
+        [i in keys(pm.ref[:gen])], basename="pg_sqr",
+        lowerbound = pm.ref[:gen][i]["pmin"]^2,
+        upperbound = pm.ref[:gen][i]["pmax"]^2
+    )
     for (i, gen) in pm.ref[:gen]
         @constraint(pm.model, norm([2*pg[i], pg_sqr[i]-1]) <= pg_sqr[i]+1)
     end
 
-    @variable(pm.model, pm.ref[:dcline][i]["pminf"]^2 <= dc_p_sqr[i in keys(pm.ref[:dcline])] <= pm.ref[:dcline][i]["pmaxf"]^2)
+    dc_p_sqr = pm.var[:dc_p_sqr] = @variable(pm.model, 
+        dc_p_sqr[i in keys(pm.ref[:dcline])], basename="dc_p_sqr",
+        lowerbound = pm.ref[:dcline][i]["pminf"]^2,
+        upperbound = pm.ref[:dcline][i]["pmaxf"]^2
+    )
     for (i, dcline) in pm.ref[:dcline]
         @constraint(pm.model, norm([2*dc_p[from_idx[i]], dc_p_sqr[i]-1]) <= dc_p_sqr[i]+1)
     end
@@ -66,7 +75,7 @@ end
 
 "Cost of building lines"
 function objective_tnep_cost(pm::GenericPowerModel)
-    line_ne = getindex(pm.model, :line_ne)
+    line_ne = pm.var[:line_ne]
     branches = pm.ref[:ne_branch]
     return @objective(pm.model, Min, sum( branches[i]["construction_cost"]*line_ne[i] for (i,branch) in branches) )
 end

src/core/solution.jl

@@ -122,7 +122,7 @@ function add_setpoint(sol, pm::GenericPowerModel, dict_name, index_name, param_n
         sol_item = sol_dict[i] = get(sol_dict, i, Dict{String,Any}())
         sol_item[param_name] = default_value(item)
         try
-            var = extract_var(getindex(pm.model, variable_symbol), idx, item)
+            var = extract_var(pm.var[variable_symbol], idx, item)
             sol_item[param_name] = scale(getvalue(var), item)
         catch
         end