Adding network expansion planning (#33)

Basic implementation of transmission network expansion planning.

This code should be revisited when #34 is closed and should be considered in the context of #35

Author: rb004f

src/PowerModels.jl

@@ -27,5 +27,6 @@ include("prob/pf.jl")
 include("prob/opf.jl")
 include("prob/ots.jl")
 include("prob/misc.jl")
+include("prob/tnep.jl")
 
 end

src/core/base.jl

@@ -31,32 +31,44 @@ type GenericPowerModel{T<:AbstractPowerFormulation} <: AbstractPowerModel
     set::PowerDataSets
     setting::Dict{AbstractString,Any}
     solution::Dict{AbstractString,Any}
+
+    # Extension dictionary
+    # Extensions should define a type to hold information particular to
+    # their functionality, and store an instance of the type in this
+    # dictionary keyed on an extension-specific symbol
+    ext::Dict{Symbol,Any}
 end
 
 # default generic constructor
-function GenericPowerModel{T}(data::Dict{AbstractString,Any}, vars::T; setting = Dict{AbstractString,Any}(), solver = JuMP.UnsetSolver())
-    data, sets = process_raw_data(data)
-
+function GenericPowerModel{T}(data::Dict{AbstractString,Any}, vars::T; setting = Dict{AbstractString,Any}(), solver = JuMP.UnsetSolver(), data_processor = process_raw_mp_data)
+    data, sets, ext = data_processor(data)
     pm = GenericPowerModel{T}(
         Model(solver = solver), # model
         data, # data
         sets, # sets
         setting, # setting
         Dict{AbstractString,Any}(), # solution
+        ext # ext
     )
 
     return pm
 end
 
-function process_raw_data(data::Dict{AbstractString,Any})
+function process_raw_mp_data(data::Dict{AbstractString,Any})
     make_per_unit(data)
-    unify_transformer_taps(data)
-    add_branch_parameters(data)
-    standardize_cost_order(data)
 
+    min_theta_delta = calc_min_phase_angle(data)
+    max_theta_delta = calc_max_phase_angle(data)
+
+    unify_transformer_taps(data["branch"])
+    add_branch_parameters(data["branch"], min_theta_delta, max_theta_delta)
+
+    standardize_cost_order(data)
     sets = build_sets(data)
 
-    return data, sets
+    ext = Dict{Symbol,Any}()
+
+    return data, sets, ext
 end
 
 #
@@ -143,7 +155,6 @@ function build_sets(data::Dict{AbstractString,Any})
     gen_idxs = collect(keys(gen_lookup))
     branch_idxs = collect(keys(branch_lookup))
 
-
     buspair_indexes = collect(Set([(i,j) for (l,i,j) in arcs_from]))
     buspairs = buspair_parameters(buspair_indexes, branch_lookup, bus_lookup)
 
@@ -232,20 +243,17 @@ function make_per_unit(mva_base::Number, data::Number)
     #println("$(parent) $(data)")
 end
 
-function unify_transformer_taps(data::Dict{AbstractString,Any})
-    for branch in data["branch"]
+function unify_transformer_taps(branches)
+    for branch in branches
         if branch["tap"] == 0.0
             branch["tap"] = 1.0
         end
     end
 end
 
 # NOTE, this function assumes all values are p.u. and angles are in radians
-function add_branch_parameters(data::Dict{AbstractString,Any})
-    min_theta_delta = calc_min_phase_angle(data)
-    max_theta_delta = calc_max_phase_angle(data)
-
-    for branch in data["branch"]
+function add_branch_parameters(branches, min_theta_delta, max_theta_delta)
+    for branch in branches
         r = branch["br_r"]
         x = branch["br_x"]
         tap_ratio = branch["tap"]
@@ -276,14 +284,18 @@ function calc_max_phase_angle(data::Dict{AbstractString,Any})
     bus_count = length(data["bus"])
     angle_max = [branch["angmax"] for branch in data["branch"]]
     sort!(angle_max, rev=true)
-
-    return sum(angle_max[1:bus_count-1])
+    if length(angle_max) > 1
+        return sum(angle_max[1:bus_count-1])
+    end
+    return angle_max[1]
 end
 
 function calc_min_phase_angle(data::Dict{AbstractString,Any})
     bus_count = length(data["bus"])
     angle_min = [branch["angmin"] for branch in data["branch"]]
     sort!(angle_min)
-
-    return sum(angle_min[1:bus_count-1])
+    if length(angle_min) > 1
+        return sum(angle_min[1:bus_count-1])
+    end
+    return angle_min[1]
 end

src/core/variable.jl

@@ -89,10 +89,7 @@ function variable_reactive_line_flow{T}(pm::GenericPowerModel{T}; bounded = true
     return q
 end
 
-function compute_voltage_product_bounds{T}(pm::GenericPowerModel{T})
-    buspairs = pm.set.buspairs
-    buspair_indexes = pm.set.buspair_indexes
-
+function compute_voltage_product_bounds(buspairs, buspair_indexes)
     wr_min = Dict([(bp, -Inf) for bp in buspair_indexes])
     wr_max = Dict([(bp,  Inf) for bp in buspair_indexes])
     wi_min = Dict([(bp, -Inf) for bp in buspair_indexes])
@@ -126,7 +123,7 @@ end
 
 function variable_complex_voltage_product{T}(pm::GenericPowerModel{T}; bounded = true)
     if bounded
-        wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm)
+        wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm.set.buspairs, pm.set.buspair_indexes)
 
         @variable(pm.model, wr_min[bp] <= wr[bp in pm.set.buspair_indexes] <= wr_max[bp], start = getstart(pm.set.buspairs, bp, "wr_start", 1.0))
         @variable(pm.model, wi_min[bp] <= wi[bp in pm.set.buspair_indexes] <= wi_max[bp], start = getstart(pm.set.buspairs, bp, "wi_start"))
