Introduce nomopyomo

.travis.yml

@@ -6,8 +6,8 @@ sudo: false # Use container-based infrastructure
 
 matrix:
   include:
-    - env:
-      - PYTHON_VERSION="2.7"
+    # - env:
+    #   - PYTHON_VERSION="2.7"
     - env:
       - PYTHON_VERSION="3.6"
     - env:
@@ -23,7 +23,7 @@ before_install:
   - export PATH="$HOME/miniconda/bin:$PATH"
   - hash -r
   - conda config --set always_yes yes --set changeps1 no
-  - conda update -q conda
+  # - conda update -q conda
   # Useful for debugging any issues with conda
   - conda info -a
   - source $HOME/miniconda/etc/profile.d/conda.sh

MANIFEST.in

@@ -1,5 +1,6 @@
 include pypsa/component_attrs/*.csv
 include pypsa/standard_types/*.csv
 include pypsa/components.csv
+include pypsa/variables.csv
 include README.rst LICENSE.txt
 include requirements.yml

environment.yaml

@@ -14,3 +14,4 @@ dependencies:
   - cartopy>=0.16
   - coincbc
   - glpk
+  - gurobi::gurobi

pypsa/__init__.py

@@ -26,7 +26,15 @@
 from __future__ import absolute_import
 
 from . import components, descriptors
-from . import pf, opf, plot, networkclustering, io, contingency, geo
+from . import pf, opf, opt, plot, networkclustering, io, contingency, geo, stats
+
+import sys
+
+#do this as long as python 2.7 should be supported
+if sys.version_info.major >= 3:
+    from . import linopf, linopt
+
+
 
 from .components import Network, SubNetwork
 

pypsa/component_attrs/generators.csv

@@ -30,3 +30,8 @@ p,series,MW,0.,active power at bus (positive if net generation),Output
 q,series,MVar,0.,reactive power (positive if net generation),Output
 p_nom_opt,float,MW,0.,Optimised nominal power.,Output
 status,series,n/a,1,"Status (1 is on, 0 is off). Only outputted if committable is True.",Output
+mu_upper,series,currency/MWh,0.,Shadow price of upper p_nom limit,Output
+mu_lower,series,currency/MWh,0.,Shadow price of lower p_nom limit,Output
+mu_p_set,series,currency/MWh,0.,Shadow price of fixed power generation p_set,Output
+mu_ramp_limit_up,series,currency/MWh,0.,Shadow price of upper ramp up limit,Output
+mu_ramp_limit_down,series,currency/MWh,0.,Shadow price of lower ramp down limit,Output

pypsa/component_attrs/links.csv

@@ -19,4 +19,5 @@ p0,series,MW,0.,Active power at bus0 (positive if branch is withdrawing power fr
 p1,series,MW,0.,Active power at bus1 (positive if branch is withdrawing power from bus1).,Output
 p_nom_opt,float,MVA,0.,Optimised capacity for active power.,Output
 mu_lower,series,currency/MVA,0.,Shadow price of lower p_nom limit  -F \leq f. Always non-negative.,Output
-mu_upper,series,currency/MVA,0.,Shadow price of upper p_nom limit f \leq F. Always non-negative.,Output
+mu_upper,series,currency/MVA,0.,Shadow price of upper p_nom limit f \leq F. Always non-negative.,Output
+mu_p_set,series,currency/MWh,0.,Shadow price of fixed power transmission p_set,Output

pypsa/component_attrs/storage_units.csv

@@ -24,7 +24,12 @@ efficiency_dispatch,float,per unit,1.,Efficiency of storage on the way out of th
 standing_loss,float,per unit,0.,Losses per hour to state of charge.,Input (optional)
 inflow,static or series,MW,0.,"Inflow to the state of charge, e.g. due to river inflow in hydro reservoir.",Input (optional)
 p,series,MW,0.,active power at bus (positive if net generation),Output
+p_dispatch,series,MW,0.,active power dispatch at bus,Output
+p_store,series,MW,0.,active power charging at bus,Output
 q,series,MVar,0.,reactive power (positive if net generation),Output
 state_of_charge,series,MWh,NaN,State of charge as calculated by the OPF.,Output
 spill,series,MW,0.,Spillage for each snapshot.,Output
 p_nom_opt,float,MW,0.,Optimised nominal power.,Output
+mu_upper,series,currency/MWh,0.,Shadow price of upper p_nom limit,Output
+mu_lower,series,currency/MWh,0.,Shadow price of lower p_nom limit,Output
+mu_state_of_charge_set,series,currency/MWh,0.,Shadow price of fixed state of charge state_of_charge_set,Output

pypsa/component_attrs/stores.csv

@@ -19,4 +19,7 @@ standing_loss,float,per unit,0.,Losses per hour to energy.,Input (optional)
 p,series,MW,0.,active power at bus (positive if net generation),Output
 q,series,MVar,0.,reactive power (positive if net generation),Output
 e,series,MWh,0.,Energy as calculated by the OPF.,Output
-e_nom_opt,float,MW,0.,Optimised nominal energy capacity outputed by OPF.,Output
+e_nom_opt,float,MW,0.,Optimised nominal energy capacity outputed by OPF.,Output
+mu_upper,series,currency/MWh,0.,Shadow price of upper e_nom limit,Output
+mu_lower,series,currency/MWh,0.,Shadow price of lower e_nom limit,Output
+mu_e_set,series,currency/MWh,0.,Shadow price of fixed energy level e_set,Output

pypsa/components.py

@@ -21,24 +21,18 @@
 
 # make the code as Python 3 compatible as possible
 from __future__ import division, absolute_import
-import six
 from six import iteritems, itervalues, iterkeys
-from six.moves import map
 from weakref import ref
 
 
 __author__ = "Tom Brown (FIAS), Jonas Hoersch (FIAS), David Schlachtberger (FIAS)"
 __copyright__ = "Copyright 2015-2017 Tom Brown (FIAS), Jonas Hoersch (FIAS), David Schlachtberger (FIAS), GNU GPL 3"
 
-import networkx as nx
 
 import numpy as np
 import pandas as pd
-import scipy as sp, scipy.sparse
 from scipy.sparse import csgraph
-from itertools import chain
 from collections import namedtuple
-from operator import itemgetter
 import os
 
 
@@ -71,10 +65,11 @@
 
 from .graph import graph, incidence_matrix, adjacency_matrix
 
-import inspect
-
 import sys
 
+if sys.version_info.major >= 3:
+    from .linopf import network_lopf as network_lopf_lowmem
+
 import logging
 logger = logging.getLogger(__name__)
 
@@ -198,7 +193,7 @@ class Network(Basic):
 
     pf = network_pf
 
-    lopf = network_lopf
+#    lopf = network_lopf
 
     opf = network_opf
 
@@ -407,6 +402,102 @@ def set_snapshots(self,snapshots):
 
         #NB: No need to rebind pnl to self, since haven't changed it
 
+    def lopf(self, snapshots=None, solver_name="glpk", extra_functionality=None,
+             solver_logfile=None, solver_options={}, keep_files=False,
+             formulation="kirchhoff", extra_postprocessing=None, pyomo=True,
+             **kwargs):
+        """
+        Linear optimal power flow for a group of snapshots.
+
+        Parameters
+        ----------
+        snapshots : list or index slice
+            A list of snapshots to optimise, must be a subset of
+            network.snapshots, defaults to network.snapshots
+        solver_name : string
+            Must be a solver name that pyomo recognises and that is
+            installed, e.g. "glpk", "gurobi"
+        pyomo : bool, default True
+            Whether to use pyomo for building and solving the model, setting
+            this to False saves a lot of memory and time.
+        solver_logfile : None|string
+            If not None, sets the logfile option of the solver.
+        solver_options : dictionary
+            A dictionary with additional options that get passed to the solver.
+            (e.g. {'threads':2} tells gurobi to use only 2 cpus)
+        keep_files : bool, default False
+            Keep the files that pyomo constructs from OPF problem
+            construction, e.g. .lp file - useful for debugging
+        formulation : string
+            Formulation of the linear power flow equations to use; must be
+            one of ["angles","cycles","kirchhoff","ptdf"]
+        extra_functionality : callable function
+            This function must take two arguments
+            `extra_functionality(network,snapshots)` and is called after
+            the model building is complete, but before it is sent to the
+            solver. It allows the user to
+            add/change constraints and add/change the objective function.
+
+        Other Parameters
+        ----------------
+
+        ptdf_tolerance : float
+            Only when pyomo is True.
+            Value below which PTDF entries are ignored
+        free_memory : set, default {'pyomo'}
+            Only when pyomo is True.
+            Any subset of {'pypsa', 'pyomo'}. Allows to stash `pypsa` time-series
+            data away while the solver runs (as a pickle to disk) and/or free
+            `pyomo` data after the solution has been extracted.
+        solver_io : string, default None
+            Only when pyomo is True.
+            Solver Input-Output option, e.g. "python" to use "gurobipy" for
+            solver_name="gurobi"
+        skip_pre : bool, default False
+            Only when pyomo is True.
+            Skip the preliminary steps of computing topology, calculating
+            dependent values and finding bus controls.
+        extra_postprocessing : callable function
+            This function must take three arguments
+            `extra_postprocessing(network,snapshots,duals)` and is called after
+            the model has solved and the results are extracted. It allows the user
+            to extract further information about the solution, such as additional
+            shadow prices.
+        warmstart : bool or string, default False
+            Only when pyomo is False.
+            Use this to warmstart the optimization. Pass a string which gives
+            the path to the basis file. If set to True, a path to
+            a basis file must be given in network.basis_fn.
+        store_basis : bool, default True
+            Only when pyomo is False.
+            Whether to store the basis of the optimization results. If True,
+            the path to the basis file is saved in network.basis_fn. Note that
+            a basis can only be stored if simplex, dual-simplex, or barrier
+            *with* crossover is used for solving.
+        keep_references : bool, default False
+            Only when pyomo is False.
+            Keep the references of variable and constraint names withing the
+            network, e.g. n.generators_t.p_varref - useful for constructing
+            extra_functionality or debugging
+        keep_shadowprices : bool or list of component names, default None
+            Only when pyomo is False.
+            Keep shadow prices for all constraints, if set to True.
+            These are stored at e.g. n.generators_t.mu_upper for upper limit
+            of p_nom. If a list of component names is passed, shadow
+            prices of variables attached to those are extracted. If set to None,
+            components default to ['Bus', 'Line', 'GlobalConstraint']
+
+        """
+        args = {'snapshots': snapshots, 'keep_files': keep_files,
+                'solver_options': solver_options, 'formulation': formulation,
+                'extra_functionality': extra_functionality,
+                'solver_name': solver_name, 'solver_logfile': solver_logfile}
+        args.update(kwargs)
+        if pyomo:
+            return network_lopf(self, **args)
+        else:
+            return network_lopf_lowmem(self, **args)
+
 
 
     def add(self, class_name, name, **kwargs):
@@ -834,6 +925,8 @@ def determine_network_topology(self):
 
         for sub in self.sub_networks.obj:
             find_cycles(sub)
+            sub.find_bus_controls()
+
 
     def iterate_components(self, components=None, skip_empty=True):
         if components is None:

pypsa/descriptors.py

@@ -22,19 +22,15 @@
 # make the code as Python 3 compatible as possible
 from __future__ import division
 from __future__ import absolute_import
-from six import iteritems, string_types
+from six import iteritems
 
 
 __author__ = "Tom Brown (FIAS), Jonas Hoersch (FIAS)"
 __copyright__ = "Copyright 2015-2017 Tom Brown (FIAS), Jonas Hoersch (FIAS), GNU GPL 3"
 
 
-
-
-
 #weak references are necessary to make sure the key-value pair are
 #destroyed if the key object goes out of scope
-from weakref import WeakKeyDictionary
 
 from collections import OrderedDict
 from itertools import repeat
@@ -303,3 +299,73 @@ def zsum(s, *args, **kwargs):
     Meant to be set as pd.Series.zsum = zsum.
     """
     return 0 if s.empty else s.sum(*args, **kwargs)
+
+#Perhaps this should rather go into components.py
+nominal_attrs = {'Generator': 'p_nom',
+                 'Line': 's_nom',
+                 'Transformer': 's_nom',
+                 'Link': 'p_nom',
+                 'Store': 'e_nom',
+                 'StorageUnit': 'p_nom'}
+
+def expand_series(ser, columns):
+    """
+    Helper function to quickly expand a series to a dataframe with according
+    column axis and every single column being the equal to the given series.
+    """
+    return ser.to_frame(columns[0]).reindex(columns=columns).ffill(axis=1)
+
+
+def get_extendable_i(n, c):
+    """
+    Getter function. Get the index of extendable elements of a given component.
+    """
+    return n.df(c)[lambda ds:
+        ds[nominal_attrs[c] + '_extendable']].index
+
+def get_non_extendable_i(n, c):
+    """
+    Getter function. Get the index of non-extendable elements of a given
+    component.
+    """
+    return n.df(c)[lambda ds:
+            ~ds[nominal_attrs[c] + '_extendable']].index
+
+def get_bounds_pu(n, c, sns, index=slice(None), attr=None):
+    """
+    Getter function to retrieve the per unit bounds of a given compoent for
+    given snapshots and possible subset of elements (e.g. non-extendables).
+    Depending on the attr you can further specify the bounds of the variable
+    you are looking at, e.g. p_store for storage units.
+
+    Parameters
+    ----------
+    n : pypsa.Network
+    c : string
+        Component name, e.g. "Generator", "Line".
+    sns : pandas.Index/pandas.DateTimeIndex
+        set of snapshots for the bounds
+    index : pd.Index, default None
+        Subset of the component elements. If None (default) bounds of all
+        elements are returned.
+    attr : string, default None
+        attribute name for the bounds, e.g. "p", "s", "p_store"
+
+    """
+    min_pu_str = nominal_attrs[c].replace('nom', 'min_pu')
+    max_pu_str = nominal_attrs[c].replace('nom', 'max_pu')
+
+    max_pu = get_switchable_as_dense(n, c, max_pu_str, sns)
+    if c in n.passive_branch_components:
+        min_pu = - max_pu
+    elif c == 'StorageUnit':
+        min_pu = pd.DataFrame(0, max_pu.index, max_pu.columns)
+        if attr == 'p_store':
+            max_pu = - get_switchable_as_dense(n, c, min_pu_str, sns)
+        if attr == 'state_of_charge':
+            max_pu = expand_series(n.df(c).max_hours, sns).T
+            min_pu = pd.DataFrame(0, *max_pu.axes)
+    else:
+        min_pu = get_switchable_as_dense(n, c, min_pu_str, sns)
+    return min_pu[index], max_pu[index]
+