Update plotting; remove unmet_demand; ad tests

* Unmet demand is now a decision variable, only triggered when
you 'ensure_feasibility' in the model-wide settings

calliope/analysis/util.py

@@ -20,9 +20,9 @@ def get_zoom(coordinate_array, width):
 
     # Keys are zoom levels, values are m/pixel at that zoom level
     zoom_dict = {0: 156412, 1: 78206, 2: 39103, 3: 19551, 4: 9776, 5: 4888,
-                6: 2444, 7: 1222, 8: 610.984, 9: 305.492, 10: 152.746,
-                11: 76.373, 12: 38.187, 13: 19.093, 14: 9.547, 15: 4.773,
-                16: 2.387, 17: 1.193, 18: 0.596, 19: 0.298}
+                 6: 2444, 7: 1222, 8: 610.984, 9: 305.492, 10: 152.746,
+                 11: 76.373, 12: 38.187, 13: 19.093, 14: 9.547, 15: 4.773,
+                 16: 2.387, 17: 1.193, 18: 0.596, 19: 0.298}
 
     bounds = [coordinate_array.max(dim='locs').values,
               coordinate_array.min(dim='locs').values]
@@ -35,7 +35,29 @@ def get_zoom(coordinate_array, width):
         if v > metres_per_pixel:
             continue
         else:
-            zoom = k-3
+            zoom = k - 3
             break
 
     return zoom
+
+
+def subset_sum_squeeze(data, subset={}, sum_dims=None, squeeze=True):
+    if subset:  # first, subset the data
+        allowed_subsets = {k: v for k, v in subset.items() if k in data.dims}
+        data = data.loc[allowed_subsets]
+
+    if sum_dims:  # second, sum along all necessary dimensions
+        data = data.sum(sum_dims)
+
+    if squeeze and len(data.techs) > 1:  # finally, squeeze out single length dimensions
+        data = data.squeeze()
+
+    return data
+
+
+def hex_to_rgba(hex_color, opacity):
+    _NUMERALS = '0123456789abcdefABCDEF'
+    _HEXDEC = {v: int(v, 16) for v in (x + y for x in _NUMERALS for y in _NUMERALS)}
+    hex_color = hex_color.lstrip('#')
+    rgb = [_HEXDEC[hex_color[0:2]], _HEXDEC[hex_color[2:4]], _HEXDEC[hex_color[4:6]]]
+    return 'rgba({1}, {2}, {3}, {0})'.format(opacity, *rgb)

calliope/backend/pyomo/constraints/energy_balance.py

