Include angle power flow formulation in MarketModel.jl

LICENSE.md

@@ -1,4 +1,4 @@
-Copyright (c) 2018: Richard Weinhold & Robert Mieth.
+Copyright (c) 2020: Richard Weinhold & Robert Mieth.
 
 The Pomato Python module and its adjacent software (i.e. code written in programming languages other than Python which the Python module utilizes) as well as provided documentation is licensed under the **[LGPL]** version 3:
 

README.md

@@ -1,4 +1,4 @@
-[![PomatoLogo](https://github.com/richard-weinhold/pomato/blob/master/docs/_static/graphics/pomato_logo_small.png "Pomato Soup")](#) POMATO - Power Market Tool
+[![PomatoLogo](https://github.com/richard-weinhold/pomato/blob/master/docs/_static/graphics/pomato_480.png "Pomato Soup")](#) POMATO - Power Market Tool
 =====================================================================================================================================
 [![Documentation Status](https://readthedocs.org/projects/pomato/badge/?version=latest)](https://pomato.readthedocs.io/en/latest/?badge=latest)
 [![Build Status](https://travis-ci.org/richard-weinhold/pomato.svg?branch=master)](https://travis-ci.org/richard-weinhold/pomato)

pomato/grid/grid_model.py

@@ -81,6 +81,7 @@ def __init__(self, wdir, grid, data, option):
                                                          multiple_slack=None,
                                                          slack_zones=None,
                                                          grid=pd.DataFrame(),
+                                                         contingency_groups={},
                                                          redispatch_grid=pd.DataFrame(),
                                                          ntc=pd.DataFrame())
 
@@ -118,19 +119,19 @@ def create_grid_representation(self):
             self.process_ntc()
         elif self.options["optimization"]["type"] == "nodal":
             self.grid_representation.grid = self.create_nodal_grid_parameters()
-            self.grid_representation.grid = self._add_zone_to_grid_representation(self.grid_representation.grid)
         elif self.options["optimization"]["type"] == "zonal":
             self.grid_representation.grid = self.create_zonal_grid_parameters()
             self.grid_representation.ntc = self.create_ntc()
         elif self.options["optimization"]["type"] == "cbco_nodal":
+            self.grid_representation.contingency_groups = self.grid.contingency_groups
             self.grid_representation.grid = self.create_cbco_nodal_grid_parameters()
-            self.grid_representation.grid = self._add_zone_to_grid_representation(self.grid_representation.grid)
         elif self.options["optimization"]["type"] == "cbco_zonal":
+            self.grid_representation.contingency_groups = self.grid.contingency_groups
             self.grid_representation.grid = self.create_cbco_zonal_grid_parameters()
             self.grid_representation.ntc = self.create_ntc()
         else:
             self.logger.info("No grid representation needed for dispatch model")
-
+        
         if self.options["optimization"]["redispatch"]["include"]:
             self.add_redispatch_grid()
 
@@ -146,20 +147,17 @@ def create_nodal_grid_parameters(self):
         grid_option = self.options["grid"]
         nodal_network = pd.DataFrame(columns=self.grid.nodes.index, data=self.grid.ptdf)
         nodal_network["ram"] = self.grid.lines.maxflow.values*self.options["grid"]["capacity_multiplier"]
+        nodal_network["cb"] = list(self.grid.lines.index)
+        nodal_network["co"] = ["basecase" for i in range(0, len(self.grid.lines.index))]
+        nodal_network = nodal_network[["cb", "co", "ram"] + list(self.grid.nodes.index)]
 
         if grid_option["cbco_option"] == "nodal_clarkson":
-            nodal_network["cb"] = list(self.grid.lines.index)
-            nodal_network["co"] = ["basecase" for i in range(0, len(self.grid.lines.index))]
-            nodal_network = nodal_network[["cb", "co", "ram"] + list(self.grid.nodes.index)]
             nodal_injection_limits = self.create_nodal_injection_limits()
 
             cbco_index = self.clarkson_algorithm(A=nodal_network.loc[:, self.grid.nodes.index].values, 
                                                  b=nodal_network.loc[:, "ram"].values, 
                                                  x_bounds=nodal_injection_limits)
             nodal_network = self.return_cbco(nodal_network, cbco_index)
-        else:
-            nodal_network.set_index(self.grid.lines.index, inplace=True)
-
 
         return nodal_network
 
@@ -170,10 +168,8 @@ def add_redispatch_grid(self, contingencies=False):
         """
         if contingencies:
             self.grid_representation.redispatch_grid = self.create_cbco_nodal_grid_parameters()
-            self.grid_representation.redispatch_grid = self._add_zone_to_grid_representation(self.grid_representation.redispatch_grid)
         else:
             self.grid_representation.redispatch_grid = self.create_nodal_grid_parameters()
-            self.grid_representation.redispatch_grid = self._add_zone_to_grid_representation(self.grid_representation.redispatch_grid)
 
     def create_zonal_grid_parameters(self):
         """Process grid information for zonal N-0 representation.
@@ -194,17 +190,15 @@ def create_zonal_grid_parameters(self):
         zonal_network = pd.DataFrame(index=self.grid.lines.index,
                                      columns=self.data.zones.index,
                                      data=np.dot(self.grid.ptdf, self.create_gsk(gsk)))
+        zonal_network["cb"] = list(self.grid.lines.index)
+        zonal_network["co"] = ["basecase" for i in range(0, len(self.grid.lines.index))]
         zonal_network["ram"] = self.grid.lines.maxflow.values*self.options["grid"]["capacity_multiplier"]
+        zonal_network = zonal_network[["cb", "co", "ram"] + list(self.data.zones.index)]
 
         if grid_option["cbco_option"] == "clarkson":
-            zonal_network["cb"] = list(self.grid.lines.index)
-            zonal_network["co"] = ["basecase" for i in range(0, len(self.grid.lines.index))]
-            zonal_network = zonal_network[["cb", "co", "ram"] + list(self.data.zones.index)]
             cbco_index = self.clarkson_algorithm(A=zonal_network.loc[:, self.data.zones.index].values, 
                                                  b=zonal_network.loc[:, "ram"].values)
             zonal_network = self.return_cbco(zonal_network, cbco_index)
-        else:
-            zonal_network.set_index(self.grid.lines.index, inplace=True)
 
         return zonal_network
 

pomato/grid/grid_topology.py

@@ -62,6 +62,9 @@ class GridTopology():
         psdf (phase shifting distribution factor) matrix :math:`(L \\times L)`.
     lodf : np.ndarray
         N-1 lodf (load outage distribution factor) matrix :math:`(L \\times L)`.
+    contingency_groups : dict
+        Dictionary that indicates which group of lines is outate, when a line is considered a 
+        contingency.
     """
     numpy_settings = np.seterr(divide="raise")
     def __init__(self):
@@ -73,7 +76,7 @@ def __init__(self):
         self.psdf = None
         self.multiple_slack = False
         self.lodf = None
-        self.combined_contingencies = None
+        self.contingency_groups = None
 
     def calculate_parameters(self, nodes, lines):        
         self.logger.info("Calculating Grid Parameters!")
@@ -88,7 +91,7 @@ def calculate_parameters(self, nodes, lines):
         self.check_grid_topology()
         self.logger.info("Calculating LODF Matrix!")
         self.lodf = self.create_n_1_lodf_matrix()
-        self.combined_contingencies = self.create_contingency_groups()
+        self.contingency_groups = self.create_contingency_groups()
         self.logger.info("Grid parameters Calculated!")
 
     def check_slack(self):
@@ -180,30 +183,30 @@ def create_contingency_groups(self, option="double_lines"):
             Rule of what outages are grouped. Options are: double_lines (default), by value or none.
 
         """
-        combined_contingencies = {line : [line] for line in self.lines.index}
+        contingency_groups = {line : [line] for line in self.lines.index}
 
         if isinstance(option, (int, float)):
-            combined_contingencies = {line : [line] for line in self.lines.index}
+            contingency_groups = {line : [line] for line in self.lines.index}
             for idx, line in enumerate(self.lines.index):
-                combined_contingencies[line].extend(list(self.lines.index[np.abs(self.lodf[idx, :]) > option]))
-                combined_contingencies[line] = list(set(combined_contingencies[line]))    
+                contingency_groups[line].extend(list(self.lines.index[np.abs(self.lodf[idx, :]) > option]))
+                contingency_groups[line] = list(set(contingency_groups[line]))    
 
         if option == "double_lines": # double lines
             if "systems" not in self.lines.columns:
                 self.add_number_of_systems()
 
-        double_lines = list(self.lines[self.lines.systems == 2].index)
-        for line in double_lines:
-            condition = (self.lines.loc[double_lines, ["node_i", "node_j"]].apply(tuple, axis=1) == tuple(self.lines.loc[line, ["node_i", "node_j"]])).values
-            double_line = list(self.lines.loc[double_lines][condition].index)
-            double_line_idx = [self.lines.index.get_loc(line) for line in double_line]
-
-            # However if the double line is radial, do not consider a combined outage
-            if not any(np.sum(np.around(np.abs(self.ptdf[double_line_idx, :]), decimals=3), axis=0) == 1):
-                combined_contingencies[line] = double_line
-            # combined_contingencies[line] = list(set(combined_contingencies[line]))
+            double_lines = list(self.lines[self.lines.systems == 2].index)
+            for line in double_lines:
+                condition = (self.lines.loc[double_lines, ["node_i", "node_j"]].apply(tuple, axis=1) == tuple(self.lines.loc[line, ["node_i", "node_j"]])).values
+                double_line = list(self.lines.loc[double_lines][condition].index)
+                double_line_idx = [self.lines.index.get_loc(line) for line in double_line]
+                
+                # However if the double line is radial, do not consider a combined outage
+                if not any(np.sum(np.around(np.abs(self.ptdf[double_line_idx, :]), decimals=3), axis=0) == 1):
+                    contingency_groups[line] = double_line
+                # contingency_groups[line] = list(set(contingency_groups[line]))
 
-        return combined_contingencies
+        return contingency_groups
 
 
     def slack_zones(self):
@@ -493,7 +496,7 @@ def create_n_1_ptdf_outage(self, outages):
 
         combines_outages = []
         for outage in outages:
-            tmp = self.combined_contingencies[self.lines.index[outage]]
+            tmp = self.contingency_groups[self.lines.index[outage]]
             combines_outages.extend([self.lines.index.get_loc(line) for line in tmp])
 
         lodf = self.create_lodf([line for line in range(0, len(self.lines))], combines_outages)
@@ -522,7 +525,7 @@ def create_n_1_ptdf_cbco(self, line, outage):
         if not isinstance(line, int):
             line = self.lines.index.get_loc(line)
 
-        outages = [self.lines.index.get_loc(line) for line in self.combined_contingencies[self.lines.index[outage]]]
+        outages = [self.lines.index.get_loc(line) for line in self.contingency_groups[self.lines.index[outage]]]
 
         n_1_ptdf_cbco = self.ptdf[line, :] + np.dot(self.create_lodf([line], outages), self.ptdf[outages, :])
         return n_1_ptdf_cbco

pomato/market_model/market_model.py

@@ -187,7 +187,7 @@ def data_to_csv(self, model_horizon):
 
         self.rolling_horizon_storage_levels(model_horizon)
 
-        for data in [d for d in self.data.model_structure if d != "lines"]:
+        for data in [d for d in self.data.model_structure]:
             cols = [col for col in self.data.model_structure[data].keys() if col != "index"]
             if "timestep" in cols:
                 getattr(self.data, data).loc[getattr(self.data, data)["timestep"].isin(model_horizon), cols] \
@@ -203,17 +203,20 @@ def data_to_csv(self, model_horizon):
         plant_types.to_csv(str(self.data_dir.joinpath('plant_types.csv')), index_label='index')
 
         if self.grid_representation.grid.empty:
-            pd.DataFrame(columns=["ram"]).to_csv(str(self.data_dir.joinpath('grid.csv')), index_label='index')
+            pd.DataFrame(columns=["cb", "co", "ram"]).to_csv(str(self.data_dir.joinpath('grid.csv')), index_label='index')
         else:
             self.grid_representation.grid \
                 .to_csv(str(self.data_dir.joinpath('grid.csv')), index_label='index')
 
         if self.grid_representation.redispatch_grid.empty:
-            pd.DataFrame(columns=["ram"]).to_csv(str(self.data_dir.joinpath('redispatch_grid.csv')), index_label='index')
+            pd.DataFrame(columns=["cb", "co", "ram"]).to_csv(str(self.data_dir.joinpath('redispatch_grid.csv')), index_label='index')
         else:
             self.grid_representation.redispatch_grid \
                 .to_csv(str(self.data_dir.joinpath('redispatch_grid.csv')), index_label='index')
 
+        with open(self.data_dir.joinpath('contingency_groups.json'), 'w') as file:
+            json.dump(self.grid_representation.contingency_groups, file, indent=2)
+
         if not self.grid_representation.ntc.empty:
             self.grid_representation.ntc.to_csv(str(self.data_dir.joinpath('ntc.csv')), index_label='index')
 

setup.py

@@ -40,7 +40,7 @@ def run(self):
         julia_instantiate(self.install_lib)
 
 setup(name='pomato',
-      version='0.2.1',
+      version='0.2.5',
       description='Power Market Tool',
       author='Richard Weinhold',
       author_email='[email protected]',

tests/test_grid_model.py

@@ -69,7 +69,8 @@ def test_zonal(self):
         self.assertRaises(AssertionError, np.testing.assert_equal, 
                           grid_representation_flat.grid.values, grid_representation_gmax.grid.values)
 
-        test_columns = list(self.grid_model.data.zones.index) + ["ram"]
+        test_columns = ["cb", "co", "ram"] + list(self.grid_model.data.zones.index)
+        print(grid_representation_flat.grid.columns)
         self.assertTrue(all(grid_representation_flat.grid.columns == test_columns))
         self.assertTrue(all(grid_representation_gmax.grid.columns == test_columns))
 
@@ -104,7 +105,6 @@ def test_cbco_nodal_no_precalc(self):
         self.grid_model.options["optimization"]["type"] = "cbco_nodal"
         self.grid_model.options["grid"]["precalc_filename"] = "random_words"
         grid = self.grid_model.create_cbco_nodal_grid_parameters()
-        grid = self.grid_model._add_zone_to_grid_representation(grid)
         c_ptdf_fallback = copy.copy(grid)
         self.grid_model.options["grid"]["precalc_filename"] = ""
         self.grid_model.options["grid"]["cbco_option"] = "full"

tests/test_grid_topology.py

@@ -110,11 +110,11 @@ def test_multiple_contingencies(self):
 
         np.testing.assert_almost_equal(ptdf_l10l9, ptdf_9_10)
 
-    def test_combined_contingencies_sensitivity(self):
+    def test_contingency_groups_sensitivity(self):
         tmp = self.grid.create_contingency_groups(0.4)
         self.assertTrue(all([outage in tmp[outage] for outage in tmp]))
 
-    def test_combined_contingencies(self):
+    def test_contingency_groups(self):
         node_in = self.grid.nodes.index.get_loc("n10")
         node_out = self.grid.nodes.index.get_loc("n112")
 
@@ -123,9 +123,9 @@ def test_combined_contingencies(self):
         inj[node_out] = -100
 
         # pre_contingency_flow = np.dot(self.grid.ptdf, inj)
-        outages = [outage for outage in self.grid.combined_contingencies if len(self.grid.combined_contingencies[outage]) > 1]
+        outages = [outage for outage in self.grid.contingency_groups if len(self.grid.contingency_groups[outage]) > 1]
         outage = outages[0]
-        outages = self.grid.combined_contingencies[outage]
+        outages = self.grid.contingency_groups[outage]
 
         c_ptdf = self.grid.create_n_1_ptdf_outage(outage)
         post_contingency_flow = np.dot(c_ptdf, inj)