Commit 34e3014e authored by Sebastian Schwarz's avatar Sebastian Schwarz
Browse files

MPI implementation for Exchange ADMM optimization algorithm.

parent de4a61cf
Pipeline #512406 failed with stages
in 4 minutes and 23 seconds
......@@ -4,6 +4,7 @@ matplotlib==3.3.4
pyomo==5.7.1
Shapely==1.7.1
pycity_base==0.3.2
mpi4py==3.0.3
pylint
sphinx
numpydoc
......
......@@ -32,9 +32,10 @@ from pycity_scheduling.classes import *
# In this example, the power schedule for a complex city district scenario is determined. The scenario is built upon the
# district setup as defined in example 'example_12_district_generator.py', but it contains more than 100 buildings and
# district setup as defined in example 'example_13_district_generator.py', but it contains more than 100 buildings and
# is hence considered more complex.
def main(do_plot=False):
print("\n\n------ Example 15: Scheduling Complex City District ------\n\n")
......
"""
The pycity_scheduling framework
Copyright (C) 2021,
Institute for Automation of Complex Power Systems (ACS),
E.ON Energy Research Center (E.ON ERC),
RWTH Aachen University
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
import numpy as np
from pycity_scheduling.classes import *
from pycity_scheduling.algorithms import *
from pycity_scheduling.util import mpi_interface
# This is a very simple power scheduling example using the distributed exchange ADMM MPI algorithm.
def main(do_plot=False):
mpi = mpi_interface.MPI_Interface()
if mpi.mpi_rank == 0:
print("\n\n------ Example 22: Algorithm Exchange-ADMM-MPI ------\n\n")
# Define timer, price, weather, and environment objects:
t = Timer(op_horizon=2, step_size=3600)
p = Prices(timer=t)
w = Weather(timer=t)
e = Environment(timer=t, weather=w, prices=p)
# City district with district operator objective "peak-shaving":
cd = CityDistrict(environment=e, objective='peak-shaving')
# Schedule two sample buildings. The buildings' objectives are defined as "price".
# Building no. one comes with fixed load, space heating, electrical heater, pv unit, thermal energy storage, and
# electrical energy storage:
bd1 = Building(environment=e, objective='price')
cd.addEntity(entity=bd1, position=[0, 0])
bes = BuildingEnergySystem(environment=e)
bd1.addEntity(bes)
ths = ThermalHeatingStorage(environment=e, e_th_max=40, soc_init=0.5)
bes.addDevice(ths)
eh = ElectricalHeater(environment=e, p_th_nom=10)
bes.addDevice(eh)
ap = Apartment(environment=e)
bd1.addEntity(ap)
load = np.array([10.0, 10.0])
fi = FixedLoad(e, method=0, demand=load)
ap.addEntity(fi)
sh = SpaceHeating(environment=e, method=0, loadcurve=load)
ap.addEntity(sh)
pv = Photovoltaic(environment=e, method=1, peak_power=4.6)
bes.addDevice(pv)
bat = Battery(environment=e, e_el_max=4.8, p_el_max_charge=3.6, p_el_max_discharge=3.6)
bes.addDevice(bat)
# Building no. two comes with deferrable load, curtailable load, space heating, chp unit, thermal energy storage and
# an electrical vehicle:
bd2 = Building(environment=e, objective='price')
cd.addEntity(entity=bd2, position=[0, 0])
bes = BuildingEnergySystem(environment=e)
bd2.addEntity(bes)
ths = ThermalHeatingStorage(environment=e, e_th_max=35, soc_init=0.5)
bes.addDevice(ths)
chp = CombinedHeatPower(environment=e, p_th_nom=20.0)
bes.addDevice(chp)
ap = Apartment(environment=e)
bd2.addEntity(ap)
load = np.array([20.0, 20.0])
dl = DeferrableLoad(environment=e, p_el_nom=2.0, e_consumption=2.0, load_time=[1, 1])
ap.addEntity(dl)
cl = CurtailableLoad(environment=e, p_el_nom=1.6, max_curtailment=0.8)
ap.addEntity(cl)
sh = SpaceHeating(environment=e, method=0, loadcurve=load)
ap.addEntity(sh)
ev = ElectricalVehicle(environment=e, e_el_max=37.0, p_el_max_charge=22.0, soc_init=0.65, charging_time=[0, 1])
ap.addEntity(ev)
# Perform the scheduling:
opt = ExchangeADMMMPI(city_district=cd, rho=2.0, eps_primal=0.001, eps_dual=0.01)
results = opt.solve()
cd.copy_schedule("admm-mpi")
# Print some ADMM results:
if mpi.mpi_rank == 0:
print("ADMM - Number of iterations:")
print(results["iterations"][-1])
print("ADMM - Norm vector 'r' over iterations:")
print(results["r_norms"])
print("ADMM - Norm vector 's' over iterations:")
print(results["s_norms"])
print("")
# Print the building's schedules:
print("Schedule building no. one:")
print(list(bd1.p_el_schedule))
print("Schedule building no. two:")
print(list(bd2.p_el_schedule))
print("Schedule of the city district:")
print(list(cd.p_el_schedule))
return
# Conclusions:
# If the distributed MPI exchange ADMM optimization algorithm is applied, the two buildings are scheduled in parallel
# such that both the local and system level objectives are satisfied. Local flexibility is used to achieve the system
# level objective. The scheduling results are close to the ones of the central algorithm, which demonstrates the
# correctness of the distributed algorithm.
if __name__ == '__main__':
# Run example:
main(do_plot=True)
"""
The pycity_scheduling framework
Copyright (C) 2021,
Institute for Automation of Complex Power Systems (ACS),
E.ON Energy Research Center (E.ON ERC),
RWTH Aachen University
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
import numpy as np
import matplotlib.pyplot as plt
import pycity_scheduling.util.factory as factory
import pycity_scheduling.util.debug as debug
from pycity_scheduling.algorithms import *
from pycity_scheduling.classes import *
from pycity_scheduling.util import mpi_interface
# In this example, the power schedule for a city district scenario is determined using the parallel MPI Exchange
# ADMM. The scenario is built upon the district setup as defined in example 'example_13_district_generator.py'.
def main(do_plot=False):
mpi = mpi_interface.MPI_Interface()
if mpi.mpi_rank == 0:
print("\n\n------ Example 23: Scheduling City District MPI ------\n\n")
# First, create an environment using the factory's "generate_standard_environment" method. The environment
# automatically encapsulates time, weather, and price data/information.
env = factory.generate_standard_environment(initial_date=(2018, 12, 6), step_size=900, op_horizon=96)
# Create 20 single-family houses:
num_sfh = 10
# 50% SFH.2002, 30% SFH.2010, 20% SFH.2016 (based on TABULA):
sfh_distribution = {
'SFH.2002': 0.5,
'SFH.2010': 0.3,
'SFH.2016': 0.2,
}
# 50% of the single-family houses are equipped with heat pump, 10% with boiler, and 40% with electrical heater:
sfh_heating_distribution = {
'HP': 0.5,
'BL': 0.1,
'EH': 0.4,
}
# All single-family houses are equipped with a fixed load, 20% have a deferrable load, and 30% have an electric
# vehicle. Moreover, 50% of all single-family houses have a battery unit and 80% have a rooftop photovoltaic unit
# installation.
# The values are rounded in case they cannot be perfectly matched to the given number of buildings.
sfh_device_probs = {
'FL': 1,
'DL': 0.2,
'EV': 0.3,
'BAT': 0.5,
'PV': 0.8,
}
# Create 0 multi-family houses (number of apartments according to TABULA):
num_mfh = 0
# 60% MFH.2002, 20% SFH.2010, 20% SFH.2016 (based on TABULA):
mfh_distribution = {
'MFH.2002': 0.6,
'MFH.2010': 0.2,
'MFH.2016': 0.2,
}
# 40% of the multi-family houses are equipped with heat pump, 20% with boiler, and 40% with electrical heater:
mfh_heating_distribution = {
'HP': 0.4,
'BL': 0.2,
'EH': 0.4,
}
# All apartments inside a multi-family houses are equipped with a fixed load, 20% have a deferrable load, and 20%
# have an electric vehicle. Moreover, 40% of all multi-family houses have a battery unit and 80% have a rooftop
# photovoltaic unit installation.
# The values are rounded in case they cannot be perfectly matched to the given number of buildings.
mfh_device_probs = {
'FL': 1,
'DL': 0.2,
'EV': 0.2,
'BAT': 0.4,
'PV': 0.8,
}
# Finally, create the desired city district using the factory's "generate_tabula_district" method. The district's/
# district operator's objective is defined as "peak-shaving" and the buildings' objectives are defined as "price".
district = factory.generate_tabula_district(env, num_sfh, num_mfh,
sfh_distribution,
sfh_heating_distribution,
sfh_device_probs,
mfh_distribution,
mfh_heating_distribution,
mfh_device_probs,
district_objective='price',
building_objective='price'
)
# Hierarchically print the district and all buildings/assets:
if mpi.mpi_rank == 0:
debug.print_district(district, 1)
# Perform the parallel city district scheduling using the exchange ADMM MPI optimization algorithm:
opt = ExchangeADMMMPI(city_district=district, rho=0.1, eps_primal=10.0, eps_dual=1.0)
results = opt.solve()
district.copy_schedule("district_schedule")
# Plot the scheduling results:
if mpi.mpi_rank == 0:
print(list(district.p_el_schedule))
plt.plot(district.p_el_schedule)
plt.ylabel("City District Power [kW]")
plt.title("Complex City District Scenario - Schedule")
plt.grid()
if do_plot:
plt.show()
return
# Conclusions:
# The parallel and distributed power scheduling for a city district scenario can be done easily using pycity_scheduling
# with MPI support and the Python module mpi4py installed.
if __name__ == '__main__':
# Run example:
main(do_plot=True)
......@@ -25,6 +25,7 @@ OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
from .stand_alone_optimization_algorithm import StandAlone
from .local_optimization_algorithm import LocalOptimization
from .exchange_admm_algorithm import ExchangeADMM
from .exchange_admm_algorithm_mpi import ExchangeADMMMPI
from .central_optimization_algorithm import CentralOptimization
from .dual_decomposition_algorithm import DualDecomposition
......@@ -33,6 +34,7 @@ __all__ = [
'StandAlone',
'LocalOptimization',
'ExchangeADMM',
'ExchangeADMMMPI',
'CentralOptimization',
'DualDecomposition',
'algorithm',
......@@ -44,6 +46,7 @@ algorithms = {
'stand-alone': StandAlone,
'local': LocalOptimization,
'exchange-admm': ExchangeADMM,
'exchange-admm-mpi': ExchangeADMMMPI,
'central': CentralOptimization,
'dual-decomposition': DualDecomposition,
}
"""
The pycity_scheduling framework
Copyright (C) 2020,
Institute for Automation of Complex Power Systems (ACS),
E.ON Energy Research Center (E.ON ERC),
RWTH Aachen University
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
import numpy as np
import pyomo.environ as pyomo
from pycity_scheduling.classes import (CityDistrict, Building, Photovoltaic, WindEnergyConverter)
from pycity_scheduling.util import extract_pyomo_values, mpi_interface
from pycity_scheduling.algorithms.algorithm import IterationAlgorithm, DistributedAlgorithm, SolverNode
from pycity_scheduling.solvers import DEFAULT_SOLVER, DEFAULT_SOLVER_OPTIONS
class ExchangeADMMMPI(IterationAlgorithm, DistributedAlgorithm):
"""Implementation of the Exchange ADMM Algorithm.
Uses the Exchange ADMM algorithm described in [1].
Parameters
----------
city_district : CityDistrict
solver : str, optional
Solver to use for solving (sub)problems.
solver_options : dict, optional
Options to pass to calls to the solver. Keys are the name of
the functions being called and are one of `__call__`, `set_instance_`,
`solve`.
`__call__` is the function being called when generating an instance
with the pyomo SolverFactory. Additionally to the options provided,
`node_ids` is passed to this call containing the IDs of the entities
being optimized.
`set_instance` is called when a pyomo Model is set as an instance of
a persistent solver. `solve` is called to perform an optimization. If
not set, `save_results` and `load_solutions` may be set to false to
provide a speedup.
mode : str, optional
Specifies which set of constraints to use.
- `convex` : Use linear constraints
- `integer` : May use non-linear constraints
eps_primal : float, optional
Primal stopping criterion for the ADMM algorithm.
eps_dual : float, optional
Dual stopping criterion for the ADMM algorithm.
rho : float, optional
Stepsize for the ADMM algorithm.
max_iterations : int, optional
Maximum number of ADMM iterations.
robustness : tuple, optional
Tuple of two floats. First entry defines how many time steps are
protected from deviations. Second entry defines the magnitude of
deviations which are considered.
References
----------
.. [1] "Alternating Direction Method of Multipliers for Decentralized
Electric Vehicle Charging Control" by Jose Rivera, Philipp Wolfrum,
Sandra Hirche, Christoph Goebel, and Hans-Arno Jacobsen
Online: https://mediatum.ub.tum.de/doc/1187583/1187583.pdf (accessed on 2020/09/28)
"""
def __init__(self, city_district, solver=DEFAULT_SOLVER, solver_options=DEFAULT_SOLVER_OPTIONS, mode="convex",
eps_primal=0.1, eps_dual=1.0, rho=2.0, max_iterations=10000, robustness=None):
super(ExchangeADMMMPI, self).__init__(city_district, solver, solver_options, mode)
self.mpi_interface = mpi_interface.MPI_Interface()
self.eps_primal = eps_primal
self.eps_dual = eps_dual
self.rho = rho
self.max_iterations = max_iterations
# create solver nodes for each entity
self.nodes = [
SolverNode(solver, solver_options, [entity], mode, robustness=robustness)
for entity in self.entities
]
# create pyomo parameters for each entity
for node, entity in zip(self.nodes, self.entities):
node.model.beta = pyomo.Param(mutable=True, initialize=1)
node.model.xs_ = pyomo.Param(entity.model.t, mutable=True, initialize=0)
node.model.us = pyomo.Param(entity.model.t, mutable=True, initialize=0)
node.model.last_p_el_schedules = pyomo.Param(entity.model.t, mutable=True, initialize=0)
self._add_objective()
def _add_objective(self):
for i, node, entity in zip(range(len(self.entities)), self.nodes, self.entities):
obj = node.model.beta * entity.get_objective()
obj += self.rho / 2 * pyomo.sum_product(entity.model.p_el_vars, entity.model.p_el_vars)
# penalty term is expanded and constant is omitted
if i == 0:
# invert sign of p_el_schedule and p_el_vars (omitted for quadratic
# term)
obj += self.rho * pyomo.sum_product(
[(-node.model.last_p_el_schedules[t] - node.model.xs_[t] - node.model.us[t])
for t in range(entity.op_horizon)],
entity.model.p_el_vars
)
else:
obj += self.rho * pyomo.sum_product(
[(-node.model.last_p_el_schedules[t] + node.model.xs_[t] + node.model.us[t])
for t in range(entity.op_horizon)],
entity.model.p_el_vars
)
node.model.o = pyomo.Objective(expr=obj)
return
def _presolve(self, full_update, beta, robustness, debug):
results, params = super()._presolve(full_update, beta, robustness, debug)
for node, entity in zip(self.nodes, self.entities):
node.model.beta = self._get_beta(params, entity)
if full_update:
node.full_update(robustness)
results["r_norms"] = []
results["s_norms"] = []
return results, params
def _is_last_iteration(self, results, params):
return results["r_norms"][-1] <= self.eps_primal and results["s_norms"][-1] <= self.eps_dual
def _iteration(self, results, params, debug):
super(ExchangeADMMMPI, self)._iteration(results, params, debug)
op_horizon = self.entities[0].op_horizon
# fill parameters if not already present
if "p_el" not in params:
params["p_el"] = np.zeros((len(self.entities), op_horizon))
if "x_" not in params:
params["x_"] = np.zeros(op_horizon)
if "u" not in params:
params["u"] = np.zeros(op_horizon)
u = params["u"]
# -----------------
# 1) optimize all entities
# -----------------
to_solve_nodes = []
variables = []
# Take actions if the number of MPI processes does not fit the total number of nodes:
if self.mpi_interface.mpi_size > len(self.nodes):
rank_id_range = np.array([i for i in range(len(self.nodes))])
if self.mpi_interface.mpi_rank > len(self.nodes):
results["r_norms"].append(0.0)
results["s_norms"].append(0.0)
return
elif self.mpi_interface.mpi_size < len(self.nodes):
if self.mpi_interface.mpi_size == 1:
rank_id_range = np.array([0 for i in range(len(self.nodes))])
else:
a, b = divmod(len(self.nodes)-1, self.mpi_interface.mpi_size-1)
rank_id_range = np.repeat(np.array([i for i in range(1, self.mpi_interface.mpi_size)]), a)
for i in range(b):
rank_id_range = np.append(rank_id_range, i+1)
rank_id_range = np.concatenate([[0], rank_id_range])
else:
rank_id_range = np.array([i for i in range(len(self.nodes))])
rank_id_range = np.sort(rank_id_range)
p_el_schedules = np.empty((len(self.entities), op_horizon))
for i, node, entity in zip(range(len(self.nodes)), self.nodes, self.entities):
if self.mpi_interface.mpi_rank == rank_id_range[i]:
if not isinstance(
entity,
(CityDistrict, Building, Photovoltaic, WindEnergyConverter)
):
continue
for t in range(op_horizon):
node.model.last_p_el_schedules[t] = params["p_el"][i][t]
node.model.xs_[t] = params["x_"][t]
node.model.us[t] = params["u"][t]
node.obj_update()
to_solve_nodes.append(node)
variables.append([entity.model.p_el_vars[t] for t in range(op_horizon)])
self._solve_nodes(results, params, to_solve_nodes, variables=variables, debug=debug)
if self.mpi_interface.mpi_rank == 0:
p_el_schedules[0] = np.array(extract_pyomo_values(self.entities[0].model.p_el_vars, float),
dtype=np.float64)
for j in range(1, len(self.nodes)):
if not isinstance(
self.entities[j],
(CityDistrict, Building, Photovoltaic, WindEnergyConverter)
):
continue
if self.mpi_interface.mpi_rank == 0:
if self.mpi_interface.mpi_size > 1:
data = np.empty(op_horizon, dtype=np.float64)
self.mpi_interface.mpi_comm.Recv(data, source=rank_id_range[j],
tag=int(results["iterations"][-1])*len(self.nodes)+j)
p_el_schedules[j] = np.array(data, dtype=np.float64)
else:
p_el_schedules[j] = np.array(extract_pyomo_values(self.entities[j].model.p_el_vars, float),
dtype=np.float64)
else:
if self.mpi_interface.mpi_rank == rank_id_range[j]:
p_el_schedules[j] = np.array(extract_pyomo_values(self.entities[j].model.p_el_vars, float),
dtype=np.float64)
if self.mpi_interface.mpi_size > 1:
self.mpi_interface.mpi_comm.Send(p_el_schedules[j], dest=0,
tag=int(results["iterations"][-1])*len(self.nodes)+j)
# --------------------------
# 2) incentive signal update
# --------------------------
if self.mpi_interface.mpi_rank == 0:
x_ = (-p_el_schedules[0] + sum(p_el_schedules[1:])) / len(self.entities)
u += x_
r_norm = np.array([np.math.sqrt(len(self.entities)) * np.linalg.norm(x_)], dtype=np.float64)
s = np.zeros_like(p_el_schedules)
s[0] = - self.rho * (-p_el_schedules[0] + params["p_el"][0] + params["x_"] - x_)
for i in range(1, len(self.entities)):
s[i] = - self.rho * (p_el_schedules[i] - params["p_el"][i] + params["x_"] - x_)
s_norm = np.array([np.linalg.norm(s.flatten())], dtype=np.float64)
else:
x_ = np.empty(op_horizon, dtype=np.float64)
u = np.empty(op_horizon, dtype=np.float64)
r_norm = np.empty(1, dtype=np.float64)
s_norm = np.empty(1, dtype=np.float64)
if self.mpi_interface.mpi_size > 1:
self.mpi_interface.mpi_comm.Bcast(x_, root=0)
self.mpi_interface.mpi_comm.Bcast(u, root=0)
self.mpi_interface.mpi_comm.Bcast(r_norm, root=0)
self.mpi_interface.mpi_comm.Bcast(s_norm, root=0)
# ------------------------------------------
# 3) Calculate parameters for stopping criteria