diff --git a/README.md b/README.md new file mode 100644 index 0000000..44fac90 --- /dev/null +++ b/README.md @@ -0,0 +1,40 @@ +

Demand Management in Smart Grids

+
+

A diploma thesis investigating the options of controlling power demand of households to reduce peaks in total power consumption in smart grids.

+

Simulation results for February 2018

+ +The problem +----------- +While the widespread adoption of electric vehicles in the future could serve to reduce our reliance on fossil fuels, provided that the energy used to charge those vehicles comes from cleaner sources, the increased electricity usage that would come as a side effect could constitute a challenge to the electrical grid. Most trips by household-owned vehicles are trips to work and back, and most people plug in their car to charge it as soon as they arrive home. As the bulk of the people arrive home around the same time, this can cause peaks in electricity consumption in the afternoon hours. + +Because traditional coal or nuclear power plants have a long ramp-up time, they are ill-suited to be used to cover these sudden peaks, and specialized peaking power plants such as gas turbines or pumped hydroelectric energy storage must be used, which is costly and not always feasible. The situation is further complicated as power plants utilizing renewable energy sources are being introduced to the grid. These plants, such as solar panel farms or wind turbines, often have an uncontrollable and intermittent power output, requiring even more energy storage to compensate. + +To combat these negative effects, there arises a need to be able to better control the electricity demand of the grid by matching the power requirements of electricity subscribers to the generation profile of the power plants. + +There are several proposed strategies to achieve that goal, but they have several key issues, one of them being their privacy implications, as most of them require the EV owners to inform the distributor or another third party about their travel schedule and charging needs, which some might not be willing to do. Another weak point is that the strategies control only electric cars, and do not take into account other appliances that could be controlled. Furthemore, most of the strategies require that subscribers are willing to adhere to them, even though it might not be directly beneficial to them. + +Our proposed solution +--------------------- +To tackle the issues with the current strategies, we decided to employ the strategy of motivating the subscribers to use energy at specific times by means of time-of-use electricity price tariffs, with different timing of low and high prices for each connected household. + +The electricity distributor decides on a total target electricity demand profile which it would like to achieve on its grid according to its needs (e.g. a flat demand profile for a grid consisting of power plants with slow demand-response characteristics, or a profile with higher demand during specific intervals for a grid containing solar plants or wind turbines when expecting a lot of sunlight or wind). + +A part of this demand does not come from subscribers whose demand can be influenced easily, but instead from other sectors like businesses and industrial customers. The part of the demand the distributor can influence, like the demand of households, can be used for balancing the demand of the rest of the grid to achieve the desired total demand. The distributor must therefore separate the total target demand profile into the uninfluenceable base demand and the influenceable target household demand. + +The electricity prices for each household are then assigned so that more households get a high price when their demand should be low, and a low price when their demand should be high. The electricity prices are finally distributed to the households, and the decision on when to turn their appliances on is left to them, expecting them to act in their own interest and optimize their appliance usage to achieve the lowest total cost of electricity. + +This achieves the goals described earlier — privacy of the subscribers is protected, since they do not need to divulge the information about the usage of their appliances; multiple types of appliances can be controlled, since they only need to be able to adjust their operation to the electricity price; and the subscribers are incentivized to use energy in coordination with the needs of the distributor, as it will also be cheaper for them. + +A more detailed description of our solution can be found in the [thesis text](thesis/thesis.pdf). + +[Simulator](simulator/) +----------------------- +For evaluating the effectiveness of the demand control algorithm, we created a simulator of the smart grid. The simulator currently simulates the Texas power grid, based on data from the [Electricity Reliability Council of Texas](http://www.ercot.com), the [Pecan Street organization](https://www.pecanstreet.org/) and the [National Household Travel Survey](https://nhts.ornl.gov/), which it can directly download and parse. + +Its structure was designed to resemble a real electrical distribution grid, with each part of the grid being a separate entity, operating independently, keeping its own clock and performing calculations exactly when needed, according to an environment-provided clock signal. This was done to allow for easier possible future extension of the simulator or its separation into multiple programs communicating over some network interface to enable simulating communication delays and other physical limitations of the distribution grid. + +The simulator implements our algorithm described in the thesis, and two simpler demand control algorithms used for comparison with the first one — an algorithm where appliances stretch their power demand across their whole usage interval, and an algorithm where they perform their operation as early as possible. + +The simulator is written in Python, and it's purpose-built for this thesis, so extending it for different uses would be quite comples, but the legwork is there. There are certainly some improvements that could be made, most prominently the use of multithreading, but that shouldn't be too complicated now with the multiprocessing improvements in Python 3.8. + +To try the simulator, download the [latest release](https://www.github.com/fnesveda/DemandManagement/releases/latest/) and unpack it. Usage instructions can be found in the included [README file](simulator/README.md). diff --git a/docs/images/.gitattributes b/docs/images/.gitattributes new file mode 100644 index 0000000..58da693 --- /dev/null +++ b/docs/images/.gitattributes @@ -0,0 +1 @@ +*.svg binary diff --git a/docs/images/screenshots/main.svg b/docs/images/screenshots/main.svg new file mode 100644 index 0000000..e0ea149 --- /dev/null +++ b/docs/images/screenshots/main.svg @@ -0,0 +1,211 @@ + + + + + + + + + + + + + + + + + + + + + Feb 19 + 00:00 + + + + + + + + + Feb 20 + 00:00 + + + + + + + + + Feb 21 + 00:00 + + + + + + + + + Feb 22 + 00:00 + + + + + + + + + 12:00 + + + + + + + + + 12:00 + + + + + + + + + 12:00 + + + + + + + + + + + 0 + + + + + + + + + 10 000 + + + + + + + + + 20 000 + + + + + + + + + 30 000 + + + + + + + + + 40 000 + + + + + + + + + 50 000 + + + + + + + + + 60 000 + + + + Demand [MW] + + + + + + + + + + + + + + + + Grid demand simulation for February 2018 + + + + + + + + + + + + + + Target demand + + + + + + + Base demand + + + + + + + Smart demand + + + + + + + Uncontrolled demand + + + + + + + + + + diff --git a/simulator/README.md b/simulator/README.md new file mode 100644 index 0000000..9645d90 --- /dev/null +++ b/simulator/README.md @@ -0,0 +1,67 @@ +Smart grid simulator +==================== + +This is a simulator of a smart grid made as a part of the masters thesis "Demand control in smart grids" at the Faculty of Mathematics and Physics of the Charles University in Prague. +It serves to simulate an electrical grid with households and appliances based on data from the National Household Travel Survey +and the Pecan Street organization research program, +while the appliances act according to several demand control algorithms described in the thesis, +which is included in the file _thesis.pdf_. + +Requirements +------------ + +The simulator can be run in an UNIX environment containing the Bash shell and a Python interpreter version 3.6 or newer. +Several third-party Python libraries are required, these are listed in the file _requirements.txt_ +and available for installation with the shell command `pip3 install -r requirements.txt`. + +Data download +------------- +Before the actual usage of the simulator, the data necessary for the simulation must be downloaded first. +This is performed in three steps: + +1. In the file _simulator/data/config.txt_, put the desired dates in the `fromdate` and `todate` fields, in the format `YYYY-MM-DD`. For a succesful simulation, +the data must be downloaded for an interval starting one week before the starting date of the simulation +and ending one week after the ending date of the simulation. + +2. In the file _simulator/data/dataport/KEY_, put the username and password +to access the Dataport database in the `username` and `password` fields, respectively. +These access credentials can be obtained at https://dataport.pecanstreet.org/. + +3. Execute the data download by using the command script `./downloadData.sh` in a terminal while being in the directory with the project. +Depending on the length of the downloaded interval, the download can take several hours and use several gigabytes of data. +Progress of the download is printed to the terminal during the run of the script. + +Running the simulator +--------------------- + +After the data is downloaded, the simulator can be executed by using the command +`./run.py startingDate simulationLength householdCount outputFolder` in a terminal while being in the directory with the project, +while replacing `startingDate` with the date from which to run the simulation in the format `YYYY-MM-DD`, +`simulationLength` with the number of days to simulate, +`householdCount` with the number of households to simulate +and `outputFolder` with the destination folder in which the simulation results should be saved, +for example: `./run.py 2018-01-01 365 10000 out/`. + +While the simulation is running, the simulator prints information about its progress to the terminal. +When the simulation finishes, the results are saved in the specified folder in two files, _desc.txt_ and _data.csv_. +The file _desc.txt_ contains information about the simulation parameters, +and the file _data.csv_ is a standard comma-separated values file containing the simulation results organized to eight columns: + +Column | Description +--------------------|------------ +Datetime | The date and time of the result +PredictedBaseDemand | The predicted grid base demand at that date and time (in kilowatts) +ActualBaseDemand | The actual grid base demand at that date and time (in kilowatts) +TargetDemand | The target power demand at that date and time (in kilowatts) +SmartDemand | The power demand of the simulated households at that date and time if the smart algorithm was used (in kilowatts) +UncontrolledDemand | The power demand of the simulated households at that date and time if the uncontrolled algorithm was used (in kilowatts) +SpreadOutDemand | The power demand of the simulated households at that date and time if the spread out algorithm was used (in kilowatts) +PriceRatio | The ratio of how many households should have a lower electricity price at that date and time + +Displaying the results +---------------------- + +To display the simulation results, one can use the included Jupyter notebook +_Results.ipynb_. After opening the notebook, first the folder with the simulation +results must be specified in the variable _resultsFolder_ in the third code cell, +and then the cells in the notebook can be executed to show the results. diff --git a/simulator/Results.ipynb b/simulator/Results.ipynb new file mode 100644 index 0000000..3f835f0 --- /dev/null +++ b/simulator/Results.ipynb @@ -0,0 +1,247 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Matplotlib preparation**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib\n", + "%matplotlib notebook\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.dates as mdates\n", + "\n", + "matplotlib.rcParams['figure.figsize'] = (9.8, 5)\n", + "matplotlib.rcParams['figure.constrained_layout.use'] = True\n", + "matplotlib.rcParams['font.size'] = 13\n", + "\n", + "# without this matplotlib emits a warning \n", + "# even when not plotting pandas data at all, because pandas messes with it\n", + "from pandas.plotting import register_matplotlib_converters\n", + "register_matplotlib_converters()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Imports necessary for loading the data**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "import pandas" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Folder with the simulation results**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "resultsFolder = \"out/\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Loading the data**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "desc = {}\n", + "with open(f\"{resultsFolder}/desc.txt\", \"r\") as descfile:\n", + "\tfor line in descfile:\n", + "\t\tkey, val = line.strip(\"\\n\").split(\"=\")\n", + "\t\tdesc[key] = val\n", + "\n", + "startingDT = datetime.datetime.strptime(desc[\"startingDatetime\"], \"%Y-%m-%d %H:%M:%S\")\n", + "simulationLength = int(desc[\"simulationLength\"])\n", + "houseCount = int(desc[\"houseCount\"])\n", + "\n", + "\n", + "data = pandas.read_csv(f\"{resultsFolder}/data.csv\", parse_dates=[0])\n", + "\n", + "datetimes = data[\"Datetime\"].values\n", + "predictedDemand = data[\"PredictedBaseDemand\"].values\n", + "actualDemand = data[\"ActualBaseDemand\"].values\n", + "targetDemand = data[\"TargetDemand\"].values\n", + "smartDemand = data[\"SmartDemand\"].values\n", + "uncontrolledDemand = data[\"UncontrolledDemand\"].values\n", + "spreadOutDemand = data[\"SpreadOutDemand\"].values\n", + "priceRatio = data[\"PriceRatio\"].values" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Plots of the simulation results**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(\"%d. %m.\\n%H:%M\"))\n", + "plt.gca().xaxis.set_minor_formatter(mdates.DateFormatter(\"%H:%M\"))\n", + "plt.gca().xaxis.set_major_locator(mdates.HourLocator(byhour=[0]))\n", + "plt.title(\"Base demand prediction vs. actual base demand\")\n", + "plt.plot(datetimes, predictedDemand, label=\"Predicted base demand\")\n", + "plt.plot(datetimes, actualDemand, label=\"Actual base demand\")\n", + "plt.plot(datetimes[0], [0])\n", + "plt.grid()\n", + "plt.legend(loc=\"upper center\", ncol=2)\n", + "plt.xlabel(\"Date and time\")\n", + "plt.ylabel(\"Demand [kW]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(\"%d. %m.\\n%H:%M\"))\n", + "plt.gca().xaxis.set_minor_formatter(mdates.DateFormatter(\"%H:%M\"))\n", + "plt.gca().xaxis.set_major_locator(mdates.HourLocator(byhour=[0]))\n", + "plt.title(\"Ideal target demand the smart homes should try to reach\")\n", + "plt.plot(datetimes, predictedDemand, label=\"Predicted base demand\")\n", + "plt.plot(datetimes, predictedDemand + targetDemand, label=\"Predicted base demand + target demand\")\n", + "plt.plot(datetimes[0], [0])\n", + "plt.grid()\n", + "plt.legend(loc=\"upper center\", ncol=2)\n", + "plt.xlabel(\"Date and time\")\n", + "plt.ylabel(\"Demand [kW]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(\"%d. %m.\\n%H:%M\"))\n", + "plt.gca().xaxis.set_minor_formatter(mdates.DateFormatter(\"%H:%M\"))\n", + "plt.gca().xaxis.set_major_locator(mdates.HourLocator(byhour=[0]))\n", + "plt.title(\"Actual target demand the smart homes will try to reach\")\n", + "plt.plot(datetimes, actualDemand, label=\"Actual base demand\")\n", + "plt.plot(datetimes, actualDemand + targetDemand, label=\"Actual base demand + target demand\")\n", + "plt.plot(datetimes[0], [0])\n", + "plt.grid()\n", + "plt.legend(loc=\"upper center\", ncol=2)\n", + "plt.xlabel(\"Date and time\")\n", + "plt.ylabel(\"Demand [kW]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(\"%d. %m.\\n%H:%M\"))\n", + "plt.gca().xaxis.set_minor_formatter(mdates.DateFormatter(\"%H:%M\"))\n", + "plt.gca().xaxis.set_major_locator(mdates.HourLocator(byhour=[0]))\n", + "plt.title(\"Target demand of the smart homes and the cheaper prices ratio generated from it\")\n", + "plt.plot(datetimes, targetDemand, label=\"Target demand\")\n", + "plt.plot(datetimes, priceRatio * houseCount, label=\"Price ratio\")\n", + "plt.plot(datetimes[0], [0])\n", + "plt.grid()\n", + "plt.legend(loc=\"upper center\", ncol=2)\n", + "plt.xlabel(\"Date and time\")\n", + "plt.ylabel(\"Demand [kW]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.gca().xaxis.set_major_formatter(mdates.DateFormatter(\"%d. %m.\\n%H:%M\"))\n", + "plt.gca().xaxis.set_minor_formatter(mdates.DateFormatter(\"%H:%M\"))\n", + "plt.gca().xaxis.set_major_locator(mdates.HourLocator(byhour=[0]))\n", + "plt.title(\"The simulated demand of the households according to different optimization algorithms\")\n", + "plt.plot(datetimes, actualDemand, label=\"Base demand\", zorder=0)\n", + "plt.plot(datetimes, actualDemand + targetDemand, label=\"Target demand\", zorder=4)\n", + "plt.plot(datetimes, actualDemand + uncontrolledDemand, label=\"Uncontrolled demand\", zorder=1)\n", + "plt.plot(datetimes, actualDemand + spreadOutDemand, label=\"Spread out demand\", zorder=2)\n", + "plt.plot(datetimes, actualDemand + smartDemand, label=\"Smart demand\", zorder=3)\n", + "plt.plot(datetimes[0], [0])\n", + "plt.grid()\n", + "plt.legend(loc=\"lower center\", ncol=3)\n", + "plt.xlabel(\"Date and time\")\n", + "plt.ylabel(\"Demand [kW]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/simulator/downloadData.sh b/simulator/downloadData.sh new file mode 100755 index 0000000..5496e6e --- /dev/null +++ b/simulator/downloadData.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env bash + +DIR=$(dirname "${BASH_SOURCE[0]}") +pushd "${DIR}" > /dev/null +simulator/data/download.sh +popd > /dev/null diff --git a/simulator/requirements.txt b/simulator/requirements.txt new file mode 100644 index 0000000..c2ffe07 --- /dev/null +++ b/simulator/requirements.txt @@ -0,0 +1,7 @@ +matplotlib>=3.1.1 +numpy>=1.16.4 +pandas>=0.24.2 +psycopg2-binary>=2.8.3 +requests>=2.22.0 +scipy>=1.3.0 +SQLAlchemy>=1.3.5 diff --git a/simulator/run.py b/simulator/run.py new file mode 100755 index 0000000..b12318e --- /dev/null +++ b/simulator/run.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python3 + +import datetime +import sys + +def usage(): + print(f"Usage: {sys.argv[0]} startingDate simulationLength houseCount outputFolder") + +if len(sys.argv) != 5: + usage() + sys.exit(1) + +try: + startingDate = datetime.datetime.strptime(sys.argv[1], "%Y-%m-%d") + simulationLength = max(0, int(sys.argv[2])) + houseCount = max(0, int(sys.argv[3])) + outputFolder = sys.argv[4] # mostly anything can be a path under POSIX, this would be too hard to validate anyway +except ValueError: + usage() + sys.exit(1) + +# importing after already running code is evil, I know +# but this takes a long time and would just delay the parameter checking +from simulator import Simulator + +Simulator.run(startingDate, simulationLength, houseCount, outputFolder) diff --git a/simulator/simulator/__init__.py b/simulator/simulator/__init__.py new file mode 100644 index 0000000..6d1b4e0 --- /dev/null +++ b/simulator/simulator/__init__.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python3 + +from .simulator import Simulator diff --git a/simulator/simulator/appliance.py b/simulator/simulator/appliance.py new file mode 100755 index 0000000..a046a60 --- /dev/null +++ b/simulator/simulator/appliance.py @@ -0,0 +1,603 @@ +#!/usr/bin/env python3 + +import datetime +import math +import random + +from abc import ABC, abstractmethod, abstractclassmethod +from types import SimpleNamespace + +import numpy + +from . import applianceStatistics +from . import utils + +from .constants import oneDay +from .profile import Profile +from .utils import minutesIn + +# abstract base class for household appliances +class Appliance(ABC): + # current date and time in the simulation + currentDT: datetime.datetime + # an object for storing information between demand calculations + memory: SimpleNamespace + # appliance usage statistics + usageStatistics: applianceStatistics.ApplianceStatistics + # the electricity price for any given minute in the simulation + priceProfile: Profile + # the electricity demand of this appliance if it was smart + smartDemand: Profile + # the electricity demand of this appliance if it would charge as early as possible + uncontrolledDemand: Profile + # the electricity demand of this appliance if it would charge evenly over their use period + spreadOutDemand: Profile + + # constructor, just prepares the variables + def __init__(self): + self.memory = SimpleNamespace() + self.priceProfile = Profile() + self.smartDemand = Profile() + self.uncontrolledDemand = Profile() + self.spreadOutDemand = Profile() + + # creates a random appliance with the right parameters + @abstractclassmethod + def random(cls): + pass + + # sets up the appliance for the simulation + def setUp(self, dt: datetime.datetime): + self.currentDT = dt + # generate usage for one day in the future + self.generateUsage(dt, dt + oneDay) + + # moves ahead one day and does all the calculations that need to be done in that day + def tick(self): + # remove past, unneeded values from the profiles to free up some memory + self.priceProfile.prune(self.currentDT - oneDay) + self.smartDemand.prune(self.currentDT - oneDay) + self.uncontrolledDemand.prune(self.currentDT - oneDay) + self.spreadOutDemand.prune(self.currentDT - oneDay) + + # generate appliance usage for one more day in the future + self.generateUsage(self.currentDT + oneDay, self.currentDT + 2 * oneDay) + # calculate the power demand for the next day + self.calculateDemand(self.currentDT, self.currentDT + oneDay) + + # move ahead one day + self.currentDT += oneDay + + # generates appliance usage for a given time interval + @abstractmethod + def generateUsage(self, fromDT: datetime.datetime, toDT: datetime.datetime): + pass + + # calculates appliance power demand for a given time interval + def calculateDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + self.calculateSmartDemand(fromDT, toDT) + self.calculateUncontrolledDemand(fromDT, toDT) + self.calculateSpreadOutDemand(fromDT, toDT) + + # calculates appliance power demand for a given time interval acting as if the appliance was smart + @abstractmethod + def calculateSmartDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + pass + + # calculates appliance power demand for a given time interval acting as if if it would charge as early as possible + @abstractmethod + def calculateUncontrolledDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + pass + + # calculates appliance power demand for a given time interval acting as if if it would charge evenly over the use period + @abstractmethod + def calculateSpreadOutDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + pass + + # sets the electricity price profile for a given time interval + def setPriceProfile(self, dt: datetime.datetime, prices: numpy.ndarray): + self.priceProfile.set(dt, prices) + +# abstract base class for battery-based household appliances (e.g. electric car) +class Battery(Appliance, ABC): + # appliance usage statistics + usageStatistics: applianceStatistics.BatteryStatistics + + # the power with which the battery charges + chargingPower: float # kW + + # constructor, just prepares the variables + def __init__(self, chargingPower: float = 0): + self.chargingPower = chargingPower + super().__init__() + # for each day keeps disconnection time, connection time and charge needed after the usage + self.memory.usages = dict() + + # creates a random appliance with the right parameters + @classmethod + def random(cls): + chargingPower = cls.usageStatistics.randomChargingPower() + return cls(chargingPower=chargingPower) + + # generates appliance usage for a given time interval + def generateUsage(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # we need to know the usage for one day ahead (at least the disconnect time) + for midnight in utils.midnightsBetween(fromDT, toDT+oneDay): + date = midnight.date() + if date not in self.memory.usages: + disconnectionTime, connectionTime = self.usageStatistics.randomUsageInterval(date) + # if the values are invalid make up an usage interval giving us as much charging time as possible + if disconnectionTime is None: + disconnectionTime = datetime.time(23, 59) + if connectionTime is None: + connectionTime = datetime.time(00, 00) + + chargeNeeded = self.usageStatistics.randomNeededCharge(date) + self.memory.usages[date] = ((disconnectionTime, connectionTime), chargeNeeded) + + # calculates appliance power demand for a given time interval acting as if the appliance was smart + def calculateSmartDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day in the interval, it charges the battery in the minutes with the cheapest electricity available + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + powerProfile = numpy.zeros(2*minutesIn(oneDay), dtype=float) + + # get the needed charge and the interval when the battery is connected + ((_, connectionTime), chargeNeeded) = self.memory.usages[date] + ((disconnectionTime, _), _) = self.memory.usages[date + oneDay] + + # get for how long the battery must be charged + chargePerSlot = self.chargingPower / 60 # kWh per minute + slotsToChargeCompletely = math.ceil(chargeNeeded / chargePerSlot) + + # if it needs to be charged, charge it + if slotsToChargeCompletely > 0: + connectionSlot = connectionTime.hour * 60 + connectionTime.minute + disconnectionSlot = minutesIn(oneDay) + disconnectionTime.hour * 60 + disconnectionTime.minute + + # if there is not enough time to charge the battery completely, just charge it all the available time + if disconnectionSlot - connectionSlot <= slotsToChargeCompletely: + powerProfile[connectionSlot:disconnectionSlot] = self.chargingPower + # otherwise pick enough of the cheapest time slots and charge the battery during those + else: + priceProfile = self.priceProfile.get(midnight, midnight+2*oneDay) + cheapestSlots = numpy.argpartition(priceProfile[connectionSlot:disconnectionSlot], slotsToChargeCompletely)[:slotsToChargeCompletely] + connectionSlot + powerProfile[cheapestSlots[:-1]] = self.chargingPower + lastSlotCharge = chargeNeeded - (chargePerSlot * (slotsToChargeCompletely - 1)) + powerProfile[cheapestSlots[-1]] = lastSlotCharge * 60 + + self.smartDemand.add(midnight, powerProfile) + + # calculates appliance power demand for a given time interval acting as if the battery wanted to charge as early as possible + def calculateUncontrolledDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day in the interval, the appliance starts charging the battery as soon as it is connected to power + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + + powerProfile = numpy.zeros(2*minutesIn(oneDay), dtype=float) + + # get the needed charge and the interval when the battery is connected + ((_, connectionTime), chargeNeeded) = self.memory.usages[date] + ((disconnectionTime, _), _) = self.memory.usages[date + oneDay] + + # get for how long the battery must be charged + chargePerSlot = self.chargingPower / 60 + slotsToChargeCompletely = math.ceil(chargeNeeded / chargePerSlot) + + # if it needs to be charged, charge it + if slotsToChargeCompletely > 0: + connectionSlot = connectionTime.hour * 60 + connectionTime.minute + disconnectionSlot = minutesIn(oneDay) + disconnectionTime.hour * 60 + disconnectionTime.minute + + # if there is not enough time to charge the battery completely, just charge it all the available time + if disconnectionSlot - connectionSlot < slotsToChargeCompletely: + powerProfile[connectionSlot:disconnectionSlot] = self.chargingPower + # otherwise start charging it as soon as it is available and charge until it's full + else: + powerProfile[connectionSlot:connectionSlot+slotsToChargeCompletely-1] = self.chargingPower + lastSlotCharge = chargeNeeded - (chargePerSlot * (slotsToChargeCompletely - 1)) + powerProfile[connectionSlot+slotsToChargeCompletely] = lastSlotCharge * 60 + + self.uncontrolledDemand.add(midnight, powerProfile) + + # calculates appliance power demand for a given time interval acting as if the battery wanted to charge as evenly as possible + def calculateSpreadOutDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day in the interval, the appliance starts charging the battery as soon as it is connected to power + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + + powerProfile = numpy.zeros(2*minutesIn(oneDay), dtype=float) + + # get the needed charge and the interval when the battery is connected + ((_, connectionTime), chargeNeeded) = self.memory.usages[date] + ((disconnectionTime, _), _) = self.memory.usages[date + oneDay] + + # get for how long the battery must be charged + chargePerSlot = self.chargingPower / 60 + slotsToChargeCompletely = math.ceil(chargeNeeded / chargePerSlot) + + # if it needs to be charged, charge it + if slotsToChargeCompletely > 0: + connectionSlot = connectionTime.hour * 60 + connectionTime.minute + disconnectionSlot = minutesIn(oneDay) + disconnectionTime.hour * 60 + disconnectionTime.minute + + # if there is not enough time to charge the battery completely, just charge it all the available time + if disconnectionSlot - connectionSlot < slotsToChargeCompletely: + powerProfile[connectionSlot:disconnectionSlot] = self.chargingPower + # otherwise charge it evenly over the whole connected period + else: + powerProfile[connectionSlot:disconnectionSlot] = chargeNeeded / ((disconnectionSlot - connectionSlot) / 60) + + self.spreadOutDemand.add(midnight, powerProfile) + +# abstract base class for accumulator-based household appliances (e.g. water heater, refrigerator) +class Accumulator(Appliance, ABC): + # appliance usage statistics + usageStatistics: applianceStatistics.AccumulatorStatistics + # the power with which the appliance charges + chargingPower: float # kW + # the charging capacity of the appliance + capacity: float # kWh + + # constructor, just prepares the variables + def __init__(self, chargingPower: float, capacity: float, dischargingProfileScale: float): + super().__init__() + self.chargingPower = chargingPower + self.capacity = capacity + + # variables for storing the state of the appliance between calculations + self.memory.dischargingProfileScale = dischargingProfileScale + self.memory.smart = SimpleNamespace() + self.memory.smart.currentCharge = random.random() * self.capacity + self.memory.uncontrolled = SimpleNamespace() + self.memory.uncontrolled.currentCharge = self.capacity + self.memory.spreadOut = SimpleNamespace() + self.memory.spreadOut.charging = random.choice([True, False]) + self.memory.spreadOut.currentCharge = random.random() * self.capacity + + # the profile of how the accumulator discharges (e.g. water heater cools down or gets used, fridge heats up) + self.memory.dischargingProfile = Profile() + + # creates a random accumulator with the right parameters + @classmethod + def random(cls): + chargingPower = cls.usageStatistics.randomChargingPower() + # in this simulator we can't have an appliance which would charge up faster than in one minute + capacity = max(cls.usageStatistics.randomCapacity(), (1.1*chargingPower/60)) + # stronger appliances are usually those which get used more, so scale the random discharging profile by the charging power of the appliance + dischargingProfileScale = cls.usageStatistics.randomDischargingProfileScale() * (chargingPower / cls.usageStatistics.averageChargingPower) + return cls(chargingPower, capacity, dischargingProfileScale) + + # moves ahead one day and does all the calculations that need to be done in that day + def tick(self): + # remove past, unneeded values from the profiles to free up some memory + self.memory.dischargingProfile.prune(self.currentDT - oneDay) + super().tick() + + # generates appliance usage for a given time interval + def generateUsage(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # generate how much will the accumulator discharge during that interval + dischargingProfile = self.usageStatistics.dischargingProfile.get(fromDT, toDT) * self.memory.dischargingProfileScale + self.memory.dischargingProfile.set(fromDT, dischargingProfile) + + # calculates appliance power demand for a given time interval acting as if the appliance was smart + def calculateSmartDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # adapted from https://ktiml.mff.cuni.cz/~fink/publication/greedy.pdf + # asymptotically, this would be faster with prefix sum trees or union-find data structures + # but in Python that is actually spreadOuter than using a quadratic algorithm with numpy + # the limit preparation, which is linear, is spreadOuter than the actual algorithm anyway, so it doesn't matter + + # get the cheapest slots in which to turn on the appliance so that it never discharges under its lower limit and never charges over its upper limit + # normally, the appliance doesn't charge more than it needs to, so at the end of the interval it would be charged barely above the lower limit + # then at the start of the next calculated interval it would need to charge more to catch up + # we calculate the charging profile with a bit of an overlap to avoid this + endMargin = oneDay + wantedSlots = utils.minutesBetween(fromDT, toDT) + totalSlots = utils.minutesBetween(fromDT, toDT+endMargin) + + # the electricity prices in the interval + priceProfile = self.priceProfile.get(fromDT, toDT+endMargin) + + # how much energy goes into the appliance each minute it's turned on + chargingRate = self.chargingPower / 60 # kWh per minute + # how much energy leaves the appliance each minute of the interval + dischargingRates = self.memory.dischargingProfile.get(fromDT, toDT+endMargin) / 60 # kWh per minute for each minute + + startingCharge = self.memory.smart.currentCharge # kWh + + # never discharge past 0 and never charge over the capacity + lowerTarget = 0 # kWh + upperTarget = self.capacity # kWh + + # convert the limits and charging rates to integer steps + dischargingSum = numpy.cumsum(dischargingRates) # total kWh cumulatively discharged for each minute + lowerLimit = numpy.ceil((lowerTarget - startingCharge + dischargingSum) / chargingRate).astype(int) + upperLimit = numpy.floor((upperTarget - startingCharge + dischargingSum) / chargingRate).astype(int) + lowerLimit = numpy.maximum(lowerLimit, 0) + upperLimit = numpy.minimum(upperLimit, totalSlots) + + # cut the limits to all the reachable charge values + for i in range(totalSlots - 1): + if lowerLimit[i+1] < lowerLimit[i]: + lowerLimit[i+1] = lowerLimit[i] + if upperLimit[i+1] < upperLimit[i]: + upperLimit[i+1] = upperLimit[i] + + for i in reversed(range(totalSlots - 1)): + if lowerLimit[i] < lowerLimit[i+1] - 1: + lowerLimit[i] = lowerLimit[i+1] - 1 + if upperLimit[i] < upperLimit[i+1] - 1: + upperLimit[i] = upperLimit[i+1] - 1 + + for i in range(totalSlots): + if lowerLimit[i] > i: + lowerLimit[i] = i + else: + break + + for i in range(totalSlots): + if upperLimit[i] > i + 1: + upperLimit[i] = i + 1 + else: + break + + # the profile of how the appliance will charge, 1 for every slot it will charge, 0 otherwise + chargingProfile = numpy.zeros(totalSlots) + + # starting point of the algorithm is 0 + lowerLimit = numpy.concatenate(([0], lowerLimit)) + upperLimit = numpy.concatenate(([0], upperLimit)) + + # the charging slots ordered from cheapest slot to most expensive slot + cheapestOrder = numpy.argsort(priceProfile) + for slot in cheapestOrder: + # if the appliance can be turned on at that slot, turn it on and update the limits + if lowerLimit[slot] < upperLimit[slot+1] and lowerLimit[slot] < lowerLimit[-1]: + chargingProfile[slot] = 1 + # update the limits to what's newly possible now + lowerSlot = (lowerLimit > lowerLimit[slot]).argmax() + upperSlot = (upperLimit == upperLimit[slot+1]).argmax() + lowerLimit[lowerSlot:] -= 1 + upperLimit[upperSlot:] -= 1 + + # save the new charge level for the appliance + self.memory.smart.currentCharge = startingCharge - dischargingSum[wantedSlots-1] + numpy.sum(chargingProfile[:wantedSlots]) * chargingRate + + # get the power profile for the charging interval and save it + powerProfile = chargingProfile[:wantedSlots] * self.chargingPower + self.smartDemand.set(fromDT, powerProfile) + + # calculates appliance power demand for a given time interval acting as if the accumulator wanted to stay as charged as possible + def calculateUncontrolledDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # simulates an uncontrolled charging algorithm, when an appliance wants to have as much energy stored as possible + + # how much energy goes into the appliance each minute it's turned on + chargingRate = self.chargingPower / 60 # kWh per minute + # how much energy leaves the appliance each minute of the interval + dischargingRates = self.memory.dischargingProfile.get(fromDT, toDT) / 60 # kWh per minute for each minute + + # never discharge past 0 and never charge over the capacity + upperLimit = self.capacity # kWh + + # get the current charge from memory + charge = self.memory.uncontrolled.currentCharge + + # the power profile for the interval + totalSlots = utils.minutesBetween(fromDT, toDT) + powerProfile = numpy.zeros(totalSlots) + + # simulate the progression of charge during the interval + for slot in range(totalSlots): + charge -= dischargingRates[slot] + if charge + chargingRate < upperLimit: + charge += chargingRate + powerProfile[slot] = self.chargingPower + + # save the new charge level to memory + self.memory.uncontrolled.currentCharge = charge + + # save the calculated demand + self.uncontrolledDemand.set(fromDT, powerProfile) + + # calculates appliance power demand for a given time interval acting as if the accumulator wanted to always charge completely and then discharge completely + def calculateSpreadOutDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # simulates a thermostat-based charging, when an appliance starts charging when it discharges past some threshhold, and stops charging when it's fully charged + + # how much energy goes into the appliance each minute it's turned on + chargingRate = self.chargingPower / 60 # kWh per minute + # how much energy leaves the appliance each minute of the interval + dischargingRates = self.memory.dischargingProfile.get(fromDT, toDT) / 60 # kWh per minute for each minute + + # never discharge past 0 and never charge over the capacity + lowerLimit = 0 # kWh + upperLimit = self.capacity # kWh + + # get the current charge from memory + charge = self.memory.spreadOut.currentCharge + # get if we're charging right now from memory + charging = self.memory.spreadOut.charging + + # the power profile for the interval + totalSlots = utils.minutesBetween(fromDT, toDT) + powerProfile = numpy.zeros(totalSlots) + # simulate the progression of charge during the interval + for slot in range(totalSlots): + charge -= dischargingRates[slot] + if charging: + if charge + chargingRate > upperLimit: + charging = False + else: + if charge <= lowerLimit: + charging = True + + if charging: + charge += chargingRate + powerProfile[slot] = self.chargingPower + + # save the new charge level and if it was charging and the end of the interval + self.memory.spreadOut.currentCharge = charge + self.memory.spreadOut.charging = charging + + # save the calculated demand + self.spreadOutDemand.set(fromDT, powerProfile) + +# abstract base class for machine-like household appliances (e.g. dishwasher, washing machine) +class Machine(Appliance, ABC): + # appliance usage statistics + usageStatistics: applianceStatistics.MachineStatistics + + # constructor, just prepares the variables + def __init__(self): + super().__init__() + # for each day keeps time the appliance should start after, time it should finish by and the usage profile of that run of the appliance + self.memory.usages = dict() + + # creates a random appliance with the right parameters + @classmethod + def random(cls): + return cls() + + # generates appliance usage for a given interval + def generateUsage(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day decide if the appliance will be used at all, + # and if so, generate the time the appliance should start after, time it should finish by + # and the usage profile of that run of the appliance + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + if date not in self.memory.usages: + if random.random() < self.usageStatistics.usageProbabilities[date]: + startAfter = self.usageStatistics.randomStartAfter() + finishBy = self.usageStatistics.randomFinishBy() + powerUsageProfile = self.usageStatistics.randomUsageProfile() + self.memory.usages[date] = ((startAfter, finishBy), powerUsageProfile) + else: + self.memory.usages[date] = None + + # calculates appliance power demand for a given time interval acting as if the appliance was NOT smart + def calculateSmartDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day in the interval, calculate the best time to start the appliance so that the run would be the cheapest + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + # the power profile of that day (plus some overlap) + powerProfile = numpy.zeros(2 * utils.minutesIn(oneDay)) + + # get the appliance usage for that day and act accordingly + usage = self.memory.usages[date] + if usage is not None: + # the price for the interval + priceProfile = self.priceProfile.get(midnight, midnight+2*oneDay) + + # the slots in which the appliance is available for being turned on + ((startAfter, finishBy), powerUsageProfile) = usage + startAfterSlot = startAfter.hour * 60 + startAfter.minute + finishBySlot = finishBy.hour * 60 + finishBy.minute + utils.minutesIn(oneDay) + + # the length of the run of the appliance + runtime = powerUsageProfile.size + cheapestSlot = startAfterSlot + + # if there is enough time to run the appliance, find the time to start at so that the run would be the cheapest + if finishBySlot - startAfterSlot > runtime: + # basically we have to try all the times between the starting and finishing slots to find the cheapest one + cheapestPrice = math.inf + for startingSlot in range(startAfterSlot, finishBySlot - runtime): + # calculate the price for if the appliance would be run starting at a given slot + slotPrice = numpy.dot(powerUsageProfile, priceProfile[startingSlot:startingSlot+runtime]) + # if it's better than what we found so far, save it + if slotPrice < cheapestPrice: + cheapestSlot = startingSlot + cheapestPrice = slotPrice + + # put the usage of the appliance at the right time in the power profile + powerProfile[cheapestSlot:cheapestSlot+runtime] = powerUsageProfile + + # save the power profile + self.smartDemand.add(midnight, powerProfile) + + # calculates appliance power demand for a given time interval acting as if the appliance wanted to be used as early as possible + def calculateUncontrolledDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day in the interval, just run the appliance as soon as possible + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + # the power profile of that day (plus some overlap) + powerProfile = numpy.zeros(2 * utils.minutesIn(oneDay)) + + # get the appliance usage for that day and act accordingly + usage = self.memory.usages[date] + if usage is not None: + ((startAfter, _), powerUsageProfile) = usage + + # the slots in which the appliance is available for being turned on + startAfterSlot = startAfter.hour * 60 + startAfter.minute + + # the length of the run of the appliance + runtime = powerUsageProfile.size + + # put the usage of the appliance at the right time in the power profile + powerProfile[startAfterSlot:startAfterSlot+runtime] = powerUsageProfile + + # save the power profile + self.uncontrolledDemand.add(midnight, powerProfile) + + # calculates appliance power demand for a given time interval acting as if the machine wanted to spread out its use across the whole possible interval + def calculateSpreadOutDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # for each day in the interval, just run the appliance in the middle of the available interval + for midnight in utils.midnightsBetween(fromDT, toDT): + date = midnight.date() + # the power profile of that day (plus some overlap) + powerProfile = numpy.zeros(2 * utils.minutesIn(oneDay)) + + # get the appliance usage for that day and act accordingly + usage = self.memory.usages[date] + if usage is not None: + # the slots in which the appliance is available for being turned on + ((startAfter, finishBy), powerUsageProfile) = usage + startAfterSlot = startAfter.hour * 60 + startAfter.minute + finishBySlot = finishBy.hour * 60 + finishBy.minute + utils.minutesIn(oneDay) + + # the length of the run of the appliance + runtime = powerUsageProfile.size + + # the slot in which the appliance should start + startingSlot = startAfterSlot + max(0, (finishBySlot - startAfterSlot - runtime) // 2) + + # put the usage of the appliance at the right time in the power profile + powerProfile[startingSlot:startingSlot+runtime] = powerUsageProfile + + # save the power profile + self.spreadOutDemand.add(midnight, powerProfile) + +# class representing an electric car +class Car(Battery): + usageStatistics = applianceStatistics.carStatistics[0] + @classmethod + def randomWithIndex(cls, index: int = 0): + usageStatistics = applianceStatistics.carStatistics[index] + chargingPower = usageStatistics.randomChargingPower() + car = cls(chargingPower=chargingPower) + car.usageStatistics = usageStatistics + return car + +# class representing air conditioning +class AirConditioning(Accumulator): + usageStatistics = applianceStatistics.airConditioningStatistics + +# class representing an electrical heating +class ElectricalHeating(Accumulator): + usageStatistics = applianceStatistics.electricalHeatingStatistics + +# class representing a refrigerator +class Fridge(Accumulator): + usageStatistics = applianceStatistics.fridgeStatistics + +# class representing a water heater +class WaterHeater(Accumulator): + usageStatistics = applianceStatistics.waterHeaterStatistics + +# class representing a dishwasher +class Dishwasher(Machine): + usageStatistics = applianceStatistics.dishwasherStatistics + +# class representing a washing machine +class WashingMachine(Machine): + usageStatistics = applianceStatistics.washingMachineStatistics diff --git a/simulator/simulator/applianceStatistics.py b/simulator/simulator/applianceStatistics.py new file mode 100755 index 0000000..2e9ff72 --- /dev/null +++ b/simulator/simulator/applianceStatistics.py @@ -0,0 +1,266 @@ +#!/usr/bin/env python3 + +import datetime +import json +import math +import random +import os + +from abc import ABC +from collections import defaultdict +from types import SimpleNamespace +from typing import Dict, List, Tuple + +import numpy +import pandas + +from .constants import oneDay +from .profile import Profile +from .utils import minutesIn + +thisDir = os.path.dirname(__file__) + +# base class for household appliance statistics +# contains parameters based on which random appliances and their random usages will be generated +class ApplianceStatistics(ABC): + pass + +# statistics of battery-based household appliances (e.g. electric car) +class BatteryStatistics(ApplianceStatistics): + # possible charging powers of the appliance type + chargingPowers: List[float] + # probabilities that the appliance will get used on a given day + usageProbabilities: Dict[datetime.date, float] + # possible charges which the appliance might need after usage on a given day + neededCharges: Dict[datetime.date, List[float]] + # average charge needed by the appliance type on a given day + averageNeededCharge: Dict[datetime.date, float] + # intervals in which the appliance might be used on a given day + usageIntervals: Dict[datetime.date, List[Tuple[datetime.time, datetime.time]]] + # profile of what ratio of appliances are available for charging at any given date and time + availabilityProfile: Profile + + # generates a random charging power for the appliance + def randomChargingPower(self) -> float: + return random.choice(self.chargingPowers) + + # generates a random needed charge for the appliance, if the appliance gets used at all + def randomNeededCharge(self, date: datetime.date) -> float: + if random.random() < self.usageProbabilities[date]: + return random.choice(self.neededCharges[date]) + else: + return 0.0 + + # generates a random usage interval for the appliance + def randomUsageInterval(self, date: datetime.date) -> Tuple[datetime.time, datetime.time]: + if len(self.usageIntervals[date]) == 0: + return (None, None) + else: + return random.choice(self.usageIntervals[date]) + + # loads charging powers from a file + # the file should have one possible charging power on each line + def loadChargingPowersFromFile(self, path): + self.chargingPowers = list(numpy.loadtxt(path)) + + # loads the usage probabilities from a CSV file + # in each row in the file there should be a date in the first column and the corresponding usage probability in the second column + def loadUsageProbabilitiesFromFile(self, path): + data = pandas.read_csv(path, parse_dates=[0]) + self.usageProbabilities = {} + for timestamp, usageprob in data.itertuples(index=False, name=None): + date = timestamp.date() + self.usageProbabilities[date] = usageprob + + # loads needed charges from a file + # on each line of the file there should be a date followed by a list of possible needed charges on that date + def loadNeededChargesFromFile(self, path): + self.neededCharges = {} + self.averageNeededCharge = {} + with open(path, "r") as chargesFile: + for line in chargesFile: + date = datetime.datetime.strptime(line.strip()[:10], "%Y-%m-%d").date() + charges = [float(x) for x in line.strip()[11:].strip(" []\n").split(",")] + self.neededCharges[date] = charges + self.averageNeededCharge[date] = numpy.mean(charges) * self.usageProbabilities[date] + + # loads the usage intervals from a file + # on each line of the file there should be a date followed by a list of possible usage intervals on that date + def loadUsageIntervalsFromFile(self, path): + def parseInterval(interval) -> Tuple[datetime.time, datetime.time]: + if len(interval) < 10: + return (None, None) + else: + start = datetime.datetime.strptime(interval[:5], "%H:%M").time() + end = datetime.datetime.strptime(interval[-5:], "%H:%M").time() + return (start, end) + + self.usageIntervals = defaultdict(list) + with open(path, "r") as usageIntervalsFile: + for line in usageIntervalsFile: + date = datetime.datetime.strptime(line.strip()[:10], "%Y-%m-%d").date() + intervals = list(map(parseInterval, line.strip()[11:].strip(" []\n").split(", "))) + self.usageIntervals[date].extend(intervals) + + # loads the availability profile from a CSV file + # in each row in the file there should be a date and time in the first column and the corresponding appliance availability in the second column + def loadAvailabilityProfileFromFile(self, path): + self.availabilityProfile = Profile.fromCSV(path) + +# statistics of accumulator-based household appliances (e.g. water heater, refrigerator) +class AccumulatorStatistics(ApplianceStatistics): + # list of possible charging powers for the appliance type + chargingPowers: List[float] + # average charging power of the appliance type + averageChargingPower: float + # mean appliance capacity and its standard deviation + capacityParameters: Tuple[float, float] + # the average discharging profile for the appliance type + dischargingProfile: Profile + # average needed charge for the appliance type for each day + averageDailyCharge: Dict[datetime.date, float] + # mean and standard deviation of the scale of discharging profiles of the appliance type + dischargingProfileScaleParameters: Tuple[float, float] + + # generates a random charging power for the appliance + def randomChargingPower(self) -> float: + return random.choice(self.chargingPowers) + + # generates a random capacity power for the appliance + def randomCapacity(self) -> float: + return random.gauss(*self.capacityParameters) + + # generates a random discharging profile scale + def randomDischargingProfileScale(self) -> float: + return random.gauss(*self.dischargingProfileScaleParameters) + + # loads charging powers from a file + # the file should have one possible charging power on each line + def loadChargingPowersFromFile(self, path: str): + self.chargingPowers = list(numpy.loadtxt(path)) + self.averageChargingPower = numpy.mean(self.chargingPowers) + + # loads the average discharging profile from a CSV file + # in each row in the file there should be a date and time in the first column and the corresponding average discharging power in the second column + def loadDischargingProfileFromFile(self, path: str): + self.dischargingProfile = Profile.fromCSV(path) + self.averageDailyCharge = {d: s * 24 for (d, s) in self.dischargingProfile.dailyAverages().items()} + +# statistics for machine-like household appliances (e.g. dishwasher, washing machine) +class MachineStatistics(ApplianceStatistics): + # mean minute of the day that the machine can start after, and its standard deviation + startAfterParameters: Tuple[int, int] + # mean minute of the day that the machine must finish by, and its standard deviation + finishByParameters: Tuple[int, int] + # probabilities that the machine will get used on a given day + usageProbabilities: Dict[datetime.date, float] + # possible power usage profiles of the machine type + usageProfiles: List[numpy.ndarray] + # average power needed by the machine type on a given day + averagePowerNeeded: Dict[datetime.date, float] + + # generates a random starting time + def randomStartAfter(self) -> datetime.time: + minutes = min(max(0, math.floor(random.gauss(*self.startAfterParameters))), minutesIn(oneDay)-1) + return datetime.time(minutes // 60, minutes % 60) + + # generates a random finishing time + def randomFinishBy(self) -> datetime.time: + minutes = min(max(0, math.floor(random.gauss(*self.finishByParameters))), minutesIn(oneDay)-1) + return datetime.time(minutes // 60, minutes % 60) + + # picks a random choosing profile of all the possible profiles + def randomUsageProfile(self) -> numpy.ndarray: + return random.choice(self.usageProfiles) + + # loads the usage probabilities from a CSV file + # in each row in the file there should be a date in the first column and the corresponding usage probability in the second column + def loadUsageProbabilitiesFromFile(self, path: str): + data = pandas.read_csv(path, parse_dates=[0]) + self.usageProbabilities = {} + for timestamp, usageprob in data.itertuples(index=False, name=None): + date = timestamp.date() + self.usageProbabilities[date] = usageprob + + # loads the power usage profiles from a file + # on each line of the file there should be a list of power draws for each minute of the operation of the machine + def loadUsageProfilesFromFile(self, path: str): + self.usageProfiles = [] + sums = [] + with open(path, "r") as sourceFile: + for line in sourceFile: + profile = numpy.fromstring(line, dtype=float, sep=",") + self.usageProfiles.append(profile) + sums.append(numpy.sum(profile)/60) + + averagePowerNeeded = numpy.mean(sums) + + self.averagePowerNeeded = {} + for date, usageProbability in self.usageProbabilities.items(): + self.averagePowerNeeded[date] = averagePowerNeeded * usageProbability + +# ownership ratios of various household appliances and vehicles +# adapted from https://www.eia.gov/consumption/residential/reports/2009/state_briefs/pdf/tx.pdf +with open(f"{thisDir}/data/manual/ownershipRatios.json", "r") as ownershipRatiosFile: + ownershipRatios = SimpleNamespace(**json.load(ownershipRatiosFile)) + +# probabilities that a house owns some number of cars +carCountProbabilities = list(pandas.read_csv(f"{thisDir}/data/nhts/cars/ownershipRatios.csv")["ratio"]) +atLeastThisManyCarsProbability = [sum(carCountProbabilities[i:]) for i in range(len(carCountProbabilities))] + +# electric car statistics +# each car in a household get used differently +carStatistics = [] +for car in range(4): + carStatistics.append(BatteryStatistics()) + carStatistics[car].loadUsageProbabilitiesFromFile(f"{thisDir}/data/nhts/cars/car{car+1}/usageRatios.csv") + carStatistics[car].loadUsageIntervalsFromFile(f"{thisDir}/data/nhts/cars/car{car+1}/trips.txt") + carStatistics[car].loadAvailabilityProfileFromFile(f"{thisDir}/data/nhts/cars/car{car+1}/availability.csv") + carStatistics[car].loadNeededChargesFromFile(f"{thisDir}/data/dataport/cars/charges.txt") + carStatistics[car].loadChargingPowersFromFile(f"{thisDir}/data/dataport/cars/maxPowers.txt") + +# load the accumulator capacities statistics from a json file +with open(f"{thisDir}/data/manual/applianceCapacities.json", "r") as capacitiesFile: + capacities = json.load(capacitiesFile, object_hook=lambda dct: SimpleNamespace(**dct)) + +# air conditioning statistics +airConditioningStatistics = AccumulatorStatistics() +airConditioningStatistics.capacityParameters = (capacities.airConditioning.mean, capacities.airConditioning.std) +airConditioningStatistics.dischargingProfileScaleParameters = (1, 0.3) +airConditioningStatistics.loadChargingPowersFromFile(f"{thisDir}/data/dataport/accumulators/airconditioning/maxPowers.txt") +airConditioningStatistics.loadDischargingProfileFromFile(f"{thisDir}/data/dataport/accumulators/airconditioning/averageUsage.csv") + +# electrical heating statistics +electricalHeatingStatistics = AccumulatorStatistics() +electricalHeatingStatistics.capacityParameters = (capacities.electricalHeating.mean, capacities.electricalHeating.std) +electricalHeatingStatistics.dischargingProfileScaleParameters = (1, 0.3) +electricalHeatingStatistics.loadChargingPowersFromFile(f"{thisDir}/data/dataport/accumulators/electricalheating/maxPowers.txt") +electricalHeatingStatistics.loadDischargingProfileFromFile(f"{thisDir}/data/dataport/accumulators/electricalheating/averageUsage.csv") + +# refrigerator statistics +fridgeStatistics = AccumulatorStatistics() +fridgeStatistics.capacityParameters = (capacities.fridge.mean, capacities.fridge.std) +fridgeStatistics.dischargingProfileScaleParameters = (1, 0.3) +fridgeStatistics.loadChargingPowersFromFile(f"{thisDir}/data/dataport/accumulators/fridge/maxPowers.txt") +fridgeStatistics.loadDischargingProfileFromFile(f"{thisDir}/data/dataport/accumulators/fridge/averageUsage.csv") + +# water heater statistics +waterHeaterStatistics = AccumulatorStatistics() +waterHeaterStatistics.capacityParameters = (capacities.waterHeater.mean, capacities.waterHeater.std) +waterHeaterStatistics.dischargingProfileScaleParameters = (1, 0.3) +waterHeaterStatistics.loadChargingPowersFromFile(f"{thisDir}/data/dataport/accumulators/waterheater/maxPowers.txt") +waterHeaterStatistics.loadDischargingProfileFromFile(f"{thisDir}/data/dataport/accumulators/waterheater/averageUsage.csv") + +# dishwasher statistics +dishwasherStatistics = MachineStatistics() +dishwasherStatistics.startAfterParameters = (21*60, 60) +dishwasherStatistics.finishByParameters = (5*60, 60) +dishwasherStatistics.loadUsageProbabilitiesFromFile(f"{thisDir}/data/dataport/machines/dishwasher/usages.csv") +dishwasherStatistics.loadUsageProfilesFromFile(f"{thisDir}/data/dataport/machines/dishwasher/profiles.txt") + +# washing machine statistics +washingMachineStatistics = MachineStatistics() +washingMachineStatistics.startAfterParameters = (21*60, 60) +washingMachineStatistics.finishByParameters = (5*60, 60) +washingMachineStatistics.loadUsageProbabilitiesFromFile(f"{thisDir}/data/dataport/machines/washingmachine/usages.csv") +washingMachineStatistics.loadUsageProfilesFromFile(f"{thisDir}/data/dataport/machines/washingmachine/profiles.txt") diff --git a/simulator/simulator/connection.py b/simulator/simulator/connection.py new file mode 100644 index 0000000..6076444 --- /dev/null +++ b/simulator/simulator/connection.py @@ -0,0 +1,150 @@ +#!/usr/bin/env python3 + +import datetime + +import numpy + +from . import utils +from . import priceConfig + +from .constants import oneDay +from .house import House +from .profile import Profile + +# class for the connections between the smart grid and houses +class Connection: + # current date and time in the simulation + currentDT: datetime.datetime + # the connected house + house: House + # probabilities with which the electricity price should be lower in a given minute + cheaperPriceRatioProfile: Profile + # times where the prices will be cheaper + cheaperMinutesProfile: Profile + # the actual price profile for the connected house + priceProfile: Profile + + # constructor, just prepares the variables + def __init__(self, house: House): + self.house = house + self.cheaperPriceRatioProfile = Profile() + self.cheaperMinutesProfile = Profile() + self.priceProfile = Profile() + + # sets everything up before the start of the simulation + def setUp(self, dt: datetime.datetime): + self.currentDT = dt + # set up the house as well + self.house.setUp(dt=dt) + + # generate cheaper price intervals and the electricity price profile ahead enough in the future + self.generateRandomCheaperIntervals(dt-1*oneDay, dt+2*oneDay) + self.generatePriceProfile(dt-1*oneDay, dt+2*oneDay) + + # send the generated price profile to the connected house + self.sendPriceProfile(dt-1*oneDay, dt+2*oneDay) + + # moves ahead one day and does all the calculations that need to be done in that day + def tick(self): + # remove past, unneeded values from the profiles to free up some memory + self.cheaperPriceRatioProfile.prune(self.currentDT - oneDay) + self.cheaperMinutesProfile.prune(self.currentDT - oneDay) + self.priceProfile.prune(self.currentDT - oneDay) + + # move one day ahead + self.currentDT += oneDay + cdt = self.currentDT + + # generate cheaper price intervals and the electricity price profile for one more day + self.generateRandomCheaperIntervals(cdt+oneDay, cdt+2*oneDay) + self.generatePriceProfile(cdt+oneDay, cdt+2*oneDay) + + # send the generated price profile to the connected house + self.sendPriceProfile(cdt+oneDay, cdt+2*oneDay) + + # generates intervals of cheaper prices based on the cheap price probabilities provided by the smart grid + def generateRandomCheaperIntervals(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # how many minutes continuously have to be cheap (60 minutes by default) + cheapIntervalLength = priceConfig.cheapIntervalLength + # at least how many minutes in total have to be cheap per day + cheapMinutesTotal = priceConfig.cheapMinutesCount + + # if it is configured to have zero minutes cheap, we act as if there was no lower or upper limit on the minutes + if cheapMinutesTotal == 0: + # get probabilities for the connected house having cheaper electricity at a given minute + probs = self.cheaperPriceRatioProfile.get(fromDT, toDT) + # generate the positions of the cheaper minutes and save them + cheaperMinutes = numpy.random.random(probs.size) < probs + self.cheaperMinutesProfile.add(fromDT, cheaperMinutes) + + # otherwise we guarantee to have the specified amount of minutes cheap + else: + # generate cheaper intervals for each day that starts in the given interval + for midnight in utils.midnightsBetween(fromDT, toDT): + # start with all prices expensive + cheaperIntervals = numpy.zeros(utils.minutesIn(oneDay)+2*cheapIntervalLength) + + # get probabilities for a cheap interval being positioned in a given spot + # we can't change prices for previous days already broadcasted to houses, so we need to make sure the intervals we set will start after midnight of the previous day + shift = datetime.timedelta(minutes=cheapIntervalLength) + probs = self.cheaperPriceRatioProfile.get(midnight+shift, midnight+shift+oneDay) + + # this distinction is just to optimize the calculation, the results are the same + if cheapIntervalLength == 1: + # generate enough cheap minutes in the given day + cheapMinutePositions = utils.randomWithRelativeProbs(probs, count=cheapMinutesTotal) + cheaperIntervals[cheapMinutePositions] = 1 + else: + # add cheap intervals until there are enough cheap minutes in the given day + # some intervals might overlap and that is okay + while numpy.sum(cheaperIntervals) < cheapMinutesTotal: + # put a cheap interval in a random position + cheapIntervalStart = utils.minutesIn(shift) + utils.randomWithRelativeProbs(probs) - cheapIntervalLength//2 + cheaperIntervals[cheapIntervalStart:cheapIntervalStart+cheapIntervalLength] = 1 + + # save the cheaper intervals + self.cheaperMinutesProfile.add(midnight, cheaperIntervals) + + # generates a price profile based on the cheaper interval locations calculated earlier + def generatePriceProfile(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # the electricity prices + lowerPrice = priceConfig.lowerPrice + higherPrice = priceConfig.higherPrice + + # get where the price should be cheaper + cheapIntervals = numpy.minimum(self.cheaperMinutesProfile.get(fromDT, toDT), 1) + + # set the prices accordingly + prices = numpy.full(cheapIntervals.size, fill_value=higherPrice, dtype=float) + prices[cheapIntervals.astype(bool)] = lowerPrice + + # add a bit of randomness so appliances don't always choose the earliest possible cheap location + prices += numpy.random.random(prices.size) * 0.01 + + # save the price profile + self.priceProfile.set(fromDT, prices) + + # sends the price profile to the connected house + def sendPriceProfile(self, fromDT: datetime.datetime, toDT: datetime.datetime): + prices = self.priceProfile.get(fromDT, toDT) + self.house.setPriceProfile(fromDT, prices) + + # sets the probabilities that the electricity will be cheaper in a given minute + # called by the grid + def setPriceRatio(self, fromDT: datetime.datetime, priceRatio: numpy.ndarray): + self.cheaperPriceRatioProfile.set(fromDT, priceRatio) + + # collect the electricity demand the house would have if it was using smart appliances + # called by the grid + def getSmartDemand(self, fromDT: datetime.datetime = None, toDT: datetime.datetime = None): + return self.house.getSmartDemand(fromDT, toDT) + + # collect the electricity demand the house would have if its appliances would charge as early as possible + # called by the grid + def getUncontrolledDemand(self, fromDT: datetime.datetime = None, toDT: datetime.datetime = None): + return self.house.getUncontrolledDemand(fromDT, toDT) + + # collect the electricity demand the house would have if its appliances would charge evenly over their use period + # called by the grid + def getSpreadOutDemand(self, fromDT: datetime.datetime = None, toDT: datetime.datetime = None): + return self.house.getSpreadOutDemand(fromDT, toDT) diff --git a/simulator/simulator/constants.py b/simulator/simulator/constants.py new file mode 100755 index 0000000..ae0adb0 --- /dev/null +++ b/simulator/simulator/constants.py @@ -0,0 +1,7 @@ +#!/usr/bin/env python3 + +import datetime + +# frequently used constants in the simulator + +oneDay = datetime.timedelta(days=1) diff --git a/simulator/simulator/data/config.txt b/simulator/simulator/data/config.txt new file mode 100644 index 0000000..9d7a5c7 --- /dev/null +++ b/simulator/simulator/data/config.txt @@ -0,0 +1,2 @@ +fromdate=yyyy-mm-dd +todate=yyyy-mm-dd diff --git a/simulator/simulator/data/dataport/KEY b/simulator/simulator/data/dataport/KEY new file mode 100644 index 0000000..4b42d3d --- /dev/null +++ b/simulator/simulator/data/dataport/KEY @@ -0,0 +1,2 @@ +username=yourDataportUsername +password=yourDataportPassword diff --git a/simulator/simulator/data/dataport/accumulators/getData.py b/simulator/simulator/data/dataport/accumulators/getData.py new file mode 100755 index 0000000..4097934 --- /dev/null +++ b/simulator/simulator/data/dataport/accumulators/getData.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 + +import os +import sys + +import numpy +import pandas +import sqlalchemy + +fromDate = sys.argv[1] +toDate = sys.argv[2] +username = sys.argv[3] +password = sys.argv[4] + +# create the postgresql connection engine +engine = sqlalchemy.create_engine(f"postgresql://{username}:{password}@dataport.pecanstreet.org:5434/dataport") + +# download the data for each appliance type +for columnID, outFolder in [("refrigerator1", "fridge"), ("furnace1", "electricalheating"), ("air1", "airconditioning"), ("waterheater1", "waterheater")]: + print(f"Downloading {outFolder} data...") + os.makedirs(outFolder, exist_ok=True) + + badIDs = [] + # these households already have a grid-controlled water heater, which would pollute our data, so we filter them out + if columnID == "waterheater1": + badIDs = [2171, 2204, 5357, 7937, 9356, 9934] + + # these households have bad data for their electrical heating, which causes errors in the simulation + if columnID == "furnace1": + badIDs = [6730, 7536, 7901, 9019] + + # these households have bad data for their air conditioning, which causes errors in the simulation + if columnID == "air1": + badIDs = [6730] + + # these households have bad data for their fridge, which causes errors in the simulation + if columnID == "refrigerator1": + badIDs = [7536, 7901, 9019] + + # download the household IDs which have some valid data + query = f''' + SELECT DISTINCT dataid + FROM electricity.eg_realpower_1min + WHERE {columnID} IS NOT NULL + AND {columnID} > 0.05 + AND localminute >= '{fromDate}' + AND localminute < '{toDate}' + ''' + dataIDs = pandas.read_sql_query(query, con=engine) + validIDs = [id for id in dataIDs["dataid"] if id not in badIDs] + idFilter = ', '.join(map(str, validIDs)) + + # download the maximum powers the appliances were drawing for each household + query = f''' + SELECT MAX({columnID}) as {columnID} + FROM electricity.eg_realpower_1min + WHERE dataid IN ({idFilter}) + AND localminute >= '{fromDate}' + AND localminute < '{toDate}' + GROUP BY dataid + ''' + maxs = pandas.read_sql_query(query, con=engine) + + # filter out unreal data (none of these appliances can realistically draw over 20kW) + maxPowers = maxs[columnID].values + maxPowers = maxPowers[maxPowers < 20] + + # download the average power draw by the appliances for each minute + query = f''' + SELECT localminute at time zone 'CST' as localminute, AVG({columnID}) as {columnID} + FROM electricity.eg_realpower_1min + WHERE dataid IN ({idFilter}) + AND localminute >= '{fromDate}' + AND localminute < '{toDate}' + GROUP BY localminute + ORDER BY localminute + ''' + avgs = pandas.read_sql_query(query, con=engine, parse_dates=["localminute"]) + + # save the downloaded data + print("Parsing...") + numpy.savetxt(f"{outFolder}/maxPowers.txt", maxPowers, fmt="%.5f") + + # for some appliances there are not many households which have that appliance, so we smooth the data out a bit + usages = avgs.set_index("localminute").resample('min').interpolate('time').rolling(window=120, center=True, min_periods=1).mean() + usages.to_csv(f"{outFolder}/averageUsage.csv", index_label="datetime", header=["usage"], float_format="%.5f") + + print("Done.") + +# dispose of the engine (sometimes the connection keeps hanging otherwise, even after the script has finished) +engine.dispose() diff --git a/simulator/simulator/data/dataport/cars/getData.py b/simulator/simulator/data/dataport/cars/getData.py new file mode 100755 index 0000000..58164b4 --- /dev/null +++ b/simulator/simulator/data/dataport/cars/getData.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 + +import sys + +from collections import defaultdict +from typing import List + +import numpy +import pandas +import sqlalchemy + +fromDate = sys.argv[1] +toDate = sys.argv[2] +username = sys.argv[3] +password = sys.argv[4] + +print("Downloading dataport car data...", flush=True) +# create the postgresql connection engine +engine = sqlalchemy.create_engine(f"postgresql://{username}:{password}@dataport.pecanstreet.org:5434/dataport") + +# download the data power draw data from the database +# we only need the sums for each charging interval, so an 1-hour resolution is enough +query = f''' + SELECT dataid, local_15min at time zone 'CST' as local_15min, car1 + FROM electricity.eg_realpower_15min + WHERE local_15min >= '{fromDate}' + AND local_15min < '{toDate}' + AND "dataid" = ANY( + SELECT DISTINCT dataid + FROM electricity.eg_realpower_15min + WHERE car1 IS NOT NULL + AND car1 > 0.1 + AND local_15min >= '{fromDate}' + AND local_15min < '{toDate}' + ) + ORDER BY dataid, local_15min +''' +data = pandas.read_sql_query(query, con=engine, parse_dates=["local_15min"]).fillna(0) + +# dispose of the engine (sometimes the connection keeps hanging otherwise, even after the script has finished) +engine.dispose() + +# parse the downloaded data +print("Parsing...") +groups = data.groupby("dataid") + +# for each car get the continuous charging intervals +# add up how much electricity did the charging consume in total in each interval +# the data is sampled in kW in 15 minute intervals, to convert it to kWh, we need to divide it by four +def getSums(df) -> List[float]: + sums = [] + startingDate = None + workingSum = 0 + for localhour, carDraw in df.itertuples(index=False): + if carDraw < 0.1: + if workingSum > 0.5: + sums.append((startingDate, workingSum)) + workingSum = 0 + startingDate = None + else: + if startingDate is None: + startingDate = localhour.date() + workingSum += carDraw / 4.0 + return sums + +# group the charging sums by date into a dictionary +charges = defaultdict(list) +for chargeSums in groups.apply(getSums).values: + for date, charge in sorted(list(chargeSums)): + charges[date].append(charge) + +# write the charging sums out into a file, each day on one line +with open("charges.txt", "w") as chargesFile: + for date, chargesOnDate in sorted(charges.items()): + chargesFile.write(f"{date}: {chargesOnDate}\n") + +# get the maximum charging power the charger was drawing for each car +maxChargingPowers = groups["car1"].max().values + +# save it into a file +numpy.savetxt("maxPowers.txt", maxChargingPowers, fmt="%.5f") + +print("Done.") diff --git a/simulator/simulator/data/dataport/download.sh b/simulator/simulator/data/dataport/download.sh new file mode 100755 index 0000000..d33fc60 --- /dev/null +++ b/simulator/simulator/data/dataport/download.sh @@ -0,0 +1,19 @@ +#!/usr/bin/env bash + +export PYTHONUNBUFFERED=1 + +while read -r line || [ -n "$line" ]; do + declare "$line"; +done < KEY + +DIRS=(accumulators machines cars household ercot) + +echo "Retreiving dataport data..." + +for dir in ${DIRS[@]}; do + pushd $dir > /dev/null + ./getData.py $1 $2 $username $password | stdbuf -oL sed 's/^/ /' + popd > /dev/null +done + +echo "Dataport data retrieved." diff --git a/simulator/simulator/data/dataport/ercot/getData.py b/simulator/simulator/data/dataport/ercot/getData.py new file mode 100755 index 0000000..c9f9b3c --- /dev/null +++ b/simulator/simulator/data/dataport/ercot/getData.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 + +import math +import os +import sys + +import pandas +import sqlalchemy + +fromDate = sys.argv[1] +toDate = sys.argv[2] +username = sys.argv[3] +password = sys.argv[4] + +print("Downloading ERCOT actual demand data...") + +# create the postgresql connection engine +engine = sqlalchemy.create_engine(f"postgresql://{username}:{password}@dataport.pecanstreet.org:5434/dataport") + +# download the actual power demand information +query = f''' + SELECT delivery_date at time zone 'CST' as deliverydate, 1000*min(demand) as demand + FROM ercot.system_wide_demand + WHERE delivery_date >= '{fromDate}' + AND delivery_date < '{toDate}' + GROUP BY delivery_date +''' +data = pandas.read_sql_query(query, con=engine, parse_dates=["deliverydate"]) + +print("Parsing...") +# resample the time series to minutes and interpolate the missing values +data = data.set_index(["deliverydate"]).resample('min').interpolate('cubic') + +# save the resampled demand to a CSV file +os.makedirs("actual", exist_ok=True) +data.to_csv("actual/systemLoad.csv") + +print("Done.") + +# download the power demand predictions +print("Downloading ERCOT predicted demand data...") +query = f''' + SELECT report_date at time zone 'CST' as reportdate, delivery_date at time zone 'CST' as deliverydate, 1000*system_total as demand + FROM ercot.load_forecast_xforecast_zone_7day + WHERE delivery_date >= '{fromDate}' + AND delivery_date < '{toDate}' +''' +data = pandas.read_sql_query(query, con=engine, parse_dates=["reportdate", "deliverydate"]) + +# dispose of the engine (sometimes the connection keeps hanging otherwise, even after the script has finished) +engine.dispose() + +print("Parsing...") + +os.makedirs("predictions", exist_ok=True) + +# from all the predictions that were made for a date, get the one which is the right time ahead +def getPrediction(group, td) -> float: + closestPrediction = 0.0 + closestPredictionDistance = math.inf + for _, row in group.iterrows(): + predictionDistance = (row["deliverydate"] - row["reportdate"]).total_seconds() // 3600 + if abs(predictionDistance - td) < abs(closestPredictionDistance - td): + closestPrediction = row["demand"] + closestPredictionDistance = predictionDistance + return closestPrediction + +# for each date that was predicted, get the predictions which were exactly 0, 1, 2, 3, 4 or 5 days ahead +groupby = data.groupby("deliverydate") +for hoursAhead in [0, 24, 48, 72, 96, 120]: + # get the predictions, smooth them out a bit, resample the resulting time series to minutes and interpolate the missing values + predictions = groupby.apply(getPrediction, hoursAhead).rolling(window=3, center=True, min_periods=1).mean().resample('min').interpolate('cubic') + # save the predictions to a CSV file + predictions.to_csv(f"predictions/{hoursAhead}.csv", header=["demand"]) + +print("Done.") diff --git a/simulator/simulator/data/dataport/household/getData.py b/simulator/simulator/data/dataport/household/getData.py new file mode 100755 index 0000000..169abda --- /dev/null +++ b/simulator/simulator/data/dataport/household/getData.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python3 + +import sys + +import pandas +import sqlalchemy + +fromDate = sys.argv[1] +toDate = sys.argv[2] +username = sys.argv[3] +password = sys.argv[4] + +print(f"Downloading household data...") + +# create the postgresql connection engine +engine = sqlalchemy.create_engine(f"postgresql://{username}:{password}@dataport.pecanstreet.org:5434/dataport") + +# download the household IDs which have some valid data +query = f''' + SELECT DISTINCT dataid + FROM electricity.eg_realpower_15min + WHERE grid IS NOT NULL + AND grid > 0.05 + AND local_15min >= '{fromDate}' + AND local_15min < '{toDate}' +''' +dataIDs = pandas.read_sql_query(query, con=engine).fillna(0) + +# this household has weird, invalid data around June 20, 2018, so we throw it out +badIDs = [6730] +validIDs = [id for id in dataIDs["dataid"] if id not in badIDs] +idFilter = ', '.join(map(str, validIDs)) + +# download the average power draw of the selected households for each minute in the interval +query = f''' + SELECT localminute at time zone 'CST' as localminute, AVG(COALESCE(solar, 0) + COALESCE(solar2, 0) + COALESCE(grid, 0)) as use + FROM electricity.eg_realpower_1min + WHERE dataid IN ({idFilter}) + AND localminute >= '{fromDate}' + AND localminute < '{toDate}' + GROUP BY localminute + ORDER BY localminute +''' +data = pandas.read_sql_query(query, con=engine, parse_dates=[0]).fillna(0) + +# smooth out the data a bit, interpolate the power draw and save it to CSV +averageDraw = data.set_index("localminute").resample('min').interpolate('time').rolling(window=120, center=True, min_periods=1).mean() +averageDraw.to_csv(f"averageDraw.csv", index_label="datetime", header=["draw"], float_format="%.5f") + +print("Done.") + +# dispose of the engine (sometimes the connection keeps hanging otherwise, even after the script has finished) +engine.dispose() diff --git a/simulator/simulator/data/dataport/machines/getData.py b/simulator/simulator/data/dataport/machines/getData.py new file mode 100755 index 0000000..6d247b5 --- /dev/null +++ b/simulator/simulator/data/dataport/machines/getData.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python3 + +import datetime +import os +import sys + +from collections import defaultdict +from typing import List, Tuple + +import pandas +import sqlalchemy + +fromDate = sys.argv[1] +toDate = sys.argv[2] +username = sys.argv[3] +password = sys.argv[4] + +# helper function to split an iterable into chunks of a given size +def batch(iterable, n=1): + l = len(iterable) + for pos in range(0, l, n): + yield iterable[pos:min(pos+n, l)] + +# create the postgresql connection engine +engine = sqlalchemy.create_engine(f"postgresql://{username}:{password}@dataport.pecanstreet.org:5434/dataport") + +# download the data for each appliance type +for columnID, outFolder in [("dishwasher1", "dishwasher"), ("clotheswasher1", "washingmachine")]: + print(f"Downloading and parsing {outFolder} data...") + allUsages = [] + + # download the household IDs which have some valid data + query = f''' + SELECT DISTINCT dataid + FROM electricity.eg_realpower_1min + WHERE {columnID} IS NOT NULL + AND {columnID} > 0.05 + AND localminute >= '{fromDate}' + AND localminute < '{toDate}' + ''' + dataIDs = list(pandas.read_sql_query(query, con=engine)["dataid"]) + + # split the household ids into batches of 10 ids + # takes a little longer but uses a lot less memory + idBatches = list(batch(dataIDs, 10)) + for idBatch in batch(dataIDs, 10): + # download the power draw from the machines for each household in the batch + idFilter = ', '.join(map(str, idBatch)) + query = f''' + SELECT dataid, localminute at time zone 'CST' as localminute, {columnID} as machinedemand + FROM electricity.eg_realpower_1min + WHERE dataid IN ({idFilter}) + AND localminute >= '{fromDate}' + AND localminute < '{toDate}' + ''' + data = pandas.read_sql_query(query, con=engine, parse_dates=["localminute"]).fillna(0) + + # for each machine, identify the intervals it was probably used and generate power profiles from those + def getUsages(df: pandas.DataFrame) -> List[Tuple[datetime.datetime, List[float]]]: + machineUsages = [] + startTime = None + workingProfile = [] + smallValues = 0 + for time, draw in df.sort_values(by="localminute")[["localminute", "machinedemand"]].values: + # if it's drawing less than 30 watts, then it's possibly done + if draw <= 0.030: + # if it's drawing less than 30 watts for more than 5 minutes, then it's probably done + if smallValues > 5: + # if it drew some power for more than 20 minutes and less than 4 hours, we call it a valid run and save the profile + if len(workingProfile) > 20 and len(workingProfile) < 240 and sum(workingProfile) > 1: + machineUsages.append((startTime, workingProfile[:-smallValues])) + startTime = None + workingProfile = [] + smallValues = 0 + # otherwise it was probably a momentary pause in the machine's program, and we continue saving the draw + else: + if len(workingProfile) > 0: + workingProfile.append(draw) + smallValues += 1 + # otherwise it's still running, we append the current draw to the power profile + else: + if startTime is None: + startTime = pandas.to_datetime(time) + workingProfile.append(draw) + smallValues = 0 + + # if the last item was not zero we do not save the profile, as it may be from an incomplete cycle + return machineUsages + + # get the power usage profiles for each machine + batchUsages = data.groupby("dataid").apply(getUsages).values + + # aggregate the downloaded usages in a list + for usages in batchUsages: + allUsages.extend(usages) + + # save the power profiles to a file, and calculate the probability that the machine will be turned on at night for each day and save that to a CSV file + os.makedirs(outFolder, exist_ok=True) + with open(f"{outFolder}/profiles.txt", "w") as profileFile, open(f"{outFolder}/usages.csv", "w") as usageRatioFile: + nightRuns = defaultdict(int) + for start, profile in allUsages: + profileFile.write(", ".join(["{:.3f}".format(draw) for draw in profile])) + profileFile.write("\n") + + # if it ran between 20:00 and 6:00 then we call it a night run + date = start.date() + if start.hour >= 20: + nightRuns[date] += 1 + if start.hour <= 5: + nightRuns[date - datetime.timedelta(days=1)] += 1 + + # save the probability of a night run for each day to a CSV file + usageRatioFile.write("date,usageRatio\n") + currentDT = datetime.datetime.strptime(fromDate, "%Y-%m-%d") + while currentDT < datetime.datetime.strptime(toDate, "%Y-%m-%d"): + date = currentDT.date() + ratio = nightRuns[date]/len(allUsages) + usageRatioFile.write(f"{date},{ratio:.5f}\n") + currentDT += datetime.timedelta(days=1) + + print("Done.") + +# dispose of the engine (sometimes the connection keeps hanging otherwise, even after the script has finished) +engine.dispose() diff --git a/simulator/simulator/data/download.sh b/simulator/simulator/data/download.sh new file mode 100755 index 0000000..bb20a72 --- /dev/null +++ b/simulator/simulator/data/download.sh @@ -0,0 +1,21 @@ +#!/usr/bin/env bash + +pushd "$(dirname ${BASH_SOURCE[0]})" > /dev/null + +while read -r line || [ -n "$line" ]; do + declare "$line"; +done < config.txt + +DIRS=(dataport nhts) + +echo "Retreiving data..." + +for dir in ${DIRS[@]}; do + pushd $dir > /dev/null + ./download.sh $fromdate $todate | stdbuf -oL sed 's/^/ /' + popd > /dev/null +done + +echo "Data retrieved." + +popd > /dev/null \ No newline at end of file diff --git a/simulator/simulator/data/manual/applianceCapacities.json b/simulator/simulator/data/manual/applianceCapacities.json new file mode 100644 index 0000000..81c8f22 --- /dev/null +++ b/simulator/simulator/data/manual/applianceCapacities.json @@ -0,0 +1,18 @@ +{ + "airConditioning": { + "mean": 0.3, + "std": 0.1 + }, + "electricalHeating": { + "mean": 3.0, + "std": 0.5 + }, + "fridge": { + "mean": 0.3, + "std": 0.1 + }, + "waterHeater": { + "mean": 2.0, + "std": 0.7 + } +} diff --git a/simulator/simulator/data/manual/ownershipRatios.json b/simulator/simulator/data/manual/ownershipRatios.json new file mode 100644 index 0000000..ca59ad1 --- /dev/null +++ b/simulator/simulator/data/manual/ownershipRatios.json @@ -0,0 +1,8 @@ +{ + "airConditioning": 0.9, + "electricalHeating": 0.5, + "fridge": 0.9, + "waterHeater": 0.9, + "washingMachine": 0.8, + "dishwasher": 0.7 +} \ No newline at end of file diff --git a/simulator/simulator/data/manual/priceConfig.json b/simulator/simulator/data/manual/priceConfig.json new file mode 100644 index 0000000..1129fa7 --- /dev/null +++ b/simulator/simulator/data/manual/priceConfig.json @@ -0,0 +1,6 @@ +{ + "cheapIntervalLength": 60, + "cheapMinutesCount": 480, + "lowerPrice": 1.0, + "higherPrice": 3.0 +} diff --git a/simulator/simulator/data/nhts/cars/getData.py b/simulator/simulator/data/nhts/cars/getData.py new file mode 100755 index 0000000..361ea17 --- /dev/null +++ b/simulator/simulator/data/nhts/cars/getData.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python3 + +import datetime +import os +import sys + +from collections import defaultdict +from io import BytesIO +from types import SimpleNamespace +from zipfile import ZipFile + +import numpy +import pandas +import requests + +fromDT = datetime.datetime.strptime(sys.argv[1], "%Y-%m-%d") +toDT = datetime.datetime.strptime(sys.argv[2], "%Y-%m-%d") + +# download the NHTS data +print("Downloading NHTS car trip data...") +r = requests.get("https://nhts.ornl.gov/assets/2016/download/Csv.zip") + +# unpack it +print("Parsing...") +with ZipFile(BytesIO(r.content)) as archive: + + # parse the household data and trip data + # for each day in the desired interval, get the probability of a car making a trip on that day + # for each minute in the desired interval write out what ratio of cars is at home + # for each day in the desired interval write out what trips were the cars taking + with archive.open("hhpub.csv") as hhpub, archive.open("trippub.csv") as trippub: + households = pandas.read_csv(hhpub) + trips = pandas.read_csv(trippub) + + # we're only interested in households from Texas + texans = households[(households["HHSTATE"] == "TX")] + + # we're only interested in trips by Texas residents with their personal cars + texanTrips = trips[(trips["VEHID"] >= 1) & (trips["VEHID"] <= 12) & (trips["HHSTATE"] == "TX")] + + # collect information about each of the households + householdInfo = defaultdict(SimpleNamespace) + + # each household is assigned some weight in the data, save that weight + # each household reports only on one day, save that day + # also save how many vehicles they have + for texan in texans.itertuples(): + householdInfo[texan.HOUSEID].weight = texan.WTHHFIN + householdInfo[texan.HOUSEID].vehicleCount = texan.HHVEHCNT + householdInfo[texan.HOUSEID].reportMonth = texan.TDAYDATE % 100 + householdInfo[texan.HOUSEID].reportWeekday = (texan.TRAVDAY - 1) % 7 + householdInfo[texan.HOUSEID].usedVehicle = [0, 0, 0, 0] + + # for each household save which cars got used on the reported day + for trip in texanTrips.itertuples(): + if trip.VEHID >= 1 and trip.VEHID <= 4: + householdInfo[trip.HOUSEID].usedVehicle[trip.VEHID - 1] = 1 + + # count how many households have how many cars + weightedCarCountOccurences = [0, 0, 0, 0, 0] + for household in householdInfo.values(): + weightedCarCountOccurences[min(household.vehicleCount, 4)] += household.weight + + carCountProbabilities = [occ / sum(weightedCarCountOccurences) for occ in weightedCarCountOccurences] + with open("ownershipRatios.csv", "w+") as ownershipRatiosFile: + ownershipRatiosFile.write("carCount,ratio\n") + for carCount, ratio in enumerate(carCountProbabilities): + ownershipRatiosFile.write(f"{carCount},{ratio:.5f}\n") + + # count the probability that a car gets used on a given day for each day and each car index + usages = [defaultdict(lambda: [0, 0]) for _ in range(4)] + for household in householdInfo.values(): + for car in range(4): + if household.vehicleCount > car: + if household.usedVehicle[car] == 1: + usages[car][(household.reportMonth, household.reportWeekday)][0] += household.weight + usages[car][(household.reportMonth, household.reportWeekday)][1] += household.weight + + # calculate all the statistics for each household car index + for car in range(4): + # select the trips by the desired vehicle index + carTrips = texanTrips[texanTrips["VEHID"] == car + 1] + + # for each vehicle and each day determines the first time it left the house and the last time it returned to the house + def getUsageInterval(df): + sortedDF = df.sort_values("TDTRPNUM") + first = sortedDF.iloc[0] + last = sortedDF.iloc[-1] + # we're interested only in the days where the first trip started from home and last trip ended at home + if (first["WHYFROM"] in [1, 2]) and (last["WHYTO"] in [1, 2]) and (first["STRTTIME"] < last["ENDTIME"]): + year = first["TDAYDATE"] // 100 + month = first["TDAYDATE"] % 100 + # Americans have Sunday as the first day of the week, fix that + weekday = first["TRAVDAY"] - 1 % 7 + departure = "{:02d}:{:02d}".format(first["STRTTIME"] // 100, first["STRTTIME"] % 100) + arrival = "{:02d}:{:02d}".format(last["ENDTIME"] // 100, last["ENDTIME"] % 100) + return (year, month, weekday, departure, arrival) + else: + return None + + # group the trips by household and get the usage intervals of the desired vehicle + times = carTrips.groupby(["HOUSEID"]).apply(getUsageInterval).values + + # filter out reports where nothing happened + times = list(times[times != None]) + daytrips = pandas.DataFrame(times, columns=["Year", "Month", "Weekday", "Departure", "Arrival"]) + + # convert data to dictionary + dateGroups = daytrips.groupby(["Month", "Weekday"])[["Departure", "Arrival"]] + intervals = dateGroups.apply(lambda group: list(group.itertuples(index=False, name=None))).to_dict() + + # create the folder for saving the data + outFolder = f"car{car+1}" + os.makedirs(outFolder, exist_ok=True) + + # write the parsed data out to files + # for each minute in the desired interval write out what ratio of cars is at home + # for each day in the desired interval write out what trips were the cars taking + with open(f"{outFolder}/usageRatios.csv", "w+") as usageRatioFile, open(f"{outFolder}/trips.txt", "w+") as tripsFile, open(f"{outFolder}/availability.csv", "w+") as availabilityFile: + # write headers + usageRatioFile.write("date,usageRatio\n") + availabilityFile.write("datetime,availability\n") + + currentDT = fromDT + while currentDT < toDT: + usage = usages[car][(currentDT.month, currentDT.weekday())] + usageRatio = usage[0] / max(1, usage[1]) + usageRatioFile.write(f"{currentDT.date()},{usageRatio:.5f}\n") + + trips = intervals.get((currentDT.month, currentDT.weekday()), []) + atHome = numpy.full(24*60, fill_value=max(len(trips), 1), dtype=float) + + for dep, arr in trips: + departureSlot = 60 * int(dep[:2]) + int(dep[-2:]) + arrivalSlot = 60 * int(arr[:2]) + int(arr[-2:]) + atHome[departureSlot:arrivalSlot] -= 1 + + availabilityForDay = 1 - (usageRatio * (1 - (atHome / max(len(trips), 1)))) + for minute, availability in enumerate(availabilityForDay): + availabilityFile.write(f"{currentDT+datetime.timedelta(minutes=minute)},{availability:.5f}\n") + + tripsFile.write(f"{currentDT.date()}: [") + tripsFile.write(", ".join([f"{dep}-{arr}" for dep, arr in trips])) + tripsFile.write("]\n") + + currentDT += datetime.timedelta(days=1) + +print("NHTS car trip data retrieved.") diff --git a/simulator/simulator/data/nhts/download.sh b/simulator/simulator/data/nhts/download.sh new file mode 100755 index 0000000..16f208b --- /dev/null +++ b/simulator/simulator/data/nhts/download.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env bash + +export PYTHONUNBUFFERED=1 + +echo "Retreiving NHTS data..." + +cd cars +./getData.py $1 $2 | stdbuf -oL sed 's/^/ /' + +echo "NHTS data retrieved." diff --git a/simulator/simulator/grid.py b/simulator/simulator/grid.py new file mode 100755 index 0000000..ac88de2 --- /dev/null +++ b/simulator/simulator/grid.py @@ -0,0 +1,242 @@ +#!/usr/bin/env python3 + +import datetime + +from typing import List + +import numpy +import scipy +import scipy.signal + +from . import applianceStatistics +from . import gridStatistics +from . import utils + +from .constants import oneDay +from .connection import Connection +from .house import House +from .profile import Profile + +# smart grid main class +class Grid: + # current date and time in the simulation + currentDT: datetime.datetime + # predicted base demand + predictedBaseDemand: Profile + # demand which the houses and their appliances should target + targetDemand: Profile + # demand from houses if they would be using smart appliances + smartDemand: Profile + # demand from houses if their appliances would charge as early as possible + uncontrolledDemand: Profile + # demand from houses if their appliances would charge evenly over their use period + spreadOutDemand: Profile + # how many households should be having cheap electricity at any given time + cheapPriceRatio: Profile + # connections to houses + connections: List[Connection] + + # constructor, just prepares all the variables + def __init__(self): + self.predictedBaseDemand = Profile() + self.targetDemand = Profile() + self.smartDemand = Profile() + self.uncontrolledDemand = Profile() + self.spreadOutDemand = Profile() + self.cheapPriceRatio = Profile() + self.connections = [] + + # connects a house to the grid + def connectHouse(self, house: House): + self.connections.append(Connection(house=house)) + + # sets everything up before the start of the simulation + def setUp(self, dt: datetime.datetime): + self.currentDT = dt + + # predict base demand, calculate target demand and price ratios enough ahead in the future + self.predictBaseDemand(fromDT=dt-3*oneDay, toDT=dt+4*oneDay) + self.calculateTargetDemand(fromDT=dt-2*oneDay, toDT=dt+3.5*oneDay) + self.calculatePriceRatio(fromDT=dt-1*oneDay, toDT=dt+2.5*oneDay) + + # pass the price ratios to the house connections + self.distributePriceRatios(fromDT=dt-1*oneDay, toDT=dt+2.5*oneDay) + + # set up all the connections to the houses + for conn in self.connections: + conn.setUp(dt) + + # moves ahead one day and does all the calculations that need to be done in that day + def tick(self): + self.currentDT += oneDay + cdt = self.currentDT + + # gather and save power demands from all the houses + self.collectDemands(fromDT=cdt-oneDay, toDT=cdt) + + # predict base demand for one more day + self.predictBaseDemand(fromDT=cdt+3*oneDay, toDT=cdt+4*oneDay) + + # calculate target demand and price ratios for one more day + self.calculateTargetDemand(fromDT=cdt+2.5*oneDay, toDT=cdt+3.5*oneDay) + self.calculatePriceRatio(fromDT=cdt+1.5*oneDay, toDT=cdt+2.5*oneDay) + + # pass the new price ratios to the house connections + self.distributePriceRatios(fromDT=cdt+1.5*oneDay, toDT=cdt+2.5*oneDay) + + # move ahead one day in all the connections to houses as well + for conn in self.connections: + conn.tick() + + # predicts the power demand on the grid without the connected houses + def predictBaseDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # right now this just takes the total demand forecast and subtracts the recorded draw from households during the specified interval + # in an actual grid this would do some fancy calculations to get the prediction + demandForecast = gridStatistics.demandForecast.demand.get(fromDT, toDT) * (len(self.connections) / gridStatistics.demandForecast.householdCount) + householdDraw = gridStatistics.averageHouseholdDraw.get(fromDT, toDT) * len(self.connections) + + baseDemandPrediction = demandForecast - householdDraw + + self.predictedBaseDemand.set(fromDT, baseDemandPrediction) + + # calculates the target demand for the connected houses, based on the base demand, household power usage estimates (and power generation predictions, if available) + def calculateTargetDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # right now this just calculates a demand that smooths out the base demand and is big enough to cover all the needs of the households + # in an actual grid this would also take into account the power generation predictions (solar and wind generation, power plant shutdowns etc) + + # calculate the expected household demand in the given interval + # get power usage estimates for all the simulated appliances for each day in the interval and sum them together + totalExpectedConsumption = 0.0 + for fraction, day in utils.dayPortionsBetween(fromDT, toDT): + expectedDayConsumption = 0.0 + for carIndex in range(4): + expectedDayConsumption += applianceStatistics.atLeastThisManyCarsProbability[carIndex+1] * applianceStatistics.carStatistics[carIndex].averageNeededCharge[day] + expectedDayConsumption += applianceStatistics.ownershipRatios.airConditioning * applianceStatistics.airConditioningStatistics.averageDailyCharge[day] + expectedDayConsumption += applianceStatistics.ownershipRatios.electricalHeating * applianceStatistics.electricalHeatingStatistics.averageDailyCharge[day] + expectedDayConsumption += applianceStatistics.ownershipRatios.fridge * applianceStatistics.fridgeStatistics.averageDailyCharge[day] + expectedDayConsumption += applianceStatistics.ownershipRatios.waterHeater * applianceStatistics.waterHeaterStatistics.averageDailyCharge[day] + expectedDayConsumption += applianceStatistics.ownershipRatios.dishwasher * applianceStatistics.dishwasherStatistics.averagePowerNeeded[day] + expectedDayConsumption += applianceStatistics.ownershipRatios.washingMachine * applianceStatistics.washingMachineStatistics.averagePowerNeeded[day] + + totalExpectedConsumption += fraction * len(self.connections) * expectedDayConsumption + + # introduce some error in the statistics + totalExpectedConsumption *= 0.9 + numpy.random.random() * 0.2 + + # find peaks in the base demand and interpolate between them to get a smooth curve + # have some margin at the ends of the desired interval to have a better interpolation + startMargin = oneDay + endMargin = 0.5*oneDay + startIndex = utils.minutesIn(startMargin) + + baseDemand = self.predictedBaseDemand.get(fromDT-startMargin, toDT+endMargin) + + peaks = list(scipy.signal.find_peaks(baseDemand, distance=18*60, width=10)[0]) + + peakLocs = [0] + peaks + [baseDemand.size-1] + peakVals = baseDemand[[peaks[0]] + peaks + [peaks[-1]]] + + if len(peakLocs) > 3: + interpolationKind = 'cubic' + else: + interpolationKind = 'quadratic' + + smoothDemand = scipy.interpolate.interp1d(peakLocs, peakVals, kind=interpolationKind)(range(baseDemand.size)) + + # set the target household demand as the difference between the smooth demand and the base demand + targetDemand = (smoothDemand - baseDemand)[startIndex:] + + # calculate the integral of the target demand + intervalLength = utils.minutesBetween(fromDT, toDT) + totalTargetIntervalConsumption = numpy.sum(targetDemand[:intervalLength]) / 60 + + if totalExpectedConsumption <= totalTargetIntervalConsumption: + # if the total target consumption is higher than the expected demand of the households, scale it down + targetDemand *= (totalExpectedConsumption / totalTargetIntervalConsumption) + else: + # otherwise shift the target demand up so it covers all the needs of the households + targetDemand += ((totalExpectedConsumption - totalTargetIntervalConsumption) / (intervalLength / 60)) + + # negative demand is not possible in this simulation (eventually could be with vehicle-to-grid systems) + targetDemand = numpy.maximum(targetDemand, 0) + + # the target demand calculations between the previous interval and this interval possibly don't join nicely + # we have to smoothly transition from the previously calculated target demand to the new one (there should be a 12-hour overlap) + self.targetDemand.transition(fromDT=fromDT, newValues=targetDemand) + + # calculates price ratios based on the target demand and appliance availability statistics + # more households should have a cheap electricity price when the target demand is higher + def calculatePriceRatio(self, fromDT: datetime.datetime, toDT: datetime.datetime): + # get the target demand with some margin at the ends of the desired interval for a better interpolation + startMargin = oneDay + endMargin = oneDay + startIndex = utils.minutesIn(startMargin) + + targetDemand = self.targetDemand.get(fromDT-startMargin, toDT+endMargin) + + # scale the target demand based on how many appliances are available to be turned on + # right now the only appliances that matter and are not available continuously are cars + # dishwashers and washing machines contribute so little to the total power consumption that scaling based on their usage intervals doesn't really make sense + + # get the ratio of how much of the target demand is expected to be used by cars + totalExpectedCarConsumption = 0.0 + for fraction, day in utils.dayPortionsBetween(fromDT-startMargin, toDT+endMargin): + for carIndex in range(4): + totalExpectedCarConsumption += fraction * len(self.connections) * applianceStatistics.atLeastThisManyCarsProbability[carIndex+1] * applianceStatistics.carStatistics[carIndex].averageNeededCharge[day] + carDemandRatio = totalExpectedCarConsumption / numpy.sum(targetDemand) + + # get the statistics for how many cars that need a charge are likely to be at home + totalNeedChargingRatio = 0 + carsAtHome = numpy.zeros(utils.minutesBetween(fromDT-startMargin, toDT+endMargin)) + for carIndex in range(4): + needChargingRatio = 0 + totalFraction = 0 + for fraction, day in utils.dayPortionsBetween(fromDT-startMargin, toDT+endMargin): + needChargingRatio += fraction * applianceStatistics.carStatistics[carIndex].usageProbabilities[day] + totalFraction += fraction + needChargingRatio /= totalFraction + carsAtHome += needChargingRatio * applianceStatistics.carStatistics[carIndex].availabilityProfile.get(fromDT-startMargin, toDT+endMargin) + totalNeedChargingRatio += needChargingRatio + carsAtHome /= totalNeedChargingRatio + + # we need to give more households cheaper prices when there are less cars at home, so those that are at home would charge and cover the target demand + availabilityScale = (1 - carDemandRatio) + carDemandRatio * carsAtHome + relativeTargetDemand = targetDemand / availabilityScale + + # we need to make sure the peaks each day correspond to cheap prices for all households + # therefore we need to scale the (already scaled) target demand so that the peaks each day scale to 1, to get nice probabilities + # get the peaks in the relative target demand, approximately one each day + peaks = list(scipy.signal.find_peaks(relativeTargetDemand, distance=18*60, width=10)[0]) + peakLocs = [0] + peaks + [relativeTargetDemand.size-1] + peakVals = relativeTargetDemand[[peaks[0]] + peaks + [peaks[-1]]] + + # interpolate between the peaks, to get a smoothed out rolling maximum curve + demandScale = utils.cosineInterpolation(peakLocs, peakVals) + # divide the relative target demand by the smoothed out maximum curve, to scale the peaks to 1 to get the price ratio + scaledDemand = relativeTargetDemand / demandScale + cheapPriceRatio = scaledDemand[startIndex:] + + # the price ratio calculations between the previous interval and this one might not join nicely + # we have to smoothly transition from the previously calculated price ratios to the new ones (there should be a 24-hour overlap) + self.cheapPriceRatio.transition(fromDT, cheapPriceRatio) + + # distributes the calculated price ratios to all the connections + def distributePriceRatios(self, fromDT: datetime.datetime, toDT: datetime.datetime): + cheapPriceRatio = self.cheapPriceRatio.get(fromDT, toDT) + for conn in self.connections: + conn.setPriceRatio(fromDT, cheapPriceRatio) + + # collects the power demands from all the connected households + def collectDemands(self, fromDT: datetime.datetime, toDT: datetime.datetime): + length = utils.minutesBetween(fromDT, toDT) + smartDemand = numpy.zeros(length) + uncontrolledDemand = numpy.zeros(length) + spreadOutDemand = numpy.zeros(length) + for conn in self.connections: + smartDemand += conn.getSmartDemand(fromDT, toDT) + uncontrolledDemand += conn.getUncontrolledDemand(fromDT, toDT) + spreadOutDemand += conn.getSpreadOutDemand(fromDT, toDT) + + self.smartDemand.set(fromDT, smartDemand) + self.uncontrolledDemand.set(fromDT, uncontrolledDemand) + self.spreadOutDemand.set(fromDT, spreadOutDemand) diff --git a/simulator/simulator/gridStatistics.py b/simulator/simulator/gridStatistics.py new file mode 100755 index 0000000..35a4fbb --- /dev/null +++ b/simulator/simulator/gridStatistics.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python3 + +import os + +from types import SimpleNamespace + +from .profile import Profile + +thisDir = os.path.dirname(__file__) + +# class for representing statistics of grid electricity demand +class GridDemandStatistics(SimpleNamespace): + # households connected to the grid + householdCount: int + # electricity demand profile of the grid + demand: Profile + +# 4-day-ahead prediction of the Texas grid electricity demand +demandForecast = GridDemandStatistics() +demandForecast.demand = Profile.fromCSV(f"{thisDir}/data/dataport/ercot/predictions/96.csv") +demandForecast.householdCount = 9500000 + + +# actual Texas grid electricity demand +actualDemand = GridDemandStatistics() +actualDemand.demand = Profile.fromCSV(f"{thisDir}/data/dataport/ercot/actual/systemLoad.csv") +actualDemand.householdCount = 9500000 + +# average power usage of a household +averageHouseholdDraw = Profile.fromCSV(f"{thisDir}/data/dataport/household/averageDraw.csv") diff --git a/simulator/simulator/house.py b/simulator/simulator/house.py new file mode 100755 index 0000000..e079473 --- /dev/null +++ b/simulator/simulator/house.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 + +import datetime +import random + +from typing import List + +import numpy + +from . import applianceStatistics, utils + +from .appliance import Appliance, Car, AirConditioning, ElectricalHeating, WaterHeater, Fridge, WashingMachine, Dishwasher +from .constants import oneDay +from .profile import Profile + +# class representing a house connected to the smart grid +class House: + # current date and time in the simulation + currentDT: datetime.datetime + # the appliances in this house + appliances: List[Appliance] + # the electricity prices for any given minute for this house + priceProfile: Profile + # the electricity demand if this house was using smart appliances + smartDemand: Profile + # the electricity demand if the appliances in this house would charge as early as possible + uncontrolledDemand: Profile + # the electricity demand if the appliances in this house would charge evenly over their use period + spreadOutDemand: Profile + + # constructor, just sets up the variables + def __init__(self): + self.appliances = [] + self.priceProfile = Profile() + self.smartDemand = Profile() + self.uncontrolledDemand = Profile() + self.spreadOutDemand = Profile() + + # creates a random house with random appliances according to appliance ownership statistics + @classmethod + def random(cls): + h = cls() + carCount = utils.randomWithRelativeProbs(applianceStatistics.carCountProbabilities) + for car in range(carCount): + h.addAppliance(Car.randomWithIndex(index=car)) + + if random.random() < applianceStatistics.ownershipRatios.airConditioning: + h.addAppliance(AirConditioning.random()) + if random.random() < applianceStatistics.ownershipRatios.electricalHeating: + h.addAppliance(ElectricalHeating.random()) + if random.random() < applianceStatistics.ownershipRatios.waterHeater: + h.addAppliance(WaterHeater.random()) + if random.random() < applianceStatistics.ownershipRatios.fridge: + h.addAppliance(Fridge.random()) + if random.random() < applianceStatistics.ownershipRatios.washingMachine: + h.addAppliance(WashingMachine.random()) + if random.random() < applianceStatistics.ownershipRatios.dishwasher: + h.addAppliance(Dishwasher.random()) + return h + + # adds an appliance to the house + def addAppliance(self, appliance: Appliance): + self.appliances.append(appliance) + + # sets up the house for the simulation + def setUp(self, dt: datetime.datetime): + self.currentDT = dt + # set up all the appliances in the house + for appliance in self.appliances: + appliance.setUp(dt) + + # moves ahead one day and does all the calculations that need to be done in that day + def tick(self): + # remove past, unneeded values from the profiles to free up some memory + self.priceProfile.prune(self.currentDT - oneDay) + self.smartDemand.prune(self.currentDT - oneDay) + self.uncontrolledDemand.prune(self.currentDT - oneDay) + self.spreadOutDemand.prune(self.currentDT - oneDay) + + # move ahead one day in all the appliances in this house as well + for appliance in self.appliances: + appliance.tick() + + # move the current time forward + self.currentDT += oneDay + + # collect the electricity demands from all the appliances + self.collectApplianceDemand(self.currentDT - oneDay, self.currentDT) + + # sets the electricity price profile for this house + # called by the connection to the grid + def setPriceProfile(self, dt: datetime.datetime, prices: numpy.ndarray): + self.priceProfile.set(dt, prices) + # pass on the price profile to all the appliances in this house + for appliance in self.appliances: + appliance.setPriceProfile(dt, prices) + + # gets the electricity demand if this house was using smart appliances + def getSmartDemand(self, fromDT: datetime.datetime = None, toDT: datetime.datetime = None): + return self.smartDemand.get(fromDT, toDT) + + # gets the electricity demand if this house was NOT using smart appliances + def getUncontrolledDemand(self, fromDT: datetime.datetime = None, toDT: datetime.datetime = None): + return self.uncontrolledDemand.get(fromDT, toDT) + + # gets the electricity demand if this house was NOT using smart appliances + def getSpreadOutDemand(self, fromDT: datetime.datetime = None, toDT: datetime.datetime = None): + return self.spreadOutDemand.get(fromDT, toDT) + + # collects the electricity demands from all the appliances in the house + def collectApplianceDemand(self, fromDT: datetime.datetime, toDT: datetime.datetime): + for appliance in self.appliances: + smartDemand = appliance.smartDemand.get(fromDT, toDT) + uncontrolledDemand = appliance.uncontrolledDemand.get(fromDT, toDT) + spreadOutDemand = appliance.spreadOutDemand.get(fromDT, toDT) + self.smartDemand.add(fromDT, smartDemand) + self.uncontrolledDemand.add(fromDT, uncontrolledDemand) + self.spreadOutDemand.add(fromDT, spreadOutDemand) diff --git a/simulator/simulator/priceConfig.py b/simulator/simulator/priceConfig.py new file mode 100644 index 0000000..53f0331 --- /dev/null +++ b/simulator/simulator/priceConfig.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 + +import json +import os + +from types import SimpleNamespace + +thisDir = os.path.dirname(__file__) + +# load the config from a file +with open(f"{thisDir}/data/manual/priceConfig.json", "r") as priceConfigFile: + priceConfig = json.load(priceConfigFile, object_hook=lambda dct: SimpleNamespace(**dct)) + +# length of each interval of cheap electricity price +cheapIntervalLength = priceConfig.cheapIntervalLength +# minimum minutes the electricity price will be cheaper each day +cheapMinutesCount = priceConfig.cheapMinutesCount + +# the electricity prices [money per kWh] +lowerPrice = priceConfig.lowerPrice +higherPrice = priceConfig.higherPrice diff --git a/simulator/simulator/profile.py b/simulator/simulator/profile.py new file mode 100755 index 0000000..6e0e868 --- /dev/null +++ b/simulator/simulator/profile.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 + +import datetime +from typing import Dict + +import numpy +import pandas + +from .utils import minutesBetween + +# helper class to deal with time series +class Profile: + # date and time of the first stored value + startingDT: datetime.datetime + # the actual stored values + values: numpy.ndarray + + # constructor + def __init__(self, startingDT: datetime.datetime = None, values: numpy.ndarray = None): + self.startingDT = startingDT + if values is None or startingDT is None: + self.values = numpy.empty(0) + else: + self.values = values.copy() + + # get values between fromDT and toDT, return zeros if the values are missing + def get(self, fromDT: datetime.datetime, toDT: datetime.datetime = None): + if self.startingDT is None: + if toDT is None: + raise IndexError() + else: + length = minutesBetween(fromDT, toDT) + return numpy.zeros(length) + + startIndex = minutesBetween(self.startingDT, fromDT) + if toDT is None: + if startIndex < 0 or startIndex >= self.values.size: + return 0 + else: + return self.values[startIndex] + else: + length = minutesBetween(fromDT, toDT) + res = numpy.zeros(length, dtype=float) + stopIndex = min(self.values.size, minutesBetween(self.startingDT, toDT)) + res[:stopIndex-startIndex] = self.values[startIndex:stopIndex] + return res + + # set values starting at fromDT to newValues + def set(self, fromDT: datetime.datetime, newValues: numpy.ndarray): + if self.startingDT is None: + self.startingDT = fromDT + self.values = newValues.copy() + else: + startIndex = minutesBetween(self.startingDT, fromDT) + newSize = startIndex + newValues.size + + if newSize > self.values.size: + self.values.resize(newSize, refcheck=False) + + self.values[startIndex:startIndex+newValues.size] = newValues + + # add up values in valuesToAdd to currently stored values starting at fromDT + def add(self, fromDT: datetime.datetime, valuesToAdd: numpy.ndarray): + if self.startingDT is None: + self.startingDT = fromDT + + startIndex = minutesBetween(self.startingDT, fromDT) + newSize = startIndex + valuesToAdd.size + + if newSize > self.values.size: + self.values.resize(newSize, refcheck=False) + self.values[startIndex:startIndex+valuesToAdd.size] += valuesToAdd + + # smoothly transition from the currently stored values to newValues starting at fromDT using a cosine interpolation + def transition(self, fromDT: datetime.datetime, newValues: numpy.ndarray): + if self.startingDT is None: + self.set(fromDT, newValues) + else: + startIndex = minutesBetween(self.startingDT, fromDT) + overlappingValues = self.values.size - startIndex + + newSize = startIndex + newValues.size + if newSize > self.values.size: + self.values.resize(newSize, refcheck=False) + + if overlappingValues == 0: + self.values = numpy.append(self.values, newValues) + elif overlappingValues < 0: + self.values[startIndex:] = newValues + else: + oldValues = self.values[startIndex:] + ratio = (numpy.cos(numpy.linspace(0, numpy.pi, num=overlappingValues, endpoint=False)) + 1) / 2 + oldMask = numpy.pad(ratio, (0, newValues.size - overlappingValues), mode='constant', constant_values=0) + newMask = numpy.pad(1-ratio, (0, newValues.size - overlappingValues), mode='constant', constant_values=1) + self.values[startIndex:] = oldValues * oldMask + newValues * newMask + + # multiply the stored values by scale + def scale(self, scale): + self.values *= scale + + # delete stored values up to, but not including, toDT + def prune(self, toDT: datetime.datetime): + if self.startingDT is None: + return + if toDT <= self.startingDT: + return + + index = minutesBetween(self.startingDT, toDT) + self.startingDT = toDT + self.values = self.values[index:].copy() + + # return a copy of this Profile + def copy(self): + return Profile(self.startingDT, self.values.copy()) + + # return a dictionary of the averate of stored values for each date in this Profile + def dailyAverages(self) -> Dict[datetime.date, float]: + if self.startingDT is None or self.values.size == 0: + return {} + + averages = {} + + currentDate = self.startingDT.date() + startIndex = minutesBetween(self.startingDT, currentDate) + endIndex = startIndex + 24*60 + + while startIndex < self.values.size: + dayValues = self.values[max(0, startIndex):min(endIndex, self.values.size)] + averages[currentDate] = numpy.mean(dayValues) + + currentDate += datetime.timedelta(days=1) + startIndex += 24*60 + endIndex += 24*60 + return averages + + @classmethod + def fromCSV(cls, path): + data = pandas.read_csv(path, parse_dates=[0]) + startingDT = data.iloc[0, 0].to_pydatetime() + values = data.iloc[:, 1].values + return cls(startingDT, values) diff --git a/simulator/simulator/simulator.py b/simulator/simulator/simulator.py new file mode 100755 index 0000000..be9f981 --- /dev/null +++ b/simulator/simulator/simulator.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 + +import datetime +import gc +import os +import time + +import pandas + +from . import gridStatistics, priceConfig + +from .constants import oneDay +from .grid import Grid +from .house import House + +# smart grid simulator main class +class Simulator: + # run the smart grid simulation + @classmethod + def run(cls, startingDT: datetime.datetime, simulationLength: int, houseCount: int, outputFolder: str = None): + # remember the starting time + st = time.time() + # create the grid + print("Creating grid...") + grid = Grid() + + # create random houses and connect them to the grid + print("Creating houses...") + houses = [] + for _ in range(houseCount): + h = House.random() + grid.connectHouse(h) + houses.append(h) + + # set up the grid with the right time + # the grid sets up the connected houses itself + print("Setting up grid...") + grid.setUp(startingDT) + + print("Preparation took {:.3f}s".format(time.time() - st)) + print() + + # for each day of the simulation, send tick signals to the grid and houses for them to tell them time has moved + # this is to make it easier to possibly change to a simulator architecture with multiple processes + # where the simulator only provides the clock signal and the grid and houses take care of everything else + endDT = startingDT + simulationLength * oneDay + currentDT = startingDT + while currentDT < endDT: + print("Calculating power draw for", currentDT.date()) + t = time.time() + for i, h in enumerate(houses): + print(f"\r{i+1}/{houseCount}... ", end="") + h.tick() + grid.tick() + # call garbage collection manually to ease memory pressure + gc.collect() + print("Calculation took {:.3f}s".format(time.time() - t)) + print() + currentDT += oneDay + + print("Simulation took {:.3f}s in total".format(time.time() - st)) + print() + + # collect the results from the grid + predictedBaseDemand = grid.predictedBaseDemand.get(startingDT, endDT) + targetDemand = grid.targetDemand.get(startingDT, endDT) + priceRatio = grid.cheapPriceRatio.get(startingDT, endDT) + smartDemand = grid.smartDemand.get(startingDT, endDT) + uncontrolledDemand = grid.uncontrolledDemand.get(startingDT, endDT) + spreadOutDemand = grid.spreadOutDemand.get(startingDT, endDT) + + # get the actual grid base demand + actualDemand = gridStatistics.actualDemand.demand.get(startingDT, endDT) * (houseCount / gridStatistics.actualDemand.householdCount) + householdDraw = gridStatistics.averageHouseholdDraw.get(startingDT, endDT) * houseCount + actualBaseDemand = actualDemand - householdDraw + + # prepare the datetime column + datetimes = [startingDT + datetime.timedelta(minutes=i) for i in range(simulationLength*24*60)] + + # gather the results in a pandas dataframe + data = pandas.DataFrame({ + "Datetime": datetimes, + "PredictedBaseDemand": predictedBaseDemand, + "ActualBaseDemand": actualBaseDemand, + "TargetDemand": targetDemand, + "SmartDemand": smartDemand, + "UncontrolledDemand": uncontrolledDemand, + "SpreadOutDemand": spreadOutDemand, + "PriceRatio": priceRatio, + }) + + # save the results to a folder if specified + if outputFolder is not None: + os.makedirs(outputFolder, exist_ok=True) + + # save the simulation parameters to a separate file + with open(f"{outputFolder}/desc.txt", "w+") as descfile: + descfile.write(f"startingDatetime={startingDT}\n") + descfile.write(f"simulationLength={simulationLength}\n") + descfile.write(f"houseCount={houseCount}\n") + descfile.write(f"lowerPrice={priceConfig.lowerPrice}\n") + descfile.write(f"higherPrice={priceConfig.higherPrice}\n") + descfile.write(f"cheapIntervalLength={priceConfig.cheapIntervalLength}\n") + descfile.write(f"cheapMinutesTotal={priceConfig.cheapMinutesCount}\n") + + # save the demand values to a csv + data.to_csv(f"{outputFolder}/data.csv", index=False, header=True, float_format="%.5f") + + # return results + return data diff --git a/simulator/simulator/utils.py b/simulator/simulator/utils.py new file mode 100755 index 0000000..403fdcb --- /dev/null +++ b/simulator/simulator/utils.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python3 + +import datetime +from typing import List, Tuple + +import numpy + +from .constants import oneDay + +# helper functions frequently used in the simulator + +# counts the minutes in a timedelta object +def minutesIn(td: datetime.timedelta) -> int: + return int(td.total_seconds() // 60) + +# counts the minutes between two datetimes +def minutesBetween(startDT: datetime.datetime, endDT: datetime.datetime) -> int: + if isinstance(startDT, datetime.date) and not isinstance(startDT, datetime.datetime): + startDT = datetime.datetime.combine(startDT, datetime.time()) + if isinstance(endDT, datetime.date) and not isinstance(endDT, datetime.datetime): + endDT = datetime.datetime.combine(endDT, datetime.time()) + return minutesIn(endDT - startDT) + +# returns the midnights between two datetimes +def midnightsBetween(startDT: datetime.datetime, endDT: datetime.datetime) -> List[datetime.datetime]: + midnights = [] + if (startDT.hour, startDT.minute, startDT.second, startDT.microsecond) == (0, 0, 0, 0): + midnights.append(startDT) + currentDT = startDT.replace(hour=0, minute=0, second=0, microsecond=0) + oneDay + while currentDT < endDT: + midnights.append(currentDT) + currentDT += oneDay + return midnights + +# returns the days between two datetimes, and for each day also a fraction which is covered in the given interval +def dayPortionsBetween(startDT: datetime.datetime, endDT: datetime.datetime) -> List[Tuple[float, datetime.date]]: + if startDT >= endDT: + return [] + elif startDT.date() == endDT.date(): + return [((endDT - startDT).total_seconds() / (24*60*60), startDT.date())] + else: + portions = [] + currentDT = startDT + nextMidnight = startDT.replace(hour=0, minute=0, second=0, microsecond=0) + oneDay + while nextMidnight <= endDT: + portions.append(((nextMidnight - currentDT).total_seconds() / (24*60*60), currentDT.date())) + currentDT = nextMidnight + nextMidnight = currentDT + oneDay + if currentDT < endDT: + portions.append(((endDT - currentDT).total_seconds() / (24*60*60), currentDT.date())) + return portions + +# randomly choose an integer with relative probabilities provided in an array +def randomWithRelativeProbs(relativeProbs: numpy.ndarray, count: int = None): + if not isinstance(relativeProbs, numpy.ndarray): + relativeProbs = numpy.array(relativeProbs) + probs = relativeProbs / numpy.sum(relativeProbs) + if count is None: + return numpy.random.choice(probs.size, p=probs) + else: + return numpy.random.choice(probs.size, size=count, replace=False, p=probs) + +# perform a cosine interpolation from a list of coordinates +def cosineInterpolation(xs: List[int], ys: List[float]) -> numpy.ndarray: + res = [] + for x1, x2, y1, y2 in zip(xs[:-1], xs[1:], ys[:-1], ys[1:]): + ratio = (numpy.cos(numpy.linspace(0, numpy.pi, num=x2-x1, endpoint=False)) + 1) / 2 + vals = y1 * ratio + y2 * (1 - ratio) + res.extend(vals) + res.append(ys[-1]) + return numpy.array(res) diff --git a/thesis/thesis.pdf b/thesis/thesis.pdf new file mode 100644 index 0000000..d6afd2c Binary files /dev/null and b/thesis/thesis.pdf differ