@@ -138,18 +135,18 @@ function variable_complex_voltage_product{T}(pm::GenericPowerModel{T}; bounded =
 end
 
 function variable_complex_voltage_product_on_off{T}(pm::GenericPowerModel{T})
-    wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm)
+    wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm.set.buspairs, pm.set.buspair_indexes)
 
     bi_bp = Dict([(i, (b["f_bus"], b["t_bus"])) for (i,b) in pm.set.branches])
 
     @variable(pm.model, min(0, wr_min[bi_bp[b]]) <= wr[b in pm.set.branch_indexes] <= max(0, wr_max[bi_bp[b]]), start = getstart(pm.set.buspairs, bi_bp[b], "wr_start", 1.0))
-    @variable(pm.model, min(0, wi_min[bi_bp[b]]) <= wi[b in pm.set.branch_indexes] <= max(0, wi_max[bi_bp[b]]), start = getstart(pm.set.buspairs, bi_bp[b], "wr_start"))
+    @variable(pm.model, min(0, wi_min[bi_bp[b]]) <= wi[b in pm.set.branch_indexes] <= max(0, wi_max[bi_bp[b]]), start = getstart(pm.set.buspairs, bi_bp[b], "wi_start"))
 
     return wr, wi
 end
 
 function variable_complex_voltage_product_matrix{T}(pm::GenericPowerModel{T})
-    wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm)
+    wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm.set.buspairs, pm.set.buspair_indexes)
 
     w_index = 1:length(pm.set.bus_indexes)
     lookup_w_index = Dict([(bi, i) for (i,bi) in enumerate(pm.set.bus_indexes)])

src/form/wr.jl

@@ -224,7 +224,7 @@ function constraint_voltage_magnitude_sqr_to_on_off{T}(pm::GenericPowerModel{T})
 end
 
 function constraint_complex_voltage_product_on_off{T}(pm::GenericPowerModel{T})
-    wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm)
+    wr_min, wr_max, wi_min, wi_max = compute_voltage_product_bounds(pm.set.buspairs, pm.set.buspair_indexes)
 
     bi_bp = Dict([(i, (b["f_bus"], b["t_bus"])) for (i,b) in pm.set.branches])