diff --git a/docs/src/index.md b/docs/src/index.md index c0df8e62..d67fd267 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -32,8 +32,6 @@ The latest stable release of PowerModelACDC can be installed using the Julia pac ```julia Pkg.add("PowerModelsACDC") ``` -The current version of PowerModelsACDC is 0.3.2 and is compatible with PowerModels v0.17.2, InfrastrucureModels v0.5.3 - !!! Important This is a research-grade optimization package. diff --git a/src/PowerModelsACDC.jl b/src/PowerModelsACDC.jl index 9e6e4827..9dc40920 100644 --- a/src/PowerModelsACDC.jl +++ b/src/PowerModelsACDC.jl @@ -25,6 +25,7 @@ __init__() = Memento.register(_LOGGER) include("prob/acdcopf.jl") include("prob/acdcpf.jl") include("prob/acdcopf_bf.jl") +include("prob/acdcopf_iv.jl") include("prob/tnepopf.jl") include("prob/acdctnepopf.jl") include("prob/tnepopf_bf.jl") @@ -52,6 +53,7 @@ include("formdcgrid/wrm.jl") include("formdcgrid/bf.jl") include("formdcgrid/lpac.jl") include("formdcgrid/shared.jl") +include("formdcgrid/iv.jl") include("formconv/acp.jl") include("formconv/dcp.jl") @@ -60,6 +62,7 @@ include("formconv/wrm.jl") include("formconv/bf.jl") include("formconv/lpac.jl") include("formconv/shared.jl") +include("formconv/iv.jl") include("core/constraint_template.jl") include("io/multinetwork.jl") diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl index b6686b39..91dafe90 100644 --- a/src/core/constraint_template.jl +++ b/src/core/constraint_template.jl @@ -20,12 +20,38 @@ function constraint_power_balance_ac(pm::_PM.AbstractPowerModel, i::Int; nw::Int constraint_power_balance_ac(pm, nw, i, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs) end +function constraint_current_balance_ac(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default) + bus = _PM.ref(pm, nw, :bus, i) + bus_arcs = _PM.ref(pm, nw, :bus_arcs, i) + bus_arcs_dc = _PM.ref(pm, nw, :bus_arcs_dc, i) + bus_gens = _PM.ref(pm, nw, :bus_gens, i) + bus_convs_ac = _PM.ref(pm, nw, :bus_convs_ac, i) + bus_loads = _PM.ref(pm, nw, :bus_loads, i) + bus_shunts = _PM.ref(pm, nw, :bus_shunts, i) + + pd = Dict(k => _PM.ref(pm, nw, :load, k, "pd") for k in bus_loads) + qd = Dict(k => _PM.ref(pm, nw, :load, k, "qd") for k in bus_loads) + + gs = Dict(k => _PM.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts) + bs = Dict(k => _PM.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts) + + constraint_current_balance_ac(pm, nw, i, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs) +end + function constraint_power_balance_dc(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default) bus_arcs_dcgrid = _PM.ref(pm, nw, :bus_arcs_dcgrid, i) bus_convs_dc = _PM.ref(pm, nw, :bus_convs_dc, i) pd = _PM.ref(pm, nw, :busdc, i)["Pdc"] constraint_power_balance_dc(pm, nw, i, bus_arcs_dcgrid, bus_convs_dc, pd) end + +function constraint_current_balance_dc(pm::_PM.AbstractIVRModel, i::Int; nw::Int=_PM.nw_id_default) + bus_arcs_dcgrid = _PM.ref(pm, nw, :bus_arcs_dcgrid, i) + bus_convs_dc = _PM.ref(pm, nw, :bus_convs_dc, i) + pd = _PM.ref(pm, nw, :busdc, i)["Pdc"] + constraint_current_balance_dc(pm, nw, bus_arcs_dcgrid, bus_convs_dc, pd) +end + # function constraint_ohms_dc_branch(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default) branch = _PM.ref(pm, nw, :branchdc, i) @@ -271,3 +297,24 @@ function constraint_conv_firing_angle_ne(pm::_PM.AbstractPowerModel, i::Int; nw: Q2 = sin(pi) * S constraint_conv_firing_angle_ne(pm, n, i, S, P1, Q1, P2, Q2) end + + +function constraint_converter_limits(pm::_PM.AbstractIVRModel, i::Int; nw::Int=_PM.nw_id_default) + bigM = 1.1; + vpu = 1; + conv = _PM.ref(pm, nw, :convdc, i) + # pmax = conv["Pacrated"] + # pmin = -conv["Pacrated"] + # qmax = conv["Qacrated"] + # qmin = -conv["Qacrated"] + # pmaxdc = conv["Pacrated"] * bigM + # pmindc = -conv["Pacrated"] * bigM + imax = conv["Pacrated"]/vpu + vmax = conv["Vmmax"] + vmin = conv["Vmmin"] + pdcmin = -conv["Pacrated"] * bigM # to account for losses + pdcmax = conv["Pacrated"] * bigM # to account for losses + b_idx = conv["busdc_i"] + + constraint_converter_limits(pm, nw, i, imax, vmax, vmin, b_idx, pdcmin, pdcmax) +end diff --git a/src/formconv/iv.jl b/src/formconv/iv.jl new file mode 100644 index 00000000..af7778eb --- /dev/null +++ b/src/formconv/iv.jl @@ -0,0 +1,461 @@ +#################### COLLECT ALL CONVERTER VARIABLES ###################### +function variable_dc_converter(pm::_PM.AbstractIVRModel; kwargs...) + variable_filter_voltage_real(pm; kwargs...) # + variable_filter_voltage_imaginary(pm; kwargs...) # + variable_converter_voltage_real(pm; kwargs...) # + variable_converter_voltage_imaginary(pm; kwargs...) # + variable_transformer_current_real_from(pm; kwargs...) # + variable_transformer_current_real_to(pm; kwargs...)# + variable_transformer_current_imaginary_from(pm; kwargs...) # + variable_transformer_current_imaginary_to(pm; kwargs...)# + variable_reactor_current_real_from(pm; kwargs...)# + variable_reactor_current_real_to(pm; kwargs...)# + variable_reactor_current_imaginary_from(pm; kwargs...)# + variable_reactor_current_imaginary_to(pm; kwargs...)# + variable_converter_current_real(pm; kwargs...)# + variable_converter_current_imaginary(pm; kwargs...)# + variable_converter_current_dc(pm; kwargs...)# + variable_converter_current_lin(pm; kwargs...)# + variable_converter_active_power(pm; kwargs...) + variable_converter_reactive_power(pm; kwargs...) + variable_dcside_power(pm; kwargs...) +end + +########### CONVERTER AC SIDE VOLTAGES ############################## +"real part of the voltage variable k at the filter bus" +function variable_filter_voltage_real(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true) + vk_r = _PM.var(pm, nw)[:vk_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vk_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "vr_start", 1.0) + ) + + if bounded + for (i, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(vk_r[i], -convdc["Vmmax"]) + JuMP.set_upper_bound(vk_r[i], convdc["Vmmax"]) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :vk_r, _PM.ids(pm, nw, :convdc), vk_r) +end + +"imaginary part of the voltage variable k at the filter bus" +function variable_filter_voltage_imaginary(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true) + vk_i = _PM.var(pm, nw)[:vk_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vk_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "vi_start", 1.0) + ) + + if bounded + for (i, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(vk_i[i], -convdc["Vmmax"]) + JuMP.set_upper_bound(vk_i[i], convdc["Vmmax"]) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :vk_i, _PM.ids(pm, nw, :convdc), vk_i) +end + +"real part of the voltage variable c at the PE converter bus" +function variable_converter_voltage_real(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true) + vc_r = _PM.var(pm, nw)[:vc_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vc_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "v_start", 1.0) + ) + + if bounded + for (i, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(vc_r[i], -convdc["Vmmax"]) + JuMP.set_upper_bound(vc_r[i], convdc["Vmmax"]) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :vc_r, _PM.ids(pm, nw, :convdc), vc_r) +end + +"imaginary part of the voltage variable c at the PE converter bus" +function variable_converter_voltage_imaginary(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool=true, report::Bool=true) + vc_i = _PM.var(pm, nw)[:vc_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vc_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "vi_start", 1.0) + ) + + if bounded + for (i, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(vc_i[i], -convdc["Vmmax"]) + JuMP.set_upper_bound(vc_i[i], convdc["Vmmax"]) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :vc_i, _PM.ids(pm, nw, :convdc), vc_i) +end +################### CONVERTER AC SIDE CURRENT VARIABLES ######################################### +function variable_transformer_current_real_from(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + iik_r = _PM.var(pm, nw)[:iik_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iik_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(iik_r[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(iik_r[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :iik_r, _PM.ids(pm, nw, :convdc), iik_r) +end + +function variable_transformer_current_imaginary_from(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + iik_i = _PM.var(pm, nw)[:iik_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iik_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(iik_i[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(iik_i[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :iik_i, _PM.ids(pm, nw, :convdc), iik_i) +end + +function variable_transformer_current_real_to(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + iki_r = _PM.var(pm, nw)[:iki_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iki_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(iki_r[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(iki_r[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :iki_r, _PM.ids(pm, nw, :convdc), iki_r) +end + +function variable_transformer_current_imaginary_to(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + iki_i = _PM.var(pm, nw)[:iki_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iki_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(iki_i[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(iki_i[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :iki_i, _PM.ids(pm, nw, :convdc), iki_i) +end + +function variable_reactor_current_real_from(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + ikc_r = _PM.var(pm, nw)[:ikc_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_ikc_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(ikc_r[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(ikc_r[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :ikc_r, _PM.ids(pm, nw, :convdc), ikc_r) +end + +function variable_reactor_current_imaginary_from(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + ikc_i = _PM.var(pm, nw)[:ikc_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_ikc_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(ikc_i[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(ikc_i[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :ikc_i, _PM.ids(pm, nw, :convdc), ikc_i) +end + +function variable_reactor_current_real_to(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + ick_r = _PM.var(pm, nw)[:ick_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_ick_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(ick_r[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(ick_r[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :ick_r, _PM.ids(pm, nw, :convdc), ick_r) +end + +function variable_reactor_current_imaginary_to(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + ick_i = _PM.var(pm, nw)[:ick_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_ick_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(ick_i[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(ick_i[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :ick_i, _PM.ids(pm, nw, :convdc), ick_i) +end + +function variable_converter_current_real(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + ic_r = _PM.var(pm, nw)[:ic_r] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_ic_r", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(ic_r[c], -convdc["Imax"] * bigM) + JuMP.set_upper_bound(ic_r[c], convdc["Imax"]* bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :ic_r, _PM.ids(pm, nw, :convdc), ic_r) +end + +function variable_converter_current_imaginary(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + ic_i = _PM.var(pm, nw)[:ic_i] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_ic_i", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(ic_i[c], -convdc["Imax"] * bigM) + JuMP.set_upper_bound(ic_i[c], convdc["Imax"] * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :ic_i, _PM.ids(pm, nw, :convdc), ic_i) +end + +function variable_converter_current_dc(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1.2; + vpu = 1; + iconv_dc = _PM.var(pm, nw)[:iconv_dc] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_dc", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(iconv_dc[c], -convdc["Pacrated"]/vpu * bigM) + JuMP.set_upper_bound(iconv_dc[c], convdc["Pacrated"]/vpu * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :iconv_dc, _PM.ids(pm, nw, :convdc), iconv_dc) +end + +function variable_converter_current_lin(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + bigM = 1; + vpu = 1; + iconv_lin = _PM.var(pm, nw)[:iconv_lin] = JuMP.@variable(pm.model, + [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_lin", + start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0) + ) + if bounded + for (c, convdc) in _PM.ref(pm, nw, :convdc) + JuMP.set_lower_bound(iconv_lin[c], 0) + JuMP.set_upper_bound(iconv_lin[c], convdc["Imax"] * bigM) + end + end + + report && _PM.sol_component_value(pm, nw, :convdc, :iconv_lin, _PM.ids(pm, nw, :convdc), iconv_lin) +end + +#################### CONVERTER CURRENT LIMITS ######################### + +function constraint_converter_limits(pm::_PM.AbstractIVRModel, n::Int, i, imax, vmax, vmin, b_idx, pdcmin, pdcmax) + iik_r = _PM.var(pm, n, :iik_r)[i] + iik_i = _PM.var(pm, n, :iik_i)[i] + iki_r = _PM.var(pm, n, :iki_r)[i] + iki_i = _PM.var(pm, n, :iki_i)[i] + + ikc_r = _PM.var(pm, n, :ikc_r)[i] + ikc_i = _PM.var(pm, n, :ikc_i)[i] + ick_r = _PM.var(pm, n, :ick_r)[i] + ick_i = _PM.var(pm, n, :ick_i)[i] + + ic_r = _PM.var(pm, n, :ic_r)[i] + ic_i = _PM.var(pm, n, :ic_i)[i] + + iconv_lin = _PM.var(pm, n, :iconv_lin)[i] + + vk_r = _PM.var(pm, n, :vk_r)[i] + vk_i = _PM.var(pm, n, :vk_i)[i] + vc_r = _PM.var(pm, n, :vc_r)[i] + vc_i = _PM.var(pm, n, :vc_i)[i] + + JuMP.@NLconstraint(pm.model, (iik_r)^2 + (iik_i)^2 <= imax^2) #(32) + JuMP.@NLconstraint(pm.model, (iki_r)^2 + (iki_i)^2 <= imax^2) #(33) + JuMP.@NLconstraint(pm.model, (ikc_r)^2 + (ikc_i)^2 <= imax^2) #(34) + JuMP.@NLconstraint(pm.model, (ick_r)^2 + (ick_i)^2 <= imax^2) #(35) + JuMP.@NLconstraint(pm.model, (ic_r)^2 + (ic_i)^2 <= imax^2) #(39) + JuMP.@NLconstraint(pm.model, (ic_r)^2 + (ic_i)^2 == (iconv_lin)^2) #(47) + ## relaxed version + #JuMP.@NLconstraint(pm.model, (ic_r)^2 + (ic_i)^2 == (iconv_lin)^2) #(49) + JuMP.@NLconstraint(pm.model, (vk_r)^2 + (vk_i)^2 <= vmax^2) #(22) + JuMP.@NLconstraint(pm.model, (vk_r)^2 + (vk_i)^2 >= vmin^2) #(22) + JuMP.@NLconstraint(pm.model, (vc_r)^2 + (vc_i)^2 <= vmax^2) #(23) + JuMP.@NLconstraint(pm.model, (vc_r)^2 + (vc_i)^2 >= vmin^2) #(23) + + + vc = _PM.var(pm, n, :vdcm)[b_idx] + pconv_dc = _PM.var(pm, n, :pconv_dc)[i] + iconv_dc = _PM.var(pm, n, :iconv_dc)[i] + + JuMP.@NLconstraint(pm.model, pconv_dc == vc * iconv_dc) +end + + +######## Constraint converter losses ###### +function constraint_converter_losses(pm::_PM.AbstractIVRModel, n::Int, i::Int, a, b, c, plmax) + iconv_lin = _PM.var(pm, n, :iconv_lin, i) + pconv_ac = _PM.var(pm, n, :pconv_ac, i) + pconv_dc = _PM.var(pm, n, :pconv_dc, i) + + JuMP.@NLconstraint(pm.model, pconv_ac + pconv_dc == a + b*iconv_lin + c*(iconv_lin)^2) +end + +function constraint_converter_current(pm::_PM.AbstractIVRModel, n::Int, i::Int, Umax, Imax) + vc_r = _PM.var(pm, n, :vc_r, i) + vc_i = _PM.var(pm, n, :vc_i, i) + ic_r = _PM.var(pm, n, :ic_r, i) + ic_i = _PM.var(pm, n, :ic_i, i) + pconv_ac = _PM.var(pm, n, :pconv_ac, i) + qconv_ac = _PM.var(pm, n, :qconv_ac, i) + + + + JuMP.@NLconstraint(pm.model, pconv_ac == vc_r * ic_r + vc_i * ic_i) + JuMP.@NLconstraint(pm.model, qconv_ac == vc_i * ic_r - vc_r * ic_i) +end + +function constraint_conv_transformer(pm::_PM.AbstractIVRModel, n::Int, i::Int, rtf, xtf, acbus, tm, transformer) + vi_r = _PM.var(pm, n, :vr, acbus) + vi_i = _PM.var(pm, n, :vi, acbus) + vk_r = _PM.var(pm, n, :vk_r, i) + vk_i = _PM.var(pm, n, :vk_i, i) + + iik_r = _PM.var(pm, n, :iik_r, i) + iik_i = _PM.var(pm, n, :iik_i, i) + iki_r = _PM.var(pm, n, :iki_r, i) + iki_i = _PM.var(pm, n, :iki_i, i) + + #TODO add transformation ratio..... + if transformer + JuMP.@constraint(pm.model, vk_r == vi_r - rtf * iik_r + xtf * iik_i) #(24) + JuMP.@constraint(pm.model, vk_i == vi_i - rtf * iik_i - xtf * iik_r) #(25) + JuMP.@constraint(pm.model, vi_r == vk_r - rtf * iki_r + xtf * iki_i) #reverse + JuMP.@constraint(pm.model, vi_i == vk_i - rtf * iki_i - xtf * iki_r) #reverse + else + JuMP.@constraint(pm.model, vk_r == vi_r) + JuMP.@constraint(pm.model, vk_i == vi_i) + JuMP.@constraint(pm.model, iik_r + iki_r == 0) + JuMP.@constraint(pm.model, iik_i + iki_i == 0) + end +end + +function constraint_conv_reactor(pm::_PM.AbstractIVRModel, n::Int, i::Int, rc, xc, reactor) + vk_r = _PM.var(pm, n, :vk_r, i) + vk_i = _PM.var(pm, n, :vk_i, i) + vc_r = _PM.var(pm, n, :vc_r, i) + vc_i = _PM.var(pm, n, :vc_i, i) + + ikc_r = _PM.var(pm, n, :ikc_r, i) + ikc_i = _PM.var(pm, n, :ikc_i, i) + ick_r = _PM.var(pm, n, :ick_r, i) + ick_i = _PM.var(pm, n, :ick_i, i) + ic_r = _PM.var(pm, n, :ic_r, i) + ic_i = _PM.var(pm, n, :ic_i, i) + + JuMP.@constraint(pm.model, ick_r + ic_r == 0) #(20) + JuMP.@constraint(pm.model, ick_i + ic_i == 0) #(21) + + if reactor + JuMP.@constraint(pm.model, vc_r == vk_r - rc * ikc_r + xc * ikc_i) #(28) + JuMP.@constraint(pm.model, vc_i == vk_i - rc * ikc_i - xc * ikc_r) #(29) + JuMP.@constraint(pm.model, vk_r == vc_r - rc * ick_r + xc * ick_i) #reverse + JuMP.@constraint(pm.model, vk_i == vc_i - rc * ick_i - xc * ick_r) #reverse + else + JuMP.@constraint(pm.model, vk_r == vc_r) + JuMP.@constraint(pm.model, vk_i == vc_i) + JuMP.@constraint(pm.model, ikc_r + ick_r == 0) + JuMP.@constraint(pm.model, ikc_i + ick_i == 0) + end +end + +function constraint_conv_filter(pm::_PM.AbstractIVRModel, n::Int, i::Int, bv, filter) + iki_r = _PM.var(pm, n, :iki_r, i) + iki_i = _PM.var(pm, n, :iki_i, i) + ikc_r = _PM.var(pm, n, :ikc_r, i) + ikc_i = _PM.var(pm, n, :ikc_i, i) + + vk_r = _PM.var(pm, n, :vk_r, i) + vk_i = _PM.var(pm, n, :vk_i, i) + + JuMP.@constraint(pm.model, iki_r + ikc_r + bv * filter * vk_i == 0) + JuMP.@constraint(pm.model, iki_i + ikc_i - bv * filter * vk_r == 0) +end + +################# Kicrchhoff's current law ############################################ + +function constraint_current_balance_ac(pm::_PM.AbstractIVRModel, n::Int, i, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, bus_pd, bus_qd, bus_gs, bus_bs) + vr = _PM.var(pm, n, :vr, i) + vi = _PM.var(pm, n, :vi, i) + + cr = _PM.var(pm, n, :cr) + ci = _PM.var(pm, n, :ci) + cidc = _PM.var(pm, n, :cidc) + + iik_r = _PM.var(pm, n, :iik_r) + iik_i = _PM.var(pm, n, :iik_i) + + crg = _PM.var(pm, n, :crg) + cig = _PM.var(pm, n, :cig) + + JuMP.@NLconstraint(pm.model, sum(cr[a] for a in bus_arcs) + sum(iik_r[c] for c in bus_convs_ac) + == + sum(crg[g] for g in bus_gens) + - (sum(pd for pd in values(bus_pd))*vr + sum(qd for qd in values(bus_qd))*vi)/(vr^2 + vi^2) + - sum(gs for gs in values(bus_gs))*vr + sum(bs for bs in values(bus_bs))*vi + ) + JuMP.@NLconstraint(pm.model, sum(ci[a] for a in bus_arcs) + sum(iik_i[c] for c in bus_convs_ac) + + sum(cidc[d] for d in bus_arcs_dc) + == + sum(cig[g] for g in bus_gens) + - (sum(pd for pd in values(bus_pd))*vi - sum(qd for qd in values(bus_qd))*vr)/(vr^2 + vi^2) + - sum(gs for gs in values(bus_gs))*vi - sum(bs for bs in values(bus_bs))*vr + ) +end \ No newline at end of file diff --git a/src/formdcgrid/bf.jl b/src/formdcgrid/bf.jl index 8d995015..afb98865 100644 --- a/src/formdcgrid/bf.jl +++ b/src/formdcgrid/bf.jl @@ -48,6 +48,12 @@ function constraint_dc_branch_current(pm::_PM.AbstractBFModel, n::Int, f_bus, f_ # do nothing end +function constraint_dc_branch_current(pm::_PM.AbstractIVRModel, n::Int, f_bus, f_idx, ccm_max, p) + variable_dcbranch_current_iv(pm; kwargs...) +end + + + ########## TNEP constraints ##################### """ diff --git a/src/formdcgrid/iv.jl b/src/formdcgrid/iv.jl new file mode 100644 index 00000000..029074e7 --- /dev/null +++ b/src/formdcgrid/iv.jl @@ -0,0 +1,44 @@ +# Explicit DC branch current variable +function variable_dcbranch_current(pm::_PM.AbstractIVRModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true) + vpu = 1; + igrid_dc = _PM.var(pm, nw)[:igrid_dc] = JuMP.@variable(pm.model, + [(l,i,j) in _PM.ref(pm, nw, :arcs_dcgrid)], base_name="$(nw)_igrid_dc", + start = (_PM.comp_start_value(_PM.ref(pm, nw, :branchdc, l), "p_start", 0.0) / vpu) + ) + if bounded + for arc in _PM.ref(pm, nw, :arcs_dcgrid) + l,i,j = arc + JuMP.set_lower_bound(igrid_dc[arc], -_PM.ref(pm, nw, :branchdc, l)["rateA"] / vpu) + JuMP.set_upper_bound(igrid_dc[arc], _PM.ref(pm, nw, :branchdc, l)["rateA"] / vpu) + end + end + report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :branchdc, :if, :it, _PM.ref(pm, nw, :arcs_dcgrid_from), _PM.ref(pm, nw, :arcs_dcgrid_to), igrid_dc) +end + +# Kirchhoff's current law for DC nodes +function constraint_current_balance_dc(pm::_PM.AbstractIVRModel, n::Int, bus_arcs_dcgrid, bus_convs_dc, pd) + igrid_dc = _PM.var(pm, n, :igrid_dc) + iconv_dc = _PM.var(pm, n, :iconv_dc) + + JuMP.@constraint(pm.model, sum(igrid_dc[a] for a in bus_arcs_dcgrid) + sum(iconv_dc[c] for c in bus_convs_dc) == 0) # deal with pd +end + +# Ohm's law for DC branches +function constraint_ohms_dc_branch(pm::_PM.AbstractIVRModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p) + i_dc_fr = _PM.var(pm, n, :igrid_dc, f_idx) + i_dc_to = _PM.var(pm, n, :igrid_dc, t_idx) + vmdc_fr = _PM.var(pm, n, :vdcm, f_bus) + vmdc_to = _PM.var(pm, n, :vdcm, t_bus) + p_fr = _PM.var(pm, n, :p_dcgrid, f_idx) + p_to = _PM.var(pm, n, :p_dcgrid, t_idx) + + if r == 0 + JuMP.@constraint(pm.model, i_dc_fr + i_dc_to == 0) + else + JuMP.@constraint(pm.model, vmdc_to == vmdc_fr - 1/p * r * i_dc_fr) + JuMP.@constraint(pm.model, vmdc_fr == vmdc_to - 1/p * r * i_dc_to) + end + + JuMP.@NLconstraint(pm.model, p_fr == vmdc_fr * i_dc_fr) + JuMP.@NLconstraint(pm.model, p_to == vmdc_to * i_dc_to) +end \ No newline at end of file diff --git a/src/prob/acdcopf_iv.jl b/src/prob/acdcopf_iv.jl new file mode 100644 index 00000000..573debf9 --- /dev/null +++ b/src/prob/acdcopf_iv.jl @@ -0,0 +1,62 @@ +"" +function run_acdcopf_iv(file::String, model_type, optimizer; kwargs...) + data = _PM.parse_file(file) + PowerModelsACDC.process_additional_data!(data) + return run_acdcopf_iv(data, model_type, optimizer; ref_extensions = [add_ref_dcgrid!], kwargs...) +end + +function run_acdcopf_iv(data::Dict{String,Any}, model_type::Type, optimizer; kwargs...) + return _PM.run_model(data, model_type, optimizer, build_acdcopf_iv; ref_extensions = [add_ref_dcgrid!], kwargs...) +end + +"" +function build_acdcopf_iv(pm::_PM.AbstractIVRModel) + _PM.variable_bus_voltage(pm) + _PM.variable_branch_current(pm) + + _PM.variable_gen_current(pm) + _PM.variable_dcline_current(pm) + + _PM.objective_min_fuel_and_flow_cost(pm) + + #DC grid variables + variable_active_dcbranch_flow(pm) + variable_dcbranch_current(pm) + variable_dcgrid_voltage_magnitude(pm) + #DC converter variables + variable_dc_converter(pm) + + for i in _PM.ids(pm, :ref_buses) + _PM.constraint_theta_ref(pm, i) + end + + for i in _PM.ids(pm, :bus) + constraint_current_balance_ac(pm, i) + end + + for i in _PM.ids(pm, :branch) + _PM.constraint_current_from(pm, i) + _PM.constraint_current_to(pm, i) + + _PM.constraint_voltage_drop(pm, i) + _PM.constraint_voltage_angle_difference(pm, i) + + _PM.constraint_thermal_limit_from(pm, i) + _PM.constraint_thermal_limit_to(pm, i) + end + + for i in _PM.ids(pm, :busdc) + constraint_current_balance_dc(pm, i) + end + for i in _PM.ids(pm, :branchdc) + constraint_ohms_dc_branch(pm, i) + end + for i in _PM.ids(pm, :convdc) + constraint_converter_limits(pm, i) + constraint_converter_losses(pm, i) + constraint_converter_current(pm, i) + constraint_conv_transformer(pm, i) + constraint_conv_reactor(pm, i) + constraint_conv_filter(pm, i) + end +end \ No newline at end of file diff --git a/test/opf.jl b/test/opf.jl index 2b0502ab..ed7d70b1 100644 --- a/test/opf.jl +++ b/test/opf.jl @@ -25,6 +25,22 @@ s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true) end end +@testset "test IVR OPF" begin + @testset "5-bus ac dc case" begin + result = PowerModelsACDC.run_acdcopf_iv("../test/data/case5_acdc.m", IVRPowerModel, ipopt_solver; setting = s) + + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["objective"], 194.16; atol = 1e0) + end + + @testset "39-bus ac dc case" begin + result = PowerModelsACDC.run_acdcopf_iv("../test/data/case39_acdc.m", IVRPowerModel, ipopt_solver; setting = s) + + @test result["termination_status"] == LOCALLY_SOLVED + @test isapprox(result["objective"], 41968.88; atol = 1e0) + end +end + @testset "test dc opf" begin @testset "3-bus case" begin result = run_acdcopf("../test/data/case3.m", DCPPowerModel, ipopt_solver; setting = s) diff --git a/test/scripts/ACDCIVRtest.jl b/test/scripts/ACDCIVRtest.jl new file mode 100644 index 00000000..75adc8e2 --- /dev/null +++ b/test/scripts/ACDCIVRtest.jl @@ -0,0 +1,74 @@ +import PowerModelsACDC +const _PMACDC = PowerModelsACDC +import PowerModels +const _PM = PowerModels +import Ipopt +import Memento +import JuMP + + +file = "./test/data/case3120sp_acdc.m" + + +data = _PM.parse_file(file) + +_PMACDC.process_additional_data!(data) + +ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6, "print_level" => 0) + +s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true) + +resultAC = _PMACDC.run_acdcopf(file, _PM.ACPPowerModel, ipopt; setting = s) +resultACPM = _PM.run_opf(file, _PM.ACPPowerModel, ipopt; setting = s) +resultIVR = _PMACDC.run_acdcopf_iv(file, _PM.IVRPowerModel, ipopt; setting = s) + +print("ACP RESULTS") +print("Objective:", resultAC["objective"],"\n") +for (c, conv) in resultAC["solution"]["convdc"] + ploss = conv["pconv"] + conv["pdc"] + ploss_tot = conv["pgrid"] + conv["pdc"] + print("Pac: ", conv["pconv"], " Pdc: ", conv["pdc"], " Ploss: ", ploss, " Plosstot: ", ploss_tot, "\n") +end + +print("IVR RESULTS") +print("Objective:", resultIVR["objective"],"\n") +for (c, conv) in resultIVR["solution"]["convdc"] + ploss = conv["pconv"] + conv["pdc"] + bus = data["convdc"][c]["busac_i"] + ploss_tot = (conv["iik_r"] * resultIVR["solution"]["bus"]["$bus"]["vr"] + conv["iik_i"] * resultIVR["solution"]["bus"]["$bus"]["vi"]) + conv["pdc"] + print("Pac: ", conv["pconv"], " Pdc: ", conv["pdc"], " Ploss: ", ploss, " Plosstot: ", ploss_tot, "\n") +end + +# il = (data["load"]["1"]["pd"] * resultIVR["solution"]["bus"]["2"]["vr"]+data["load"]["1"]["qd"] * resultIVR["solution"]["bus"]["2"]["vi"]) / (resultIVR["solution"]["bus"]["2"]["vr"]^2 + resultIVR["solution"]["bus"]["2"]["vi"]^2) + + +# so = resultIVR["solution"] +# ir_balance = (so["gen"]["2"]["crg"] - il) - (so["convdc"]["1"]["iik_r"] + so["branch"]["1"]["cr_to"] + so["branch"]["3"]["cr_fr"] + so["branch"]["4"]["cr_fr"] + so["branch"]["5"]["cr_fr"]) + + +# for (b, bus) in resultIVR["solution"]["bus"] +# vm = sqrt(bus["vi"]^2 + bus["vr"]^2) +# print(b, ": ", vm, "\n") +# end + +# for (c, conv) in resultIVR["solution"]["convdc"] +# iik = sqrt(conv["iik_r"]^2 + conv["iik_i"]^2) +# ikc = sqrt(conv["ikc_r"]^2 + conv["ikc_i"]^2) +# ic = sqrt(conv["ic_r"]^2 + conv["ic_i"]^2) +# print(c, "-iik: ", iik, "\n") +# print(c, "-ikc: ", ikc, "\n") +# print(c, "-ic: ", ic, "\n") + +# vk = sqrt(conv["vk_r"]^2 + conv["vk_i"]^2) +# vc = sqrt(conv["vc_r"]^2 + conv["vc_i"]^2) +# print(c, "-vk: ", vk, "\n") +# print(c, "-vc: ", vc, "\n") +# end +# pg_acp = resultAC["solution"]["gen"]["1"]["pg"] + resultAC["solution"]["gen"]["2"]["pg"] +# pg_ivr = resultIVR["solution"]["gen"]["1"]["pg"] + resultIVR["solution"]["gen"]["2"]["pg"] + +# print("pac: ", pg_acp, "\n") +# print("pivr: ", pg_ivr, "\n") + + +