diff --git a/CHANGELOG.md b/CHANGELOG.md index 88497f6ec..8a065028f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,9 @@ 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 @@ -10,13 +12,13 @@ PowerModels.jl Change Log ### 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 diff --git a/docs/src/specifications.md b/docs/src/specifications.md index eed4c8209..cbd5721e1 100644 --- a/docs/src/specifications.md +++ b/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 ``` diff --git a/src/PowerModels.jl b/src/PowerModels.jl index cd6b7fa17..02f697ab9 100644 --- a/src/PowerModels.jl +++ b/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") diff --git a/src/core/base.jl b/src/core/base.jl index b73468bc0..fe2c41676 100644 --- a/src/core/base.jl +++ b/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,6 +66,7 @@ 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 ) @@ -69,19 +74,6 @@ function GenericPowerModel(data::Dict{String,Any}, T::DataType; setting = Dict{S 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) diff --git a/src/core/constraint.jl b/src/core/constraint.jl index 6d70d6eba..2a9c7d91a 100644 --- a/src/core/constraint.jl +++ b/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) diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 103c36ac8..865a00f3a 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -91,23 +91,9 @@ 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"] @@ -115,7 +101,7 @@ function constraint_dcline_voltage(pm::GenericPowerModel, dcline; epsilon = 0.0) 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 diff --git a/src/core/objective.jl b/src/core/objective.jl index ec6b2a29c..63a9d3b98 100644 --- a/src/core/objective.jl +++ b/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 diff --git a/src/core/solution.jl b/src/core/solution.jl index a00950248..5dd4926d5 100644 --- a/src/core/solution.jl +++ b/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 diff --git a/src/core/variable.jl b/src/core/variable.jl index 7b456bafc..7c35908d7 100644 --- a/src/core/variable.jl +++ b/src/core/variable.jl @@ -11,33 +11,52 @@ end "variable: `t[i]` for `i` in `bus`es" function variable_phase_angle(pm::GenericPowerModel; bounded::Bool = true) - @variable(pm.model, t[i in keys(pm.ref[:bus])], start = getstart(pm.ref[:bus], i, "t_start")) - return t + pm.var[:t] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="t", + start = getstart(pm.ref[:bus], i, "t_start") + ) + return pm.var[:t] end "variable: `v[i]` for `i` in `bus`es" function variable_voltage_magnitude(pm::GenericPowerModel; bounded = true) if bounded - @variable(pm.model, pm.ref[:bus][i]["vmin"] <= v[i in keys(pm.ref[:bus])] <= pm.ref[:bus][i]["vmax"], start = getstart(pm.ref[:bus], i, "v_start", 1.0)) + pm.var[:v] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="v", + lowerbound = pm.ref[:bus][i]["vmin"], + upperbound = pm.ref[:bus][i]["vmax"], + start = getstart(pm.ref[:bus], i, "v_start", 1.0) + ) else - @variable(pm.model, v[i in keys(pm.ref[:bus])] >= 0, start = getstart(pm.ref[:bus], i, "v_start", 1.0)) + pm.var[:v] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="v", + lowerbound = 0, + start = getstart(pm.ref[:bus], i, "v_start", 1.0)) end - return v + return pm.var[:v] end "real part of the voltage variable `i` in `bus`es" function variable_voltage_real(pm::GenericPowerModel; bounded::Bool = true) - vm_ub = Dict(i => bus["vmax"] for (i,bus) in pm.ref[:bus]) - @variable(pm.model, -vm_ub[i] <= vr[i in keys(pm.ref[:bus])] <= vm_ub[i], start = getstart(pm.ref[:bus], i, "vr_start", 1.0)) - return vr + pm.var[:vr] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="vr", + lowerbound = -pm.ref[:bus][i]["vmax"], + upperbound = pm.ref[:bus][i]["vmax"], + start = getstart(pm.ref[:bus], i, "vr_start", 1.0) + ) + return pm.var[:vr] end "real part of the voltage variable `i` in `bus`es" function variable_voltage_imaginary(pm::GenericPowerModel; bounded::Bool = true) - vm_ub = Dict(i => bus["vmax"] for (i,bus) in pm.ref[:bus]) - @variable(pm.model, -vm_ub[i] <= vi[i in keys(pm.ref[:bus])] <= vm_ub[i], start = getstart(pm.ref[:bus], i, "vi_start")) - return vi + pm.var[:vi] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="vi", + lowerbound = -pm.ref[:bus][i]["vmax"], + upperbound = pm.ref[:bus][i]["vmax"], + start = getstart(pm.ref[:bus], i, "vi_start") + ) + return pm.var[:vi] end @@ -47,9 +66,14 @@ function variable_voltage_magnitude_from_on_off(pm::GenericPowerModel) buses = pm.ref[:bus] branches = pm.ref[:branch] - @variable(pm.model, 0 <= v_from[i in keys(pm.ref[:branch])] <= buses[branches[i]["f_bus"]]["vmax"], start = getstart(pm.ref[:bus], i, "v_from_start", 1.0)) + pm.var[:v_from] = @variable(pm.model, + [i in keys(pm.ref[:branch])], basename="v_from", + lowerbound = 0, + upperbound = buses[branches[i]["f_bus"]]["vmax"], + start = getstart(pm.ref[:bus], i, "v_from_start", 1.0) + ) - return v_from + return pm.var[:v_from] end "variable: `0 <= v_to[l] <= buses[branches[l][\"t_bus\"]][\"vmax\"]` for `l` in `branch`es" @@ -57,20 +81,34 @@ function variable_voltage_magnitude_to_on_off(pm::GenericPowerModel) buses = pm.ref[:bus] branches = pm.ref[:branch] - @variable(pm.model, 0 <= v_to[i in keys(pm.ref[:branch])] <= buses[branches[i]["t_bus"]]["vmax"], start = getstart(pm.ref[:bus], i, "v_to_start", 1.0)) + pm.var[:v_to] = @variable(pm.model, + [i in keys(pm.ref[:branch])], basename="v_to", + lowerbound = 0, + upperbound = buses[branches[i]["t_bus"]]["vmax"], + start = getstart(pm.ref[:bus], i, "v_to_start", 1.0) + ) - return v_to + return pm.var[:v_to] end "variable: `w[i] >= 0` for `i` in `bus`es" function variable_voltage_magnitude_sqr(pm::GenericPowerModel; bounded = true) if bounded - @variable(pm.model, pm.ref[:bus][i]["vmin"]^2 <= w[i in keys(pm.ref[:bus])] <= pm.ref[:bus][i]["vmax"]^2, start = getstart(pm.ref[:bus], i, "w_start", 1.001)) + pm.var[:w] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="w", + lowerbound = pm.ref[:bus][i]["vmin"]^2, + upperbound = pm.ref[:bus][i]["vmax"]^2, + start = getstart(pm.ref[:bus], i, "w_start", 1.001) + ) else - @variable(pm.model, w[i in keys(pm.ref[:bus])] >= 0, start = getstart(pm.ref[:bus], i, "w_start", 1.001)) + pm.var[:w] = @variable(pm.model, + [i in keys(pm.ref[:bus])], basename="w", + lowerbound = 0, + start = getstart(pm.ref[:bus], i, "w_start", 1.001) + ) end - return w + return pm.var[:w] end "variable: `0 <= w_from[l] <= buses[branches[l][\"f_bus\"]][\"vmax\"]^2` for `l` in `branch`es" @@ -78,9 +116,14 @@ function variable_voltage_magnitude_sqr_from_on_off(pm::GenericPowerModel) buses = pm.ref[:bus] branches = pm.ref[:branch] - @variable(pm.model, 0 <= w_from[i in keys(pm.ref[:branch])] <= buses[branches[i]["f_bus"]]["vmax"]^2, start = getstart(pm.ref[:bus], i, "w_from_start", 1.001)) + pm.var[:w_from] = @variable(pm.model, + [i in keys(pm.ref[:branch])], basename="w_from", + lowerbound = 0, + upperbound = buses[branches[i]["f_bus"]]["vmax"]^2, + start = getstart(pm.ref[:bus], i, "w_from_start", 1.001) + ) - return w_from + return pm.var[:w_from] end "variable: `0 <= w_to[l] <= buses[branches[l][\"t_bus\"]][\"vmax\"]^2` for `l` in `branch`es" @@ -88,9 +131,14 @@ function variable_voltage_magnitude_sqr_to_on_off(pm::GenericPowerModel) buses = pm.ref[:bus] branches = pm.ref[:branch] - @variable(pm.model, 0 <= w_to[i in keys(pm.ref[:branch])] <= buses[branches[i]["t_bus"]]["vmax"]^2, start = getstart(pm.ref[:bus], i, "w_to_start", 1.001)) + pm.var[:w_to] = @variable(pm.model, + [i in keys(pm.ref[:branch])], basename="w_to", + lowerbound = 0, + upperbound = buses[branches[i]["t_bus"]]["vmax"]^2, + start = getstart(pm.ref[:bus], i, "w_to_start", 1.001) + ) - return w_to + return pm.var[:w_to] end @@ -99,25 +147,50 @@ function variable_voltage_product(pm::GenericPowerModel; bounded = true) if bounded wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:buspairs]) - @variable(pm.model, wr_min[bp] <= wr[bp in keys(pm.ref[:buspairs])] <= wr_max[bp], start = getstart(pm.ref[:buspairs], bp, "wr_start", 1.0)) - @variable(pm.model, wi_min[bp] <= wi[bp in keys(pm.ref[:buspairs])] <= wi_max[bp], start = getstart(pm.ref[:buspairs], bp, "wi_start")) + pm.var[:wr] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="wr", + lowerbound = wr_min[bp], + upperbound = wr_max[bp], + start = getstart(pm.ref[:buspairs], bp, "wr_start", 1.0) + ) + pm.var[:wi] = @variable(pm.model, + wi[bp in keys(pm.ref[:buspairs])], basename="wi", + lowerbound = wi_min[bp], + upperbound = wi_max[bp], + start = getstart(pm.ref[:buspairs], bp, "wi_start") + ) else - @variable(pm.model, wr[bp in keys(pm.ref[:buspairs])], start = getstart(pm.ref[:buspairs], bp, "wr_start", 1.0)) - @variable(pm.model, wi[bp in keys(pm.ref[:buspairs])], start = getstart(pm.ref[:buspairs], bp, "wi_start")) + pm.var[:wr] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="wr", + start = getstart(pm.ref[:buspairs], bp, "wr_start", 1.0) + ) + pm.var[:wi] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="wi", + start = getstart(pm.ref[:buspairs], bp, "wi_start") + ) end - return wr, wi + return pm.var[:wr], pm.var[:wi] end "" function variable_voltage_product_on_off(pm::GenericPowerModel) wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:buspairs]) - bi_bp = Dict([(i, (b["f_bus"], b["t_bus"])) for (i,b) in pm.ref[:branch]]) - @variable(pm.model, min(0, wr_min[bi_bp[b]]) <= wr[b in keys(pm.ref[:branch])] <= max(0, wr_max[bi_bp[b]]), start = getstart(pm.ref[:buspairs], bi_bp[b], "wr_start", 1.0)) - @variable(pm.model, min(0, wi_min[bi_bp[b]]) <= wi[b in keys(pm.ref[:branch])] <= max(0, wi_max[bi_bp[b]]), start = getstart(pm.ref[:buspairs], bi_bp[b], "wi_start")) - - return wr, wi + pm.var[:wr] = @variable(pm.model, + wr[b in keys(pm.ref[:branch])], basename="wr", + lowerbound = min(0, wr_min[bi_bp[b]]), + upperbound = max(0, wr_max[bi_bp[b]]), + start = getstart(pm.ref[:buspairs], bi_bp[b], "wr_start", 1.0) + ) + pm.var[:wi] = @variable(pm.model, + wi[b in keys(pm.ref[:branch])], basename="wi", + lowerbound = min(0, wi_min[bi_bp[b]]), + upperbound = max(0, wi_max[bi_bp[b]]), + start = getstart(pm.ref[:buspairs], bi_bp[b], "wi_start") + ) + + return pm.var[:wr], pm.var[:wi] end @@ -127,24 +200,41 @@ function variable_generation(pm::GenericPowerModel; kwargs...) variable_reactive_generation(pm; kwargs...) end + "variable: `pg[j]` for `j` in `gen`" function variable_active_generation(pm::GenericPowerModel; bounded = true) if bounded - @variable(pm.model, pm.ref[:gen][i]["pmin"] <= pg[i in keys(pm.ref[:gen])] <= pm.ref[:gen][i]["pmax"], start = getstart(pm.ref[:gen], i, "pg_start")) + pm.var[:pg] = @variable(pm.model, + [i in keys(pm.ref[:gen])], basename="pg", + lowerbound = pm.ref[:gen][i]["pmin"], + upperbound = pm.ref[:gen][i]["pmax"], + start = getstart(pm.ref[:gen], i, "pg_start") + ) else - @variable(pm.model, pg[i in keys(pm.ref[:gen])], start = getstart(pm.ref[:gen], i, "pg_start")) + pm.var[:pg] = @variable(pm.model, + [i in keys(pm.ref[:gen])], basename="pg", + start = getstart(pm.ref[:gen], i, "pg_start") + ) end - return pg + return pm.var[:pg] end "variable: `qq[j]` for `j` in `gen`" function variable_reactive_generation(pm::GenericPowerModel; bounded = true) if bounded - @variable(pm.model, pm.ref[:gen][i]["qmin"] <= qg[i in keys(pm.ref[:gen])] <= pm.ref[:gen][i]["qmax"], start = getstart(pm.ref[:gen], i, "qg_start")) + pm.var[:qg] = @variable(pm.model, + [i in keys(pm.ref[:gen])], basename="qg", + lowerbound = pm.ref[:gen][i]["qmin"], + upperbound = pm.ref[:gen][i]["qmax"], + start = getstart(pm.ref[:gen], i, "qg_start") + ) else - @variable(pm.model, qg[i in keys(pm.ref[:gen])], start = getstart(pm.ref[:gen], i, "qg_start")) + pm.var[:qg] = @variable(pm.model, + [i in keys(pm.ref[:gen])], basename="qg", + start = getstart(pm.ref[:gen], i, "qg_start") + ) end - return qg + return pm.var[:qg] end "" @@ -157,21 +247,37 @@ end "variable: `p[l,i,j]` for `(l,i,j)` in `arcs`" function variable_active_line_flow(pm::GenericPowerModel; bounded = true) if bounded - @variable(pm.model, -pm.ref[:branch][l]["rate_a"] <= p[(l,i,j) in pm.ref[:arcs]] <= pm.ref[:branch][l]["rate_a"], start = getstart(pm.ref[:branch], l, "p_start")) + pm.var[:p] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs]], basename="p", + lowerbound = -pm.ref[:branch][l]["rate_a"], + upperbound = pm.ref[:branch][l]["rate_a"], + start = getstart(pm.ref[:branch], l, "p_start") + ) else - @variable(pm.model, p[(l,i,j) in pm.ref[:arcs]], start = getstart(pm.ref[:branch], l, "p_start")) + pm.var[:p] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs]], basename="p", + start = getstart(pm.ref[:branch], l, "p_start") + ) end - return p + return pm.var[:p] end "variable: `q[l,i,j]` for `(l,i,j)` in `arcs`" function variable_reactive_line_flow(pm::GenericPowerModel; bounded = true) if bounded - @variable(pm.model, -pm.ref[:branch][l]["rate_a"] <= q[(l,i,j) in pm.ref[:arcs]] <= pm.ref[:branch][l]["rate_a"], start = getstart(pm.ref[:branch], l, "q_start")) + pm.var[:q] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs]], basename="q", + lowerbound = -pm.ref[:branch][l]["rate_a"], + upperbound = pm.ref[:branch][l]["rate_a"], + start = getstart(pm.ref[:branch], l, "q_start") + ) else - @variable(pm.model, q[(l,i,j) in pm.ref[:arcs]], start = getstart(pm.ref[:branch], l, "q_start")) + pm.var[:q] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs]], basename="q", + start = getstart(pm.ref[:branch], l, "q_start") + ) end - return q + return pm.var[:q] end function variable_dcline_flow(pm::GenericPowerModel; kwargs...) @@ -181,48 +287,68 @@ end "variable: `p_dc[l,i,j]` for `(l,i,j)` in `arcs_dc`" function variable_active_dcline_flow(pm::GenericPowerModel; bounded = true) - pmin = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - pref = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - pmax = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - loss0 = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - for (l,i,j) in pm.ref[:arcs_from_dc] - pmin[(l,i,j)] = pm.ref[:dcline][l]["pminf"] - pmax[(l,i,j)] = pm.ref[:dcline][l]["pmaxf"] - pmin[(l,j,i)] = pm.ref[:dcline][l]["pmint"] - pmax[(l,j,i)] = pm.ref[:dcline][l]["pmaxt"] - pref[(l,i,j)] = pm.ref[:dcline][l]["pf"] - pref[(l,j,i)] = pm.ref[:dcline][l]["pt"] - loss0[(l,i,j)] = 0 #loss completely assigned to to side as per matpower - loss0[(l,j,i)] = pm.ref[:dcline][l]["loss0"] #loss completely assigned to to side as per matpower - end - if bounded - @variable(pm.model, pmin[(l,i,j)] <= p_dc[(l,i,j) in pm.ref[:arcs_dc]] <= pmax[(l,i,j)], start = pref[(l,i,j)]) - else - @variable(pm.model, p_dc[(l,i,j) in pm.ref[:arcs_dc]], start = pref[(l,i,j)]) - end - return p_dc + # build bounds lookups based over arcs set + pmin = Dict() + pref = Dict() + pmax = Dict() + for (l,i,j) in pm.ref[:arcs_from_dc] + pmin[(l,i,j)] = pm.ref[:dcline][l]["pminf"] + pmax[(l,i,j)] = pm.ref[:dcline][l]["pmaxf"] + pref[(l,i,j)] = pm.ref[:dcline][l]["pf"] + end + for (l,i,j) in pm.ref[:arcs_to_dc] + pmin[(l,i,j)] = pm.ref[:dcline][l]["pmint"] + pmax[(l,i,j)] = pm.ref[:dcline][l]["pmaxt"] + pref[(l,i,j)] = pm.ref[:dcline][l]["pt"] + end + + if bounded + pm.var[:p_dc] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs_dc]], basename="p_dc", + lowerbound = pmin[(l,i,j)], + upperbound = pmax[(l,i,j)], + start = pref[(l,i,j)] + ) + else + pm.var[:p_dc] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs_dc]], basename="p_dc", + start = pref[(l,i,j)] + ) + end + return pm.var[:p_dc] end "variable: `q_dc[l,i,j]` for `(l,i,j)` in `arcs_dc`" function variable_reactive_dcline_flow(pm::GenericPowerModel; bounded = true) - qmin = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - qref = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - qmax = Dict([(a, 0.0) for a in pm.ref[:arcs_dc]]) - for (l,i,j) in pm.ref[:arcs_from_dc] - qmin[(l,i,j)] = pm.ref[:dcline][l]["qminf"] - qmax[(l,i,j)] = pm.ref[:dcline][l]["qmaxf"] - qmin[(l,j,i)] = pm.ref[:dcline][l]["qmint"] - qmax[(l,j,i)] = pm.ref[:dcline][l]["qmaxt"] - qref[(l,i,j)] = pm.ref[:dcline][l]["qf"] - qref[(l,j,i)] = pm.ref[:dcline][l]["qt"] + # build bounds lookups based over arcs set + qmin = Dict() + qref = Dict() + qmax = Dict() + for (l,i,j) in pm.ref[:arcs_from_dc] + qmin[(l,i,j)] = pm.ref[:dcline][l]["qminf"] + qmax[(l,i,j)] = pm.ref[:dcline][l]["qmaxf"] + qref[(l,i,j)] = pm.ref[:dcline][l]["qf"] + end + for (l,i,j) in pm.ref[:arcs_to_dc] + qmin[(l,i,j)] = pm.ref[:dcline][l]["qmint"] + qmax[(l,i,j)] = pm.ref[:dcline][l]["qmaxt"] + qref[(l,i,j)] = pm.ref[:dcline][l]["qt"] end - if bounded - @variable(pm.model, qmin[(l,i,j)] <= q_dc[(l,i,j) in pm.ref[:arcs_dc]] <= qmax[(l,i,j)], start = qref[(l,i,j)]) + if bounded + pm.var[:q_dc] = @variable(pm.model, + q_dc[(l,i,j) in pm.ref[:arcs_dc]], basename="q_dc", + lowerbound = qmin[(l,i,j)], + upperbound = qmax[(l,i,j)], + start = qref[(l,i,j)] + ) else - @variable(pm.model, q_dc[(l,i,j) in pm.ref[:arcs_dc]], start = qref[(l,i,j)]) + pm.var[:q_dc] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs_dc]], basename="q_dc", + start = qref[(l,i,j)] + ) end - return q_dc + return pm.var[:q_dc] end ################################################################## @@ -235,25 +361,47 @@ end "variable: `-ne_branch[l][\"rate_a\"] <= p_ne[l,i,j] <= ne_branch[l][\"rate_a\"]` for `(l,i,j)` in `ne_arcs`" function variable_active_line_flow_ne(pm::GenericPowerModel) - @variable(pm.model, -pm.ref[:ne_branch][l]["rate_a"] <= p_ne[(l,i,j) in pm.ref[:ne_arcs]] <= pm.ref[:ne_branch][l]["rate_a"], start = getstart(pm.ref[:ne_branch], l, "p_start")) - return p_ne + pm.var[:p_ne] = @variable(pm.model, + [(l,i,j) in pm.ref[:ne_arcs]], basename="p_ne", + lowerbound = -pm.ref[:ne_branch][l]["rate_a"], + upperbound = pm.ref[:ne_branch][l]["rate_a"], + start = getstart(pm.ref[:ne_branch], l, "p_start") + ) + return pm.var[:p_ne] end "variable: `-ne_branch[l][\"rate_a\"] <= q_ne[l,i,j] <= ne_branch[l][\"rate_a\"]` for `(l,i,j)` in `ne_arcs`" function variable_reactive_line_flow_ne(pm::GenericPowerModel) - @variable(pm.model, -pm.ref[:ne_branch][l]["rate_a"] <= q_ne[(l,i,j) in pm.ref[:ne_arcs]] <= pm.ref[:ne_branch][l]["rate_a"], start = getstart(pm.ref[:ne_branch], l, "q_start")) - return q_ne + pm.var[:q_ne] = @variable(pm.model, + q_ne[(l,i,j) in pm.ref[:ne_arcs]], basename="q_ne", + lowerbound = -pm.ref[:ne_branch][l]["rate_a"], + upperbound = pm.ref[:ne_branch][l]["rate_a"], + start = getstart(pm.ref[:ne_branch], l, "q_start") + ) + return pm.var[:q_ne] end "variable: `0 <= line_z[l] <= 1` for `l` in `branch`es" function variable_line_indicator(pm::GenericPowerModel) - @variable(pm.model, 0 <= line_z[l in keys(pm.ref[:branch])] <= 1, Int, start = getstart(pm.ref[:branch], l, "line_z_start", 1.0)) - return line_z + pm.var[:line_z] = @variable(pm.model, + [l in keys(pm.ref[:branch])], basename="line_z", + lowerbound = 0, + upperbound = 1, + category = :Int, + start = getstart(pm.ref[:branch], l, "line_z_start", 1.0) + ) + return pm.var[:line_z] end "variable: `0 <= line_ne[l] <= 1` for `l` in `branch`es" function variable_line_ne(pm::GenericPowerModel) branches = pm.ref[:ne_branch] - @variable(pm.model, 0 <= line_ne[l in keys(branches)] <= 1, Int, start = getstart(branches, l, "line_tnep_start", 1.0)) - return line_ne + pm.var[:line_ne] = @variable(pm.model, + [l in keys(branches)], basename="line_ne", + lowerbound = 0, + upperbound = 1, + category = :Int, + start = getstart(branches, l, "line_tnep_start", 1.0) + ) + return pm.var[:line_ne] end diff --git a/src/form/acp.jl b/src/form/acp.jl index c2ba1b200..71c90d9a3 100644 --- a/src/form/acp.jl +++ b/src/form/acp.jl @@ -1,6 +1,7 @@ +### polar form of the non-convex AC equations + export ACPPowerModel, StandardACPForm, - ACRPowerModel, StandardACRForm, APIACPPowerModel, APIACPForm "" @@ -31,13 +32,9 @@ constraint_voltage{T <: AbstractACPForm}(pm::GenericPowerModel{T}) = Set() "do nothing, this model does not have complex voltage constraints" constraint_voltage_ne{T <: AbstractACPForm}(pm::GenericPowerModel{T}) = nothing -"`t[ref_bus] == 0`" -constraint_theta_ref{T <: AbstractACPForm}(pm::GenericPowerModel{T}, ref_bus::Int) = - Set([@constraint(pm.model, getindex(pm.model, :t)[ref_bus] == 0)]) - "`vm - epsilon <= v[i] <= vm + epsilon`" function constraint_voltage_magnitude_setpoint{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, vm, epsilon) - v = getindex(pm.model, :v)[i] + v = pm.var[:v][i] if epsilon == 0.0 c = @constraint(pm.model, v == vm) @@ -55,9 +52,9 @@ v_from - epsilon <= v[i] <= v_from + epsilon v_to - epsilon <= v[i] <= v_to + epsilon ''' """ -function constraint_dcline_voltage{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) - v_f = getindex(pm.model, :v)[f_bus] - v_t = getindex(pm.model, :v)[t_bus] +function constraint_voltage_dcline_setpoint{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) + v_f = pm.var[:v][f_bus] + v_t = pm.var[:v][t_bus] if epsilon == 0.0 c1 = @constraint(pm.model, v_f == vf) @@ -79,13 +76,13 @@ sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[ ``` """ function constraint_kcl_shunt{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) - v = getindex(pm.model, :v)[i] - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) - p_dc = getindex(pm.model, :p_dc) - q_dc = getindex(pm.model, :q_dc) + v = pm.var[:v][i] + p = pm.var[:p] + q = pm.var[:q] + pg = pm.var[:pg] + qg = pm.var[:qg] + p_dc = pm.var[:p_dc] + q_dc = pm.var[:q_dc] c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*v^2) c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*v^2) @@ -99,15 +96,15 @@ sum(q[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) + sum(q_ne ``` """ function constraint_kcl_shunt_ne{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, pd, qd, gs, bs) - v = getindex(pm.model, :v)[i] - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - p_ne = getindex(pm.model, :p_ne) - q_ne = getindex(pm.model, :q_ne) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) - p_dc = getindex(pm.model, :p_dc) - q_dc = getindex(pm.model, :q_dc) + v = pm.var[:v][i] + p = pm.var[:p] + q = pm.var[:q] + p_ne = pm.var[:p_ne] + q_ne = pm.var[:q_ne] + pg = pm.var[:pg] + qg = pm.var[:qg] + p_dc = pm.var[:p_dc] + q_dc = pm.var[:q_dc] c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) + sum(p_ne[a] for a in bus_arcs_ne) == sum(pg[g] for g in bus_gens) - pd - gs*v^2) c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) + sum(q_ne[a] for a in bus_arcs_ne) == sum(qg[g] for g in bus_gens) - qd + bs*v^2) @@ -123,12 +120,12 @@ q[f_idx] == -(b+c/2)/tm*v[f_bus]^2 - (-b*tr-g*ti)/tm*(v[f_bus]*v[t_bus]*cos(t[f_ ``` """ function constraint_ohms_yt_from{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] c1 = @NLconstraint(pm.model, p_fr == g/tm*v_fr^2 + (-g*tr+b*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-b*tr-g*ti)/tm*(v_fr*v_to*sin(t_fr-t_to)) ) c2 = @NLconstraint(pm.model, q_fr == -(b+c/2)/tm*v_fr^2 - (-b*tr-g*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-g*tr+b*ti)/tm*(v_fr*v_to*sin(t_fr-t_to)) ) @@ -144,12 +141,12 @@ q[t_idx] == -(b+c/2)*v[t_bus]^2 - (-b*tr+g*ti)/tm*(v[t_bus]*v[f_bus]*cos(t[f_bus ``` """ function constraint_ohms_yt_to{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - p_to = getindex(pm.model, :p)[t_idx] - q_to = getindex(pm.model, :q)[t_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] + p_to = pm.var[:p][t_idx] + q_to = pm.var[:q][t_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] c1 = @NLconstraint(pm.model, p_to == g*v_to^2 + (-g*tr-b*ti)/tm*(v_to*v_fr*cos(t_to-t_fr)) + (-b*tr+g*ti)/tm*(v_to*v_fr*sin(t_to-t_fr)) ) c2 = @NLconstraint(pm.model, q_to == -(b+c/2)*v_to^2 - (-b*tr+g*ti)/tm*(v_to*v_fr*cos(t_fr-t_to)) + (-g*tr-b*ti)/tm*(v_to*v_fr*sin(t_to-t_fr)) ) @@ -165,12 +162,12 @@ q[f_idx] == -(b+c/2)*(v[f_bus]/tr)^2 + b*v[f_bus]/tr*v[t_bus]*cos(t[f_bus]-t[t_b ``` """ function constraint_ohms_y_from{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, as) - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] c1 = @NLconstraint(pm.model, p_fr == g*(v_fr/tr)^2 + -g*v_fr/tr*v_to*cos(t_fr-t_to-as) + -b*v_fr/tr*v_to*sin(t_fr-t_to-as) ) c2 = @NLconstraint(pm.model, q_fr == -(b+c/2)*(v_fr/tr)^2 + b*v_fr/tr*v_to*cos(t_fr-t_to-as) + -g*v_fr/tr*v_to*sin(t_fr-t_to-as) ) @@ -186,32 +183,18 @@ q_to == -(b+c/2)*v[t_bus]^2 + b*v[t_bus]*v[f_bus]/tr*cos(t[f_bus]-t[t_bus]+as) + ``` """ function constraint_ohms_y_to{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, as) - p_to = getindex(pm.model, :p)[t_idx] - q_to = getindex(pm.model, :q)[t_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] + p_to = pm.var[:p][t_idx] + q_to = pm.var[:q][t_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] c1 = @NLconstraint(pm.model, p_to == g*v_to^2 + -g*v_to*v_fr/tr*cos(t_to-t_fr+as) + -b*v_to*v_fr/tr*sin(t_to-t_fr+as) ) c2 = @NLconstraint(pm.model, q_to == -(b+c/2)*v_to^2 + b*v_to*v_fr/tr*cos(t_fr-t_to+as) + -g*v_to*v_fr/tr*sin(t_to-t_fr+as) ) return Set([c1, c2]) end -""" -``` -t[f_bus] - t[t_bus] <= angmax -t[f_bus] - t[t_bus] >= angmin -``` -""" -function constraint_phase_angle_difference{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - - c1 = @constraint(pm.model, t_fr - t_to <= angmax) - c2 = @constraint(pm.model, t_fr - t_to >= angmin) - return Set([c1, c2]) -end "" function variable_voltage_on_off{T <: AbstractACPForm}(pm::GenericPowerModel{T}; kwargs...) @@ -229,13 +212,13 @@ q[f_idx] == z*(-(b+c/2)/tm*v[f_bus]^2 - (-b*tr-g*ti)/tm*(v[f_bus]*v[t_bus]*cos(t ``` """ function constraint_ohms_yt_from_on_off{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] c1 = @NLconstraint(pm.model, p_fr == z*(g/tm*v_fr^2 + (-g*tr+b*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-b*tr-g*ti)/tm*(v_fr*v_to*sin(t_fr-t_to))) ) c2 = @NLconstraint(pm.model, q_fr == z*(-(b+c/2)/tm*v_fr^2 - (-b*tr-g*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-g*tr+b*ti)/tm*(v_fr*v_to*sin(t_fr-t_to))) ) @@ -249,13 +232,13 @@ q[t_idx] == z*(-(b+c/2)*v[t_bus]^2 - (-b*tr+g*ti)/tm*(v[t_bus]*v[f_bus]*cos(t[f_ ``` """ function constraint_ohms_yt_to_on_off{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_to = getindex(pm.model, :p)[t_idx] - q_to = getindex(pm.model, :q)[t_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + p_to = pm.var[:p][t_idx] + q_to = pm.var[:q][t_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] c1 = @NLconstraint(pm.model, p_to == z*(g*v_to^2 + (-g*tr-b*ti)/tm*(v_to*v_fr*cos(t_to-t_fr)) + (-b*tr+g*ti)/tm*(v_to*v_fr*sin(t_to-t_fr))) ) c2 = @NLconstraint(pm.model, q_to == z*(-(b+c/2)*v_to^2 - (-b*tr+g*ti)/tm*(v_to*v_fr*cos(t_fr-t_to)) + (-g*tr-b*ti)/tm*(v_to*v_fr*sin(t_to-t_fr))) ) @@ -269,13 +252,13 @@ q_ne[f_idx] == z*(-(b+c/2)/tm*v[f_bus]^2 - (-b*tr-g*ti)/tm*(v[f_bus]*v[t_bus]*co ``` """ function constraint_ohms_yt_from_ne{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p_ne)[f_idx] - q_fr = getindex(pm.model, :q_ne)[f_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + p_fr = pm.var[:p_ne][f_idx] + q_fr = pm.var[:q_ne][f_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] c1 = @NLconstraint(pm.model, p_fr == z*(g/tm*v_fr^2 + (-g*tr+b*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-b*tr-g*ti)/tm*(v_fr*v_to*sin(t_fr-t_to))) ) c2 = @NLconstraint(pm.model, q_fr == z*(-(b+c/2)/tm*v_fr^2 - (-b*tr-g*ti)/tm*(v_fr*v_to*cos(t_fr-t_to)) + (-g*tr+b*ti)/tm*(v_fr*v_to*sin(t_fr-t_to))) ) @@ -289,13 +272,13 @@ q_ne[t_idx] == z*(-(b+c/2)*v[t_bus]^2 - (-b*tr+g*ti)/tm*(v[t_bus]*v[f_bus]*cos(t ``` """ function constraint_ohms_yt_to_ne{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_to = getindex(pm.model, :p_ne)[t_idx] - q_to = getindex(pm.model, :q_ne)[t_idx] - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + p_to = pm.var[:p_ne][t_idx] + q_to = pm.var[:q_ne][t_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] c1 = @NLconstraint(pm.model, p_to == z*(g*v_to^2 + (-g*tr-b*ti)/tm*(v_to*v_fr*cos(t_to-t_fr)) + (-b*tr+g*ti)/tm*(v_to*v_fr*sin(t_to-t_fr))) ) c2 = @NLconstraint(pm.model, q_to == z*(-(b+c/2)*v_to^2 - (-b*tr+g*ti)/tm*(v_to*v_fr*cos(t_fr-t_to)) + (-g*tr-b*ti)/tm*(v_to*v_fr*sin(t_to-t_fr))) ) @@ -304,9 +287,9 @@ end "`angmin <= line_z[i]*(t[f_bus] - t[t_bus]) <= angmax`" function constraint_phase_angle_difference_on_off{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, angmin, angmax, t_min, t_max) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] c1 = @constraint(pm.model, z*(t_fr - t_to) <= angmax) c2 = @constraint(pm.model, z*(t_fr - t_to) >= angmin) @@ -315,9 +298,9 @@ end "`angmin <= line_ne[i]*(t[f_bus] - t[t_bus]) <= angmax`" function constraint_phase_angle_difference_ne{T <: AbstractACPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, angmin, angmax, t_min, t_max) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] c1 = @constraint(pm.model, z*(t_fr - t_to) <= angmax) c2 = @constraint(pm.model, z*(t_fr - t_to) >= angmin) @@ -331,12 +314,12 @@ q[f_idx] + q[t_idx] >= -c/2*(v[f_bus]^2/tr^2 + v[t_bus]^2) ``` """ function constraint_loss_lb{T <: AbstractACPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, c, tr) - v_fr = getindex(pm.model, :v)[f_bus] - v_to = getindex(pm.model, :v)[t_bus] - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - p_to = getindex(pm.model, :p)[t_idx] - q_to = getindex(pm.model, :q)[t_idx] + v_fr = pm.var[:v][f_bus] + v_to = pm.var[:v][t_bus] + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + p_to = pm.var[:p][t_idx] + q_to = pm.var[:q][t_idx] c1 = @constraint(m, p_fr + p_to >= 0) c2 = @constraint(m, q_fr + q_to >= -c/2*(v_fr^2/tr^2 + v_to^2)) @@ -345,167 +328,6 @@ end - -"" -@compat abstract type AbstractACRForm <: AbstractPowerFormulation end - -"" -@compat abstract type StandardACRForm <: AbstractACRForm end - -"" -const ACRPowerModel = GenericPowerModel{StandardACRForm} - -"default rectangular AC constructor" -ACRPowerModel(data::Dict{String,Any}; kwargs...) = - GenericPowerModel(data, StandardACRForm; kwargs...) - - -"" -function variable_voltage{T <: AbstractACRForm}(pm::GenericPowerModel{T}; kwargs...) - variable_voltage_real(pm; kwargs...) - variable_voltage_imaginary(pm; kwargs...) -end - - -"add constraints for voltage magnitude" -function constraint_voltage{T <: AbstractACRForm}(pm::GenericPowerModel{T}; kwargs...) - vr = getindex(pm.model, :vr) - vi = getindex(pm.model, :vi) - - cs = Set([]) - for (i,bus) in pm.ref[:bus] - c1 = @constraint(pm.model, bus["vmin"]^2 <= (vr[i]^2 + vi[i]^2)) - c2 = @constraint(pm.model, bus["vmax"]^2 >= (vr[i]^2 + vi[i]^2)) - push!(cs, Set([c1, c2])) - end - - # does not seem to improve convergence - #wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:buspairs]) - #for bp in keys(pm.ref[:buspairs]) - # i,j = bp - # c1 = @constraint(pm.model, wr_min[bp] <= vr[i]*vr[j] + vi[i]*vi[j]) - # c2 = @constraint(pm.model, wr_max[bp] >= vr[i]*vr[j] + vi[i]*vi[j]) - # - # c3 = @constraint(pm.model, wi_min[bp] <= vi[i]*vr[j] - vr[i]*vi[j]) - # c4 = @constraint(pm.model, wi_max[bp] >= vi[i]*vr[j] - vr[i]*vi[j]) - # - # push!(cs, Set([c1, c2, c3, c4])) - #end - - return cs -end - - -"reference bus angle constraint" -function constraint_theta_ref{T <: AbstractACRForm}(pm::GenericPowerModel{T}, ref_bus::Int) - vi = getindex(pm.model, :vi) - c = @constraint(pm.model, vi[ref_bus] == 0) - return Set([c]) -end - - -function constraint_kcl_shunt{T <: AbstractACRForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) - vr = getindex(pm.model, :vr)[i] - vi = getindex(pm.model, :vi)[i] - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) - p_dc = getindex(pm.model, :p_dc) - q_dc = getindex(pm.model, :q_dc) - - c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*(vr^2 + vi^2)) - c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*(vr^2 + vi^2)) - return Set([c1, c2]) -end - - -""" -Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form) -""" -function constraint_ohms_yt_from{T <: AbstractACRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - vr_fr = getindex(pm.model, :vr)[f_bus] - vr_to = getindex(pm.model, :vr)[t_bus] - vi_fr = getindex(pm.model, :vi)[f_bus] - vi_to = getindex(pm.model, :vi)[t_bus] - - c1 = @NLconstraint(pm.model, p_fr == g/tm*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm*(vi_fr*vr_to - vr_fr*vi_to) ) - c2 = @NLconstraint(pm.model, q_fr == -(b+c/2)/tm*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm*(vi_fr*vr_to - vr_fr*vi_to) ) - return Set([c1, c2]) -end - -""" -Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form) -""" -function constraint_ohms_yt_to{T <: AbstractACRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - p_to = getindex(pm.model, :p)[t_idx] - q_to = getindex(pm.model, :q)[t_idx] - vr_fr = getindex(pm.model, :vr)[f_bus] - vr_to = getindex(pm.model, :vr)[t_bus] - vi_fr = getindex(pm.model, :vi)[f_bus] - vi_to = getindex(pm.model, :vi)[t_bus] - - c1 = @NLconstraint(pm.model, p_to == g*(vr_to^2 + vi_to^2) + (-g*tr-b*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr+g*ti)/tm*(-(vi_fr*vr_to - vr_fr*vi_to)) ) - c2 = @NLconstraint(pm.model, q_to == -(b+c/2)*(vr_to^2 + vi_to^2) - (-b*tr+g*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr-b*ti)/tm*(-(vi_fr*vr_to - vr_fr*vi_to)) ) - return Set([c1, c2]) -end - - -""" -branch phase angle difference bounds -""" -function constraint_phase_angle_difference{T <: AbstractACRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) - vr_fr = getindex(pm.model, :vr)[f_bus] - vr_to = getindex(pm.model, :vr)[t_bus] - vi_fr = getindex(pm.model, :vi)[f_bus] - vi_to = getindex(pm.model, :vi)[t_bus] - - # this form appears to be more numerically stable than the one below - c1 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to)/(vr_fr*vr_to + vi_fr*vi_to) <= tan(angmax)) - c2 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to)/(vr_fr*vr_to + vi_fr*vi_to) >= tan(angmin)) - - #c1 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to) <= tan(angmax)*(vr_fr*vr_to + vi_fr*vi_to)) - #c2 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to) >= tan(angmin)*(vr_fr*vr_to + vi_fr*vi_to)) - - return Set([c1, c2]) -end - - -"extracts voltage set points from rectangular voltage form and converts into polar voltage form" -function add_bus_voltage_setpoint{T <: AbstractACRForm}(sol, pm::GenericPowerModel{T}) - sol_dict = sol["bus"] = get(sol, "bus", Dict{String,Any}()) - for (i,item) in pm.data["bus"] - idx = Int(item["bus_i"]) - sol_item = sol_dict[i] = get(sol_dict, i, Dict{String,Any}()) - sol_item["vm"] = NaN - sol_item["va"] = NaN - try - vr = getvalue(getindex(pm.model, :vr)[idx]) - vi = getvalue(getindex(pm.model, :vi)[idx]) - - vm = sqrt(vr^2 + vi^2) - sol_item["vm"] = vm - - if vr == 0.0 - if vi >= 0 - va = pi/2 - else - va = 3*pi/2 - end - else - va = atan(vi/vr) - end - sol_item["va"] = va - catch - end - end -end - - - - "" @compat abstract type APIACPForm <: AbstractACPForm end @@ -517,32 +339,38 @@ APIACPPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, APIACPForm; kwargs...) "variable: load_factor >= 1.0" -variable_load_factor(pm::GenericPowerModel) = - @variable(pm.model, load_factor >= 1.0, start = 1.0) +function variable_load_factor(pm::GenericPowerModel) + pm.var[:load_factor] = @variable(pm.model, + basename="load_factor", + lowerbound=1.0, + start = 1.0 + ) + return pm.var[:load_factor] +end "objective: Max. load_factor" objective_max_loading(pm::GenericPowerModel) = - @objective(pm.model, Max, getindex(pm.model, :load_factor)) + @objective(pm.model, Max, pm.var[:load_factor]) "" function objective_max_loading_voltage_norm(pm::GenericPowerModel) # Seems to create too much reactive power and makes even small models hard to converge - load_factor = getindex(pm.model, :load_factor) + load_factor = pm.var[:load_factor] scale = length(pm.ref[:bus]) - v = getindex(pm.model, :v) + v = pm.var[:v] - return @objective(pm.model, Max, 10*scale*load_factor - sum(((bus["vmin"] + bus["vmax"])/2 - v[i])^2 for (i,bus) in pm.ref[:bus] )) + return @objective(pm.model, Max, 10*scale*load_factor - sum(((bus["vmin"] + bus["vmax"])/2 - v[i])^2 for (i,bus) in pm.ref[:bus])) end "" function objective_max_loading_gen_output(pm::GenericPowerModel) # Works but adds unnecessary runtime - load_factor = getindex(pm.model, :load_factor) + load_factor = pm.var[:load_factor] scale = length(pm.ref[:gen]) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) + pg = pm.var[:pg] + qg = pm.var[:qg] return @NLobjective(pm.model, Max, 100*scale*load_factor - sum( (pg[i]^2 - (2*qg[i])^2)^2 for (i,gen) in pm.ref[:gen] )) end @@ -550,7 +378,7 @@ end "" function bounds_tighten_voltage(pm::APIACPPowerModel; epsilon = 0.001) for (i,bus) in pm.ref[:bus] - v = getindex(pm.model, :v)[i] + v = pm.var[:v][i] setupperbound(v, bus["vmax"]*(1.0-epsilon)) setlowerbound(v, bus["vmin"]*(1.0+epsilon)) end @@ -560,7 +388,7 @@ end function upperbound_negative_active_generation(pm::APIACPPowerModel) for (i,gen) in pm.ref[:gen] if gen["pmax"] <= 0 - pg = getindex(pm.model, :pg)[i] + pg = pm.var[:pg][i] setupperbound(pg, gen["pmax"]) end end @@ -572,12 +400,12 @@ function constraint_kcl_shunt_scaled(pm::APIACPPowerModel, bus) bus_arcs = pm.ref[:bus_arcs][i] bus_gens = pm.ref[:bus_gens][i] - load_factor = getindex(pm.model, :load_factor) - v = getindex(pm.model, :v) - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) + load_factor = pm.var[:load_factor] + v = pm.var[:v] + p = pm.var[:p] + q = pm.var[:q] + pg = pm.var[:pg] + qg = pm.var[:qg] if bus["pd"] > 0 && bus["qd"] > 0 c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) == sum(pg[g] for g in bus_gens) - bus["pd"]*load_factor - bus["gs"]*v[i]^2) diff --git a/src/form/acr.jl b/src/form/acr.jl new file mode 100644 index 000000000..28a29b016 --- /dev/null +++ b/src/form/acr.jl @@ -0,0 +1,161 @@ +### rectangular form of the non-convex AC equations + +export + ACRPowerModel, StandardACRForm + +"" +@compat abstract type AbstractACRForm <: AbstractPowerFormulation end + +"" +@compat abstract type StandardACRForm <: AbstractACRForm end + +"" +const ACRPowerModel = GenericPowerModel{StandardACRForm} + +"default rectangular AC constructor" +ACRPowerModel(data::Dict{String,Any}; kwargs...) = + GenericPowerModel(data, StandardACRForm; kwargs...) + + +"" +function variable_voltage{T <: AbstractACRForm}(pm::GenericPowerModel{T}; kwargs...) + variable_voltage_real(pm; kwargs...) + variable_voltage_imaginary(pm; kwargs...) +end + + +"add constraints for voltage magnitude" +function constraint_voltage{T <: AbstractACRForm}(pm::GenericPowerModel{T}; kwargs...) + vr = pm.var[:vr] + vi = pm.var[:vi] + + cs = Set([]) + for (i,bus) in pm.ref[:bus] + c1 = @constraint(pm.model, bus["vmin"]^2 <= (vr[i]^2 + vi[i]^2)) + c2 = @constraint(pm.model, bus["vmax"]^2 >= (vr[i]^2 + vi[i]^2)) + push!(cs, Set([c1, c2])) + end + + # does not seem to improve convergence + #wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:buspairs]) + #for bp in keys(pm.ref[:buspairs]) + # i,j = bp + # c1 = @constraint(pm.model, wr_min[bp] <= vr[i]*vr[j] + vi[i]*vi[j]) + # c2 = @constraint(pm.model, wr_max[bp] >= vr[i]*vr[j] + vi[i]*vi[j]) + # + # c3 = @constraint(pm.model, wi_min[bp] <= vi[i]*vr[j] - vr[i]*vi[j]) + # c4 = @constraint(pm.model, wi_max[bp] >= vi[i]*vr[j] - vr[i]*vi[j]) + # + # push!(cs, Set([c1, c2, c3, c4])) + #end + + return cs +end + + +"reference bus angle constraint" +function constraint_theta_ref{T <: AbstractACRForm}(pm::GenericPowerModel{T}, ref_bus::Int) + c = @constraint(pm.model, pm.var[:vi][ref_bus] == 0) + return Set([c]) +end + + +function constraint_kcl_shunt{T <: AbstractACRForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) + vr = pm.var[:vr][i] + vi = pm.var[:vi][i] + p = pm.var[:p] + q = pm.var[:q] + pg = pm.var[:pg] + qg = pm.var[:qg] + p_dc = pm.var[:p_dc] + q_dc = pm.var[:q_dc] + + c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*(vr^2 + vi^2)) + c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*(vr^2 + vi^2)) + return Set([c1, c2]) +end + + +""" +Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form) +""" +function constraint_ohms_yt_from{T <: AbstractACRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + vr_fr = pm.var[:vr][f_bus] + vr_to = pm.var[:vr][t_bus] + vi_fr = pm.var[:vi][f_bus] + vi_to = pm.var[:vi][t_bus] + + c1 = @NLconstraint(pm.model, p_fr == g/tm*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm*(vi_fr*vr_to - vr_fr*vi_to) ) + c2 = @NLconstraint(pm.model, q_fr == -(b+c/2)/tm*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm*(vi_fr*vr_to - vr_fr*vi_to) ) + return Set([c1, c2]) +end + +""" +Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form) +""" +function constraint_ohms_yt_to{T <: AbstractACRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) + p_to = pm.var[:p][t_idx] + q_to = pm.var[:q][t_idx] + vr_fr = pm.var[:vr][f_bus] + vr_to = pm.var[:vr][t_bus] + vi_fr = pm.var[:vi][f_bus] + vi_to = pm.var[:vi][t_bus] + + c1 = @NLconstraint(pm.model, p_to == g*(vr_to^2 + vi_to^2) + (-g*tr-b*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr+g*ti)/tm*(-(vi_fr*vr_to - vr_fr*vi_to)) ) + c2 = @NLconstraint(pm.model, q_to == -(b+c/2)*(vr_to^2 + vi_to^2) - (-b*tr+g*ti)/tm*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr-b*ti)/tm*(-(vi_fr*vr_to - vr_fr*vi_to)) ) + return Set([c1, c2]) +end + + +""" +branch phase angle difference bounds +""" +function constraint_phase_angle_difference{T <: AbstractACRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) + vr_fr = pm.var[:vr][f_bus] + vr_to = pm.var[:vr][t_bus] + vi_fr = pm.var[:vi][f_bus] + vi_to = pm.var[:vi][t_bus] + + # this form appears to be more numerically stable than the one below + c1 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to)/(vr_fr*vr_to + vi_fr*vi_to) <= tan(angmax)) + c2 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to)/(vr_fr*vr_to + vi_fr*vi_to) >= tan(angmin)) + + #c1 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to) <= tan(angmax)*(vr_fr*vr_to + vi_fr*vi_to)) + #c2 = @NLconstraint(pm.model, (vi_fr*vr_to - vr_fr*vi_to) >= tan(angmin)*(vr_fr*vr_to + vi_fr*vi_to)) + + return Set([c1, c2]) +end + + +"extracts voltage set points from rectangular voltage form and converts into polar voltage form" +function add_bus_voltage_setpoint{T <: AbstractACRForm}(sol, pm::GenericPowerModel{T}) + sol_dict = sol["bus"] = get(sol, "bus", Dict{String,Any}()) + for (i,item) in pm.data["bus"] + idx = Int(item["bus_i"]) + sol_item = sol_dict[i] = get(sol_dict, i, Dict{String,Any}()) + sol_item["vm"] = NaN + sol_item["va"] = NaN + try + vr = getvalue(pm.var[:vr][idx]) + vi = getvalue(pm.var[:vi][idx]) + + vm = sqrt(vr^2 + vi^2) + sol_item["vm"] = vm + + if vr == 0.0 + if vi >= 0 + va = pi/2 + else + va = 3*pi/2 + end + else + va = atan(vi/vr) + end + sol_item["va"] = va + catch + end + end +end + diff --git a/src/form/act.jl b/src/form/act.jl index 76cd4200a..9cb8ff73c 100644 --- a/src/form/act.jl +++ b/src/form/act.jl @@ -1,4 +1,4 @@ -### W-Theta form non-convex AC equations +### w-theta form of the non-convex AC equations export ACTPowerModel, StandardACTForm @@ -24,10 +24,10 @@ function variable_voltage{T <: AbstractACTForm}(pm::GenericPowerModel{T}; kwargs end function constraint_voltage{T <: StandardACTForm}(pm::GenericPowerModel{T}) - t = getindex(pm.model, :t) - w = getindex(pm.model, :w) - wr = getindex(pm.model, :wr) - wi = getindex(pm.model, :wi) + t = pm.var[:t] + w = pm.var[:w] + wr = pm.var[:wr] + wi = pm.var[:wi] for (i,j) in keys(pm.ref[:buspairs]) @NLconstraint(pm.model, wr[(i,j)]^2 + wi[(i,j)]^2 == w[i]*w[j]) @@ -35,25 +35,6 @@ function constraint_voltage{T <: StandardACTForm}(pm::GenericPowerModel{T}) end end -"`t[ref_bus] == 0`" -constraint_theta_ref{T <: AbstractACTForm}(pm::GenericPowerModel{T}, ref_bus::Int) = - Set([@constraint(pm.model, getindex(pm.model, :t)[ref_bus] == 0)]) - - -""" -``` -t[f_bus] - t[t_bus] <= angmax -t[f_bus] - t[t_bus] >= angmin -``` -""" -function constraint_phase_angle_difference{T <: AbstractACTForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - - c1 = @constraint(pm.model, t_fr - t_to <= angmax) - c2 = @constraint(pm.model, t_fr - t_to >= angmin) - return Set([c1, c2]) -end "" function add_bus_voltage_setpoint{T <: AbstractACTForm}(sol, pm::GenericPowerModel{T}) diff --git a/src/form/dcp.jl b/src/form/dcp.jl index 6e87939a9..c17c2e3be 100644 --- a/src/form/dcp.jl +++ b/src/form/dcp.jl @@ -19,49 +19,58 @@ DCPPowerModel(data::Dict{String,Any}; kwargs...) = variable_voltage{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) = variable_phase_angle(pm; kwargs...) -"" -variable_voltage_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) = nothing +"nothing to add, there are no voltage variables on branches" +variable_voltage_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) = Set([]) -"" -function variable_generation{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) - variable_active_generation(pm; kwargs...) - # omit reactive variables -end +"dc models ignore reactive power flows" +variable_reactive_generation{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; bounded = true) = Set() -"" -function variable_line_flow{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) - variable_active_line_flow(pm; kwargs...) - # omit reactive variables -end +"dc models ignore reactive power flows" +variable_reactive_line_flow{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; bounded = true) = Set() + +"dc models ignore reactive power flows" +variable_reactive_line_flow_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}) = Set([]) -"" -function variable_line_flow_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) - # do nothing, this model does not have reactive variables - variable_active_line_flow_ne(pm; kwargs...) -end "" function variable_active_line_flow{T <: StandardDCPForm}(pm::GenericPowerModel{T}; bounded = true) if bounded - @variable(pm.model, -pm.ref[:branch][l]["rate_a"] <= p[(l,i,j) in pm.ref[:arcs_from]] <= pm.ref[:branch][l]["rate_a"], start = getstart(pm.ref[:branch], l, "p_start")) + pm.var[:p] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs_from]], basename="p", + lowerbound = -pm.ref[:branch][l]["rate_a"], + upperbound = pm.ref[:branch][l]["rate_a"], + start = getstart(pm.ref[:branch], l, "p_start") + ) else - @variable(pm.model, p[(l,i,j) in pm.ref[:arcs_from]], start = getstart(pm.ref[:branch], l, "p_start")) + pm.var[:p] = @variable(pm.model, + [(l,i,j) in pm.ref[:arcs_from]], basename="p", + start = getstart(pm.ref[:branch], l, "p_start") + ) end - p_expr = Dict([((l,i,j), 1.0*p[(l,i,j)]) for (l,i,j) in pm.ref[:arcs_from]]) - p_expr = merge(p_expr, Dict([((l,j,i), -1.0*p[(l,i,j)]) for (l,i,j) in pm.ref[:arcs_from]])) + # this explicit type erasure is necessary + p_expr = Dict{Any,Any}([((l,i,j), pm.var[:p][(l,i,j)]) for (l,i,j) in pm.ref[:arcs_from]]) + p_expr = merge(p_expr, Dict([((l,j,i), -1.0*pm.var[:p][(l,i,j)]) for (l,i,j) in pm.ref[:arcs_from]])) + pm.var[:p] = p_expr - pm.model.ext[:p_expr] = p_expr + return pm.var[:p] end "" function variable_active_line_flow_ne{T <: StandardDCPForm}(pm::GenericPowerModel{T}) - @variable(pm.model, -pm.ref[:ne_branch][l]["rate_a"] <= p_ne[(l,i,j) in pm.ref[:ne_arcs_from]] <= pm.ref[:ne_branch][l]["rate_a"], start = getstart(pm.ref[:ne_branch], l, "p_start")) - - p_ne_expr = Dict([((l,i,j), 1.0*p_ne[(l,i,j)]) for (l,i,j) in pm.ref[:ne_arcs_from]]) - p_ne_expr = merge(p_ne_expr, Dict([((l,j,i), -1.0*p_ne[(l,i,j)]) for (l,i,j) in pm.ref[:ne_arcs_from]])) - - pm.model.ext[:p_ne_expr] = p_ne_expr + pm.var[:p_ne] = @variable(pm.model, + [(l,i,j) in pm.ref[:ne_arcs_from]], basename="p_ne", + lowerbound = -pm.ref[:ne_branch][l]["rate_a"], + upperbound = pm.ref[:ne_branch][l]["rate_a"], + start = getstart(pm.ref[:ne_branch], l, "p_start") + ) + + # this explicit type erasure is necessary + p_ne_expr = Dict{Any,Any}([((l,i,j), 1.0*pm.var[:p_ne][(l,i,j)]) for (l,i,j) in pm.ref[:ne_arcs_from]]) + p_ne_expr = merge(p_ne_expr, Dict([((l,j,i), -1.0*pm.var[:p_ne][(l,i,j)]) for (l,i,j) in pm.ref[:ne_arcs_from]])) + pm.var[:p_ne] = p_ne_expr + + return pm.var[:p_ne] end "do nothing, this model does not have complex voltage variables" @@ -70,10 +79,6 @@ constraint_voltage{T <: AbstractDCPForm}(pm::GenericPowerModel{T}) = nothing "do nothing, this model does not have complex voltage variables" constraint_voltage_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}) = nothing -"" -constraint_theta_ref{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, ref_bus::Int) = - Set([@constraint(pm.model, getindex(pm.model, :t)[ref_bus] == 0)]) - "do nothing, this model does not have voltage variables" constraint_voltage_magnitude_setpoint{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, vm, epsilon) = Set() @@ -81,26 +86,27 @@ constraint_voltage_magnitude_setpoint{T <: AbstractDCPForm}(pm::GenericPowerMode constraint_reactive_gen_setpoint{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, qg) = Set() "do nothing, this model does not have voltage variables" -constraint_dcline_voltage{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) = Set() +constraint_voltage_dcline_setpoint{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) = Set() "" function constraint_kcl_shunt{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) - pg = getindex(pm.model, :pg) - p_expr = pm.model.ext[:p_expr] - p_dc = getindex(pm.model, :p_dc) + pg = pm.var[:pg] + p = pm.var[:p] + p_dc = pm.var[:p_dc] - c = @constraint(pm.model, sum(p_expr[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2) + c = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2) # omit reactive constraint return Set([c]) end "" function constraint_kcl_shunt_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, pd, qd, gs, bs) - pg = getindex(pm.model, :pg) - p_expr = pm.model.ext[:p_expr] - p_ne_expr = pm.model.ext[:p_ne_expr] + pg = pm.var[:pg] + p = pm.var[:p] + p_ne = pm.var[:p_ne] + p_dc = pm.var[:p_dc] - c = @constraint(pm.model, sum(p_expr[a] for a in bus_arcs) + sum(p_ne_expr[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2) + c = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2) return Set([c]) end @@ -112,9 +118,9 @@ p[f_idx] == -b*(t[f_bus] - t[t_bus]) ``` """ function constraint_ohms_yt_from{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - p_fr = getindex(pm.model, :p)[f_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] + p_fr = pm.var[:p][f_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] c = @constraint(pm.model, p_fr == -b*(t_fr - t_to)) # omit reactive constraint @@ -125,10 +131,10 @@ end constraint_ohms_yt_to{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) = Set() function constraint_ohms_yt_from_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p_ne)[f_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + p_fr = pm.var[:p_ne][f_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] c1 = @constraint(pm.model, p_fr <= -b*(t_fr - t_to + t_max*(1-z)) ) c2 = @constraint(pm.model, p_fr >= -b*(t_fr - t_to + t_min*(1-z)) ) @@ -138,19 +144,10 @@ end "Do nothing, this model is symmetric" constraint_ohms_yt_to_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) = Set() -"`angmin <= t[f_bus] - t[t_bus] <= angmax`" -function constraint_phase_angle_difference{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - - c1 = @constraint(pm.model, t_fr - t_to <= angmax) - c2 = @constraint(pm.model, t_fr - t_to >= angmin) - return Set([c1, c2]) -end "`-rate_a <= p[f_idx] <= rate_a`" function constraint_thermal_limit_from{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, f_idx, rate_a) - p_fr = getindex(pm.model, :p)[f_idx] + p_fr = pm.var[:p][f_idx] if getlowerbound(p_fr) < -rate_a setlowerbound(p_fr, -rate_a) @@ -172,17 +169,6 @@ function add_bus_voltage_setpoint{T <: AbstractDCPForm}(sol, pm::GenericPowerMod add_setpoint(sol, pm, "bus", "bus_i", "va", :t) end -"" -function add_branch_flow_setpoint{T <: AbstractDCPForm}(sol, pm::GenericPowerModel{T}) - # check the line flows were requested - if haskey(pm.setting, "output") && haskey(pm.setting["output"], "line_flows") && pm.setting["output"]["line_flows"] == true - add_setpoint(sol, pm, "branch", "index", "pf", :p; extract_var = (var,idx,item) -> var[(idx, item["f_bus"], item["t_bus"])]) - add_setpoint(sol, pm, "branch", "index", "qf", :q; extract_var = (var,idx,item) -> var[(idx, item["f_bus"], item["t_bus"])]) - add_setpoint(sol, pm, "branch", "index", "pt", :p; scale = (x,item) -> -x, extract_var = (var,idx,item) -> var[(idx, item["f_bus"], item["t_bus"])]) - add_setpoint(sol, pm, "branch", "index", "qt", :q; scale = (x,item) -> -x, extract_var = (var,idx,item) -> var[(idx, item["f_bus"], item["t_bus"])]) - end -end - "" variable_voltage_on_off{T <: AbstractDCPForm}(pm::GenericPowerModel{T}; kwargs...) = variable_phase_angle(pm; kwargs...) @@ -191,10 +177,10 @@ constraint_voltage_on_off{T <: AbstractDCPForm}(pm::GenericPowerModel{T}) = noth "`-b*(t[f_bus] - t[t_bus] + t_min*(1-line_z[i])) <= p[f_idx] <= -b*(t[f_bus] - t[t_bus] + t_max*(1-line_z[i]))`" function constraint_ohms_yt_from_on_off{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p)[f_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + p_fr = pm.var[:p][f_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] c1 = @constraint(pm.model, p_fr <= -b*(t_fr - t_to + t_max*(1-z)) ) c2 = @constraint(pm.model, p_fr >= -b*(t_fr - t_to + t_min*(1-z)) ) @@ -212,8 +198,8 @@ Generic on/off thermal limit constraint ``` """ function constraint_thermal_limit_from_on_off{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_idx, rate_a) - p_fr = getindex(pm.model, :p)[f_idx] - z = getindex(pm.model, :line_z)[i] + p_fr = pm.var[:p][f_idx] + z = pm.var[:line_z][i] c1 = @constraint(pm.model, p_fr <= rate_a*z) c2 = @constraint(pm.model, p_fr >= -rate_a*z) @@ -228,8 +214,8 @@ Generic on/off thermal limit constraint ``` """ function constraint_thermal_limit_from_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_idx, rate_a) - p_fr = getindex(pm.model, :p_ne)[f_idx] - z = getindex(pm.model, :line_ne)[i] + p_fr = pm.var[:p_ne][f_idx] + z = pm.var[:line_ne][i] c1 = @constraint(pm.model, p_fr <= rate_a*z) c2 = @constraint(pm.model, p_fr >= -rate_a*z) @@ -244,9 +230,9 @@ constraint_thermal_limit_to_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i "`angmin*line_z[i] + t_min*(1-line_z[i]) <= t[f_bus] - t[t_bus] <= angmax*line_z[i] + t_max*(1-line_z[i])`" function constraint_phase_angle_difference_on_off{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, angmin, angmax, t_min, t_max) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] c1 = @constraint(pm.model, t_fr - t_to <= angmax*z + t_max*(1-z)) c2 = @constraint(pm.model, t_fr - t_to >= angmin*z + t_min*(1-z)) @@ -255,9 +241,9 @@ end "`angmin*line_ne[i] + t_min*(1-line_ne[i]) <= t[f_bus] - t[t_bus] <= angmax*line_ne[i] + t_max*(1-line_ne[i])`" function constraint_phase_angle_difference_ne{T <: AbstractDCPForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, angmin, angmax, t_min, t_max) - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] c1 = @constraint(pm.model, t_fr - t_to <= angmax*z + t_max*(1-z)) c2 = @constraint(pm.model, t_fr - t_to >= angmin*z + t_min*(1-z)) @@ -278,9 +264,9 @@ DCPLLPowerModel(data::Dict{String,Any}; kwargs...) = GenericPowerModel(data, Sta "`sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc)== sum(pg[g] for g in bus_gens) - pd - gs*1.0^2`" function constraint_kcl_shunt{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) - pg = getindex(pm.model, :pg) - p = getindex(pm.model, :p) - p_dc = getindex(pm.model, :p_dc) + pg = pm.var[:pg] + p = pm.var[:p] + p_dc = pm.var[:p_dc] c = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2) return Set([c]) @@ -288,10 +274,10 @@ end "`sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2`" function constraint_kcl_shunt_ne{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, pd, qd, gs, bs) - p = getindex(pm.model, :p) - p_ne = getindex(pm.model, :p_ne) - p_dc = getindex(pm.model, :p_dc) - pg = getindex(pm.model, :pg) + p = pm.var[:p] + p_ne = pm.var[:p_ne] + p_dc = pm.var[:p_dc] + pg = pm.var[:pg] c = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*1.0^2) return Set([c]) @@ -305,11 +291,11 @@ Creates Ohms constraints (yt post fix indicates that Y and T values are in recta ``` """ function constraint_ohms_yt_from_on_off{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p)[f_idx] - p_to = getindex(pm.model, :p)[t_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + p_fr = pm.var[:p][f_idx] + p_to = pm.var[:p][t_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] c1 = @constraint(pm.model, p_fr <= -b*(t_fr - t_to + t_max*(1-z)) ) c2 = @constraint(pm.model, p_fr >= -b*(t_fr - t_to + t_min*(1-z)) ) @@ -326,11 +312,11 @@ p[f_idx] + p[t_idx] >= r*( (-b*(t[f_bus] - t[t_bus]))^2 - (-b*(t_m))^2*(1-line_z where `r = g/(g^2 + b^2)` and `t_m = max(|t_min|, |t_max|)` """ function constraint_ohms_yt_to_on_off{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p)[f_idx] - p_to = getindex(pm.model, :p)[t_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_z)[i] + p_fr = pm.var[:p][f_idx] + p_to = pm.var[:p][t_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_z][i] r = g/(g^2 + b^2) t_m = max(abs(t_min),abs(t_max)) @@ -346,11 +332,11 @@ Creates Ohms constraints (yt post fix indicates that Y and T values are in recta ``` """ function constraint_ohms_yt_from_ne{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p_ne)[f_idx] - p_to = getindex(pm.model, :p_ne)[t_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + p_fr = pm.var[:p_ne][f_idx] + p_to = pm.var[:p_ne][t_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] c1 = @constraint(pm.model, p_fr <= -b*(t_fr - t_to + t_max*(1-z)) ) c2 = @constraint(pm.model, p_fr >= -b*(t_fr - t_to + t_min*(1-z)) ) @@ -367,11 +353,11 @@ p_ne[f_idx] + p_ne[t_idx] >= r*( (-b*(t[f_bus] - t[t_bus]))^2 - (-b*(t_m))^2*(1- where `r = g/(g^2 + b^2)` and `t_m = max(|t_min|, |t_max|)` """ function constraint_ohms_yt_to_ne{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p_ne)[f_idx] - p_to = getindex(pm.model, :p_ne)[t_idx] - t_fr = getindex(pm.model, :t)[f_bus] - t_to = getindex(pm.model, :t)[t_bus] - z = getindex(pm.model, :line_ne)[i] + p_fr = pm.var[:p_ne][f_idx] + p_to = pm.var[:p_ne][t_idx] + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + z = pm.var[:line_ne][i] r = g/(g^2 + b^2) t_m = max(abs(t_min),abs(t_max)) @@ -381,8 +367,8 @@ end "`-rate_a*line_z[i] <= p[t_idx] <= rate_a*line_z[i]`" function constraint_thermal_limit_to_on_off{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, t_idx, rate_a) - p_to = getindex(pm.model, :p)[t_idx] - z = getindex(pm.model, :line_z)[i] + p_to = pm.var[:p][t_idx] + z = pm.var[:line_z][i] c1 = @constraint(pm.model, p_to <= rate_a*z) c2 = @constraint(pm.model, p_to >= -rate_a*z) @@ -391,8 +377,8 @@ end "`-rate_a*line_ne[i] <= p_ne[t_idx] <= rate_a*line_ne[i]`" function constraint_thermal_limit_to_ne{T <: AbstractDCPLLForm}(pm::GenericPowerModel{T}, i, t_idx, rate_a) - p_to = getindex(pm.model, :p_ne)[t_idx] - z = getindex(pm.model, :line_ne)[i] + p_to = pm.var[:p_ne][t_idx] + z = pm.var[:line_ne][i] c1 = @constraint(pm.model, p_to <= rate_a*z) c2 = @constraint(pm.model, p_to >= -rate_a*z) diff --git a/src/form/shared.jl b/src/form/shared.jl index 00d989044..61ee67a51 100644 --- a/src/form/shared.jl +++ b/src/form/shared.jl @@ -19,10 +19,30 @@ # AbstractWRForms = Union{AbstractACTForm, AbstractWRForm} +AbstractPForms = Union{AbstractACPForm, AbstractACTForm, AbstractDCPForm} + +"`t[ref_bus] == 0`" +constraint_theta_ref{T <: AbstractPForms}(pm::GenericPowerModel{T}, ref_bus::Int) = + Set([@constraint(pm.model, pm.var[:t][ref_bus] == 0)]) + +""" +``` +t[f_bus] - t[t_bus] <= angmax +t[f_bus] - t[t_bus] >= angmin +``` +""" +function constraint_phase_angle_difference{T <: AbstractPForms}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) + t_fr = pm.var[:t][f_bus] + t_to = pm.var[:t][t_bus] + + c1 = @constraint(pm.model, t_fr - t_to <= angmax) + c2 = @constraint(pm.model, t_fr - t_to >= angmin) + return Set([c1, c2]) +end function constraint_voltage_magnitude_setpoint{T <: AbstractWRForms}(pm::GenericPowerModel{T}, i, vm, epsilon) - w = getindex(pm.model, :w)[i] + w = pm.var[:w][i] if epsilon == 0.0 c = @constraint(pm.model, w == vm^2) @@ -39,9 +59,9 @@ end """ enforces pv-like buses on both sides of a dcline """ -function constraint_dcline_voltage{T <: AbstractWRForms}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) - w_f = getindex(pm.model, :w)[f_bus] - w_t = getindex(pm.model, :w)[t_bus] +function constraint_voltage_dcline_setpoint{T <: AbstractWRForms}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) + w_f = pm.var[:w][f_bus] + w_t = pm.var[:w][t_bus] if epsilon == 0.0 c1 = @constraint(pm.model, w_f == vf^2) c2 = @constraint(pm.model, w_t == vt^2) @@ -63,13 +83,13 @@ sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[ ``` """ function constraint_kcl_shunt{T <: AbstractWRForms}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) - w = getindex(pm.model, :w)[i] - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) - p_dc = getindex(pm.model, :p_dc) - q_dc = getindex(pm.model, :q_dc) + w = pm.var[:w][i] + p = pm.var[:p] + q = pm.var[:q] + pg = pm.var[:pg] + qg = pm.var[:qg] + p_dc = pm.var[:p_dc] + q_dc = pm.var[:q_dc] c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w) c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w) @@ -81,11 +101,11 @@ end Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form) """ function constraint_ohms_yt_from{T <: AbstractWRForms}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - w_fr = getindex(pm.model, :w)[f_bus] - wr = getindex(pm.model, :wr)[(f_bus, t_bus)] - wi = getindex(pm.model, :wi)[(f_bus, t_bus)] + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + w_fr = pm.var[:w][f_bus] + wr = pm.var[:wr][(f_bus, t_bus)] + wi = pm.var[:wi][(f_bus, t_bus)] c1 = @constraint(pm.model, p_fr == g/tm*w_fr + (-g*tr+b*ti)/tm*(wr) + (-b*tr-g*ti)/tm*( wi) ) c2 = @constraint(pm.model, q_fr == -(b+c/2)/tm*w_fr - (-b*tr-g*ti)/tm*(wr) + (-g*tr+b*ti)/tm*( wi) ) @@ -97,11 +117,11 @@ end Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form) """ function constraint_ohms_yt_to{T <: AbstractWRForms}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - q_to = getindex(pm.model, :q)[t_idx] - p_to = getindex(pm.model, :p)[t_idx] - w_to = getindex(pm.model, :w)[t_bus] - wr = getindex(pm.model, :wr)[(f_bus, t_bus)] - wi = getindex(pm.model, :wi)[(f_bus, t_bus)] + q_to = pm.var[:q][t_idx] + p_to = pm.var[:p][t_idx] + w_to = pm.var[:w][t_bus] + wr = pm.var[:wr][(f_bus, t_bus)] + wi = pm.var[:wi][(f_bus, t_bus)] c1 = @constraint(pm.model, p_to == g*w_to + (-g*tr-b*ti)/tm*(wr) + (-b*tr+g*ti)/tm*(-wi) ) c2 = @constraint(pm.model, q_to == -(b+c/2)*w_to - (-b*tr+g*ti)/tm*(wr) + (-g*tr-b*ti)/tm*(-wi) ) diff --git a/src/form/wr.jl b/src/form/wr.jl index 2d0a76795..da6753bb5 100644 --- a/src/form/wr.jl +++ b/src/form/wr.jl @@ -22,9 +22,9 @@ end "" function constraint_voltage{T <: AbstractWRForm}(pm::GenericPowerModel{T}) - w = getindex(pm.model, :w) - wr = getindex(pm.model, :wr) - wi = getindex(pm.model, :wi) + w = pm.var[:w] + wr = pm.var[:wr] + wi = pm.var[:wi] for (i,j) in keys(pm.ref[:buspairs]) relaxation_complex_product(pm.model, w[i], w[j], wr[(i,j)], wi[(i,j)]) @@ -42,13 +42,13 @@ sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(q_dc[a_dc] ``` """ function constraint_kcl_shunt_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_arcs_ne, bus_gens, pd, qd, gs, bs) - w = getindex(pm.model, :w)[i] - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - p_ne = getindex(pm.model, :p_ne) - q_ne = getindex(pm.model, :q_ne) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) + w = pm.var[:w][i] + p = pm.var[:p] + q = pm.var[:q] + p_ne = pm.var[:p_ne] + q_ne = pm.var[:q_ne] + pg = pm.var[:pg] + qg = pm.var[:qg] c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w) c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w) @@ -66,11 +66,11 @@ q[f_idx] == -(b+c/2)/tm*w_from_ne[i] - (-b*tr-g*ti)/tm*(wr_ne[i]) + (-g*tr+b*ti) ``` """ function constraint_ohms_yt_from_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p_ne)[f_idx] - q_fr = getindex(pm.model, :q_ne)[f_idx] - w_fr = getindex(pm.model, :w_from_ne)[i] - wr = getindex(pm.model, :wr_ne)[i] - wi = getindex(pm.model, :wi_ne)[i] + p_fr = pm.var[:p_ne][f_idx] + q_fr = pm.var[:q_ne][f_idx] + w_fr = pm.var[:w_from_ne][i] + wr = pm.var[:wr_ne][i] + wi = pm.var[:wi_ne][i] c1 = @constraint(pm.model, p_fr == g/tm*w_fr + (-g*tr+b*ti)/tm*(wr) + (-b*tr-g*ti)/tm*( wi) ) c2 = @constraint(pm.model, q_fr == -(b+c/2)/tm*w_fr - (-b*tr-g*ti)/tm*(wr) + (-g*tr+b*ti)/tm*( wi) ) @@ -86,11 +86,11 @@ q[t_idx] == -(b+c/2)*w_to_ne[i] - (-b*tr+g*ti)/tm*(wr_ne[i]) + (-g*tr-b*ti)/tm*( ``` """ function constraint_ohms_yt_to_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_to = getindex(pm.model, :p_ne)[t_idx] - q_to = getindex(pm.model, :q_ne)[t_idx] - w_to = getindex(pm.model, :w_to_ne)[i] - wr = getindex(pm.model, :wr_ne)[i] - wi = getindex(pm.model, :wi_ne)[i] + p_to = pm.var[:p_ne][t_idx] + q_to = pm.var[:q_ne][t_idx] + w_to = pm.var[:w_to_ne][i] + wr = pm.var[:wr_ne][i] + wi = pm.var[:wi_ne][i] c1 = @constraint(pm.model, p_to == g*w_to + (-g*tr-b*ti)/tm*(wr) + (-b*tr+g*ti)/tm*(-wi) ) c2 = @constraint(pm.model, q_to == -(b+c/2)*w_to - (-b*tr+g*ti)/tm*(wr) + (-g*tr-b*ti)/tm*(-wi) ) @@ -99,10 +99,10 @@ end "" function constraint_phase_angle_difference{T <: AbstractWRForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) - w_fr = getindex(pm.model, :w)[f_bus] - w_to = getindex(pm.model, :w)[t_bus] - wr = getindex(pm.model, :wr)[(f_bus, t_bus)] - wi = getindex(pm.model, :wi)[(f_bus, t_bus)] + w_fr = pm.var[:w][f_bus] + w_to = pm.var[:w][t_bus] + wr = pm.var[:wr][(f_bus, t_bus)] + wi = pm.var[:wi][(f_bus, t_bus)] c1 = @constraint(pm.model, wi <= tan(angmax)*wr) c2 = @constraint(pm.model, wi >= tan(angmin)*wr) @@ -129,13 +129,13 @@ end "" function constraint_voltage_on_off{T <: AbstractWRForm}(pm::GenericPowerModel{T}) - w = getindex(pm.model, :w) - wr = getindex(pm.model, :wr) - wi = getindex(pm.model, :wi) - z = getindex(pm.model, :line_z) + w = pm.var[:w] + wr = pm.var[:wr] + wi = pm.var[:wi] + z = pm.var[:line_z] - w_from = getindex(pm.model, :w_from) - w_to = getindex(pm.model, :w_to) + w_from = pm.var[:w_from] + w_to = pm.var[:w_to] cs = Set() cs1 = constraint_voltage_magnitude_sqr_from_on_off(pm) @@ -161,13 +161,13 @@ function constraint_voltage_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}) wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:ne_buspairs]) bi_bp = Dict([(i, (b["f_bus"], b["t_bus"])) for (i,b) in branches]) - w = getindex(pm.model, :w) - wr = getindex(pm.model, :wr_ne) - wi = getindex(pm.model, :wi_ne) - z = getindex(pm.model, :line_ne) + w = pm.var[:w] + wr = pm.var[:wr_ne] + wi = pm.var[:wi_ne] + z = pm.var[:line_ne] - w_from = getindex(pm.model, :w_from_ne) - w_to = getindex(pm.model, :w_to_ne) + w_from = pm.var[:w_from_ne] + w_to = pm.var[:w_to_ne] cs = Set() for (l,i,j) in pm.ref[:ne_arcs_from] @@ -196,8 +196,8 @@ function constraint_voltage_magnitude_from_on_off{T <: AbstractWRForm}(pm::Gener buses = pm.ref[:bus] branches = pm.ref[:branch] - v_from = getindex(pm.model, :v_from) - z = getindex(pm.model, :line_z) + v_from = pm.var[:v_from] + z = pm.var[:line_z] cs = Set() for (i, branch) in pm.ref[:branch] @@ -214,8 +214,8 @@ function constraint_voltage_magnitude_to_on_off{T <: AbstractWRForm}(pm::Generic buses = pm.ref[:bus] branches = pm.ref[:branch] - v_to = getindex(pm.model, :v_to) - z = getindex(pm.model, :line_z) + v_to = pm.var[:v_to] + z = pm.var[:line_z] cs = Set() for (i, branch) in pm.ref[:branch] @@ -233,8 +233,8 @@ function constraint_voltage_magnitude_sqr_from_on_off{T <: AbstractWRForm}(pm::G buses = pm.ref[:bus] branches = pm.ref[:branch] - w_from = getindex(pm.model, :w_from) - z = getindex(pm.model, :line_z) + w_from = pm.var[:w_from] + z = pm.var[:line_z] cs = Set() for (i, branch) in pm.ref[:branch] @@ -251,8 +251,8 @@ function constraint_voltage_magnitude_sqr_to_on_off{T <: AbstractWRForm}(pm::Gen buses = pm.ref[:bus] branches = pm.ref[:branch] - w_to = getindex(pm.model, :w_to) - z = getindex(pm.model, :line_z) + w_to = pm.var[:w_to] + z = pm.var[:line_z] cs = Set() for (i, branch) in pm.ref[:branch] @@ -270,9 +270,9 @@ function constraint_voltage_product_on_off{T <: AbstractWRForm}(pm::GenericPower bi_bp = Dict([(i, (b["f_bus"], b["t_bus"])) for (i,b) in pm.ref[:branch]]) - wr = getindex(pm.model, :wr) - wi = getindex(pm.model, :wi) - z = getindex(pm.model, :line_z) + wr = pm.var[:wr] + wi = pm.var[:wi] + z = pm.var[:line_z] cs = Set() for b in keys(pm.ref[:branch]) @@ -297,11 +297,11 @@ q[f_idx] == -(b+c/2)/tm*w_from[i] - (-b*tr-g*ti)/tm*(wr[i]) + (-g*tr+b*ti)/tm*(w ``` """ function constraint_ohms_yt_from_on_off{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_fr = getindex(pm.model, :p)[f_idx] - q_fr = getindex(pm.model, :q)[f_idx] - w_fr = getindex(pm.model, :w_from)[i] - wr = getindex(pm.model, :wr)[i] - wi = getindex(pm.model, :wi)[i] + p_fr = pm.var[:p][f_idx] + q_fr = pm.var[:q][f_idx] + w_fr = pm.var[:w_from][i] + wr = pm.var[:wr][i] + wi = pm.var[:wi][i] c1 = @constraint(pm.model, p_fr == g/tm*w_fr + (-g*tr+b*ti)/tm*(wr) + (-b*tr-g*ti)/tm*( wi) ) c2 = @constraint(pm.model, q_fr == -(b+c/2)/tm*w_fr - (-b*tr-g*ti)/tm*(wr) + (-g*tr+b*ti)/tm*( wi) ) @@ -317,11 +317,11 @@ q[t_idx] == -(b+c/2)*w_to[i] - (-b*tr+g*ti)/tm*(wr[i]) + (-g*tr-b*ti)/tm*(-wi[i] ``` """ function constraint_ohms_yt_to_on_off{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm, t_min, t_max) - p_to = getindex(pm.model, :p)[t_idx] - q_to = getindex(pm.model, :q)[t_idx] - w_to = getindex(pm.model, :w_to)[i] - wr = getindex(pm.model, :wr)[i] - wi = getindex(pm.model, :wi)[i] + p_to = pm.var[:p][t_idx] + q_to = pm.var[:q][t_idx] + w_to = pm.var[:w_to][i] + wr = pm.var[:wr][i] + wi = pm.var[:wi][i] c1 = @constraint(pm.model, p_to == g*w_to + (-g*tr-b*ti)/tm*(wr) + (-b*tr+g*ti)/tm*(-wi) ) c2 = @constraint(pm.model, q_to == -(b+c/2)*w_to - (-b*tr+g*ti)/tm*(wr) + (-g*tr-b*ti)/tm*(-wi) ) @@ -330,8 +330,8 @@ end "`angmin*wr[i] <= wi[i] <= angmax*wr[i]`" function constraint_phase_angle_difference_on_off{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, angmin, angmax, t_min, t_max) - wr = getindex(pm.model, :wr)[i] - wi = getindex(pm.model, :wi)[i] + wr = pm.var[:wr][i] + wi = pm.var[:wi][i] c1 = @constraint(pm.model, wi <= tan(angmax)*wr) c2 = @constraint(pm.model, wi >= tan(angmin)*wr) @@ -340,8 +340,8 @@ end "`angmin*wr_ne[i] <= wi_ne[i] <= angmax*wr_ne[i]`" function constraint_phase_angle_difference_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}, i, f_bus, t_bus, angmin, angmax, t_min, t_max) - wr = getindex(pm.model, :wr_ne)[i] - wi = getindex(pm.model, :wi_ne)[i] + wr = pm.var[:wr_ne][i] + wi = pm.var[:wi_ne][i] c1 = @constraint(pm.model, wi <= tan(angmax)*wr) c2 = @constraint(pm.model, wi >= tan(angmin)*wr) @@ -359,25 +359,52 @@ end function variable_voltage_magnitude_sqr_from_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}) buses = pm.ref[:bus] branches = pm.ref[:ne_branch] - @variable(pm.model, 0 <= w_from_ne[i in keys(pm.ref[:ne_branch])] <= buses[branches[i]["f_bus"]]["vmax"]^2, start = getstart(pm.ref[:bus], i, "w_from_start", 1.001)) - return w_from_ne + + pm.var[:w_from_ne] = @variable(pm.model, + [i in keys(pm.ref[:ne_branch])], basename="w_from_ne", + lowerbound = 0, + upperbound = buses[branches[i]["f_bus"]]["vmax"]^2, + start = getstart(pm.ref[:bus], i, "w_from_start", 1.001) + ) + + return pm.var[:w_from_ne] end "" function variable_voltage_magnitude_sqr_to_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}) buses = pm.ref[:bus] branches = pm.ref[:ne_branch] - @variable(pm.model, 0 <= w_to_ne[i in keys(pm.ref[:ne_branch])] <= buses[branches[i]["t_bus"]]["vmax"]^2, start = getstart(pm.ref[:bus], i, "w_to", 1.001)) - return w_to_ne + + pm.var[:w_to_ne] = @variable(pm.model, + [i in keys(pm.ref[:ne_branch])], basename="w_to_ne", + lowerbound = 0, + upperbound = buses[branches[i]["t_bus"]]["vmax"]^2, + start = getstart(pm.ref[:bus], i, "w_to", 1.001) + ) + + return pm.var[:w_to_ne] end "" function variable_voltage_product_ne{T <: AbstractWRForm}(pm::GenericPowerModel{T}) wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:ne_buspairs]) bi_bp = Dict([(i, (b["f_bus"], b["t_bus"])) for (i,b) in pm.ref[:ne_branch]]) - @variable(pm.model, min(0, wr_min[bi_bp[b]]) <= wr_ne[b in keys(pm.ref[:ne_branch])] <= max(0, wr_max[bi_bp[b]]), start = getstart(pm.ref[:ne_buspairs], bi_bp[b], "wr_start", 1.0)) - @variable(pm.model, min(0, wi_min[bi_bp[b]]) <= wi_ne[b in keys(pm.ref[:ne_branch])] <= max(0, wi_max[bi_bp[b]]), start = getstart(pm.ref[:ne_buspairs], bi_bp[b], "wi_start")) - return wr_ne, wi_ne + + pm.var[:wr_ne] = @variable(pm.model, + [b in keys(pm.ref[:ne_branch])], basename="wr_ne", + lowerbound = min(0, wr_min[bi_bp[b]]), + upperbound = max(0, wr_max[bi_bp[b]]), + start = getstart(pm.ref[:ne_buspairs], bi_bp[b], "wr_start", 1.0) + ) + + pm.var[:wi_ne] = @variable(pm.model, + [b in keys(pm.ref[:ne_branch])], basename="wi_ne", + lowerbound = min(0, wi_min[bi_bp[b]]), + upperbound = max(0, wi_max[bi_bp[b]]), + start = getstart(pm.ref[:ne_buspairs], bi_bp[b], "wi_start") + ) + + return pm.var[:wr_ne], pm.var[:wi_ne] end "" @@ -391,19 +418,30 @@ function QCWRPowerModel(data::Dict{String,Any}; kwargs...) return GenericPowerModel(data, QCWRForm; kwargs...) end + + + "Creates variables associated with differences in phase angles" function variable_phase_angle_difference{T}(pm::GenericPowerModel{T}) - @variable(pm.model, pm.ref[:buspairs][bp]["angmin"] <= td[bp in keys(pm.ref[:buspairs])] <= pm.ref[:buspairs][bp]["angmax"], start = getstart(pm.ref[:buspairs], bp, "td_start")) - return td + pm.var[:td] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="td", + lowerbound = pm.ref[:buspairs][bp]["angmin"], + upperbound = pm.ref[:buspairs][bp]["angmax"], + start = getstart(pm.ref[:buspairs], bp, "td_start") + ) + return pm.var[:td] end "Creates the voltage magnitude product variables" function variable_voltage_magnitude_product{T}(pm::GenericPowerModel{T}) - vv_min = Dict([(bp, buspair["v_from_min"]*buspair["v_to_min"]) for (bp, buspair) in pm.ref[:buspairs]]) - vv_max = Dict([(bp, buspair["v_from_max"]*buspair["v_to_max"]) for (bp, buspair) in pm.ref[:buspairs]]) - - @variable(pm.model, vv_min[bp] <= vv[bp in keys(pm.ref[:buspairs])] <= vv_max[bp], start = getstart(pm.ref[:buspairs], bp, "vv_start", 1.0)) - return vv + buspairs = pm.ref[:buspairs] + pm.var[:vv] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="vv", + lowerbound = buspairs[bp]["v_from_min"]*buspairs[bp]["v_to_min"], + upperbound = buspairs[bp]["v_from_max"]*buspairs[bp]["v_to_max"], + start = getstart(pm.ref[:buspairs], bp, "vv_start", 1.0) + ) + return pm.var[:vv] end "" @@ -426,22 +464,36 @@ function variable_cosine{T}(pm::GenericPowerModel{T}) end end - @variable(pm.model, cos_min[bp] <= cs[bp in keys(pm.ref[:buspairs])] <= cos_max[bp], start = getstart(pm.ref[:buspairs], bp, "cs_start", 1.0)) - return cs + pm.var[:cs] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="cs", + lowerbound = cos_min[bp], + upperbound = cos_max[bp], + start = getstart(pm.ref[:buspairs], bp, "cs_start", 1.0) + ) + return pm.var[:cs] end "" -variable_sine(pm::GenericPowerModel) = - @variable(pm.model, sin(pm.ref[:buspairs][bp]["angmin"]) <= si[bp in keys(pm.ref[:buspairs])] <= sin(pm.ref[:buspairs][bp]["angmax"]), start = getstart(pm.ref[:buspairs], bp, "si_start")) +function variable_sine(pm::GenericPowerModel) + pm.var[:si] = @variable(pm.model, + [bp in keys(pm.ref[:buspairs])], basename="si", + lowerbound = sin(pm.ref[:buspairs][bp]["angmin"]), + upperbound = sin(pm.ref[:buspairs][bp]["angmax"]), + start = getstart(pm.ref[:buspairs], bp, "si_start") + ) + return pm.var[:si] +end "" function variable_current_magnitude_sqr{T}(pm::GenericPowerModel{T}) buspairs = pm.ref[:buspairs] - cm_min = Dict([(bp, 0) for bp in keys(pm.ref[:buspairs])]) - cm_max = Dict([(bp, (buspair["rate_a"]*buspair["tap"]/buspair["v_from_min"])^2) for (bp, buspair) in pm.ref[:buspairs]]) - - @variable(pm.model, cm_min[bp] <= cm[bp in keys(pm.ref[:buspairs])] <= cm_max[bp], start = getstart(pm.ref[:buspairs], bp, "cm_start")) - return cm + pm.var[:cm] = @variable(pm.model, + cm[bp in keys(pm.ref[:buspairs])], basename="cm", + lowerbound = 0, + upperbound = (buspairs[bp]["rate_a"]*buspairs[bp]["tap"]/buspairs[bp]["v_from_min"])^2, + start = getstart(pm.ref[:buspairs], bp, "cm_start") + ) + return pm.var[:cm] end "" @@ -461,17 +513,17 @@ end "" function constraint_voltage(pm::QCWRPowerModel) - v = getindex(pm.model, :v) - t = getindex(pm.model, :t) + v = pm.var[:v] + t = pm.var[:t] - td = getindex(pm.model, :td) - si = getindex(pm.model, :si) - cs = getindex(pm.model, :cs) - vv = getindex(pm.model, :vv) + td = pm.var[:td] + si = pm.var[:si] + cs = pm.var[:cs] + vv = pm.var[:vv] - w = getindex(pm.model, :w) - wr = getindex(pm.model, :wr) - wi = getindex(pm.model, :wi) + w = pm.var[:w] + wr = pm.var[:wr] + wi = pm.var[:wi] const_set = Set() for (i,b) in pm.ref[:bus] @@ -512,10 +564,10 @@ end "`p[f_idx]^2 + q[f_idx]^2 <= w[f_bus]/tm*cm[f_bus,t_bus]`" function constraint_power_magnitude_sqr(pm::QCWRPowerModel, f_bus, t_bus, arc_from, tm) - w_i = getindex(pm.model, :w)[f_bus] - p_fr = getindex(pm.model, :p)[arc_from] - q_fr = getindex(pm.model, :q)[arc_from] - cm = getindex(pm.model, :cm)[(f_bus, t_bus)] + w_i = pm.var[:w][f_bus] + p_fr = pm.var[:p][arc_from] + q_fr = pm.var[:q][arc_from] + cm = pm.var[:cm][(f_bus, t_bus)] c = @constraint(pm.model, p_fr^2 + q_fr^2 <= w_i/tm*cm) return Set([c]) @@ -523,12 +575,12 @@ end "`cm[f_bus,t_bus] == (g^2 + b^2)*(w[f_bus]/tm + w[t_bus] - 2*(tr*wr[f_bus,t_bus] + ti*wi[f_bus,t_bus])/tm) - c*q[f_idx] - ((c/2)/tm)^2*w[f_bus]`" function constraint_power_magnitude_link(pm::QCWRPowerModel, f_bus, t_bus, arc_from, g, b, c, tr, ti, tm) - w_fr = getindex(pm.model, :w)[f_bus] - w_to = getindex(pm.model, :w)[t_bus] - q_fr = getindex(pm.model, :q)[arc_from] - wr = getindex(pm.model, :wr)[(f_bus, t_bus)] - wi = getindex(pm.model, :wi)[(f_bus, t_bus)] - cm = getindex(pm.model, :cm)[(f_bus, t_bus)] + w_fr = pm.var[:w][f_bus] + w_to = pm.var[:w][t_bus] + q_fr = pm.var[:q][arc_from] + wr = pm.var[:wr][(f_bus, t_bus)] + wi = pm.var[:wi][(f_bus, t_bus)] + cm = pm.var[:cm][(f_bus, t_bus)] c = @constraint(pm.model, cm == (g^2 + b^2)*(w_fr/tm + w_to - 2*(tr*wr + ti*wi)/tm) - c*q_fr - ((c/2)/tm)^2*w_fr) return Set([c]) @@ -536,11 +588,11 @@ end "`t[ref_bus] == 0`" constraint_theta_ref(pm::QCWRPowerModel, ref_bus::Int) = - @constraint(pm.model, getindex(pm.model, :t)[ref_bus] == 0) + Set([@constraint(pm.model, pm.var[:t][ref_bus] == 0)]) "" function constraint_phase_angle_difference(pm::QCWRPowerModel, f_bus, t_bus, angmin, angmax) - td = getindex(pm.model, :td)[(f_bus, t_bus)] + td = pm.var[:td][(f_bus, t_bus)] if getlowerbound(td) < angmin setlowerbound(td, angmin) @@ -550,10 +602,10 @@ function constraint_phase_angle_difference(pm::QCWRPowerModel, f_bus, t_bus, ang setupperbound(td, angmax) end - w_fr = getindex(pm.model, :w)[f_bus] - w_to = getindex(pm.model, :w)[t_bus] - wr = getindex(pm.model, :wr)[(f_bus, t_bus)] - wi = getindex(pm.model, :wi)[(f_bus, t_bus)] + w_fr = pm.var[:w][f_bus] + w_to = pm.var[:w][t_bus] + wr = pm.var[:wr][(f_bus, t_bus)] + wi = pm.var[:wi][(f_bus, t_bus)] c1 = @constraint(pm.model, wi <= tan(angmax)*wr) c2 = @constraint(pm.model, wi >= tan(angmin)*wr) @@ -594,8 +646,13 @@ end "" function variable_phase_angle_difference_on_off{T}(pm::GenericPowerModel{T}) - @variable(pm.model, min(0, pm.ref[:branch][l]["angmin"]) <= td[l in keys(pm.ref[:branch])] <= max(0, pm.ref[:branch][l]["angmax"]), start = getstart(pm.ref[:branch], l, "td_start")) - return td + pm.var[:td] = @variable(pm.model, + td[l in keys(pm.ref[:branch])], basename="td", + lowerbound = min(0, pm.ref[:branch][l]["angmin"]), + upperbound = max(0, pm.ref[:branch][l]["angmax"]), + start = getstart(pm.ref[:branch], l, "td_start") + ) + return pm.var[:td] end "" @@ -603,8 +660,14 @@ function variable_voltage_magnitude_product_on_off{T}(pm::GenericPowerModel{T}) vv_min = Dict([(l, pm.ref[:bus][branch["f_bus"]]["vmin"]*pm.ref[:bus][branch["t_bus"]]["vmin"]) for (l, branch) in pm.ref[:branch]]) vv_max = Dict([(l, pm.ref[:bus][branch["f_bus"]]["vmax"]*pm.ref[:bus][branch["t_bus"]]["vmax"]) for (l, branch) in pm.ref[:branch]]) - @variable(pm.model, min(0, vv_min[l]) <= vv[l in keys(pm.ref[:branch])] <= max(0, vv_max[l]), start = getstart(pm.ref[:branch], l, "vv_start", 1.0)) - return vv + pm.var[:vv] = @variable(pm.model, + [l in keys(pm.ref[:branch])], basename="vv", + lowerbound = min(0, vv_min[l]), + upperbound = max(0, vv_max[l]), + start = getstart(pm.ref[:branch], l, "vv_start", 1.0) + ) + + return pm.var[:vv] end @@ -628,14 +691,25 @@ function variable_cosine_on_off{T}(pm::GenericPowerModel{T}) end end - @variable(pm.model, min(0, cos_min[l]) <= cs[l in keys(pm.ref[:branch])] <= max(0, cos_max[l]), start = getstart(pm.ref[:branch], l, "cs_start", 1.0)) - return cs + pm.var[:cs] = @variable(pm.model, + [l in keys(pm.ref[:branch])], basename="cs", + lowerbound = min(0, cos_min[l]), + upperbound = max(0, cos_max[l]), + start = getstart(pm.ref[:branch], l, "cs_start", 1.0) + ) + + return pm.var[:cs] end "" function variable_sine_on_off(pm::GenericPowerModel) - @variable(pm.model, min(0, sin(pm.ref[:branch][l]["angmin"])) <= si[l in keys(pm.ref[:branch])] <= max(0, sin(pm.ref[:branch][l]["angmax"])), start = getstart(pm.ref[:branch], l, "si_start")) - return si + pm.var[:si] = @variable(pm.model, + [l in keys(pm.ref[:branch])], basename="si", + lowerbound = min(0, sin(pm.ref[:branch][l]["angmin"])), + upperbound = max(0, sin(pm.ref[:branch][l]["angmax"])), + start = getstart(pm.ref[:branch], l, "si_start") + ) + return pm.var[:si] end @@ -644,31 +718,37 @@ function variable_current_magnitude_sqr_on_off{T}(pm::GenericPowerModel{T}) cm_min = Dict([(l, 0) for l in keys(pm.ref[:branch])]) cm_max = Dict([(l, (branch["rate_a"]*branch["tap"]/pm.ref[:bus][branch["f_bus"]]["vmin"])^2) for (l, branch) in pm.ref[:branch]]) - @variable(pm.model, cm_min[l] <= cm[l in keys(pm.ref[:branch])] <= cm_max[l], start = getstart(pm.ref[:branch], l, "cm_start")) - return cm + pm.var[:cm] = @variable(pm.model, + [l in keys(pm.ref[:branch])], basename="cm", + lowerbound = cm_min[l], + upperbound = cm_max[l], + start = getstart(pm.ref[:branch], l, "cm_start") + ) + + return pm.var[:cm] end "" function constraint_voltage_on_off(pm::QCWRPowerModel) - v = getindex(pm.model, :v) - t = getindex(pm.model, :t) - v_from = getindex(pm.model, :v_from) - v_to = getindex(pm.model, :v_to) + v = pm.var[:v] + t = pm.var[:t] + v_from = pm.var[:v_from] + v_to = pm.var[:v_to] - td = getindex(pm.model, :td) - si = getindex(pm.model, :si) - cs = getindex(pm.model, :cs) - vv = getindex(pm.model, :vv) + td = pm.var[:td] + si = pm.var[:si] + cs = pm.var[:cs] + vv = pm.var[:vv] - w = getindex(pm.model, :w) - w_from = getindex(pm.model, :w_from) - w_to = getindex(pm.model, :w_to) + w = pm.var[:w] + w_from = pm.var[:w_from] + w_to = pm.var[:w_to] - wr = getindex(pm.model, :wr) - wi = getindex(pm.model, :wi) + wr = pm.var[:wr] + wi = pm.var[:wi] - z = getindex(pm.model, :line_z) + z = pm.var[:line_z] td_lb = pm.ref[:off_angmin] td_ub = pm.ref[:off_angmax] @@ -725,11 +805,11 @@ end "`p[arc_from]^2 + q[arc_from]^2 <= w[f_bus]/tm*cm[i]`" function constraint_power_magnitude_sqr_on_off(pm::QCWRPowerModel, i, f_bus, arc_from, tm) - w = getindex(pm.model, :w)[f_bus] - p_fr = getindex(pm.model, :p)[arc_from] - q_fr = getindex(pm.model, :q)[arc_from] - cm = getindex(pm.model, :cm)[i] - z = getindex(pm.model, :line_z)[i] + w = pm.var[:w][f_bus] + p_fr = pm.var[:p][arc_from] + q_fr = pm.var[:q][arc_from] + cm = pm.var[:cm][i] + z = pm.var[:line_z][i] # TODO see if there is a way to leverage relaxation_complex_product_on_off here w_ub = getupperbound(w) @@ -745,12 +825,12 @@ end "`cm[f_bus,t_bus] == (g^2 + b^2)*(w[f_bus]/tm + w[t_bus] - 2*(tr*wr[f_bus,t_bus] + ti*wi[f_bus,t_bus])/tm) - c*q[f_idx] - ((c/2)/tm)^2*w[f_bus]`" function constraint_power_magnitude_link_on_off(pm::QCWRPowerModel, i, arc_from, g, b, c, tr, ti, tm) - w_fr = getindex(pm.model, :w_from)[i] - w_to = getindex(pm.model, :w_to)[i] - q_fr = getindex(pm.model, :q)[arc_from] - wr = getindex(pm.model, :wr)[i] - wi = getindex(pm.model, :wi)[i] - cm = getindex(pm.model, :cm)[i] + w_fr = pm.var[:w_from][i] + w_to = pm.var[:w_to][i] + q_fr = pm.var[:q][arc_from] + wr = pm.var[:wr][i] + wi = pm.var[:wi][i] + cm = pm.var[:cm][i] c = @constraint(pm.model, cm == (g^2 + b^2)*(w_fr/tm + w_to - 2*(tr*wr + ti*wi)/tm) - c*q_fr - ((c/2)/tm)^2*w_fr) return Set([c]) diff --git a/src/form/wrm.jl b/src/form/wrm.jl index a1680d7b0..5b21f3349 100644 --- a/src/form/wrm.jl +++ b/src/form/wrm.jl @@ -21,10 +21,14 @@ function variable_voltage_product_matrix{T <: AbstractWRMForm}(pm::GenericPowerM wr_min, wr_max, wi_min, wi_max = calc_voltage_product_bounds(pm.ref[:buspairs]) w_index = 1:length(keys(pm.ref[:bus])) - lookup_w_index = Dict([(bi, i) for (i,bi) in enumerate(keys(pm.ref[:bus]))]) + lookup_w_index = Dict([(bi,i) for (i,bi) in enumerate(keys(pm.ref[:bus]))]) - @variable(pm.model, WR[1:length(keys(pm.ref[:bus])), 1:length(keys(pm.ref[:bus]))], Symmetric) - @variable(pm.model, WI[1:length(keys(pm.ref[:bus])), 1:length(keys(pm.ref[:bus]))]) + WR = pm.var[:WR] = @variable(pm.model, + [1:length(keys(pm.ref[:bus])), 1:length(keys(pm.ref[:bus]))], Symmetric, basename="WR" + ) + WI = pm.var[:WI] = @variable(pm.model, + [1:length(keys(pm.ref[:bus])), 1:length(keys(pm.ref[:bus]))], basename="WI" + ) # bounds on diagonal for (i, bus) in pm.ref[:bus] @@ -52,14 +56,14 @@ function variable_voltage_product_matrix{T <: AbstractWRMForm}(pm::GenericPowerM setlowerbound(WI[wi_idx, wj_idx], wi_min[(i,j)]) end - pm.model.ext[:lookup_w_index] = lookup_w_index + pm.ext[:lookup_w_index] = lookup_w_index return WR, WI end "" function constraint_voltage{T <: AbstractWRMForm}(pm::GenericPowerModel{T}) - WR = getindex(pm.model, :WR) - WI = getindex(pm.model, :WI) + WR = pm.var[:WR] + WI = pm.var[:WI] c = @SDconstraint(pm.model, [WR WI; -WI WR] >= 0) @@ -75,16 +79,15 @@ constraint_theta_ref{T <: AbstractWRMForm}(pm::GenericPowerModel{T}, ref_bus::In "" function constraint_kcl_shunt{T <: AbstractWRMForm}(pm::GenericPowerModel{T}, i, bus_arcs, bus_arcs_dc, bus_gens, pd, qd, gs, bs) - WR = getindex(pm.model, :WR) - w_index = pm.model.ext[:lookup_w_index][i] - w = WR[w_index, w_index] + w_index = pm.ext[:lookup_w_index][i] + w = pm.var[:WR][w_index, w_index] - p = getindex(pm.model, :p) - q = getindex(pm.model, :q) - pg = getindex(pm.model, :pg) - qg = getindex(pm.model, :qg) - p_dc = getindex(pm.model, :p_dc) - q_dc = getindex(pm.model, :q_dc) + p = pm.var[:p] + q = pm.var[:q] + pg = pm.var[:pg] + qg = pm.var[:qg] + p_dc = pm.var[:p_dc] + q_dc = pm.var[:q_dc] c1 = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - pd - gs*w) c2 = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - qd + bs*w) @@ -93,18 +96,16 @@ end "Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)" function constraint_ohms_yt_from{T <: AbstractWRMForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - 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] - WR = getindex(pm.model, :WR) - WI = getindex(pm.model, :WI) - w_fr_index = pm.model.ext[:lookup_w_index][f_bus] - w_to_index = pm.model.ext[:lookup_w_index][t_bus] + w_fr_index = pm.ext[:lookup_w_index][f_bus] + w_to_index = pm.ext[:lookup_w_index][t_bus] - w_fr = WR[w_fr_index, w_fr_index] - w_to = WR[w_to_index, w_to_index] - wr = WR[w_fr_index, w_to_index] - wi = WI[w_fr_index, w_to_index] + w_fr = pm.var[:WR][w_fr_index, w_fr_index] + w_to = pm.var[:WR][w_to_index, w_to_index] + wr = pm.var[:WR][w_fr_index, w_to_index] + wi = pm.var[:WI][w_fr_index, w_to_index] c1 = @constraint(pm.model, p_fr == g/tm*w_fr + (-g*tr+b*ti)/tm*(wr) + (-b*tr-g*ti)/tm*( wi) ) c2 = @constraint(pm.model, q_fr == -(b+c/2)/tm*w_fr - (-b*tr-g*ti)/tm*(wr) + (-g*tr+b*ti)/tm*( wi) ) @@ -113,18 +114,16 @@ end "" function constraint_ohms_yt_to{T <: AbstractWRMForm}(pm::GenericPowerModel{T}, f_bus, t_bus, f_idx, t_idx, g, b, c, tr, ti, tm) - q_to = getindex(pm.model, :q)[t_idx] - p_to = getindex(pm.model, :p)[t_idx] + q_to = pm.var[:q][t_idx] + p_to = pm.var[:p][t_idx] - WR = getindex(pm.model, :WR) - WI = getindex(pm.model, :WI) - w_fr_index = pm.model.ext[:lookup_w_index][f_bus] - w_to_index = pm.model.ext[:lookup_w_index][t_bus] + w_fr_index = pm.ext[:lookup_w_index][f_bus] + w_to_index = pm.ext[:lookup_w_index][t_bus] - w_fr = WR[w_fr_index, w_fr_index] - w_to = WR[w_to_index, w_to_index] - wr = WR[w_fr_index, w_to_index] - wi = WI[w_fr_index, w_to_index] + w_fr = pm.var[:WR][w_fr_index, w_fr_index] + w_to = pm.var[:WR][w_to_index, w_to_index] + wr = pm.var[:WR][w_fr_index, w_to_index] + wi = pm.var[:WI][w_fr_index, w_to_index] c1 = @constraint(pm.model, p_to == g*w_to + (-g*tr-b*ti)/tm*(wr) + (-b*tr+g*ti)/tm*(-wi) ) c2 = @constraint(pm.model, q_to == -(b+c/2)*w_to - (-b*tr+g*ti)/tm*(wr) + (-g*tr-b*ti)/tm*(-wi) ) @@ -133,15 +132,13 @@ end "" function constraint_phase_angle_difference{T <: AbstractWRMForm}(pm::GenericPowerModel{T}, f_bus, t_bus, angmin, angmax) - WR = getindex(pm.model, :WR) - WI = getindex(pm.model, :WI) - w_fr_index = pm.model.ext[:lookup_w_index][f_bus] - w_to_index = pm.model.ext[:lookup_w_index][t_bus] + w_fr_index = pm.ext[:lookup_w_index][f_bus] + w_to_index = pm.ext[:lookup_w_index][t_bus] - w_fr = WR[w_fr_index, w_fr_index] - w_to = WR[w_to_index, w_to_index] - wr = WR[w_fr_index, w_to_index] - wi = WI[w_fr_index, w_to_index] + w_fr = pm.var[:WR][w_fr_index, w_fr_index] + w_to = pm.var[:WR][w_to_index, w_to_index] + wr = pm.var[:WR][w_fr_index, w_to_index] + wi = pm.var[:WI][w_fr_index, w_to_index] c1 = @constraint(pm.model, wi <= tan(angmax)*wr) c2 = @constraint(pm.model, wi >= tan(angmin)*wr) @@ -153,11 +150,9 @@ end "" function add_bus_voltage_setpoint{T <: AbstractWRMForm}(sol, pm::GenericPowerModel{T}) - add_setpoint(sol, pm, "bus", "bus_i", "vm", :WR; scale = (x,item) -> sqrt(x), extract_var = (var,idx,item) -> var[pm.model.ext[:lookup_w_index][idx], pm.model.ext[:lookup_w_index][idx]]) + add_setpoint(sol, pm, "bus", "bus_i", "vm", :WR; scale = (x,item) -> sqrt(x), extract_var = (var,idx,item) -> var[pm.ext[:lookup_w_index][idx], pm.ext[:lookup_w_index][idx]]) # What should the default value be? #add_setpoint(sol, pm, "bus", "bus_i", "va", :t; default_value = 0) end -"DC Line voltage constraint not supported" -constraint_dcline_voltage{T <: AbstractWRMForm}(pm::GenericPowerModel{T}, f_bus, t_bus, vf, vt, epsilon) = Set() diff --git a/src/prob/misc.jl b/src/prob/misc.jl index b304d15eb..674273cf1 100644 --- a/src/prob/misc.jl +++ b/src/prob/misc.jl @@ -32,7 +32,7 @@ function post_api_opf(pm::GenericPowerModel) end for (i,gen) in pm.ref[:gen] - pg = getindex(pm.model, :pg)[i] + pg = pm.var[:pg][i] @constraint(pm.model, pg >= gen["pmin"]) end @@ -87,8 +87,8 @@ function post_sad_opf{T <: Union{AbstractACPForm, AbstractDCPForm}}(pm::GenericP constraint_ohms_yt_to(pm, branch) constraint_phase_angle_difference(pm, branch) - theta_fr = getindex(pm.model, :t)[branch["f_bus"]] - theta_to = getindex(pm.model, :t)[branch["t_bus"]] + theta_fr = pm.var[:t][branch["f_bus"]] + theta_to = pm.var[:t][branch["t_bus"]] @constraint(pm.model, theta_fr - theta_to <= theta_delta_bound) @constraint(pm.model, theta_fr - theta_to >= -theta_delta_bound) diff --git a/src/prob/pf.jl b/src/prob/pf.jl index d69f57fcc..de575659f 100644 --- a/src/prob/pf.jl +++ b/src/prob/pf.jl @@ -53,6 +53,6 @@ function post_pf(pm::GenericPowerModel) for (i,dcline) in pm.ref[:dcline] #constraint_dcline(pm, dcline) not needed, active power flow fully defined by dc line setpoints constraint_active_dcline_setpoint(pm, dcline) - constraint_dcline_voltage(pm, dcline; epsilon = 0.00001) + constraint_voltage_dcline_setpoint(pm, dcline; epsilon = 0.00001) end end diff --git a/test/ots.jl b/test/ots.jl index cf8d83120..143e44ea7 100644 --- a/test/ots.jl +++ b/test/ots.jl @@ -141,6 +141,14 @@ end @test result["status"] == :Optimal @test isapprox(result["objective"], 15001.6; atol = 1e0) end + @testset "5-bus asymmetric case" begin + result = run_ots("../test/data/case5_asym.m", QCWRPowerModel, pajarito_solver) + + check_br_status(result["solution"]) + + @test result["status"] == :Optimal + @test isapprox(result["objective"], 15009.4; atol = 1e0) + end @testset "6-bus case" begin result = run_ots("../test/data/case6.m", QCWRPowerModel, pajarito_solver)