Skip to content

Commit c6ba865

Browse files
CarlosEpiaCarlosEpia
committed
Merge branch 'feature/14a-fixes' into project/411_LoMa_14aOptimization_with_virtual_generators
2 parents b678f9b + d377dbc commit c6ba865

9 files changed

Lines changed: 928 additions & 371 deletions

File tree

‎edisgo/io/powermodels_io.py‎

Lines changed: 421 additions & 121 deletions
Large diffs are not rendered by default.

‎edisgo/opf/eDisGo_OPF.jl/src/core/constraint_cp_14a.jl‎

Lines changed: 54 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,14 @@ end
3434
Ensures that the net electrical load (charging point load - virtual generator support)
3535
stays above the §14a minimum power level (typically 4.2 kW = 0.0042 MW).
3636
37+
For flexible CPs (in electromobility dict): uses optimization variable `pcp`.
38+
For fixed CPs (in load dict): uses fixed load parameter `load["pd"]`.
39+
40+
Big-M formulation:
41+
p_cp - p_cp14a >= p_min_14a - M * (1 - z_cp14a)
42+
When z=0 (inactive): constraint relaxed
43+
When z=1 (active): p_cp - p_cp14a >= p_min_14a (net load >= 4.2 kW)
44+
3745
# Arguments
3846
- `pm::AbstractBFModelEdisgo`: PowerModels model
3947
- `i::Int`: Virtual generator index
@@ -42,75 +50,69 @@ stays above the §14a minimum power level (typically 4.2 kW = 0.0042 MW).
4250
function constraint_cp_14a_min_net_load(pm::AbstractBFModelEdisgo, i::Int, nw::Int=nw_id_default)
4351
gen_cp14a = PowerModels.ref(pm, nw, :gen_cp_14a, i)
4452
cp_idx = gen_cp14a["cp_index"]
45-
53+
54+
# Virtual generator support variable
55+
p_cp14a = PowerModels.var(pm, nw, :p_cp14a, i)
56+
57+
# §14a minimum power (per unit)
58+
p_min_14a = gen_cp14a["p_min_14a"]
59+
60+
# Maximum support capacity
61+
p_max_support = gen_cp14a["pmax"]
62+
63+
if p_max_support < 1e-6
64+
# Charging point too small for §14a curtailment, disable virtual generator
65+
JuMP.@constraint(pm.model, p_cp14a == 0.0)
66+
return
67+
end
68+
69+
# Binary variable for Big-M formulation
70+
z_cp14a = PowerModels.var(pm, nw, :z_cp14a, i)
71+
M = p_max_support + p_min_14a
72+
4673
# Check if CP is flexible (in electromobility dict) or simple load
47-
p_cp_load = nothing
4874
if haskey(PowerModels.ref(pm, nw), :electromobility) && haskey(PowerModels.ref(pm, nw, :electromobility), cp_idx)
49-
# Flexible CP: use electromobility variable
50-
cp = PowerModels.ref(pm, nw, :electromobility, cp_idx)
51-
p_cp_load = cp["pcp"] # Charging power variable (optimization variable)
75+
# Flexible CP: use optimization VARIABLE
76+
pcp = PowerModels.var(pm, nw, :pcp, cp_idx)
77+
JuMP.@constraint(pm.model, pcp - p_cp14a >= p_min_14a - M * (1 - z_cp14a))
5278
else
5379
# Non-flexible CP: use fixed load timeseries
54-
# Find the load by CP name
80+
# ============================================
81+
# NOTE: This is an O(n) linear search through all loads
82+
#
83+
# Why this search is necessary:
84+
# - cp_index may reference electromobility dict (flexible CPs with optimization vars)
85+
# - OR it may reference load dict (fixed CPs with constant power draw)
86+
# - No unified index mapping exists between these two data structures
87+
# - Load name matching is the only reliable way to find the correct load
88+
#
89+
# Performance impact:
90+
# - Acceptable for <1,000 loads per network
91+
# - For larger networks, consider pre-indexing loads by name in Python
92+
# before serializing to PowerModels dict
93+
# ============================================
5594
cp_name = gen_cp14a["cp_name"]
56-
# Search for load with matching name
5795
load_found = false
5896
for (load_id, load) in PowerModels.ref(pm, nw, :load)
5997
if haskey(load, "name") && load["name"] == cp_name
60-
p_cp_load = load["pd"] # Fixed load value (parameter, not variable)
98+
p_cp_load = load["pd"]
6199
load_found = true
100+
if p_cp_load > 1e-6
101+
# Fixed load with Big-M: when z=1, enforce min net load
102+
JuMP.@constraint(pm.model, p_cp_load - p_cp14a >= p_min_14a - M * (1 - z_cp14a))
103+
else
104+
# CP is off (p_cp_load ≈ 0), no support needed
105+
JuMP.@constraint(pm.model, p_cp14a == 0.0)
106+
end
62107
break
63108
end
64109
end
65-
110+
66111
if !load_found
67112
@warn "Could not find load for charging point $(cp_name), skipping constraint"
68113
return
69114
end
70115
end
71-
72-
# Virtual generator support
73-
p_cp14a = PowerModels.var(pm, nw, :p_cp14a, i)
74-
75-
# §14a minimum power (per unit)
76-
p_min_14a = gen_cp14a["p_min_14a"]
77-
78-
# Maximum support capacity (matches Python field name "pmax")
79-
p_max_support = gen_cp14a["pmax"]
80-
81-
if i == 3
82-
println(" [DEBUG NW $nw HP $i]")
83-
println(" > p_cp_load: $(round(p_cp_load, digits=6))")
84-
println(" > p_cp14a (Var): $(p_cp14a)") # Zeigt die JuMP-Variable
85-
println(" > p_min_14a: $(round(p_min_14a, digits=6))")
86-
println(" > p_max_support: $(round(p_max_support, digits=6))")
87-
88-
# Hilfswert für die Logik unten berechnen
89-
p_min_net_debug = min(p_cp_load, p_min_14a)
90-
println(" > p_min_net: $(round(p_min_net_debug, digits=6))")
91-
println(" " * "-"^20)
92-
end
93-
94-
# Net load must stay ≥ minimum net load allowed
95-
# The minimum is the LOWER of: current load or §14a limit
96-
# This handles cases where CP draws less than 4.2 kW (e.g., 3 kW due to low charging demand)
97-
# p_cp_load - p_cp14a ≥ min(p_cp_load, p_min_14a)
98-
#
99-
# Special cases:
100-
# - If p_max_support ≈ 0 (CP too small), force virtual gen to zero
101-
# - If CP is off (p_cp_load ≈ 0), no support needed
102-
if p_max_support < 1e-6
103-
# Charging point too small for §14a curtailment, disable virtual generator
104-
JuMP.@constraint(pm.model, p_cp14a == 0.0)
105-
elseif p_cp_load > 1e-6
106-
# Normal case: enforce minimum net load
107-
# Net load cannot go below current load or §14a minimum, whichever is lower
108-
p_min_net = min(p_cp_load, p_min_14a)
109-
JuMP.@constraint(pm.model, p_cp14a <= p_cp_load - p_min_net)
110-
else
111-
# Charging point is off, no support needed
112-
JuMP.@constraint(pm.model, p_cp14a == 0.0)
113-
end
114116
end
115117