@@ -91,10 +91,15 @@ def system_balance_constraint_rule(backend_model, loc_carrier, timestep):
 
     """
     prod, con, export = get_loc_tech_carriers(backend_model, loc_carrier)
+    if hasattr(backend_model, 'unmet_demand'):
+        unmet_demand = backend_model.unmet_demand[loc_carrier, timestep]
+    else:
+        unmet_demand = 0
 
     backend_model.system_balance[loc_carrier, timestep].expr = (
         sum(backend_model.carrier_prod[loc_tech_carrier, timestep] for loc_tech_carrier in prod) +
-        sum(backend_model.carrier_con[loc_tech_carrier, timestep] for loc_tech_carrier in con)
+        sum(backend_model.carrier_con[loc_tech_carrier, timestep] for loc_tech_carrier in con) +
+        unmet_demand
     )
 
     return backend_model.system_balance[loc_carrier, timestep] == 0

calliope/backend/pyomo/initialize.py

@@ -172,7 +172,7 @@ def run_iterative(model_data, timings):
         # We always save the window data. Until the last window(s) this will crop
         # the window_to_horizon timesteps. In the last window(s), optimistion will
         # only be occurring over a window length anyway
-        results = results.loc[dict(timesteps=slice(None,window_ends.index[i]))]
+        results = results.loc[dict(timesteps=slice(None, window_ends.index[i]))]
         result_array.append(results)
 
         # Set up initial storage for the next iteration

calliope/backend/pyomo/model.py

@@ -55,7 +55,8 @@ def generate_model(model_data):
             model_data[k].to_series().dropna().replace('inf', np.inf).to_dict()
             for k in model_data.data_vars},
         'dims': {k: model_data[k].dims for k in model_data.data_vars},
-        'sets': list(model_data.coords)
+        'sets': list(model_data.coords),
+        'attrs': {k: v for k, v in model_data.attrs.items() if k is not 'defaults'}
     }
     # Dims in the dict's keys are ordered as in model_data, which is enforced
     # in model_data generation such that timesteps are always last and the
@@ -191,4 +192,7 @@ def get_result_array(backend_model):
         i.name: get_var(backend_model, i.name) for i in backend_model.component_objects()
         if isinstance(i, po.base.var.IndexedVar)
     }
+    # if unmet_demand was unused, delete it before it reaches the user
+    if 'unmet_demand' in all_variables.keys() and not sum(backend_model.unmet_demand.get_values().values()):
+        del all_variables['unmet_demand']
     return reorganise_dataset_dimensions(xr.Dataset(all_variables))

calliope/backend/pyomo/objective.py

@@ -21,10 +21,22 @@ def cost_minimization(backend_model):
         min: z = \sum_{loc::tech_{cost}} cost(loc::tech, cost=cost_{monetary}))
 
     """
+
     def obj_rule(backend_model):
