import time
from typing import Optional, Union
import gurobipy as gp
from gurobipy import GRB
from apem.unit_based_model.allocation.allocation import Allocation
from apem.unit_based_model.solver_configuration import SolverConfiguration
from apem.unit_based_model.error import Error
from apem.unit_based_model.data.parsing.scenario import Scenario
from apem.unit_based_model.pricing.algorithms.fbmc_support import (
add_fbmc_price_coupling_constraints,
get_fbmc_pricing_data,
)
from apem.unit_based_model.pricing.algorithms.pricing_algorithm import PricingAlgorithm
from apem.unit_based_model.pricing.analysis.pricing import MWPS, Pricing
from apem.unit_based_model.pricing.analysis.write_prices import write_prices, write_prices_failure
from apem.unit_based_model.utils.extraction import preprocess_as_dict
[docs]
class MinMWP(PricingAlgorithm):
"""Implementation of Minimum Make-Whole Payments Pricing.
"""
[docs]
def compute_prices(self, allocation: Allocation, scenario: Scenario, configuration: SolverConfiguration,
file_prices: Optional[str] = None,
fixed_prices: Optional[Pricing] = None) -> Union[Pricing, Error]:
"""
Formulates and solves a Min-MWP problem similar to the one from
"Pricing Optimal Outcomes in Coupled and Non-Convex Markets: Theory and Applications to Electricity
Markets" (Appendix E, https://arxiv.org/abs/2209.07386). The method can also be used to compute the
MWPs for an allocation-prices pair.
:param allocation: allocation for which supporting prices are computed
:param scenario: scenario for which prices are computed
:param configuration: configuration object containing the parameters for the pricing algorithm
:param file_prices: name of the file in which results are written
:param fixed_prices: prices for which MWPs should be computed
:return: Pricing object if prices could be computed or Error object otherwise
"""
if allocation.status != 1:
if file_prices:
write_prices_failure(file_prices, str(self), -1)
print(f'{self} pricing error with code -1')
return Error(-1) # Allocation computation failed
start = time.time()
model = gp.Model('Min-MWP-Pricing')
# apply Gurobi configuration parameters
configuration.apply_to_model(model)
df_buyers = scenario.df_buyers
df_sellers = scenario.df_sellers
network = scenario.network
periods = scenario.periods
blocks_buyers = scenario.blocks_buyers
blocks_sellers = scenario.blocks_sellers
r_star = scenario.r_star
nodes = network.nodes
buyers = df_buyers['buyer'].unique().tolist()
sellers = df_sellers['seller'].unique().tolist()
# precompute dictionaries for fast access
buyer_val_dict, seller_cost_dict = {}, {}
buyer_node_dict = preprocess_as_dict(df_buyers, ['buyer', 'period'], 'node')
seller_no_load_cost_dict = preprocess_as_dict(df_sellers, ['seller', 'period'], 'no_load_cost')
seller_node_dict = preprocess_as_dict(df_sellers, ['seller', 'period'], 'node')
for block in blocks_buyers:
buyer_val_dict[block] = preprocess_as_dict(df_buyers, ['buyer', 'period'], 'val', block)
for block in blocks_sellers:
seller_cost_dict[block] = preprocess_as_dict(df_sellers, ['seller', 'period'], 'cost', block)
y_st = allocation.SellersAllocation.y_st
x_btl = allocation.BuyersAllocation.x_btl
y_stl = allocation.SellersAllocation.y_stl
f_vwt = allocation.TransmissionNetworkAllocation.f_vwt
f_vwkt = getattr(allocation.TransmissionNetworkAllocation, "f_vwkt", None)
u_st = allocation.SellersAllocation.u_st
fbmc_data = get_fbmc_pricing_data(allocation)
use_fbmc_network = fbmc_data is not None
p_vt = model.addVars(nodes, periods, lb=-GRB.INFINITY, ub=GRB.INFINITY, name='p_vt')
if use_fbmc_network:
fb_constraint_ids = fbmc_data["constraint_ids"]
flow_lt = fbmc_data["flow"]
gamma_lt = model.addVars(fb_constraint_ids, periods, lb=-GRB.INFINITY, ub=GRB.INFINITY, name='gamma_l_t')
else:
is_multigraph = network.is_multigraph()
undirected_edges = list(network.edges(keys=True, data=True)) if is_multigraph else \
[(u, v, None, data) for u, v, data in network.edges(data=True)]
if is_multigraph:
if not f_vwkt:
if file_prices:
write_prices_failure(file_prices, str(self), -2)
print(f'{self} pricing error with code -2: missing multigraph per-edge flows')
return Error(-2)
missing_flow_key = next(
(
(u, v, k, t)
for (u, v, k, _) in undirected_edges
for t in periods
if (u, v, k, t) not in f_vwkt
),
None,
)
if missing_flow_key is not None:
if file_prices:
write_prices_failure(file_prices, str(self), -2)
print(f'{self} pricing error with code -2: missing multigraph flow key {missing_flow_key}')
return Error(-2)
directed_edges = []
for idx, (u, v, k, data) in enumerate(undirected_edges):
directed_edges.append((idx, u, v, k, data))
directed_edges.append((idx, v, u, k, data))
flow_et = {}
for idx, (u, v, k, data) in enumerate(undirected_edges):
for t in periods:
if is_multigraph:
base = f_vwkt[(u, v, k, t)]
else:
base = f_vwt[(u, v, t)]
flow_et[(idx, u, v, t)] = base
flow_et[(idx, v, u, t)] = -base
gamma_et = model.addVars([(e, v, w, t) for (e, v, w, _, _) in directed_edges for t in periods],
lb=-GRB.INFINITY, ub=GRB.INFINITY, name='gamma_e_t')
if fixed_prices:
model.addConstrs(p_vt[v, t] == fixed_prices.node_prices[v, t] for v in nodes for t in periods)
r_t = model.addVars(periods, lb=-GRB.INFINITY, ub=GRB.INFINITY, name='r_t')
lambda_b = model.addVars(buyers, ub=GRB.INFINITY, name='lambda_b')
lambda_s = model.addVars(sellers, ub=GRB.INFINITY, name='lambda_s')
if use_fbmc_network:
lambda_lt = model.addVars(fb_constraint_ids, periods, ub=GRB.INFINITY, name='lambda_l_t')
else:
lambda_et = model.addVars([(e, v, w, t) for (e, v, w, _, _) in directed_edges for t in periods],
ub=GRB.INFINITY, name='lambda_e_t')
model.update()
model.setObjective(
gp.quicksum(lambda_b[b] for b in buyers) +
gp.quicksum(lambda_s[s] for s in sellers) +
(
gp.quicksum(lambda_lt[line_id, t] for line_id in fb_constraint_ids for t in periods)
if use_fbmc_network
else gp.quicksum(lambda_et[e, v, w, t] for (e, v, w, _, _) in directed_edges for t in periods)
),
GRB.MINIMIZE
)
# 1
model.addConstrs(
- gp.quicksum(
buyer_val_dict[lb][b, t] * x_btl[b, t, lb]
for t in periods
for lb in blocks_buyers
)
+ gp.quicksum(
p_vt[buyer_node_dict[b, t], t] * x_btl[b, t, lb]
for t in periods
for lb in blocks_buyers
)
- lambda_b[b]
<= 0
for b in buyers
)
# 2
model.addConstrs(
-gp.quicksum(
p_vt[seller_node_dict[s, t], t] * y_st[s, t]
for t in periods
)
+ gp.quicksum(
seller_cost_dict[ls][s, t] * y_stl[s, t, ls]
for t in periods
for ls in blocks_sellers
)
+ gp.quicksum(
seller_no_load_cost_dict[s, t] * u_st[s, t]
for t in periods
)
- lambda_s[s]
<= 0
for s in sellers
)
if use_fbmc_network:
for line_id in fb_constraint_ids:
for t in periods:
model.addConstr(-gamma_lt[line_id, t] * flow_lt[(line_id, t)] - lambda_lt[line_id, t] <= 0)
add_fbmc_price_coupling_constraints(model, p_vt, r_t, gamma_lt, nodes, periods, fbmc_data)
else:
for (e, v, w, k, data) in directed_edges:
for t in periods:
model.addConstr(-gamma_et[e, v, w, t] * flow_et[(e, v, w, t)] - lambda_et[e, v, w, t] <= 0)
for v in nodes:
if v == r_star:
continue
for t in periods:
inflow = gp.quicksum(
data['B'] * (p_vt[w, t] + gamma_et[e, w, v, t])
for (e, w, v2, k, data) in directed_edges if v2 == v
)
outflow = gp.quicksum(
data['B'] * (p_vt[v, t] + gamma_et[e, v, w, t])
for (e, v2, w, k, data) in directed_edges if v2 == v
)
model.addConstr(inflow - outflow == 0)
for t in periods:
inflow = gp.quicksum(
data['B'] * (p_vt[w, t] + gamma_et[e, w, r_star, t])
for (e, w, v2, k, data) in directed_edges if v2 == r_star
)
outflow = gp.quicksum(
data['B'] * (p_vt[r_star, t] + gamma_et[e, r_star, w, t])
for (e, v2, w, k, data) in directed_edges if v2 == r_star
)
model.addConstr(r_t[t] + inflow - outflow == 0)
model.optimize()
end = time.time()
runtime = end - start
num_vars = model.NumVars
num_constrs = model.NumConstrs
status = model.getAttr('Status')
if status == GRB.OPTIMAL:
total_mwps = round(model.getObjective().getValue(), 2)
mwps_buyers = round(sum(lambda_b[b].X for b in buyers), 2)
mwps_sellers = round(sum(lambda_s[s].X for s in sellers), 2)
if use_fbmc_network:
mwps_network = round(sum(lambda_lt[line_id, t].X for line_id in fb_constraint_ids for t in periods), 2)
else:
mwps_network = round(
sum(lambda_et[e, v, w, t].X for (e, v, w, _, _) in directed_edges for t in periods), 2)
mwps_per_buyer = {b: round(lambda_b[b].X, 2) for b in buyers}
mwps_per_seller = {s: round(lambda_s[s].X, 2) for s in sellers}
if use_fbmc_network:
mwps_per_line = {
line_id: round(sum(lambda_lt[line_id, t].X for t in periods), 2)
for line_id in fb_constraint_ids
}
else:
mwps_per_line = {}
for (e, v, w, _, _) in directed_edges:
mwps_per_line[(v, w, e)] = round(sum(lambda_et[e, v, w, t].X for t in periods), 2)
p_vt = {(v, t): p_vt[v, t].X for v in nodes for t in periods}
if use_fbmc_network:
gamma_vwt = {(line_id, t): gamma_lt[line_id, t].X for line_id in fb_constraint_ids for t in periods}
gamma_vwkt = {}
else:
gamma_vwt = {}
gamma_vwkt = {}
for (e, v, w, k, _) in directed_edges:
for t in periods:
gamma_val = gamma_et[e, v, w, t].X
gamma_vwt[(v, w, t)] = gamma_vwt.get((v, w, t), 0) + gamma_val
gamma_vwkt[(v, w, k, t)] = gamma_val
pricing = Pricing(p_vt, gamma_vwt, str(self), runtime, num_vars, num_constrs,
mwps=MWPS(total_mwps, mwps_buyers, mwps_sellers, mwps_network,
mwps_per_buyer, mwps_per_seller, mwps_per_line),
line_congestion_prices_per_edge=gamma_vwkt)
if file_prices:
write_prices(file_prices, pricing, scenario)
return pricing
if file_prices:
write_prices_failure(file_prices, str(self), status)
print(f'{self} pricing error with code {status}')
return Error(status)
def __str__(self):
return 'Min_MWP'