116118

‎edisgo/opf/eDisGo_OPF.jl/src/core/constraint_hp_14a.jl‎

Lines changed: 47 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,18 @@ end
3434
Ensures that the net electrical load (heat pump load - virtual generator support)
3535
stays above the §14a minimum power level (typically 4.2 kW = 0.0042 MW).
3636
37+
Uses the HP optimization variable `php` (not the fixed parameter) so the constraint
38+
correctly tracks the actual HP electrical draw after heat storage optimization.
39+
40+
Big-M formulation:
41+
php - p_hp14a >= p_min_14a - M * (1 - z_hp14a)
42+
When z=0 (inactive): constraint relaxed (always satisfied)
43+
When z=1 (active): php - p_hp14a >= p_min_14a (net load >= 4.2 kW)
44+
45+
If php < p_min_14a at a timestep, z is forced to 0 (no curtailment possible),
46+
which in turn forces p_hp14a = 0 via binary coupling. This correctly handles
47+
cases where the HP draws less than 4.2 kW.
48+
3749
# Arguments
3850
- `pm::AbstractBFModelEdisgo`: PowerModels model
3951
- `i::Int`: Virtual generator index
@@ -42,55 +54,54 @@ stays above the §14a minimum power level (typically 4.2 kW = 0.0042 MW).
4254
function constraint_hp_14a_min_net_load(pm::AbstractBFModelEdisgo, i::Int, nw::Int=nw_id_default)
4355
gen_hp14a = PowerModels.ref(pm, nw, :gen_hp_14a, i)
4456
hp_idx = gen_hp14a["hp_index"]
45-
hp = PowerModels.ref(pm, nw, :heatpumps, hp_idx)
46-
47-
# Electrical power demand of heat pump (thermal demand / COP)
48-
p_hp_load = hp["pd"] / hp["cop"]
49-
50-
# Virtual generator support
57+
58+
# Get the actual HP electrical power VARIABLE (not the fixed parameter)
59+
php = PowerModels.var(pm, nw, :php, hp_idx)
60+
61+
# Virtual generator support variable
5162
p_hp14a = PowerModels.var(pm, nw, :p_hp14a, i)
52-
63+
5364
# §14a minimum power (per unit)
5465
p_min_14a = gen_hp14a["p_min_14a"]
55-
56-
# Maximum support capacity (matches Python field name "pmax")
66+
67+
# Maximum support capacity
5768
p_max_support = gen_hp14a["pmax"]
5869

