diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 43e8d86..4eb1a8d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -17,6 +17,10 @@ jobs: runs-on: ubuntu-latest + defaults: + run: + shell: bash -l {0} + steps: - uses: actions/checkout@v3 - name: Set up conda @@ -29,10 +33,6 @@ jobs: use-only-tar-bz2: true - run: | conda config --env --set pip_interop_enabled True - - name: Set up Python 3.10 - uses: actions/setup-python@v3 - with: - python-version: "3.10" - name: Install dependencies run: | python -m pip install --upgrade pip diff --git a/docs/source/conf.py b/docs/source/conf.py index 8a35f6f..d4e0f7a 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -32,14 +32,15 @@ # ones. extensions = [ "myst_parser", - "numpydoc", "sphinx.ext.napoleon", "sphinx.ext.autodoc", "sphinx.ext.autosummary", "sphinx.ext.intersphinx", "sphinx.ext.viewcode", "sphinx_autodoc_typehints", - "sphinx_design" + "sphinx_design", + "sphinx.ext.mathjax", + "sphinx.ext.coverage", ] # Add any paths that contain templates here, relative to this directory. @@ -50,6 +51,13 @@ # This pattern also affects html_static_path and html_extra_path. exclude_patterns = [] +intersphinx_mapping = { + 'python': ('https://docs.python.org/3', None), + 'numpy': ('https://numpy.org/doc/stable/', None), + 'astropy': ('https://docs.astropy.org/en/stable/', None), + 'jinja2': ('https://jinja.palletsprojects.com/en/3.0.x/', None), + 'unyt': ('https://unyt.readthedocs.io/en/stable/', None), +} # -- Options for HTML output ------------------------------------------------- diff --git a/docs/source/getting-started/index.md b/docs/source/getting-started/index.md new file mode 100644 index 0000000..0a93ea1 --- /dev/null +++ b/docs/source/getting-started/index.md @@ -0,0 +1,19 @@ +# Getting Started + +`osier` is intended to be extremely user friendly and most problems +should be solvable without an external solver. However, the current implementation +of the `osier.DispatchModel` as a linear program requires an external solver +such as CPLEX, CBC, GLPK, or Gurobi (since these are all supported by `pyomo`). + +CPLEX and Gurobi are commercial solvers, however it is quite simple to obtain a free +academic license (instructions forthcoming). GLPK will work on Windows, but requires +installing external binaries. CBC is a good open-source solver for a unix operating +system such as Mac or Linux. + +In order to use CBC on the latter two operating systems you must have a version +of [Anaconda/conda](https://www.anaconda.com/products/distribution) installed. +Then you can install it using + +```bash +$ conda install -c conda-forge coincbc +``` diff --git a/docs/source/index.md b/docs/source/index.md index 1a941c6..f06cccd 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -17,6 +17,7 @@ This package is in active development. contrib reference/index +getting-started/index ``` ## Indices and tables diff --git a/docs/source/reference/index.md b/docs/source/reference/index.md index edea6a5..03411cc 100644 --- a/docs/source/reference/index.md +++ b/docs/source/reference/index.md @@ -7,5 +7,13 @@ :template: my_class.rst osier.Technology + osier.DispatchModel -``` \ No newline at end of file +``` + + \ No newline at end of file diff --git a/osier/__init__.py b/osier/__init__.py index 44224b1..8c0c5d3 100644 --- a/osier/__init__.py +++ b/osier/__init__.py @@ -1 +1,2 @@ from .technology import * +from .models.dispatch import * \ No newline at end of file diff --git a/osier/models/dispatch.py b/osier/models/dispatch.py new file mode 100644 index 0000000..858a418 --- /dev/null +++ b/osier/models/dispatch.py @@ -0,0 +1,203 @@ +import pandas as pd +import numpy as np +import pyomo.environ as pe +import pyomo.opt as po +from pyomo.environ import ConcreteModel +from unyt import unyt_array +import itertools as it + + +class DispatchModel(): + """ + The :class:`DispatchModel` class creates and solves a basic dispatch + model from the perspective of a "grid operator." The model uses + `pyomo `_ + to create and solve a linear programming model. + The mathematical formulation for this problem is: + + Minimize + + + .. math:: + \\text{C}_{\\text{total}} = \\sum_t^T \\sum_u^U + \\left[c^{\\text{fuel}}_{u,t} + c^{\\text{om,var}}_{u,t}\\right]x_{u,t} + + Such that, + + 1. The generation meets demand within a user-specified tolerance + (undersupply and oversupply) + + .. math:: + \\sum_u^Ux_{u,t} &\\geq \\left(1-\\text{undersupply}\\right)\\text{D}_t + \\ \\forall \\ t \\in T + + \\sum_u^Ux_{u,t} &\\leq \\left(1+\\text{oversupply}\\right)\\text{D}_t + \\ \\forall \\ t \\in T + + 2. A technology's generation (:math:`x_u`) does not exceed its capacity to + generate at any time, :math:`t`. + + .. math:: + x_{u,t} \\leq \\textbf{CAP}_{u} \\ \\forall \\ u,t \\in U,T + + + Parameters + ---------- + technology_list : list of :class:`osier.Technology` + The list of :class:`Technology` objects to dispatch -- i.e. decide + how much energy each technology should produce. + net_demand : List, :class:`numpy.ndarray`, :class:`unyt.array.unyt_array` + The remaining energy demand to be fulfilled by the technologies in + :attr:`technology_list`. For example: + + >>> from osier import DispatchModel + >>> from osier import Technology + >>> wind = Technology("wind") + >>> nuclear = Technology("nuclear") + >>> import numpy as np + >>> hours = np.arange(24) + >>> demand = np.cos(hours*np.pi/180) + >>> model = DispatchModel([wind, nuclear], demand) + >>> model.solve() + + solver : str + Indicates which solver to use. May require separate installation. + Accepts: ['cplex', 'cbc']. Other solvers will be added in the future. + lower_bound : float + The minimum amount of energy each technology can produce per time + period. Default is 0.0. + oversupply : float + The amount of allowed oversupply as a percentage of demand. + Default is 0.0 (no oversupply allowed). + undersupply : float + The amount of allowed undersupply as a percentage of demand. + Default is 0.0 (no undersupply allowed). + + Attributes + ---------- + model : :class:`pyomo.environ.ConcreteModel` + The :mod:`pyomo` model class that converts python code into a + set of linear equations. + upperbound : float + The upper bound for all decision variables. Chosen to be equal + to the maximum capacity of all technologies in :attr:`tech_set`. + """ + + def __init__(self, + technology_list, + net_demand, + solver='cplex', + lower_bound=0.0, + oversupply=0.0, + undersupply=0.0): + self.net_demand = net_demand + self.technology_list = technology_list + self.lower_bound = lower_bound + self.oversupply = oversupply + self.undersupply = undersupply + self.model = ConcreteModel() + self.solver = solver + self.results = None + self.objective = None + + @property + def n_hours(self): + return len(self.net_demand) + + @property + def tech_set(self): + tech_names = [t.technology_name for t in self.technology_list] + return tech_names + + @property + def capacity_dict(self): + capacity_set = unyt_array([t.capacity for t in self.technology_list]) + return dict(zip(self.tech_set, capacity_set)) + + @property + def time_set(self): + return range(self.n_hours) + + @property + def indices(self): + return list(it.product(self.tech_set, self.time_set)) + + @property + def cost_params(self): + v_costs = np.array([ + (t.variable_cost_ts(self.n_hours)) + for t in self.technology_list + ]).flatten() + return dict(zip(self.indices, v_costs)) + + @property + def upper_bound(self): + caps = unyt_array([t.capacity for t in self.technology_list]) + return caps.max().to_value() + + def _create_model_indices(self): + self.model.U = pe.Set(initialize=self.tech_set, ordered=True) + self.model.T = pe.Set(initialize=self.time_set, ordered=True) + + def _create_demand_param(self): + self.model.D = pe.Param(self.model.T, initialize=dict( + zip(self.model.T, self.net_demand))) + + def _create_cost_param(self): + self.model.C = pe.Param( + self.model.U, + self.model.T, + initialize=self.cost_params) + + def _create_model_variables(self): + self.model.x = pe.Var(self.model.U, self.model.T, + domain=pe.Reals, + bounds=(self.lower_bound, self.upper_bound)) + + def _objective_function(self): + expr = sum(self.model.C[u, t] * self.model.x[u, t] + for u in self.model.U for t in self.model.T) + self.model.objective = pe.Objective(sense=pe.minimize, expr=expr) + + def _supply_constraints(self): + self.model.oversupply = pe.ConstraintList() + self.model.undersupply = pe.ConstraintList() + for t in self.model.T: + generation = sum(self.model.x[u, t] for u in self.model.U) + over_demand = self.model.D[t] * (1 + self.oversupply) + under_demand = self.model.D[t] * (1 - self.undersupply) + self.model.oversupply.add(generation <= over_demand) + self.model.undersupply.add(generation >= under_demand) + + def _generation_constraint(self): + self.model.gen_limit = pe.ConstraintList() + for t in self.model.T: + for u in self.model.U: + unit_generation = self.model.x[u, t] + unit_capacity = self.capacity_dict[u].to_value() + self.model.gen_limit.add(unit_generation <= unit_capacity) + + def _write_model_equations(self): + + self._create_model_indices() + self._create_demand_param() + self._create_cost_param() + self._create_model_variables() + self._objective_function() + self._supply_constraints() + self._generation_constraint() + + def _format_results(self): + df = pd.DataFrame(index=self.model.T) + for u in self.model.U: + df[u] = [self.model.x[u, t].value for t in self.model.T] + + return df + + def solve(self): + self._write_model_equations() + solver = po.SolverFactory(self.solver) + results = solver.solve(self.model, tee=True) + self.objective = self.model.objective() + + self.results = self._format_results() diff --git a/osier/technology.py b/osier/technology.py index 35f0a78..ea3aebf 100644 --- a/osier/technology.py +++ b/osier/technology.py @@ -3,6 +3,8 @@ from unyt import unyt_quantity from unyt.exceptions import UnitParseError, UnitConversionError +import numpy as np + _dim_opts = {'time': hr, 'power': MW, @@ -64,10 +66,10 @@ def _validate_quantity(value, dimension): ---------- value : string, float, int, or :class:`unyt.unyt_quantity` The value being tested. Should be something like - + >>> _validate_quantity("10 MW", dimension='power') unyt_quantity(10., 'MW') - + dimension : string The expected dimensions of `value`. Currently accepts: ['time', 'energy', 'power', 'spec_power', 'spec_energy']. @@ -111,12 +113,27 @@ def _validate_quantity(value, dimension): class Technology(object): """ + The :class:`Technology` base class contains the minimum required + data to solve an energy systems problem. All other technologies in + :mod:`osier` inherit from this class. + Parameters ---------- - technology_name : string + technology_name : str The name identifier of the technology. - technology_type : string + technology_type : str The string identifier for the type of technology. + Two common types are: ["production", "storage"]. + technology_category : str + The string identifier the the technology category. + For example: "renewable," "fossil," or "nuclear." + dispatchable : bool + Indicates whether the technology can be dispatched by a + grid operator, or if it produces variable electricity + that must be used or stored the moment it is produced. + For example, solar panels and wind turbines are not + dispatchable, but nuclear and biopower are dispatchable. + Default value is true. capital_cost : float or :class:`unyt.array.unyt_quantity` Specifies the capital cost. If float, the default unit is $/MW. @@ -168,7 +185,9 @@ class Technology(object): def __init__(self, technology_name, - technology_type='base', + technology_type='production', + technology_category='base', + dispatchable=True, capital_cost=0.0, om_cost_fixed=0.0, om_cost_variable=0.0, @@ -180,6 +199,8 @@ def __init__(self, self.technology_name = technology_name self.technology_type = technology_type + self.technology_category = technology_category + self.dispatchable = dispatchable self.unit_power = default_power_units self.unit_time = default_time_units @@ -209,7 +230,7 @@ def unit_time(self, value): @property def unit_energy(self): - return self._unit_power*self._unit_time + return self._unit_power * self._unit_time @unit_energy.setter def unit_energy(self, value): @@ -256,3 +277,41 @@ def fuel_cost(self): @fuel_cost.setter def fuel_cost(self, value): self._fuel_cost = _validate_quantity(value, dimension="spec_energy") + + @property + def total_capital_cost(self): + return self.capacity * self.capital_cost + + @property + def annual_fixed_cost(self): + return self.capacity * self.om_cost_fixed + + @property + def variable_cost(self): + return self.fuel_cost + self.om_cost_variable + + def variable_cost_ts(self, size): + """ + Returns the total variable cost as a time series of + length :attr:`size`. + + .. warning:: + The current implementation assumes a single constant cost + for the variable cost. In the future, users will be able to + pass their own time series data. + + Parameters + ---------- + size : int + The number of periods, i.e. length, of the + time series. + + Returns + ------- + var_cost_ts : :class:`numpy.ndarray` + The variable cost time series. + """ + + # assumes single cost + var_cost_ts = np.ones(size) * self.variable_cost + return var_cost_ts diff --git a/setup.py b/setup.py index ed9804f..a6b2033 100644 --- a/setup.py +++ b/setup.py @@ -58,7 +58,8 @@ 'dill', 'openpyxl', 'nrelpy', - 'unyt'] + 'unyt', + 'pyomo'] EXTRAS_REQUIRE = { 'doc': [ 'sphinx>=5.1', @@ -66,7 +67,8 @@ "sphinx_design", "sphinx-autodoc-typehints", 'numpydoc', - 'pydata_sphinx_theme', ]} + 'pydata_sphinx_theme' + ]} PYTHON_REQUIRES = ">= 3.6" PACKAGES = find_packages() diff --git a/tests/test_models.py b/tests/test_models.py new file mode 100644 index 0000000..4b219fd --- /dev/null +++ b/tests/test_models.py @@ -0,0 +1,77 @@ +from osier import DispatchModel +from osier import Technology +from unyt import unyt_array +import numpy as np +import pytest +import sys + +if "win32" in sys.platform: + solver = 'cplex' +elif "linux" in sys.platform: + solver = "cbc" +else: + solver = "cbc" + +TOL = 1e-5 + +@pytest.fixture +def technology_set(): + nuclear = Technology(technology_name='Nuclear', + technology_type='production', + capacity=1e3, + capital_cost=6, + om_cost_variable=20, + om_cost_fixed=50, + fuel_cost=5 + ) + natural_gas = Technology(technology_name='NaturalGas', + technology_type='production', + capacity=1e3, + capital_cost=1, + om_cost_variable=12, + om_cost_fixed=30, + fuel_cost=20 + ) + + return [nuclear, natural_gas] + + +@pytest.fixture +def net_demand(): + n_hours = 24 + np.random.seed(123) + x = np.arange(0, n_hours, 1) + y = np.sin(8 * x * np.pi / 180) + \ + np.random.normal(loc=0, scale=0.1, size=n_hours) + y[y < 0] = 0 + return y + + +def test_dispatch_model_initialize(technology_set, net_demand): + """ + Tests that the dispatch model is properly initialized. + """ + model = DispatchModel(technology_set, + net_demand=net_demand, + solver=solver) + assert model.technology_list == technology_set + assert model.tech_set == [t.technology_name for t in technology_set] + assert model.solver == solver + assert len(model.capacity_dict) == len(technology_set) + assert len(model.indices) == len(net_demand) * len(technology_set) + + +def test_dispatch_model_solve(technology_set, net_demand): + """ + Tests that the dispatch model produces expected results. + """ + model = DispatchModel(technology_set, + net_demand=net_demand, + solver=solver) + model.solve() + cheapest_tech = unyt_array([t.variable_cost for t in technology_set]).min() + expected_result = cheapest_tech * net_demand.sum() + + assert model.objective == pytest.approx(expected_result, TOL) + assert model.results['Nuclear'].sum() == pytest.approx(net_demand.sum(), TOL) + assert model.results['NaturalGas'].sum() == pytest.approx(0.0, TOL) diff --git a/tests/test_technology.py b/tests/test_technology.py index ec8cb24..363f4fe 100644 --- a/tests/test_technology.py +++ b/tests/test_technology.py @@ -25,11 +25,13 @@ dict_type = {"value": 10, "unit": MW} + @pytest.fixture def advanced_tech(): PLANET_EXPRESS = Technology(TECH_NAME) return PLANET_EXPRESS + def test_validate_unit(): assert _validate_unit("MW", 'power').same_dimensions_as(Horsepower) assert _validate_unit("BTU", 'energy').same_dimensions_as(MW * hr) @@ -70,7 +72,7 @@ def test_validate_quantity(): def test_initialize(advanced_tech): assert advanced_tech.technology_name == TECH_NAME - assert advanced_tech.technology_type == 'base' + assert advanced_tech.technology_type == 'production' assert advanced_tech.capacity == 0.0 assert advanced_tech.capital_cost == 0.0 assert advanced_tech.om_cost_fixed == 0.0 @@ -79,6 +81,23 @@ def test_initialize(advanced_tech): assert advanced_tech.unit_power == MW assert advanced_tech.unit_time == hr assert advanced_tech.unit_energy == MW * hr + assert advanced_tech.annual_fixed_cost == 0.0 + assert advanced_tech.total_capital_cost == 0.0 + + +def test_total_capital_cost(advanced_tech): + advanced_tech.capital_cost = spec_power_unyt + advanced_tech.capacity = power_unyt + assert advanced_tech.total_capital_cost == 100 + + advanced_tech.capacity = 0.5 * power_unyt + assert advanced_tech.total_capital_cost == 50 + + +def test_variable_cost(advanced_tech): + advanced_tech.fuel_cost = spec_energy_unyt + advanced_tech.om_cost_variable = spec_energy_unyt + assert advanced_tech.variable_cost == 2 * spec_energy_unyt def test_attribute_types(advanced_tech): @@ -219,7 +238,7 @@ def test_om_cost_variable(advanced_tech): advanced_tech.unit_power = "kW" advanced_tech.unit_time = "day" - assert advanced_tech.om_cost_variable.units == (kW*day)**-1 + assert advanced_tech.om_cost_variable.units == (kW * day)**-1 def test_fuel_cost(advanced_tech): @@ -248,10 +267,10 @@ def test_fuel_cost(advanced_tech): assert advanced_tech.fuel_cost.value == pytest.approx( 3412141.47989694, 0.5) assert advanced_tech.fuel_cost.units == (MW * hr)**-1 - + advanced_tech.unit_power = "kW" advanced_tech.unit_time = "day" - assert advanced_tech.fuel_cost.units == (kW*day)**-1 + assert advanced_tech.fuel_cost.units == (kW * day)**-1 def test_unit_power(advanced_tech): @@ -288,14 +307,14 @@ def test_unit_time(advanced_tech): def test_unit_energy(advanced_tech): advanced_tech.unit_energy = "darkmatter" - assert advanced_tech.unit_energy == MW*hr + assert advanced_tech.unit_energy == MW * hr advanced_tech.unit_energy = BTU - assert advanced_tech.unit_energy == MW*hr + assert advanced_tech.unit_energy == MW * hr advanced_tech.unit_energy = "MW" - assert advanced_tech.unit_energy == MW*hr + assert advanced_tech.unit_energy == MW * hr advanced_tech.unit_energy = 10 - assert advanced_tech.unit_energy == MW*hr + assert advanced_tech.unit_energy == MW * hr advanced_tech.unit_energy = Horsepower * day - assert advanced_tech.unit_energy == MW*hr + assert advanced_tech.unit_energy == MW * hr advanced_tech.unit_energy = "Horsepower*day" - assert advanced_tech.unit_energy == MW*hr + assert advanced_tech.unit_energy == MW * hr