From 481cc61606373bb19ac1130fdc9065e683c437da Mon Sep 17 00:00:00 2001
From: Christian Doh Dinga <cdohdinga@tudelft.nl>
Date: Tue, 25 Jun 2024 16:36:59 +0200
Subject: [PATCH] add centralized optimization solution approach
---
.../centralized_optimization.py | 109 ++++++++++++++++++
1 file changed, 109 insertions(+)
create mode 100644 src/demoses_distibuted_optimization/centralized_optimization.py
diff --git a/src/demoses_distibuted_optimization/centralized_optimization.py b/src/demoses_distibuted_optimization/centralized_optimization.py
new file mode 100644
index 0000000..44eff32
--- /dev/null
+++ b/src/demoses_distibuted_optimization/centralized_optimization.py
@@ -0,0 +1,109 @@
+import yaml
+import numpy as np
+import pandas as pd
+from typing import Dict
+import pyomo.environ as pyo
+
+
+def create_centralized_optimization_problem(data: Dict, ts: pd.DataFrame) -> pyo.ConcreteModel:
+ """Create optimization problem to solve the equilibrium problem in a centralized manner."""
+ model = pyo.ConcreteModel(name='Optimization-problem')
+
+ # Define sets
+ number_of_timesteps = data["General"]["nTimesteps"]
+ generator_list = data["Generators"].keys()
+ model.time = pyo.Set(initialize=list(range(number_of_timesteps)), name='timesteps')
+ model.generators = pyo.Set(initialize=list(generator_list), name='generators')
+
+ # Generator-specific parameters
+ generator_cost_param_a_value = {gen: data['Generators'][gen]['a'] for gen in data['Generators']}
+ generator_cost_param_b_value = {gen: data['Generators'][gen]['b'] for gen in data['Generators']}
+ generator_capacity_value = {gen: data['Generators'][gen]['C'] for gen in data['Generators']}
+ model.generator_cost_param_a = pyo.Param(model.generators, name='generator_cost_param_a', initialize=generator_cost_param_a_value, mutable=False)
+ model.generator_cost_param_b = pyo.Param(model.generators, name='generator_cost_param_b', initialize=generator_cost_param_b_value, mutable=False)
+ model.generator_capacity = pyo.Param(model.generators, name='generator_capacity', initialize=generator_capacity_value, mutable=False)
+
+ # Consumer-specific parameters
+ demand_timeseries = dict(enumerate(ts.loc[:,'LOAD'].values)) # combined demand from all consumers
+ model.demand_profile = pyo.Param(model.time, name='demand_profile', initialize=demand_timeseries, mutable=False)
+
+ # For each consumer, we want to get their pv generation timeseries, then combined all together
+ consumer_list = data["Consumers"].keys()
+ total_consumers = data['General']['totConsumers']
+ all_pv_profile = np.zeros(data["General"]["nTimesteps"])
+ for con in consumer_list:
+ pv_availability_factor = ts.loc[:, data['Consumers'][con]['PV_AF']].values
+ pv_capacity = data['Consumers'][con]['PV_cap']
+ share_of_agent = data['Consumers'][con]['Share']
+ pv_profile = total_consumers * share_of_agent * pv_capacity * pv_availability_factor
+ try:
+ pv_profile.reshape(1, 24)
+ all_pv_profile = np.add(all_pv_profile, pv_profile)
+ except ValueError:
+ print("Error: trying to add np.zeros(24) array with a different shape array")
+ return 0
+
+ model.pv_profile = pyo.Param(model.time, name='demand_profile', initialize=dict(enumerate(all_pv_profile)), mutable=False)
+
+ # Decision variables for generators
+ model.var_g = pyo.Var(model.generators, model.time, name='generation', domain=pyo.NonNegativeReals, initialize=0)
+
+ # Define generator capacity limit constraint
+ def available_capacity_value(model, gen, t):
+ if 'AF' in data['Generators'][gen].keys():
+ return model.generator_capacity[gen] * ts.loc[t, data['Generators'][gen]["AF"]]
+ else:
+ return model.generator_capacity[gen]
+ model.available_capacity = pyo.Param(model.generators, model.time, name='available_capacity', initialize=available_capacity_value, mutable=False)
+
+ # Define objective function (OPEX are parametrized as a/2 * g**2 + b * g)
+ def generation_cost_rule(model):
+ return sum(
+ model.generator_cost_param_a[gen]/2 * model.var_g[gen, t]**2
+ + model.generator_cost_param_b[gen] * model.var_g[gen, t]
+ for gen in model.generators for t in model.time
+ )
+ model.objective_function = pyo.Objective(rule=generation_cost_rule, sense=pyo.minimize)
+
+ # Define generator capacity limit constraint
+ def capacity_limit_rule(model, gen, t):
+ return model.var_g[gen, t] <= model.available_capacity[gen, t]
+ model.capacity_limit = pyo.Constraint(model.generators, model.time, rule=capacity_limit_rule)
+
+ # Define energy balance constraint
+ def energy_balance_rule(model, t):
+ return (sum(model.var_g[gen, t] for gen in model.generators) >= model.demand_profile[t] - model.pv_profile[t])
+ model.energy_balance = pyo.Constraint(model.time, rule=energy_balance_rule)
+
+ return model
+
+
+def solve_model():
+ """Solve model and return results."""
+ def read_config(config_file):
+ with open(config_file, 'r') as file:
+ config = yaml.safe_load(file)
+ return config
+ data = read_config('config.yaml')
+ ts = pd.read_csv('timeseries.csv', delimiter=';')
+ model = create_centralized_optimization_problem(data, ts)
+ solver = pyo.SolverFactory('gurobi')
+ model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
+ solver.solve(model)
+ print('############## Objective function value ###########')
+ model.objective_function.display()
+ print()
+ print('############## Generator output ###########')
+ model.var_g.display()
+ print()
+ dual_values = []
+ for t in model.time:
+ dual_values.append(model.dual[model.energy_balance[t]])
+ market_prices = pd.DataFrame(data={"market-prices": dual_values}, index=model.time)
+ market_prices.index.name = 'timesteps'
+ print('############## Market prices ###########')
+ print(market_prices)
+
+
+if __name__ == "__main__":
+ solve_model()
\ No newline at end of file
--
GitLab