5970
# --- DEBUG PRINT START ---
60-
# Ändere die '12' in die ID, die du genauer prüfen möchtest
61-
if i == 12
62-
println(" [DEBUG NW $nw HP $i]")
63-
println(" > p_hp_load: $(round(p_hp_load, digits=6))")
64-
println(" > p_hp14a (Var): $(p_hp14a)") # Zeigt die JuMP-Variable
65-
println(" > p_min_14a: $(round(p_min_14a, digits=6))")
66-
println(" > p_max_support: $(round(p_max_support, digits=6))")
71+
if i == 12 # ID bei Bedarf anpassen
72+
println("\n[DEBUG §14a Min Net Load | NW $nw | HP $i]")
73+
println(" > hp_idx: $hp_idx")
74+
println(" > php (Var): $(php)")
75+
println(" > p_hp14a (Var): $(p_hp14a)")
76+
println(" > p_min_14a: $p_min_14a")
77+
println(" > p_max_support: $p_max_support")
78+
println(keys(gen_hp14a))
79+
#println(" > Aktuelle Last (pd p.u.): ", gen_hp14a["pd"])
80+
#println(" > Last: $(gen_hp14a["pd"])")
6781

68-
# Hilfswert für die Logik unten berechnen
69-
p_min_net_debug = min(p_hp_load, p_min_14a)
70-
println(" > p_min_net: $(round(p_min_net_debug, digits=6))")
71-
println(" " * "-"^20)
82+
if p_max_support >= 1e-6
83+
z_hp14a = PowerModels.var(pm, nw, :z_hp14a, i)
84+
M = p_max_support + p_min_14a
85+
println(" > z_hp14a (Var): $(z_hp14a)")
86+
println(" > Big-M: $M")
87+
# Zeigt die mathematische Formel der Constraint vor der Lösung
88+
println(" > Constraint: php - p_hp14a >= $p_min_14a - $M * (1 - z_hp14a)")
89+
else
90+
println(" > Status: Small HP (Constraint: p_hp14a == 0)")
91+
end
92+
println("-"^40)
7293
end
7394
# --- DEBUG PRINT ENDE ---
74-
75-
# Net load must stay ≥ minimum net load allowed
76-
# The minimum is the LOWER of: current load or §14a limit
77-
# This handles cases where HP draws less than 4.2 kW (e.g., 3 kW due to low thermal demand)
78-
# p_hp_load - p_hp14a ≥ min(p_hp_load, p_min_14a)
79-
#
80-
# Special cases:
81-
# - If p_max_support ≈ 0 (HP too small), force virtual gen to zero
82-
# - If HP is off (p_hp_load ≈ 0), no support needed
95+
8396
if p_max_support < 1e-6
8497
# Heat pump too small for §14a curtailment, disable virtual generator
8598
JuMP.@constraint(pm.model, p_hp14a == 0.0)
86-
elseif p_hp_load > 1e-6
87-
# Normal case: enforce minimum net load
88-
# Net load cannot go below current load or §14a minimum, whichever is lower
89-
p_min_net = min(p_hp_load, p_min_14a)
90-
JuMP.@constraint(pm.model, p_hp_load - p_hp14a >= p_min_net)
9199
else
92-
# Heat pump is off, no support needed
93-
JuMP.@constraint(pm.model, p_hp14a == 0.0)
100+
# Big-M formulation: when z=1 (curtailment active), enforce min net load
101+
# when z=0 (curtailment inactive), constraint is relaxed
102+
z_hp14a = PowerModels.var(pm, nw, :z_hp14a, i)
103+
M = p_max_support + p_min_14a
104+
JuMP.@constraint(pm.model, php - p_hp14a >= p_min_14a - M * (1 - z_hp14a))
94105
end
95106
end
96107