-        return sum(
-            backend_model.cost['monetary', loc_tech]
-            for loc_tech in backend_model.loc_techs_cost
+        if hasattr(backend_model, 'unmet_demand'):
+            unmet_demand = sum(
+                backend_model.unmet_demand[loc_carrier, timestep]
+                for loc_carrier in backend_model.loc_carriers
+                for timestep in backend_model.timesteps
+            ) * backend_model.bigM
+        else:
+            unmet_demand = 0
+
+        return (
+            sum(
+                backend_model.cost['monetary', loc_tech]
+                for loc_tech in backend_model.loc_techs_cost
+            ) + unmet_demand
         )
 
     backend_model.obj = po.Objective(sense=po.minimize, rule=obj_rule)

calliope/backend/pyomo/variables.py

@@ -7,6 +7,7 @@
 import pyomo.core as po  # pylint: disable=import-error
 import numpy as np
 
+
 def initialize_decision_variables(backend_model):
     """
     Defines variables
@@ -88,3 +89,7 @@ def cap_initializer(cap_var_name, within):
         if backend_model.mode == 'operate':
             for k, v in backend_model.units.items():
                 backend_model.energy_cap[k] = v * backend_model.energy_cap_per_unit[k]
+
+    if model_data_dict['attrs'].get('model.ensure_feasibility', False):
+        backend_model.unmet_demand = po.Var(backend_model.loc_carriers, backend_model.timesteps, within=po.NonNegativeReals)
+        backend_model.bigM = model_data_dict['attrs'].get('model.bigM')

calliope/config/model.yaml

@@ -28,6 +28,8 @@ model:
     operation:
         horizon: null
         window: null
+    bigM: 1e9 # Used for unmet demand, but should be of a similar order of magnitude as the largest cost that the model could achieve. Too high and the model will not converge
+    ensure_feasibility: false # If true, unmet_demand will be a decision variable, to account for an ability to meet demand with the available supply. If False and a mismatch occurs, the optimisation will fail due to infeasibility
 
 ##
 # Base technology groups
@@ -65,21 +67,6 @@ tech_groups:
             resource_unit: power
             force_resource: true
             energy_con: true
-    unmet_demand:
-        required_constraints: []
-        allowed_constraints: ['energy_prod', 'energy_cap_max', 'energy_prod', 'resource', 'resource_unit']
-        allowed_costs: ['om_prod']
-        essentials:
-            parent: null
-            color: '#666666'
-        constraints:
-            resource: inf
-            resource_unit: power
-            energy_cap_max: inf
-            energy_prod: true
-        costs:
-            monetary:
-                om_prod: 1.0e+9
     storage:
         required_constraints: [['energy_cap_max', 'energy_cap_equals', 'storage_cap_max', 'storage_cap_equals', 'energy_cap_per_unit'], 'charge_rate']
         allowed_constraints: ['energy_prod', 'energy_con', 'lifetime', 'resource_cap_max', 'resource_cap_equals', 'resource_cap_equals_energy_cap', 'resource_scale', 'resource_scale_to_peak', 'energy_eff', 'energy_cap_min', 'energy_cap_max', 'energy_cap_equals', 'energy_cap_max_systemwide', 'energy_cap_equals_systemwide', 'energy_cap_scale', 'energy_cap_min_use', 'energy_cap_per_unit', 'energy_ramping', 'storage_initial', 'storage_cap_min', 'storage_cap_max', 'storage_cap_equals', 'storage_cap_per_unit', 'charge_rate', 'storage_time_max', 'storage_loss', 'units_min', 'units_max', 'units_equals']

calliope/core/preprocess/model_run.py

@@ -279,7 +279,7 @@ def process_techs(config_model):
         # CHECK: If necessary, populate carrier_in and carrier_out in essentials, but
         # also break on missing carrier data
         if 'carrier_in' not in tech_result.essentials:
-            if tech_result.inheritance[-1] in ['supply', 'supply_plus', 'unmet_demand']:
+            if tech_result.inheritance[-1] in ['supply', 'supply_plus']:
                 tech_result.essentials.carrier_in = 'resource'
             elif tech_result.inheritance[-1] in ['demand', 'transmission',
                                                  'storage']:
@@ -300,7 +300,7 @@ def process_techs(config_model):
         if 'carrier_out' not in tech_result.essentials:
             if tech_result.inheritance[-1] == 'demand':
                 tech_result.essentials.carrier_out = 'resource'
-            elif tech_result.inheritance[-1] in ['supply', 'supply_plus', 'unmet_demand',
+            elif tech_result.inheritance[-1] in ['supply', 'supply_plus',
                                                  'transmission', 'storage']:
                 try:
                     tech_result.essentials.carrier_out = \

calliope/core/preprocess/sets.py

@@ -40,7 +40,6 @@
 * loc_techs_supply_plus
 * loc_techs_conversion
 * loc_techs_conversion_plus
-* loc_techs_unmet_demand
 
 Subsets based on active constraints
 
@@ -229,7 +228,7 @@ def generate_loc_tech_sets(model_run, simple_sets):
 
     for group in [
             'storage', 'demand', 'supply', 'supply_plus',
-            'unmet_demand', 'conversion', 'conversion_plus']:
+            'conversion', 'conversion_plus']:
         tech_set = set(
             k for k in sets.loc_techs_non_transmission
             if model_run.techs[k.split('::')[1]].inheritance[-1] == group
@@ -245,8 +244,7 @@ def generate_loc_tech_sets(model_run, simple_sets):
     # Techs that introduce energy into the system
     sets.loc_techs_supply_all = (
         sets.loc_techs_supply |
-        sets.loc_techs_supply_plus |
-        sets.loc_techs_unmet_demand
+        sets.loc_techs_supply_plus
     )
 
     ##

calliope/core/util/dataset.py

@@ -74,6 +74,8 @@ def split_loc_techs(data_var, as_='DataArray'):
 
     # Separately find the loc_techs(_carriers) dimension and all other dimensions
     loc_tech_dim = [i for i in data_var.dims if 'loc_tech' in i]
+    if not loc_tech_dim:
+        loc_tech_dim = [i for i in data_var.dims if 'loc_carrier' in i]
     non_loc_tech_dims = list(set(data_var.dims).difference(loc_tech_dim))
 
     if not loc_tech_dim:
@@ -97,8 +99,8 @@ def split_loc_techs(data_var, as_='DataArray'):
 
     # carrier_prod, carrier_con, and carrier_export will return an index_list
     # of size 3, all others will be an index list of size 2
-    possible_names = ['locs', 'techs', 'carriers']
-    names = [possible_names[i] for i in range(len(index_list[0]))]
+    possible_names = ['loc', 'tech', 'carrier']
+    names = [i + 's' for i in possible_names if i in loc_tech_dim]
 
     data_var_df.index = pd.MultiIndex.from_tuples(index_list, names=names)
 

calliope/example_models/national_scale/model.yaml

@@ -19,5 +19,9 @@ model:
 
     subset_time: ['2005-01-01', '2005-01-05']  # Subset of timesteps
 
+    ensure_feasibility: true # Switching on unmet demand
+
+    bigM: 1e6 # setting the scale of unmet demand, which cannot be too high, otherwise the optimisation will not converge
+
 run:
     solver: glpk

calliope/example_models/national_scale/model_config/locations.yaml

@@ -6,7 +6,6 @@ locations:
     region1:
         coordinates: {lat: 40, lon: -2}
         techs:
-            unmet_demand_power:
             demand_power:
                 constraints:
                     resource: file=demand-1.csv:demand
@@ -17,7 +16,6 @@ locations:
     region2:
         coordinates: {lat: 40, lon: -8}
         techs:
-            unmet_demand_power:
             demand_power:
                 constraints:
                     resource: file=demand-2.csv:demand

calliope/example_models/national_scale/model_config/techs.yaml

@@ -78,11 +78,6 @@ techs:
             name: 'Power demand'
             parent: demand
             carrier: power
-    unmet_demand_power:
-        essentials:
-            name: 'Unmet power demand'
-            parent: unmet_demand
-            carrier: power
 
     ##
     # Transmission

calliope/example_models/urban_scale/model.yaml

@@ -15,5 +15,9 @@ model:
 
     subset_time: ['2005-07-01', '2005-07-02']  # Subset of timesteps
 
+    ensure_feasibility: true # Switching on unmet demand
+
+    bigM: 1e6 # setting the scale of unmet demand, which cannot be too high, otherwise the optimisation will not converge
+
 run:
     solver: glpk

calliope/example_models/urban_scale/model_config/locations.yaml

@@ -10,8 +10,6 @@ locations:
                 constraints.resource: file=demand_power.csv
             demand_heat:
                 constraints.resource: file=demand_heat.csv
-            unmet_demand_electricity:
-            unmet_demand_heat:
         available_area: 500
         coordinates: {x: 2, y: 7}
 
@@ -28,8 +26,6 @@ locations:
                 constraints.resource: file=demand_power.csv
             demand_heat:
                 constraints.resource: file=demand_heat.csv
-            unmet_demand_electricity:
-            unmet_demand_heat:
         available_area: 1300
         coordinates: {x: 8, y: 7}
 
@@ -47,8 +43,6 @@ locations:
                 constraints.resource: file=demand_power.csv
             demand_heat:
                 constraints.resource: file=demand_heat.csv
-            unmet_demand_electricity:
-            unmet_demand_heat:
         available_area: 900
         coordinates: {x: 5, y: 3}
 

calliope/example_models/urban_scale/model_config/techs.yaml

@@ -110,24 +110,12 @@ techs:
             parent: demand
             carrier: electricity
 
-    unmet_demand_electricity:
-        essentials:
-            name: 'Unmet electrical demand'
-            parent: unmet_demand
-            carrier: electricity
-
     demand_heat:
         essentials:
             name: 'Heat demand'
             parent: demand
             carrier: heat
 
-    unmet_demand_heat:
-        essentials:
-            name: 'Unmet heat demand'
-            parent: unmet_demand
-            carrier: heat
-
 ##-DISTRIBUTION-##
 
     power_lines: