Merge pull request #79 from ghulam41/ACR-Formulation

Add ACR Power Model formulation

docs/src/formulations.md

@@ -7,6 +7,7 @@ For details on `GenericPowerModel`, see _PM.jl [documentation](https://lanl-ansi
 
 Extending PowerModels,  formulations for balanced  OPF in DC grids have been implemented and mapped to the following AC grid formulations:
 - ACPPowerModel
+- ACRPowerModel
 - DCPPowerModel
 - LPACPowerModel
 - SOCWRPowerModel
@@ -45,12 +46,26 @@ Note that generally, $a \geq 0, b \geq 0, c \geq 0$ as physical losses are posit
 
 ### ACDC converters
 - Power balance: $P^{conv, ac}_{ij} + P^{conv, dc}_{ji}$ = $a + b \cdot I^{conv, ac} + c \cdot (I^{conv, ac})^2$.
-- Current variable model: $(P^{conv,ac}_{ij})^2$ + $(Q^{conv,ac}_{ij})^2$ = $U_i^2 \cdot  (I^{conv, ac})^2$.
+- Converter current variable model: $(P^{conv,ac}_{ij})^2$ + $(Q^{conv,ac}_{ij})^2$ = $U_i^2 \cdot  (I^{conv, ac})^2$.
 - LCC converters, active /reactive power:
 
 $P^{conv, ac} = \cos\varphi_{c} \cdot S^{conv,ac,rated}$
 $Q^{conv, ac} = \sin\varphi_{c} \cdot S^{conv,ac,rated}$
 
+
+## ACRPowerModel (BIM)
+### DC lines
+- Same model as ACP formulation
+
+
+### ACDC converters
+Two separate current variables, $I^{conv, ac}$ and $i^{conv, ac, sq}$ are defined, the nonconvex relation $i^{conv, ac, sq} = (I^{conv, ac})^2$ is convexified.
+- Linking both current variables: $(I^{conv, ac})^2$ $\leq$ $i^{conv, ac, sq}$
+- Power balance: $P^{conv, ac}_{ij} + P^{conv, dc}_{ji}$ = $a + b\cdot I^{conv, ac} + c\cdot i^{conv, ac, sq}$.
+- Converter current variable model: $(P^{conv,ac}_{ij})^2$ + $(Q^{conv,ac}_{ij})^2$ = $(U_{ri}^2+U_{ii}^2) \cdot  i^{conv, ac, sq}$.
+- LCC converters, active /reactive power: Same model as ACP formulation
+
+
 ## DCPPowerModel (NF)
 Due to the absence of voltage angles in DC grids, the DC power flow model reduces to network flow (NF) under the 'DC' assumptions
 ### DC lines

docs/src/index.md

@@ -14,6 +14,7 @@ PowerModelsACDC.jl adds new problem types:
 
 PowerModelsACDC.jl extends the formulation hierarchy developed for AC grids, with equivalent DC grid and converter station formulations:
 - ACPPowerModel
+- ACRPowerModel
 - DCPPowerModel
 - SOCWRPowerModel
 - SDPWRMPowerModel
@@ -22,7 +23,7 @@ PowerModelsACDC.jl extends the formulation hierarchy developed for AC grids, wit
 - LPACPowerModel
 
 Developed by:
-- Hakan Ergun, Jay Dave KU Leuven / EnergyVille
+- Hakan Ergun, Jay Dave, KU Leuven / EnergyVille
 - Frederik Geth, CSIRO
 
 ## Installation of PowerModelACDC

src/PowerModelsACDC.jl

@@ -43,6 +43,7 @@ include("core/objective.jl")
 include("core/relaxation_scheme.jl")
 
 include("formdcgrid/acp.jl")
+include("formdcgrid/acr.jl")
 include("formdcgrid/dcp.jl")
 include("formdcgrid/wr.jl")
 include("formdcgrid/wrm.jl")
@@ -52,6 +53,7 @@ include("formdcgrid/shared.jl")
 include("formdcgrid/iv.jl")
 
 include("formconv/acp.jl")
+include("formconv/acr.jl")
 include("formconv/dcp.jl")
 include("formconv/wr.jl")
 include("formconv/wrm.jl")

src/core/variableconv.jl

@@ -237,6 +237,29 @@ function variable_acside_current(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_defa
     report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :iconv_ac_sq, _PM.ids(pm, nw, :convdc), icsq)
 end
 
+"variable: `iconv_ac[j]` and `iconv_ac_sq[j]` for `j` in `convdc` for ACR formulation"
+function variable_acside_current(pm::_PM.AbstractACRModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
+    ic = _PM.var(pm, nw)[:iconv_ac] = JuMP.@variable(pm.model,
+    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_ac",
+    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
+    )
+    icsq = _PM.var(pm, nw)[:iconv_ac_sq] = JuMP.@variable(pm.model,
+    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_ac_sq",
+    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[c],  0)
+            JuMP.set_upper_bound(ic[c],  convdc["Imax"])
+            JuMP.set_lower_bound(icsq[c],  0)
+            JuMP.set_upper_bound(icsq[c],  convdc["Imax"]^2)
+        end
+    end
+
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :iconv_ac, _PM.ids(pm, nw, :convdc), ic)
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :iconv_ac_sq, _PM.ids(pm, nw, :convdc), icsq)
+end
+
 "variable: `itf_sq[j]` for `j` in `convdc`"
 function variable_conv_transformer_current_sqr(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
     bigM = 2; #TODO derive exact bound
@@ -311,6 +334,45 @@ function variable_converter_filter_voltage_angle(pm::_PM.AbstractPowerModel; nw:
 end
 
 
+function variable_converter_filter_voltage(pm::_PM.AbstractACRModel; kwargs...)
+    variable_converter_filter_voltage_real(pm; kwargs...)
+    variable_converter_filter_voltage_imaginary(pm; kwargs...)
+end
+
+
+"real part of the voltage variable `vrf[j]` for `j` in `convdc`"
+function variable_converter_filter_voltage_real(pm::_PM.AbstractACRModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
+    bigM = 1.2; # only internal converter voltage is strictly regulated
+    vrf = _PM.var(pm, nw)[:vrf] = JuMP.@variable(pm.model,
+    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vrf",
+    start = _PM.ref(pm, nw, :convdc, i, "Vtar")
+    )
+    if bounded
+        for (c, convdc) in _PM.ref(pm, nw, :convdc)
+            JuMP.set_lower_bound(vrf[c],  -convdc["Vmmax"] * bigM)   
+            JuMP.set_upper_bound(vrf[c],  convdc["Vmmax"] * bigM)   
+        end
+    end
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vrfilt, _PM.ids(pm, nw, :convdc), vrf)
+end
+
+"imaginary part of the voltage variable `vif[j]` for `j` in `convdc`"
+function variable_converter_filter_voltage_imaginary(pm::_PM.AbstractACRModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
+    bigM = 1.2; 
+    vif = _PM.var(pm, nw)[:vif] = JuMP.@variable(pm.model,
+    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vif",
+    start = 0
+    )
+    if bounded
+        for (c, convdc) in _PM.ref(pm, nw, :convdc)
+            JuMP.set_lower_bound(vif[c], -convdc["Vmmax"] * bigM)
+            JuMP.set_upper_bound(vif[c],  convdc["Vmmax"] * bigM)
+        end
+    end
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vifilt, _PM.ids(pm, nw, :convdc), vif)
+end
+
+
 function variable_converter_internal_voltage(pm::_PM.AbstractPowerModel; kwargs...)
     variable_converter_internal_voltage_magnitude(pm; kwargs...)
     variable_converter_internal_voltage_angle(pm; kwargs...)
@@ -350,6 +412,43 @@ end
 
 
 
+function variable_converter_internal_voltage(pm::_PM.AbstractACRModel; kwargs...)
+    variable_converter_internal_voltage_real(pm; kwargs...)
+    variable_converter_internal_voltage_imaginary(pm; kwargs...)
+end
+
+"real part of the voltage variable `vrc[j]` for `j` in `convdc`"
+function variable_converter_internal_voltage_real(pm::_PM.AbstractACRModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
+    vrc = _PM.var(pm, nw)[:vrc] = JuMP.@variable(pm.model,
+    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vrc",
+    start = _PM.ref(pm, nw, :convdc, i, "Vtar")
+    )
+    if bounded
+        for (c, convdc) in _PM.ref(pm, nw, :convdc)
+            JuMP.set_lower_bound(vrc[c], -convdc["Vmmax"])
+            JuMP.set_upper_bound(vrc[c],  convdc["Vmmax"])
+        end
+    end
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vrconv, _PM.ids(pm, nw, :convdc), vrc)
+end
+
+"imaginary part of the voltage variable `vic[j]` for `j` in `convdc`"
+function variable_converter_internal_voltage_imaginary(pm::_PM.AbstractACRModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
+    vic = _PM.var(pm, nw)[:vic] = JuMP.@variable(pm.model,
+    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vic",
+    start = 0
+    )
+    if bounded
+        for (c, convdc) in _PM.ref(pm, nw, :convdc)
+            JuMP.set_lower_bound(vic[c], -convdc["Vmmax"])
+            JuMP.set_upper_bound(vic[c],  convdc["Vmmax"])
+        end
+    end
+    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :viconv, _PM.ids(pm, nw, :convdc), vic)
+end
+
+
+
 "variable: `wrf_ac[j]` and `wif_ac`  for `j` in `convdc`"
 function variable_converter_filter_voltage_cross_products(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
     bigM = 1.2; # only internal converter voltage is strictly regulated

src/formconv/acr.jl

@@ -0,0 +1,157 @@
+"""
+Creates lossy converter model between AC and DC grid using a lifted converter current magnitude variable
+```
+pconv_ac[i] + pconv_dc[i] == a + b * iconv[i] + c * iconv_sq[i]
+```
+Links the converter current magnitude variable with the squared converter current magnitude variable
+```
+iconv_sq[i] == iconv[i]^2 
+```
+"""
+function constraint_converter_losses(pm::_PM.AbstractACRModel, n::Int, i::Int, a, b, c, plmax)
+    pconv_ac = _PM.var(pm, n, :pconv_ac, i)
+    pconv_dc = _PM.var(pm, n, :pconv_dc, i)
+    iconv = _PM.var(pm, n, :iconv_ac, i)
+    iconv_sq = _PM.var(pm, n, :iconv_ac_sq, i)
+
+    JuMP.@constraint(pm.model, iconv_sq == iconv^2)
+    JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a + b*iconv + c* iconv_sq)
+end
+"""
+Links converter power & current
+```
+pconv_ac[i]^2 + pconv_dc[i]^2 == (vrc[i]^2 + vic[i]^2) * iconv_sq[i]
+```
+"""
+function constraint_converter_current(pm::_PM.AbstractACRModel, n::Int, i::Int, Umax, Imax)
+    vrc = _PM.var(pm, n, :vrc, i)
+    vic = _PM.var(pm, n, :vic, i)
+    pconv_ac = _PM.var(pm, n, :pconv_ac, i)
+    qconv_ac = _PM.var(pm, n, :qconv_ac, i)
+    iconv_sq = _PM.var(pm, n, :iconv_ac_sq, i)
+
+    JuMP.@NLconstraint(pm.model, pconv_ac^2 + qconv_ac^2 == (vrc^2 + vic^2) * iconv_sq)         
+end
+"""
+Converter transformer constraints
+```
+p_tf_fr ==  g/(tm^2)*(vr_fr^2+vi_fr^2) + -g/(tm)*(vr_fr*vr_to + vi_fr*vi_to) + -b/(tm)*(vi_fr*vr_to-vr_fr*vi*to)
+q_tf_fr == -b/(tm^2)*(vr_fr^2+vi_fr^2) +  b/(tm)*(vr_fr*vr_to + vi_fr*vi_to) + -g/(tm)*(vi_fr*vr_to-vr_fr*vi*to)
+p_tf_to ==  g*(vr_to^2+vi_to^2)        + -g/(tm)*(vr_fr*vr_to + vi_fr*vi_to) + -b/(tm)*(-(vi_fr*vr_to-vr_fr*vi*to))
+q_tf_to == -b*(vr_to^2+vi_to^2)        +  b/(tm)*(vr_fr*vr_to + vi_fr*vi_to) + -g/(tm)*(-(vi_fr*vr_to-vr_fr*vi*to))
+```
+"""
+function constraint_conv_transformer(pm::_PM.AbstractACRModel, n::Int, i::Int, rtf, xtf, acbus, tm, transformer)
+    ptf_fr = _PM.var(pm, n, :pconv_tf_fr, i)
+    qtf_fr = _PM.var(pm, n, :qconv_tf_fr, i)
+    ptf_to = _PM.var(pm, n, :pconv_tf_to, i)
+    qtf_to = _PM.var(pm, n, :qconv_tf_to, i)
+
+    vr = _PM.var(pm, n, :vr, acbus)
+    vi = _PM.var(pm, n, :vi, acbus)
+    vrf = _PM.var(pm, n, :vrf, i)      
+    vif = _PM.var(pm, n, :vif, i) 
+
+    ztf = rtf + im*xtf
+    if transformer
+        ytf = 1/(rtf + im*xtf)
+        gtf = real(ytf)
+        btf = imag(ytf)
+
+        JuMP.@NLconstraint(pm.model, ptf_fr ==  gtf / tm^2 * (vr^2 + vi^2)  + -gtf / tm * (vr * vrf + vi * vif) + -btf / tm * (vi * vrf - vr * vif) )
+        JuMP.@NLconstraint(pm.model, qtf_fr == -btf / tm^2 * (vr^2 + vi^2)  - -btf / tm * (vr * vrf + vi * vif) + -gtf / tm * (vi * vrf - vr * vif) )
+        JuMP.@NLconstraint(pm.model, ptf_to ==  gtf * (vrf^2 + vif^2)       + -gtf / tm * (vr * vrf + vi * vif) + -btf / tm * (-(vi * vrf - vr * vif)) )
+        JuMP.@NLconstraint(pm.model, qtf_to == -btf * (vrf^2 + vif^2)       - -btf / tm * (vr * vrf + vi * vif) + -gtf / tm * (-(vi * vrf - vr * vif)) )
+
+    else
+
+        JuMP.@constraint(pm.model, ptf_fr + ptf_to == 0)
+        JuMP.@constraint(pm.model, qtf_fr + qtf_to == 0)
+        JuMP.@constraint(pm.model, vr == vrf)
+        JuMP.@constraint(pm.model, vi == vif)
+    end
+end
+"""
+Converter reactor constraints
+```
+-pconv_ac == gc*(vrc^2 + vic^2) + -gc*(vrc * vrf + vic * vif) + -bc*(vic * vrf - vrc * vif)
+-qconv_ac ==-bc*(vrc^2 + vic^2) +  bc*(vrc * vrf + vic * vif) + -gc*(vic * vrf - vrc * vif)
+p_pr_fr ==  gc *(vrf^2 + vif^2) + -gc *(vrc * vrf + vic * vif) + -bc *(-(vic * vrf - vrc * vif))
+q_pr_fr == -bc *(vrf^2 + vif^2) +  bc *(vrc * vrf + vic * vif) + -gc *(-(vic * vrf - vrc * vif))
+```
+"""
+function constraint_conv_reactor(pm::_PM.AbstractACRModel, n::Int, i::Int, rc, xc, reactor)
+    pconv_ac = _PM.var(pm, n,  :pconv_ac, i)
+    qconv_ac = _PM.var(pm, n,  :qconv_ac, i)
+    ppr_to = - pconv_ac
+    qpr_to = - qconv_ac
+    ppr_fr = _PM.var(pm, n,  :pconv_pr_fr, i)
+    qpr_fr = _PM.var(pm, n,  :qconv_pr_fr, i)
+
+    vrf = _PM.var(pm, n, :vrf, i) 
+    vif = _PM.var(pm, n, :vif, i) 
+    vrc = _PM.var(pm, n, :vrc, i) 
+    vic = _PM.var(pm, n, :vic, i) 
+
+    zc = rc + im*xc
+    if reactor
+        yc = 1/(zc)
+        gc = real(yc)
+        bc = imag(yc)                                      
+        JuMP.@NLconstraint(pm.model, - pconv_ac ==  gc * (vrc^2 + vic^2) + -gc * (vrc * vrf + vic * vif) + -bc * (vic * vrf - vrc * vif))  
+        JuMP.@NLconstraint(pm.model, - qconv_ac == -bc * (vrc^2 + vic^2) +  bc * (vrc * vrf + vic * vif) + -gc * (vic * vrf - vrc * vif)) 
+        JuMP.@NLconstraint(pm.model, ppr_fr ==  gc * (vrf^2 + vif^2) + -gc * (vrc * vrf + vic * vif) + -bc * (-(vic * vrf - vrc * vif)))
+        JuMP.@NLconstraint(pm.model, qpr_fr == -bc * (vrf^2 + vif^2) +  bc * (vrc * vrf + vic * vif) + -gc * (-(vic * vrf - vrc * vif)))
+
+
+    else
+        JuMP.@constraint(pm.model, ppr_fr + ppr_to == 0)
+        JuMP.@constraint(pm.model, qpr_fr + qpr_to == 0)
+        JuMP.@constraint(pm.model, vrc == vrf)
+        JuMP.@constraint(pm.model, vic == vif)
+    end
+end
+"""
+Converter filter constraints
+```
+ppr_fr + ptf_to == 0
+qpr_fr + qtf_to +  (-bv) * filter *(vrf^2 + vif^2) == 0
+```
+"""
+function constraint_conv_filter(pm::_PM.AbstractACRModel, n::Int, i::Int, bv, filter)
+    ppr_fr = _PM.var(pm, n, :pconv_pr_fr, i)
+    qpr_fr = _PM.var(pm, n, :qconv_pr_fr, i)
+    ptf_to = _PM.var(pm, n, :pconv_tf_to, i)
+    qtf_to = _PM.var(pm, n, :qconv_tf_to, i)
+
+    vrf = _PM.var(pm, n, :vrf, i)
+    vif = _PM.var(pm, n, :vif, i)
+
+    JuMP.@constraint(pm.model,   ppr_fr + ptf_to == 0 )
+   JuMP.@constraint(pm.model, qpr_fr + qtf_to +  (-bv) * filter *(vrf^2 + vif^2) == 0)
+
+end
+"""
+LCC firing angle constraints
+```
+pconv_ac == cos(phi) * Srated
+qconv_ac == sin(phi) * Srated
+```
+"""
+function constraint_conv_firing_angle(pm::_PM.AbstractACRModel, n::Int, i::Int, S, P1, Q1, P2, Q2)
+    p = _PM.var(pm, n, :pconv_ac, i)
+    q = _PM.var(pm, n, :qconv_ac, i)
+    phi = _PM.var(pm, n, :phiconv, i)
+
+    JuMP.@NLconstraint(pm.model,   p == cos(phi) * S)
+    JuMP.@NLconstraint(pm.model,   q == sin(phi) * S)
+end
+
+function constraint_dc_droop_control(pm::_PM.AbstractACRModel, n::Int, i::Int, busdc_i, vref_dc, pref_dc, k_droop)
+    pconv_dc = _PM.var(pm, n, :pconv_dc, i)
+    vdc = _PM.var(pm, n, :vdcm, busdc_i)
+
+    JuMP.@constraint(pm.model, pconv_dc == pref_dc - sign(pref_dc) * 1 / k_droop * (vdc - vref_dc))
+end
+
+#################### TNEP Constraints #########################