‎edisgo/opf/eDisGo_OPF.jl/src/core/objective.jl‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,43 @@ function objective_min_losses_slacks_OG(pm::AbstractBFModelEdisgo)
109109
)
110110
end
111111

112+
# OPF Version 5: Minimize line losses, use ONLY §14a curtailment as flexibility
113+
# Feasibility slacks exist but are penalized at 1e8 to ensure model remains feasible
114+
function objective_min_losses_14a_only(pm::AbstractBFModelEdisgo)
115+
nws = PowerModels.nw_ids(pm)
116+
ccm = Dict(n => PowerModels.var(pm, n, :ccm) for n in nws)
117+
r = Dict(n => Dict(i => get(branch, "br_r", 1.0) for (i,branch) in PowerModels.ref(pm, n, :branch)) for n in nws)
118+
pgc = Dict(n => PowerModels.var(pm, n, :pgc) for n in nws)
119+
pgens = Dict(n => PowerModels.var(pm, n, :pgens) for n in nws)
120+
pds = Dict(n => PowerModels.var(pm, n, :pds) for n in nws)
121+
pcps = Dict(n => PowerModels.var(pm, n, :pcps) for n in nws)
122+
phps = Dict(n => PowerModels.var(pm, n, :phps) for n in nws)
123+
phps2 = Dict(n => PowerModels.var(pm, n, :phps2) for n in nws)
124+
phss = Dict(n => PowerModels.var(pm, n, :phss) for n in nws)
125+
126+
# §14a virtual generators for HPs and CPs
127+
p_hp14a = Dict(n => get(PowerModels.var(pm, n), :p_hp14a, Dict()) for n in nws)
128+
p_cp14a = Dict(n => get(PowerModels.var(pm, n), :p_cp14a, Dict()) for n in nws)
129+
130+
factor_14a = 0.5 # Weight for §14a curtailment
131+
factor_feasibility = 1e8 # Extreme penalty - slacks should be zero in normal operation
132+
133+
return JuMP.@objective(pm.model, Min,
134+
# Primary: minimize line losses
135+
0.4 * sum(sum(ccm[n][b] * r[n][b] for (b,i,j) in PowerModels.ref(pm, n, :arcs_from)) for n in nws)
136+
# Secondary: minimize §14a curtailment usage
137+
+ factor_14a * sum(sum(p_hp14a[n][i] for i in keys(p_hp14a[n])) for n in nws) # §14a HP curtailment
138+
+ factor_14a * sum(sum(p_cp14a[n][i] for i in keys(p_cp14a[n])) for n in nws) # §14a CP curtailment
139+
# Feasibility slacks (extreme penalty - should never be used if §14a is sufficient)
140+
+ factor_feasibility * sum(sum(pgc[n][i] for i in keys(PowerModels.ref(pm, 1, :gen_nd))) for n in nws)
141+
+ factor_feasibility * sum(sum(pgens[n][i] for i in keys(PowerModels.ref(pm, 1, :gen))) for n in nws)
142+
+ factor_feasibility * sum(sum(pds[n][i] for i in keys(PowerModels.ref(pm, 1, :load))) for n in nws)
143+
+ factor_feasibility * sum(sum(pcps[n][i] for i in keys(PowerModels.ref(pm, 1, :electromobility))) for n in nws)
144+
+ factor_feasibility * sum(sum(phps[n][i] for i in keys(PowerModels.ref(pm, 1, :heatpumps))) for n in nws)
145+
+ 1e4 * sum(sum(phss[n][i] + phps2[n][i] for i in keys(PowerModels.ref(pm, 1, :heatpumps))) for n in nws)
146+
)
147+
end
148+
112149
# OPF Version 3 (alternative): Minimize line losses, maximal line loading and HV slacks (with overlying grid)
113150
function objective_min_line_loading_max_OG(pm::AbstractBFModelEdisgo)
114151
nws = PowerModels.nw_ids(pm)

