popt.cost_functions.ecalc_npv

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