diff --git a/CHANGELOG.md b/CHANGELOG.md index 2b1b84756..e94ffdaa0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## staged +- Add wye-connected CapControl for IVR and FOT (polar) formulations - Fixed indexing issue for single-phase delta load models in linear formulations (LinDist3Flow, FOTP, FOTR, FBS) - Added ZIP load model - Updated documentation in `make_multiconductor!` to better indicate its unsupported nature diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index 37342e800..583f23364 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -395,6 +395,26 @@ function constraint_mc_current_balance(pm::AbstractUnbalancedPowerModel, i::Int; end +""" + constraint_mc_current_balance_capc(pm::AbstractUnbalancedPowerModel, i::Int; nw::Int=nw_id_default)::Nothing + +Template function for KCL constraints in current-voltage variable space with capacitor control variables. +""" +function constraint_mc_current_balance_capc(pm::AbstractUnbalancedPowerModel, i::Int; nw::Int=nw_id_default)::Nothing + bus = ref(pm, nw, :bus, i) + bus_arcs = ref(pm, nw, :bus_arcs_conns_branch, i) + bus_arcs_sw = ref(pm, nw, :bus_arcs_conns_switch, i) + bus_arcs_trans = ref(pm, nw, :bus_arcs_conns_transformer, i) + bus_gens = ref(pm, nw, :bus_conns_gen, i) + bus_storage = ref(pm, nw, :bus_conns_storage, i) + bus_loads = ref(pm, nw, :bus_conns_load, i) + bus_shunts = ref(pm, nw, :bus_conns_shunt, i) + + constraint_mc_current_balance_capc(pm, nw, i, bus["terminals"], bus["grounded"], bus_arcs, bus_arcs_sw, bus_arcs_trans, bus_gens, bus_storage, bus_loads, bus_shunts) + nothing +end + + """ constraint_mc_network_power_balance(pm::AbstractUnbalancedPowerModel, i::Int; nw::Int=nw_id_default)::Nothing diff --git a/src/core/variable.jl b/src/core/variable.jl index 0876fc5b5..a5350ae89 100644 --- a/src/core/variable.jl +++ b/src/core/variable.jl @@ -729,6 +729,16 @@ function variable_mc_switch_state(pm::AbstractUnbalancedPowerModel; nw::Int=nw_i end +""" + variable_mc_capcontrol(pm::AbstractUnbalancedPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) + +Capacitor switching variables. +""" +function variable_mc_capcontrol(pm::AbstractUnbalancedPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) + variable_mc_capacitor_switch_state(pm; nw=nw, relax=relax, report=report) +end + + """ variable_mc_capacitor_switch_state(pm::AbstractUnbalancedPowerModel, relax::Bool; nw::Int=nw_id_default, report::Bool=true) diff --git a/src/form/acp.jl b/src/form/acp.jl index f6d5a6627..cfbde1b5a 100644 --- a/src/form/acp.jl +++ b/src/form/acp.jl @@ -68,16 +68,6 @@ function variable_mc_bus_voltage_on_off(pm::AbstractUnbalancedACPModel; nw::Int= end -""" - variable_mc_capcontrol(pm::AbstractUnbalancedACPModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - -Capacitor switching variables. -""" -function variable_mc_capcontrol(pm::AbstractUnbalancedACPModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - variable_mc_capacitor_switch_state(pm; nw=nw, relax=relax, report=report) -end - - "" function constraint_mc_switch_state_closed(pm::AbstractUnbalancedACPModel, nw::Int, f_bus::Int, t_bus::Int, f_connections::Vector{Int}, t_connections::Vector{Int}) vm_fr = var(pm, nw, :vm, f_bus) diff --git a/src/form/acr.jl b/src/form/acr.jl index 3fd1a20b8..3981a78a5 100644 --- a/src/form/acr.jl +++ b/src/form/acr.jl @@ -107,16 +107,6 @@ function variable_mc_bus_voltage_on_off(pm::AbstractUnbalancedACRModel; nw::Int= end -""" - variable_mc_capcontrol(pm::AbstractUnbalancedACRModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - -Capacitor switching variables. -""" -function variable_mc_capcontrol(pm::AbstractUnbalancedACRModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - variable_mc_capacitor_switch_state(pm; nw=nw, relax=relax, report=report) -end - - "" function constraint_mc_switch_state_closed(pm::AbstractUnbalancedACRModel, nw::Int, f_bus::Int, t_bus::Int, f_connections::Vector{Int}, t_connections::Vector{Int}) vr_fr = var(pm, nw, :vr, f_bus) diff --git a/src/form/bf_fbs.jl b/src/form/bf_fbs.jl index 86c332037..7d8b7ab9a 100644 --- a/src/form/bf_fbs.jl +++ b/src/form/bf_fbs.jl @@ -89,16 +89,6 @@ function variable_mc_bus_voltage(pm::FBSUBFPowerModel; nw::Int=nw_id_default, bo end -""" - variable_mc_capcontrol(pm::FBSUBFPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - -Capacitor switching variables. -""" -function variable_mc_capcontrol(pm::FBSUBFPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - variable_mc_capacitor_switch_state(pm; nw=nw, relax=relax, report=report) -end - - "Creates variables for both `active` and `reactive` power flow at each transformer." function variable_mc_transformer_power(pm::FBSUBFPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true) variable_mc_transformer_power_real(pm; nw=nw, bounded=bounded, report=report) diff --git a/src/form/fotp.jl b/src/form/fotp.jl index 29f0facdf..409bb2550 100644 --- a/src/form/fotp.jl +++ b/src/form/fotp.jl @@ -200,6 +200,238 @@ function constraint_mc_power_balance(pm::FOTPUPowerModel, nw::Int, i::Int, termi end +@doc raw""" + constraint_mc_power_balance_capc(pm::FOTPUPowerModel, nw::Int, i::Int, terminals::Vector{Int}, grounded::Vector{Bool}, bus_arcs::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_sw::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_trans::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_gens::Vector{Tuple{Int,Vector{Int}}}, bus_storage::Vector{Tuple{Int,Vector{Int}}}, bus_loads::Vector{Tuple{Int,Vector{Int}}}, bus_shunts::Vector{Tuple{Int,Vector{Int}}}) + +Power balance constraints with capacitor control with shunt current calculated using initial operating point. + +```math +\begin{align} + & B_s = b_s ⋅ z,~~ cq_{sh} = B_s ⋅ v, \\ + & B_s \cdot v_m^t \cdot v_m^u \cdot \cos(v_a^t-v_a^u) \Rightarrow B_{s0} \cdot v_{m0}^t \cdot v_{m0}^u \cdot \cos(v_{a0}^t-v_{a0}^u) + +\begin{bmatrix} +B_{s0} \cdot v_{m0}^u \cdot \cos(v_{a0}^t-v_{a0}^u) \\ +B_{s0} \cdot v_{m0}^t \cdot \cos(v_{a0}^t-v_{a0}^u) \\ +-B_{s0} \cdot v_{m0}^t \cdot v_{m0}^u \cdot \sin(v_{a0}^t-v_{a0}^u) \\ +B_{s0} \cdot v_{m0}^t \cdot v_{m0}^u \cdot \sin(v_{a0}^t-v_{a0}^u) \\ +v_{m0}^t \cdot v_{m0}^u \cdot \cos(v_{a0}^t-v_{a0}^u) +\end{bmatrix}^\top +\begin{bmatrix} +v_m^t-v_{m0}^t \\ +v_m^u-v_{m0}^u \\ +v_a^t-v_{a0}^t \\ +v_a^u-v_{a0}^u \\ +B_{s} -B_{s0} +\end{bmatrix} \\ +& B_s \cdot v_m^t \cdot v_m^u \cdot \sin(v_a^t-v_a^u) \Rightarrow B_{s0} \cdot v_{m0}^t \cdot v_{m0}^u \cdot \sin(v_{a0}^t-v_{a0}^u) + +\begin{bmatrix} + B_{s0} \cdot v_{m0}^u \cdot \sin(v_{a0}^t-v_{a0}^u) \\ + B_{s0} \cdot v_{m0}^t \cdot \sin(v_{a0}^t-v_{a0}^u) \\ + B_{s0} \cdot v_{m0}^t \cdot v_{m0}^u \cdot \cos(v_{a0}^t-v_{a0}^u) \\ + -B_{s0} \cdot v_{m0}^t \cdot v_{m0}^u \cdot \cos(v_{a0}^t-v_{a0}^u) \\ + v_{m0}^t \cdot v_{m0}^u \cdot \sin(v_{a0}^t-v_{a0}^u) +\end{bmatrix}^\top +\begin{bmatrix} +v_m^t-v_{m0}^t \\ +v_m^u-v_{m0}^u \\ +v_a^t-v_{a0}^t \\ +v_a^u-v_{a0}^u \\ +B_{s} -B_{s0} +\end{bmatrix} + +\end{align} +``` +""" +function constraint_mc_power_balance_capc(pm::FOTPUPowerModel, nw::Int, i::Int, terminals::Vector{Int}, grounded::Vector{Bool}, bus_arcs::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_sw::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_trans::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_gens::Vector{Tuple{Int,Vector{Int}}}, bus_storage::Vector{Tuple{Int,Vector{Int}}}, bus_loads::Vector{Tuple{Int,Vector{Int}}}, bus_shunts::Vector{Tuple{Int,Vector{Int}}}) + vm = var(pm, nw, :vm, i) + va = var(pm, nw, :va, i) + vm0 = var(pm, nw, :vm0, i) + va0 = var(pm, nw, :va0, i) + p = get(var(pm, nw), :p, Dict()); _check_var_keys( p, bus_arcs, "active power", "branch") + q = get(var(pm, nw), :q, Dict()); _check_var_keys( q, bus_arcs, "reactive power", "branch") + pg = get(var(pm, nw), :pg_bus, Dict()); _check_var_keys( pg, bus_gens, "active power", "generator") + qg = get(var(pm, nw), :qg_bus, Dict()); _check_var_keys( qg, bus_gens, "reactive power", "generator") + ps = get(var(pm, nw), :ps, Dict()); _check_var_keys( ps, bus_storage, "active power", "storage") + qs = get(var(pm, nw), :qs, Dict()); _check_var_keys( qs, bus_storage, "reactive power", "storage") + psw = get(var(pm, nw), :psw, Dict()); _check_var_keys(psw, bus_arcs_sw, "active power", "switch") + qsw = get(var(pm, nw), :qsw, Dict()); _check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch") + pt = get(var(pm, nw), :pt, Dict()); _check_var_keys( pt, bus_arcs_trans, "active power", "transformer") + qt = get(var(pm, nw), :qt, Dict()); _check_var_keys( qt, bus_arcs_trans, "reactive power", "transformer") + pd = get(var(pm, nw), :pd_bus, Dict()); _check_var_keys( pd, bus_loads, "active power", "load") + qd = get(var(pm, nw), :qd_bus, Dict()); _check_var_keys( pd, bus_loads, "reactive power", "load") + + # calculate Gs, Bs + ncnds = length(terminals) + Gs = fill(0.0, ncnds, ncnds) + Bs0 = fill(0.0, ncnds, ncnds) + Bs = convert(Matrix{JuMP.AffExpr}, JuMP.@expression(pm.model, [idx=1:ncnds, jdx=1:ncnds], 0.0)) + for (val, connections) in bus_shunts + shunt = ref(pm,nw,:shunt,val) + for (idx,c) in enumerate(connections) + cap_state = haskey(shunt,"controls") ? var(pm, nw, :capacitor_state, val)[c] : 1.0 + for (jdx,d) in enumerate(connections) + Gs[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] += shunt["gs"][idx,jdx] + Bs0[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] += shunt["bs"][idx,jdx] + Bs[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] = JuMP.@expression(pm.model, Bs[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] + shunt["bs"][idx,jdx]*cap_state) + end + end + end + + cstr_p = [] + cstr_q = [] + ungrounded_terminals = [(idx,t) for (idx,t) in enumerate(terminals) if !grounded[idx]] + + # add constraints to model capacitor switching + if !isempty(bus_shunts) && haskey(ref(pm, nw, :shunt, bus_shunts[1][1]), "controls") + constraint_capacitor_on_off(pm, nw, i, bus_shunts) + + for (idx,t) in ungrounded_terminals + if any(Bs[idx,jdx] != 0 for (jdx, u) in ungrounded_terminals if idx != jdx) || any(Gs[idx,jdx] != 0 for (jdx, u) in ungrounded_terminals if idx != jdx) + cp = JuMP.@constraint(pm.model, + sum( p[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(psw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( pt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( pg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( ps[s][t] for (s, conns) in bus_storage if t in conns) + + sum( pd[l][t] for (l, conns) in bus_loads if t in conns) + + ( Gs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + +sum( Gs[idx,jdx] * vm0[idx]*vm0[jdx] * cos(va0[idx]-va0[jdx]) + +Bs0[idx,jdx] * vm0[idx]*vm0[jdx] * sin(va0[idx]-va0[jdx]) + +[Gs[idx,jdx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) Gs[idx,jdx]*vm0[idx]*cos(va0[idx]-va0[jdx]) -Gs[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) Gs[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx]] + +[Bs0[idx,jdx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) Bs0[idx,jdx]*vm0[idx]*sin(va0[idx]-va0[jdx]) Bs0[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) -Bs0[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) vm0[t]*vm0[u]*sin(va0[t]-va0[u])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx];Bs[idx,jdx]-Bs0[idx,jdx]] + for (jdx,u) in ungrounded_terminals if idx != jdx) + ) + == + 0.0 + ) + push!(cstr_p, cp) + + cq = JuMP.@constraint(pm.model, + sum( q[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(qsw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( qt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( qg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( qs[s][t] for (s, conns) in bus_storage if t in conns) + + sum( qd[l][t] for (l, conns) in bus_loads if t in conns) + + ( -Bs0[idx,idx]*vm0[idx]^2 - Bs0[idx,idx]*2*vm0[idx]*(vm[t]-vm0[idx]) - vm0[idx]^2*(Bs[idx,idx]-Bs0[idx,idx]) + -sum( Bs0[idx,jdx] * vm0[idx]*vm0[jdx] * cos(va0[idx]-va0[jdx]) + -Gs[idx,jdx] * vm0[idx]*vm0[jdx] * sin(va0[idx]-va0[jdx]) + +[Bs0[idx,jdx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) Bs0[idx,jdx]*vm0[idx]*cos(va0[idx]-va0[jdx]) -Bs0[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) Bs0[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) vm0[t]*vm0[u]*cos(va0[t]-va0[u])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx];Bs[idx,jdx]-Bs0[idx,jdx]] + +[-Gs[idx,jdx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) -Gs[idx,jdx]*vm0[idx]*sin(va0[idx]-va0[jdx]) -Gs[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) Gs[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx]] + for (jdx,u) in ungrounded_terminals if idx != jdx) + ) + == + 0.0 + ) + push!(cstr_q, cq) + else + cp = JuMP.@constraint(pm.model, + sum( p[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(psw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( pt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( pg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( ps[s][t] for (s, conns) in bus_storage if t in conns) + + sum( pd[l][t] for (l, conns) in bus_loads if t in conns) + + Gs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + == + 0.0 + ) + push!(cstr_p, cp) + + cq = JuMP.@constraint(pm.model, + sum( q[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(qsw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( qt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( qg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( qs[s][t] for (s, conns) in bus_storage if t in conns) + + sum( qd[l][t] for (l, conns) in bus_loads if t in conns) + - Bs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + == + 0.0 + ) + push!(cstr_q, cq) + end + end + else + for (idx,t) in ungrounded_terminals + if any(Bs[idx,jdx] != 0 for (jdx, u) in ungrounded_terminals if idx != jdx) || any(Gs[idx,jdx] != 0 for (jdx, u) in ungrounded_terminals if idx != jdx) + cp = JuMP.@constraint(pm.model, + sum( p[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(psw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( pt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( pg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( ps[s][t] for (s, conns) in bus_storage if t in conns) + + sum( pd[l][t] for (l, conns) in bus_loads if t in conns) + + ( Gs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + +sum( Gs[idx,jdx] * vm0[idx]*vm0[jdx] * cos(va0[idx]-va0[jdx]) + +Bs[idx,jdx] * vm0[idx]*vm0[jdx] * sin(va0[idx]-va0[jdx]) + +[Gs[idx,jdx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) Gs[idx,jdx]*vm0[idx]*cos(va0[idx]-va0[jdx]) -Gs[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) Gs[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx]] + +[Bs[idx,jdx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) Bs[idx,jdx]*vm0[idx]*sin(va0[idx]-va0[jdx]) Bs[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) -Bs[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx]] + for (jdx,u) in ungrounded_terminals if idx != jdx) + ) + == + 0.0 + ) + push!(cstr_p, cp) + + cq = JuMP.@constraint(pm.model, + sum( q[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(qsw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( qt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( qg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( qs[s][t] for (s, conns) in bus_storage if t in conns) + + sum( qd[l][t] for (l, conns) in bus_loads if t in conns) + + ( -Bs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + -sum( Bs[idx,jdx] * vm0[idx]*vm0[jdx] * cos(va0[idx]-va0[jdx]) + -Gs[idx,jdx] * vm0[idx]*vm0[jdx] * sin(va0[idx]-va0[jdx]) + +[Bs[idx,jdx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) Bs[idx,jdx]*vm0[idx]*cos(va0[idx]-va0[jdx]) -Bs[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) Bs[idx,jdx]*vm0[idx]*vm0[jdx]*sin(va0[idx]-va0[jdx])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx]] + +[-Gs[idx,jdx]*vm0[jdx]*sin(va0[idx]-va0[jdx]) -Gs[idx,jdx]*vm0[idx]*sin(va0[idx]-va0[jdx]) -Gs[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx]) Gs[idx,jdx]*vm0[idx]*vm0[jdx]*cos(va0[idx]-va0[jdx])]*[vm[t]-vm0[idx];vm[u]-vm0[jdx];va[t]-va0[idx];va[u]-va0[jdx]] + for (jdx,u) in ungrounded_terminals if idx != jdx) + ) + == + 0.0 + ) + push!(cstr_q, cq) + else + cp = JuMP.@constraint(pm.model, + sum( p[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(psw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( pt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( pg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( ps[s][t] for (s, conns) in bus_storage if t in conns) + + sum( pd[l][t] for (l, conns) in bus_loads if t in conns) + + Gs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + == + 0.0 + ) + push!(cstr_p, cp) + + cq = JuMP.@constraint(pm.model, + sum( q[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(qsw[a][t] for (a, conns) in bus_arcs_sw if t in conns) + + sum( qt[a][t] for (a, conns) in bus_arcs_trans if t in conns) + - sum( qg[g][t] for (g, conns) in bus_gens if t in conns) + + sum( qs[s][t] for (s, conns) in bus_storage if t in conns) + + sum( qd[l][t] for (l, conns) in bus_loads if t in conns) + - Bs[idx,idx]*(vm0[idx]^2+2*vm0[idx]*(vm[t]-vm0[idx])) + == + 0.0 + ) + push!(cstr_q, cq) + end + end + end + + con(pm, nw, :lam_kcl_r)[i] = cstr_p + con(pm, nw, :lam_kcl_i)[i] = cstr_q + + if _IM.report_duals(pm) + sol(pm, nw, :bus, i)[:lam_kcl_r] = cstr_p + sol(pm, nw, :bus, i)[:lam_kcl_i] = cstr_q + end +end + + """ constraint_mc_ohms_yt_from(pm::FOTPUPowerModel, nw::Int, f_bus::Int, t_bus::Int, f_idx::Tuple{Int,Int,Int}, t_idx::Tuple{Int,Int,Int}, f_connections::Vector{Int}, t_connections::Vector{Int}, G::Matrix{<:Real}, B::Matrix{<:Real}, G_fr::Matrix{<:Real}, B_fr::Matrix{<:Real}) diff --git a/src/form/fotr.jl b/src/form/fotr.jl index 0eaae0f4e..f5a216f1a 100644 --- a/src/form/fotr.jl +++ b/src/form/fotr.jl @@ -57,16 +57,6 @@ function variable_mc_bus_voltage(pm::FOTRUPowerModel; nw::Int=nw_id_default, bou end -""" - variable_mc_capcontrol(pm::FOTRUPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - -Capacitor switching variables. -""" -function variable_mc_capcontrol(pm::FOTRUPowerModel; nw::Int=nw_id_default, relax::Bool=false, report::Bool=true) - variable_mc_capacitor_switch_state(pm; nw=nw, relax=relax, report=report) -end - - """ constraint_mc_voltage_magnitude_bounds(pm::FOTRUPowerModel, nw::Int, i::Int, vmin::Vector{<:Real}, vmax::Vector{<:Real}) diff --git a/src/form/ivr.jl b/src/form/ivr.jl index 293dd94c0..4d3bbd20c 100644 --- a/src/form/ivr.jl +++ b/src/form/ivr.jl @@ -308,6 +308,138 @@ function constraint_mc_current_balance(pm::AbstractUnbalancedIVRModel, nw::Int, end +@doc raw""" + constraint_mc_current_balance_capc(pm::AbstractUnbalancedIVRModel, nw::Int, i::Int, terminals::Vector{Int}, grounded::Vector{Bool}, bus_arcs::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_sw::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_trans::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_gens::Vector{Tuple{Int,Vector{Int}}}, bus_storage::Vector{Tuple{Int,Vector{Int}}}, bus_loads::Vector{Tuple{Int,Vector{Int}}}, bus_shunts::Vector{Tuple{Int,Vector{Int}}}) + +Current balance constraints with capacitor control. +""" +function constraint_mc_current_balance_capc(pm::AbstractUnbalancedIVRModel, nw::Int, i::Int, terminals::Vector{Int}, grounded::Vector{Bool}, bus_arcs::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_sw::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_arcs_trans::Vector{Tuple{Tuple{Int,Int,Int},Vector{Int}}}, bus_gens::Vector{Tuple{Int,Vector{Int}}}, bus_storage::Vector{Tuple{Int,Vector{Int}}}, bus_loads::Vector{Tuple{Int,Vector{Int}}}, bus_shunts::Vector{Tuple{Int,Vector{Int}}}) + vr = var(pm, nw, :vr, i) + vi = var(pm, nw, :vi, i) + + cr = get(var(pm, nw), :cr, Dict()); _check_var_keys(cr, bus_arcs, "real current", "branch") + ci = get(var(pm, nw), :ci, Dict()); _check_var_keys(ci, bus_arcs, "imaginary current", "branch") + crd = get(var(pm, nw), :crd_bus, Dict()); _check_var_keys(crd, bus_loads, "real current", "load") + cid = get(var(pm, nw), :cid_bus, Dict()); _check_var_keys(cid, bus_loads, "imaginary current", "load") + crg = get(var(pm, nw), :crg_bus, Dict()); _check_var_keys(crg, bus_gens, "real current", "generator") + cig = get(var(pm, nw), :cig_bus, Dict()); _check_var_keys(cig, bus_gens, "imaginary current", "generator") + crs = get(var(pm, nw), :crs, Dict()); _check_var_keys(crs, bus_storage, "real currentr", "storage") + cis = get(var(pm, nw), :cis, Dict()); _check_var_keys(cis, bus_storage, "imaginary current", "storage") + crsw = get(var(pm, nw), :crsw, Dict()); _check_var_keys(crsw, bus_arcs_sw, "real current", "switch") + cisw = get(var(pm, nw), :cisw, Dict()); _check_var_keys(cisw, bus_arcs_sw, "imaginary current", "switch") + crt = get(var(pm, nw), :crt, Dict()); _check_var_keys(crt, bus_arcs_trans, "real current", "transformer") + cit = get(var(pm, nw), :cit, Dict()); _check_var_keys(cit, bus_arcs_trans, "imaginary current", "transformer") + + # add constraints to model capacitor switching + if !isempty(bus_shunts) && haskey(ref(pm, nw, :shunt, bus_shunts[1][1]), "controls") + constraint_capacitor_on_off(pm, nw, i, bus_shunts) + end + + # calculate Gs, Bs + ncnds = length(terminals) + Gt = fill(0.0, ncnds, ncnds) + Bt = convert(Matrix{JuMP.NonlinearExpression}, JuMP.@NLexpression(pm.model, [idx=1:ncnds, jdx=1:ncnds], 0.0)) + for (val, connections) in bus_shunts + shunt = ref(pm,nw,:shunt,val) + for (idx,c) in enumerate(connections) + cap_state = haskey(shunt,"controls") ? var(pm, nw, :capacitor_state, val)[c] : 1.0 + for (jdx,d) in enumerate(connections) + Gt[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] += shunt["gs"][idx,jdx] + Bt[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] = JuMP.@NLexpression(pm.model, Bt[findfirst(isequal(c), terminals),findfirst(isequal(d), terminals)] + shunt["bs"][idx,jdx]*cap_state) + end + end + end + + ungrounded_terminals = [(idx,t) for (idx,t) in enumerate(terminals) if !grounded[idx]] + + for (idx, t) in ungrounded_terminals + @smart_constraint(pm.model, [cr, crd, crg, crs, crsw, crt, vr, vi, Bt], + sum(cr[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(crsw[a_sw][t] for (a_sw, conns) in bus_arcs_sw if t in conns) + + sum(crt[a_trans][t] for (a_trans, conns) in bus_arcs_trans if t in conns) + == + sum(crg[g][t] for (g, conns) in bus_gens if t in conns) + - sum(crs[s][t] for (s, conns) in bus_storage if t in conns) + - sum(crd[d][t] for (d, conns) in bus_loads if t in conns) + - sum( Gt[idx,jdx]*vr[u] -Bt[idx,jdx]*vi[u] for (jdx,u) in ungrounded_terminals) # shunts + ) + @smart_constraint(pm.model, [ci, cid, cig, cis, cisw, cit, vr, vi, Bt], + sum(ci[a][t] for (a, conns) in bus_arcs if t in conns) + + sum(cisw[a_sw][t] for (a_sw, conns) in bus_arcs_sw if t in conns) + + sum(cit[a_trans][t] for (a_trans, conns) in bus_arcs_trans if t in conns) + == + sum(cig[g][t] for (g, conns) in bus_gens if t in conns) + - sum(cis[s][t] for (s, conns) in bus_storage if t in conns) + - sum(cid[d][t] for (d, conns) in bus_loads if t in conns) + - sum( Gt[idx,jdx]*vi[u] +Bt[idx,jdx]*vr[u] for (jdx,u) in ungrounded_terminals) # shunts + ) + end +end + + +@doc raw""" + constraint_capacitor_on_off(pm::AbstractUnbalancedIVRModel, i::Int; nw::Int=nw_id_default) + +Add constraints to model capacitor switching + +```math +\begin{align} +&\text{kvar control: } s = (vr_fr+im*vi_fr).*(cr_fr-im*ci_fr),\\ +&\text{kvar control (ON): } \Im{s}-q_\text{on} ≤ M_q ⋅ z - ϵ ⋅ (1-z), \\ +&\text{kvar control (OFF): } \Im{s}-q_\text{off} ≥ -M_q ⋅ (1-z) - ϵ ⋅ z, \\ +&\text{voltage control (ON): } v_r^2 + v_i^2 - v_\text{min}^2 ≥ -M_v ⋅ z + ϵ ⋅ (1-z), \\ +&\text{voltage control (OFF): } v_r^2 + v_i^2 - v_\text{max}^2 ≤ M_v ⋅ (1-z) - ϵ ⋅ z. +\end{align} +``` +""" +function constraint_capacitor_on_off(pm::AbstractUnbalancedIVRModel, nw::Int, i::Int, bus_shunts::Vector{Tuple{Int,Vector{Int}}}) + cap_state = var(pm, nw, :capacitor_state, bus_shunts[1][1]) + shunt = ref(pm, nw, :shunt, bus_shunts[1][1]) + ϵ = 1e-5 + M_q = 1e5 + M_v = 2 + elem_type = shunt["controls"]["element"]["type"] + if shunt["controls"]["type"] == CAP_REACTIVE_POWER + bus_idx = shunt["controls"]["terminal"] == 1 ? (shunt["controls"]["element"]["index"], shunt["controls"]["element"]["f_bus"], shunt["controls"]["element"]["t_bus"]) : (shunt["controls"]["element"]["index"], shunt["controls"]["element"]["t_bus"], shunt["controls"]["element"]["f_bus"]) + vr_fr = var(pm, nw, :vr, bus_idx[2]) + vi_fr = var(pm, nw, :vi, bus_idx[2]) + cr_fr = elem_type == "branch" ? var(pm, nw, :cr)[bus_idx] : elem_type == "switch" ? var(pm, nw, :crsw) : var(pm, nw, :crt, bus_idx) + ci_fr = elem_type == "branch" ? var(pm, nw, :ci)[bus_idx] : elem_type == "switch" ? var(pm, nw, :cisw) : var(pm, nw, :cit, bus_idx) + q_fr = JuMP.@expression(pm.model, vi_fr.*cr_fr .- vr_fr.*ci_fr) + JuMP.@constraint(pm.model, sum(q_fr) - shunt["controls"]["onsetting"] ≤ M_q*cap_state[shunt["connections"][1]] - ϵ*(1-cap_state[shunt["connections"][1]])) + JuMP.@constraint(pm.model, sum(q_fr) - shunt["controls"]["offsetting"] ≥ -M_q*(1-cap_state[shunt["connections"][1]]) - ϵ*cap_state[shunt["connections"][1]]) + JuMP.@constraint(pm.model, cap_state .== cap_state[shunt["connections"][1]]) + if shunt["controls"]["voltoverride"] + for (idx,val) in enumerate(shunt["connections"]) + vr_cap = var(pm, nw, :vr, i)[val] + vi_cap = var(pm, nw, :vi, i)[val] + JuMP.@constraint(pm.model, vr_cap^2 + vi_cap^2 - shunt["controls"]["vmin"]^2 ≥ -M_v*cap_state[val] + ϵ*(1-cap_state[val])) + JuMP.@constraint(pm.model, vr_cap^2 + vi_cap^2 - shunt["controls"]["vmax"]^2 ≤ M_v*(1-cap_state[val]) - ϵ*cap_state[val]) + end + end + else + for (idx,val) in enumerate(shunt["connections"]) + if shunt["controls"]["type"][idx] == CAP_VOLTAGE + bus_idx = shunt["controls"]["terminal"][idx] == 1 ? shunt["controls"]["element"]["f_bus"] : shunt["controls"]["element"]["t_bus"] + vr_cap = var(pm, nw, :vr, i)[val] + vi_cap = var(pm, nw, :vi, i)[val] + JuMP.@constraint(pm.model, vr_cap^2 + vi_cap^2 - shunt["controls"]["onsetting"][idx]^2 ≤ M_v*cap_state[val] - ϵ*(1-cap_state[val])) + JuMP.@constraint(pm.model, vr_cap^2 + vi_cap^2 - shunt["controls"]["offsetting"][idx]^2 ≥ -M_v*(1-cap_state[val]) - ϵ*cap_state[val]) + end + if shunt["controls"]["voltoverride"][idx] + vr_cap = var(pm, nw, :vr, i)[val] + vi_cap = var(pm, nw, :vi, i)[val] + JuMP.@constraint(pm.model, vr_cap^2 + vi_cap^2 - shunt["controls"]["vmin"][idx]^2 ≥ -M_v*cap_state[val] + ϵ*(1-cap_state[val])) + JuMP.@constraint(pm.model, vr_cap^2 + vi_cap^2 - shunt["controls"]["vmax"][idx]^2 ≤ M_v*(1-cap_state[val]) - ϵ*cap_state[val]) + end + if shunt["controls"]["type"][idx] == CAP_DISABLED + JuMP.@constraint(pm.model, cap_state[val] == 1 ) + end + end + end +end + + "`p[f_idx]^2 + q[f_idx]^2 <= rate_a^2`" function constraint_mc_thermal_limit_from(pm::AbstractUnbalancedIVRModel, nw::Int, f_idx::Tuple{Int,Int,Int}, f_connections::Vector{Int}, rate_a::Vector{<:Real}) (l, f_bus, t_bus) = f_idx diff --git a/src/prob/opf_capc.jl b/src/prob/opf_capc.jl index 784a103b5..a97a41538 100644 --- a/src/prob/opf_capc.jl +++ b/src/prob/opf_capc.jl @@ -67,7 +67,7 @@ function build_mc_opf_capc(pm::AbstractUnbalancedPowerModel) end -"constructor for branch flow opf" +"constructor for branch flow opf with capcontrol" function build_mc_opf_capc(pm::AbstractUBFModels) # Variables variable_mc_bus_voltage(pm) @@ -124,3 +124,60 @@ function build_mc_opf_capc(pm::AbstractUBFModels) # Objective objective_mc_min_fuel_cost(pm) end + + +"constructor for capcontrol OPF in current-voltage variable space" +function build_mc_opf_capc(pm::AbstractUnbalancedIVRModel) + # Variables + variable_mc_bus_voltage(pm) + variable_mc_branch_current(pm) + variable_mc_switch_current(pm) + variable_mc_transformer_current(pm) + variable_mc_generator_current(pm) + variable_mc_load_current(pm) + + variable_mc_capcontrol(pm; relax=true) + + # Constraints + for i in ids(pm, :ref_buses) + constraint_mc_theta_ref(pm, i) + end + + # gens should be constrained before KCL, or Pd/Qd undefined + for id in ids(pm, :gen) + constraint_mc_generator_power(pm, id) + end + + # loads should be constrained before KCL, or Pd/Qd undefined + for id in ids(pm, :load) + constraint_mc_load_power(pm, id) + end + + for i in ids(pm, :bus) + constraint_mc_current_balance_capc(pm, i) + end + + for i in ids(pm, :branch) + constraint_mc_current_from(pm, i) + constraint_mc_current_to(pm, i) + + constraint_mc_bus_voltage_drop(pm, i) + + constraint_mc_voltage_angle_difference(pm, i) + + constraint_mc_thermal_limit_from(pm, i) + constraint_mc_thermal_limit_to(pm, i) + end + + for i in ids(pm, :switch) + constraint_mc_switch_state(pm, i) + constraint_mc_switch_current_limit(pm, i) + end + + for i in ids(pm, :transformer) + constraint_mc_transformer_power(pm, i) + end + + # Objective + objective_mc_min_fuel_cost(pm) +end diff --git a/test/capacitor.jl b/test/capacitor.jl index 158a9747e..1d320a734 100644 --- a/test/capacitor.jl +++ b/test/capacitor.jl @@ -31,6 +31,21 @@ @test all(isapprox.(result["solution"]["shunt"]["c1"]["cap_state"], [1.0]; atol=2e-1)) end + @testset "capcontrol_ivr" begin + result = solve_mc_opf_capc(IEEE13_CapControl, IVRUPowerModel, ipopt_solver; solution_processors=[sol_data_model!]) + + @test result["termination_status"] == LOCALLY_SOLVED + + @test isapprox(sum(result["solution"]["voltage_source"]["source"]["pg"]), 405.556; atol=5) + @test isapprox(sum(result["solution"]["voltage_source"]["source"]["qg"]), -527.15; atol=300) + + vbase,_ = calc_voltage_bases(IEEE13_CapControl, IEEE13_CapControl["settings"]["vbases_default"]) + @test all(isapprox.(result["solution"]["bus"]["646"]["vm"] ./ vbase["646"], [1.03533, 1.06987]; atol=2e-2)) + @test all(isapprox.(result["solution"]["bus"]["646"]["va"], [-90.1768, 148.498]; atol=3e-1)) + + @test all(isapprox.(result["solution"]["shunt"]["c1"]["cap_state"], [1.0]; atol=2e-1)) + end + @testset "capcontrol_fbs" begin result = solve_mc_opf_capc(IEEE13_CapControl, FBSUBFPowerModel, ipopt_solver; solution_processors=[sol_data_model!]) @@ -74,5 +89,20 @@ @test all(isapprox.(result["solution"]["shunt"]["c1"]["cap_state"], [1.0]; atol=6e-1)) end + + @testset "capcontrol_fotp" begin + result = solve_mc_opf_capc(IEEE13_CapControl, FOTPUPowerModel, ipopt_solver) + + @test result["termination_status"] == LOCALLY_SOLVED + + @test isapprox(sum(result["solution"]["voltage_source"]["source"]["pg"]), 404.784; atol=5) + @test isapprox(sum(result["solution"]["voltage_source"]["source"]["qg"]), -328.146; atol=400) + + vbase,_ = calc_voltage_bases(IEEE13_CapControl, IEEE13_CapControl["settings"]["vbases_default"]) + @test all(isapprox.(result["solution"]["bus"]["646"]["vm"] ./ vbase["646"], [1.03928, 1.05688]; atol=2e-1)) + @test all(isapprox.(result["solution"]["bus"]["646"]["va"], [-89.997, 148.711]; atol=1e0)) + + @test all(isapprox.(result["solution"]["shunt"]["c1"]["cap_state"], [1.0]; atol=6e-1)) + end end