Merge pull request #77 from EnergyAssetOptimization/develop

added fill level output

.gitignore

@@ -30,3 +30,4 @@ tests/.~lock.plant_test_data.csv\#
 doc/source/samples/benchmark_problems/res_cross_commodity_benchmark.xlsx
 doc/source/samples/benchmark_problems/DK1_input_data.xlsx
 .vscode/launch.json
+doc/source/samples/grid and battery/~$test.xlsx

doc/source/samples/order_book/order_book_battery.ipynb

@@ -171,7 +171,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -187,7 +187,8 @@
     "                                        start_level = 20,\n",
     "                                        end_level   = 20,\n",
     "                                        eff_in      = efficiency,\n",
-    "                                        size        = 40)\n",
+    "                                        size        = 40,\n",
+    "                                        no_simult_in_out = True) # at negative prices, we want to ensure the battery does not charge & discharge at the same time to \"burn\" power\n",
     "# last resort - battery end level. May allow battery not to be completely full, \"borrowing\" in last hours\n",
     "extra_power = eao.assets.SimpleContract('fill_level_adjust', node,\n",
     "                                max_cap   = 10,\n",
@@ -246,7 +247,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [
     {
@@ -274,10 +275,12 @@
     "ax.plot([S, S], [0,0],'r-',linewidth = 5, alpha = 1, label = 'executed buy - we sell')\n",
     "ax.set_xlim(S, E)\n",
     "ax2 = ax.twinx() \n",
-    "# calculate fill level from dispatch\n",
-    "fill_level = -out['dispatch'].loc[:,'battery']\n",
-    "fill_level[fill_level>0] *= efficiency\n",
-    "fill_level = fill_level.cumsum()+20 \n",
+    "### show how to manually calculate fill level from dispatch\n",
+    "# fill_level = -out['dispatch'].loc[:,'battery']\n",
+    "# fill_level[fill_level>0] *= efficiency\n",
+    "# fill_level = fill_level.cumsum()+20 \n",
+    "### this is the automatic output\n",
+    "fill_level = out['internal_variables']['battery_fill_level']\n",
     "ax2.plot(fill_level, label = 'battery fill level')\n",
     "ax2.set_ylabel('fill level battery in MWh')\n",
     "ax.legend(loc = 'upper left')\n",

eaopack/assets.py

@@ -109,7 +109,7 @@ def dcf(self, optim_problem:OptimProblem, results:Results) -> np.array:
         dcf = np.zeros(self.timegrid.T)
         # filter for right asset in case larger problem is given
         my_mapping =  optim_problem.mapping.loc[optim_problem.mapping['asset']==self.name].copy()
-        # drop duplicate index - since mapping may contain several rows per varaible (indexes enumerate variables)
+        # drop duplicate index - since mapping may contain several rows per variable (indexes enumerate variables)
         my_mapping = pd.DataFrame(my_mapping[~my_mapping.index.duplicated(keep = 'first')])
 
         for i, r in my_mapping.iterrows():
@@ -355,7 +355,9 @@ def setup_optim_problem(self, prices: dict, timegrid:Timegrid = None, costs_only
                 price = np.asarray(myprice)
             else: # simply restrict prices to  asset time window
                 price           = price[self.timegrid.restricted.I]
-
+            # warning if there are neg. prices and no_simult_in_out is not True
+            if (not self.no_simult_in_out) and (self.eff_in <= 1) and (not (price >= 0).all()):
+                print('Storage --'+self.name+'--: no_simult_in_out is set to False, but there are neg. prices. Likely storage will load & unload simultaneously. You want that?')
         # separation into in/out needed?  Only one or two dispatch variables per time step
         # new separation reason: separate nodes in and out
         sep_needed =  (self.eff_in != 1) or (self.cost_in !=0) or (self.cost_out !=0) or (len(self.nodes)==2)
@@ -553,7 +555,46 @@ def setup_optim_problem(self, prices: dict, timegrid:Timegrid = None, costs_only
                                 periodic_period_length = self.periodicity,
                                 periodic_duration      = self.periodicity_duration,
                                 timegrid               = self.timegrid)
+
+    def fill_level(self, optim_problem:OptimProblem, results:Results) -> np.array:
+        """ Calculate discounted cash flow for the asset given the optimization results
+
+        Args:
+            optim_problem (OptimProblem): optimization problem created by this asset
+            results (Results): Results given by optimizer
 
+        Returns:
+            np.array: array with DCF per time step as per timegrid of asset
+        """
+        # for this asset simply from cost vector and optimal dispatch
+        # mapped to the full timegrid
+
+        ######### missing: mapping in optim problem
+        fill_level = np.zeros(self.timegrid.T)
+        # filter for right asset in case larger problem is given
+        my_mapping =  optim_problem.mapping.loc[(optim_problem.mapping['asset']==self.name) & (optim_problem.mapping['type']=='d')].copy()
+        # drop duplicate index - since mapping may contain several rows per variable (indexes enumerate variables)
+        my_mapping = pd.DataFrame(my_mapping[~my_mapping.index.duplicated(keep = 'first')])
+
+        # # for only one variable
+        # I_d  = my_mapping['var_name']=="disp"
+        # if sum(I_d) != 0:
+        #     fill_level =   np.maximum(0,-results.x[my_mapping.loc[I_d].index]*self.eff_in) \
+        #                  + np.minimum(0,-results.x[my_mapping.loc[I_d].index])
+        # else:
+        #     I_in  = my_mapping['var_name']=="disp_in"
+        #     I_out = my_mapping['var_name']=="disp_out"
+        #     fill_level = -results.x[my_mapping.loc[I_in].index]*self.eff_in \
+        #                 - results.x[my_mapping.loc[I_out].index]
+        # fill_level = fill_level.cumsum() + self.start_level
+
+        fill_level = np.zeros(self.timegrid.T)
+        for i, r in my_mapping.iterrows():
+            fill_level[r['time_step']] +=  max(0,-results.x[i])*self.eff_in \
+                                         + min(0,-results.x[i])
+        fill_level = fill_level.cumsum() + self.start_level            
+        return fill_level
+
 class SimpleContract(Asset):
     """ Contract Class """
     def __init__(self,

eaopack/io.py

@@ -2,10 +2,12 @@
 import pandas as pd
 import datetime as dt
 import copy
+from typing import Union, List, Dict
 
 from eaopack.portfolio import Portfolio
 from eaopack.optimization import Results, OptimProblem
 from eaopack import serialization 
+from eaopack.assets import Storage
 
 
 def extract_output(portf: Portfolio, op: OptimProblem, res:Results, prices: dict = None) -> dict:
@@ -58,6 +60,7 @@ def extract_output(portf: Portfolio, op: OptimProblem, res:Results, prices: dict
         # in case an asset links nodes, dispatch should be separate per node
         for a in portf.assets:
             for i_n, n in enumerate(a.nodes):
+                # sum up dispatch
                 I =   (op.mapping['asset'] == a.name) \
                     & (op.mapping['type'] == 'd')      \
                     & (op.mapping['node'] == n.name)
@@ -68,10 +71,10 @@ def extract_output(portf: Portfolio, op: OptimProblem, res:Results, prices: dict
                     myCol = (a.name +' ('+  n.name + ')')
                 disp[myCol] = 0.
                 for i,r in my_mapping.iterrows():
-                    disp.loc[times[r.time_step], myCol] += res.x[i]*r.disp_factor
+                    disp.loc[times[r.time_step], myCol] += res.x[i]*r.disp_factor                 
         # extract internal variables per asset
         for a in portf.assets:
-            variable_names = op.mapping[(op.mapping['asset'] == a.name)& (op.mapping['type'] == 'i')]['var_name'].unique()
+            variable_names = op.mapping[(op.mapping['asset'] == a.name) & (op.mapping['type'] == 'i')]['var_name'].unique()
             for v in variable_names:
                 I =   (op.mapping['asset'] == a.name) \
                     & (op.mapping['type'] == 'i')      \
@@ -80,8 +83,35 @@ def extract_output(portf: Portfolio, op: OptimProblem, res:Results, prices: dict
                 myCol = (a.name +' ('+  v + ')')
                 internal_variables[myCol] = None
                 for i,r in my_mapping.iterrows():
-                    pass
                     internal_variables.loc[times[r.time_step], myCol] = res.x[i]
+            # specific case: Storage; also extract disp_in,  disp_out and fill level separately
+            if isinstance(a, Storage): 
+                I =   (op.mapping['asset'] == a.name) \
+                    & (op.mapping['type'] == 'd')      \
+                    & (op.mapping['node'] == n.name)
+                my_mapping = op.mapping.loc[I,:]
+                ### extract ... disp in 
+                what = 'charge'
+                if len(portf.nodes)==1:
+                    myCol = a.name+'_'+what
+                else: # add node information
+                    myCol = (a.name +' ('+  n.name +'_'+ what + ')')
+                internal_variables[myCol] = 0.
+                for i,r in my_mapping.iterrows():
+                    internal_variables.loc[times[r.time_step], myCol] += max(0,-res.x[i])*r.disp_factor
+                ### extract ... disp out
+                what = 'discharge'
+                if len(portf.nodes)==1:
+                    myCol = a.name+'_'+what
+                else: # add node information
+                    myCol = (a.name +' ('+  n.name +'_'+ what + ')')
+                internal_variables[myCol] = 0.
+                for i,r in my_mapping.iterrows():
+                    internal_variables.loc[times[r.time_step], myCol] += min(0,-res.x[i])*r.disp_factor                                                                           
+                ### extract ... fill level 
+                myCol = a.name+'_fill_level'
+                internal_variables[myCol] = 0.
+                internal_variables.loc[:, myCol] = a.fill_level(op, res)
         # extract duals from nodal restrictions
         # looping through nodes and their recorded nodal restrictions and extract dual
         if not res.duals is None and not res.duals['N'] is None:
@@ -181,7 +211,7 @@ def output_to_file(output, file_name:str, format_output:str = 'xlsx',csv_ger:boo
 
 #### easy access to object parameters e.g. for assets & portfolio
 ## get tree, get parameter, set parameter
-def get_params_tree(obj) -> (list, dict):
+def get_params_tree(obj) -> Union[List, Dict]:
     """ get parameters of object - typically asset or portfolio
 
     Args:
@@ -191,7 +221,7 @@ def get_params_tree(obj) -> (list, dict):
         list of parameter names (nested)
     """
 
-    def make_dict(dd) -> (list, dict):
+    def make_dict(dd) -> Union[List, Dict]:
         if isinstance(dd, list): dd_keys = range(0,len(dd))
         elif isinstance(dd, dict): dd_keys = list(dd)
         else: 

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = eaopack
-version = 2.1.3
+version = 2.1.4
 author = The EAO development Team
 description = A Framework for Optimizing Decentralized Portfolios and Green Supply
 long_description = file: README.md