popt.cost_functions.ecalc_pareto_npv

  1"""Net present value."""
  2import numpy
  3import numpy as np
  4import csv
  5from pathlib import Path
  6import pandas as pd
  7import sys
  8import os
  9
 10HERE = Path().cwd()  # fallback for ipynb's
 11HERE = HERE.resolve()
 12
 13
 14def ecalc_pareto_npv(pred_data, kwargs):
 15    """
 16    Net present value cost function using eCalc to calculate emmisions
 17
 18    Parameters
 19    ----------
 20    pred_data : array_like
 21        Ensemble of predicted data.
 22
 23    **kwargs : dict
 24        Other arguments sent to the npv function
 25
 26        keys_opt : list
 27            Keys with economic data.
 28
 29        report : list
 30            Report dates.
 31
 32    Returns
 33    -------
 34    objective_values : array_like
 35        Objective function values (NPV) for all ensemble members.
 36    """
 37
 38    from libecalc.application.energy_calculator import EnergyCalculator
 39    from libecalc.common.time_utils import Frequency
 40    from libecalc.presentation.yaml.model import YamlModel
 41
 42    # Get the necessary input
 43    keys_opt = kwargs.get('input_dict', {})
 44    report = kwargs.get('true_order', [])
 45
 46    # Economic values
 47    npv_const = {}
 48    for name, value in keys_opt['npv_const']:
 49        npv_const[name] = value
 50
 51    # Collect production data
 52    Qop = []
 53    Qgp = []
 54    Qwp = []
 55    Qwi = []
 56    Dd = []
 57    for i in np.arange(1, len(pred_data)):
 58
 59        Qop.append(np.squeeze(pred_data[i]['fopt']) - np.squeeze(pred_data[i - 1]['fopt']))
 60        Qgp.append(np.squeeze(pred_data[i]['fgpt']) - np.squeeze(pred_data[i - 1]['fgpt']))
 61        Qwp.append(np.squeeze(pred_data[i]['fwpt']) - np.squeeze(pred_data[i - 1]['fwpt']))
 62        Qwi.append(np.squeeze(pred_data[i]['fwit']) - np.squeeze(pred_data[i - 1]['fwit']))
 63        Dd.append((report[1][i] - report[1][i - 1]).days)
 64
 65    # Write production data to .csv file for eCalc input, for each ensemble member
 66    Qop = np.array(Qop).T
 67    Qwp = np.array(Qwp).T
 68    Qgp = np.array(Qgp).T
 69    Qwi = np.array(Qwi).T
 70    Dd = np.array(Dd)
 71    if len(Qop.shape) == 1:
 72        Qop = np.expand_dims(Qop,0)
 73        Qwp = np.expand_dims(Qwp, 0)
 74        Qgp = np.expand_dims(Qgp, 0)
 75        Qwi = np.expand_dims(Qwi, 0)
 76
 77    N = Qop.shape[0]
 78    T = Qop.shape[1]
 79    values = []
 80    pareto_values = []
 81    for n in range(N):
 82        with open('ecalc_input.csv', 'w') as csvfile:
 83            writer = csv.writer(csvfile, delimiter=',')
 84            writer.writerow(['dd/mm/yyyy', 'GAS_PROD', 'OIL_PROD', 'WATER_INJ'])
 85            for t in range(T):
 86                D = report[1][t]
 87                writer.writerow([D.strftime("%d/%m/%Y"), Qgp[n, t]/Dd[t], Qop[n, t]/Dd[t], Qwi[n, t]/Dd[t]])
 88
 89        # Config
 90        model_path = HERE / "ecalc_config.yaml"  # "drogn.yaml"
 91        yaml_model = YamlModel(path=model_path, output_frequency=Frequency.NONE)
 92        # comps = {c.name: id_hash for (id_hash, c) in yaml_model.graph.components.items()}
 93
 94        # Compute energy, emissions
 95        model = EnergyCalculator(graph=yaml_model.graph)
 96        consumer_results = model.evaluate_energy_usage(yaml_model.variables)
 97        emission_results = model.evaluate_emissions(yaml_model.variables, consumer_results)
 98
 99        # Extract
100        # energy = results_as_df(yaml_model, consumer_results, lambda r: r.component_result.energy_usage)
101        emissions = results_as_df(yaml_model, emission_results, lambda r: r['co2_fuel_gas'].rate)
102        emissions_total = emissions.sum(1).rename("emissions_total")
103        emissions_total.to_csv(HERE / "emissions.csv")
104        Qem = emissions_total.values * Dd  # total number of tons
105
106        value1 = (Qop[n, :] * npv_const['wop'] + Qgp[n, :] * npv_const['wgp'] - Qwp[n, :] * npv_const['wwp'] -
107                  Qwi[n, :] * npv_const['wwi'] - Qem * npv_const['wem']) / (
108                         (1 + npv_const['disc']) ** (Dd / 365))
109        value1 = np.sum(value1)
110        if 'obj_scaling' in npv_const:
111            value1 /= npv_const['obj_scaling']
112
113        value2 = np.array([])
114        if 'wemc' in npv_const:  # multi-opjective with co2 cost correction
115            value2 = - Qem * npv_const['wemc'] / ((1 + npv_const['disc']) ** (Dd / 365))
116            value2 = np.sum(value2)
117            if 'obj_scaling' in npv_const:
118                value2 /= npv_const['obj_scaling']
119        elif 'w' in npv_const:  # multi-objective with emission intensity
120            rho_o = 840.0  # oil density
121            rho_g = 1  # gas density
122            conv = 1000  # convert from kilo to tonne
123            value2 = np.sum(Qem*conv) / (np.sum(Qop[n, :]*rho_o + Qgp[n, :]*rho_g)/conv) # kg/toe
124
125        pareto_values.append([np.sum(Qem), value1, value2])
126        if 'w' in npv_const:
127            values.append((1-npv_const['w'])*value1 + npv_const['w']*value2)
128        else:
129            values.append(value1 + value2)  # total objective function
130
131    # Save emissions and both objective functions for later analysis
132    pareto_file = 'pareto_values.npz'
133    if os.path.exists(pareto_file):
134        pareto_iterations = np.load(pareto_file,allow_pickle=True)['pareto_iterations'][()]
135    else:
136        pareto_iterations = {}
137    num_eval = len(pareto_iterations)
138    pareto_iterations[str(num_eval)] = pareto_values
139    np.savez('pareto_values.npz', pareto_iterations=pareto_iterations)
140
141    return np.array(values)
142
143def results_as_df(yaml_model, results, getter) -> pd.DataFrame:
144    """Extract relevant values, as well as some meta (`attrs`)."""
145    df = {}
146    attrs = {}
147    res = None
148    for id_hash in results:
149        res = results[id_hash]
150        res = getter(res)
151        component = yaml_model.graph.get_node(id_hash)
152        df[component.name] = res.values
153        attrs[component.name] = {'id_hash': id_hash,
154                                 'kind': type(component).__name__,
155                                 'unit': res.unit}
156    if res is None:
157        sys.exit('No emission results from eCalc!')
158    df = pd.DataFrame(df, index=res.timesteps)
159    df.index.name = "dates"
160    df.attrs = attrs
161    return df