‎edisgo/opf/eDisGo_OPF.jl/src/core/variables.jl‎

Lines changed: 66 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -596,4 +596,69 @@ function variable_slack_HV_requirements_imaginary(pm::AbstractPowerModel; nw::In
596596

597597
end
598598

599-
""
599+
600+
### Version 5: Fixed variables (no flexibility except §14a)
601+
602+
"""
603+
variable_heat_pump_power_fixed(pm; nw, report)
604+
605+
For OPF version 5: Fix heat pump power to pd/cop (current demand).
606+
No optimization flexibility - only §14a can reduce load.
607+
"""
608+
function variable_heat_pump_power_fixed(pm::AbstractPowerModel; nw::Int=nw_id_default, report::Bool=true)
609+
# Create php variable fixed to pd/cop for each heat pump
610+
php = PowerModels.var(pm, nw)[:php] = JuMP.@variable(pm.model,
611+
[i in PowerModels.ids(pm, nw, :heatpumps)], base_name="$(nw)_php"
612+
)
613+
614+
# Fix to current demand: php = pd/cop
615+
for (i, hp) in PowerModels.ref(pm, nw, :heatpumps)
616+
p_demand = hp["pd"] / hp["cop"]
617+
JuMP.fix(php[i], p_demand; force=true)
618+
end
619+
620+
report && PowerModels.sol_component_value(pm, nw, :heatpumps, :php, PowerModels.ids(pm, nw, :heatpumps), php)
621+
end
622+
623+
"""
624+
variable_cp_power_fixed(pm; nw, report)
625+
626+
For OPF version 5: Fix charging point power to pd (current demand from load).
627+
No optimization flexibility - only §14a can reduce load.
628+
Also creates cpe (energy) variable fixed at mid-range.
629+
"""
630+
function variable_cp_power_fixed(pm::AbstractPowerModel; nw::Int=nw_id_default, report::Bool=true)
631+
# Check if electromobility dict exists and has entries
632+
if !haskey(PowerModels.ref(pm, nw), :electromobility) || isempty(PowerModels.ref(pm, nw, :electromobility))
633+
# Create empty variables to avoid errors in power balance
634+
PowerModels.var(pm, nw)[:pcp] = Dict{Int, JuMP.VariableRef}()
635+
PowerModels.var(pm, nw)[:cpe] = Dict{Int, JuMP.VariableRef}()
636+
return
637+
end
638+
639+
# Create pcp variable fixed to pd for each charging point
640+
pcp = PowerModels.var(pm, nw)[:pcp] = JuMP.@variable(pm.model,
641+
[i in PowerModels.ids(pm, nw, :electromobility)], base_name="$(nw)_pcp"
642+
)
643+
644+
# Fix to current demand: pcp = pd
645+
for (i, cp) in PowerModels.ref(pm, nw, :electromobility)
646+
p_demand = cp["pd"]
647+
JuMP.fix(pcp[i], p_demand; force=true)
648+
end
649+
650+
report && PowerModels.sol_component_value(pm, nw, :electromobility, :pcp, PowerModels.ids(pm, nw, :electromobility), pcp)
651+
652+
# Also need cpe (energy) variable for constraints - fix at midpoint
653+
cpe = PowerModels.var(pm, nw)[:cpe] = JuMP.@variable(pm.model,
654+
[i in PowerModels.ids(pm, nw, :electromobility)], base_name="$(nw)_cpe"
655+
)
656+
657+
for (i, cp) in PowerModels.ref(pm, nw, :electromobility)
658+
e_mid = 0.5 * (cp["e_min"] + cp["e_max"])
659+
JuMP.fix(cpe[i], e_mid; force=true)
660+
end
661+
662+
report && PowerModels.sol_component_value(pm, nw, :electromobility, :cpe, PowerModels.ids(pm, nw, :electromobility), cpe)
663+
end
664+

‎edisgo/opf/eDisGo_OPF.jl/src/form/bf.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -150,7 +150,7 @@ function constraint_power_balance(pm::AbstractBFModelEdisgo, n::Int, i, bus_gens
150150
p_hp14a = get(PowerModels.var(pm, n), :p_hp14a, Dict()) # §14a virtual generators for heat pumps
151151
p_cp14a = get(PowerModels.var(pm, n), :p_cp14a, Dict()) # §14a virtual generators for charging points
152152

153-
if PowerModels.ref(pm, 1, :opf_version) in(2, 4) # Eq. (3.3iii), (3.4iii)
153+
if PowerModels.ref(pm, 1, :opf_version) in(2, 4, 5) # Eq. (3.3iii), (3.4iii)
154154
pgens = get(PowerModels.var(pm, n), :pgens, Dict()); PowerModels._check_var_keys(pgens, bus_gens, "active power slack", "curtailment")
155155
pds = get(PowerModels.var(pm, n), :pds, Dict()); PowerModels._check_var_keys(pds, bus_loads, "active power slack", "load")
156156
pcps = get(PowerModels.var(pm, n), :pcps, Dict()); PowerModels._check_var_keys(pcps, bus_cps, "active power slack", "charging point")
@@ -191,7 +191,7 @@ function constraint_power_balance(pm::AbstractBFModelEdisgo, n::Int, i, bus_gens
191191
+ sum(pgens[g] * bus_gen_d_pf[g] for g in bus_gens)
192192
+ sum(pdsm[dsm] * bus_dsm_pf[dsm] for dsm in bus_dsm)
193193
+ sum((php[hp] - phps[hp]) * bus_hps_pf[hp] for hp in bus_hps)
194-
+ sum((pcp[cp] - pcps[cp]) * bus_cps_pf[cp] for hp in bus_cps)
194+
+ sum((pcp[cp] - pcps[cp]) * bus_cps_pf[cp] for cp in bus_cps)
195195
# §14a generators have pf=1, q=0
196196
)
197197
else # Eq. (3.3ii), (3.4ii)

0 commit comments

Comments
 (0)
⚡