add centralized optimization solution approach

src/demoses_distibuted_optimization/centralized_optimization.py

+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