Skip to content

Commit 586b98e

Browse files
rolfjlMathias Methlie NilsenMathiasMNilsenRolf Johan Lorentzen
authored
Merge Mathias code (#93)
* corrected AdaMax optimizer (it was wrong!) * updated eCalc import * Added line-search method from scipy * fixed linesearch * made LineSearch very efficient! * added max step-size * changed co2 npv * included more in co2 calculation * corrected npv * fixed some bugs * some changes * small fix * made it much more efficient * re-implemented GenOpt * small change * added loop over Keras optimizers * added f0 and g0 * some modification * GenOpt hessian (not finished) * some update * added Newton method * added resampling and cleaned up some code * minor changes * Added resampling and some minor changes * updated LineSearch * Added BFGS and re-designed LineSearch * added interpolation of step-size * Made LineSearch a function (instead of class) * some additional changes * some additional changes * fixed NPV * minor changes * Cleaned up code * Some changes * big cleanup * fixed bug * Added Newton method and made some small changes to logic * cleaned up some code * Some changes * Added TrustRegion * re-designed TrustRegion * improved a transformation * fixed a bug * fixed a bug * minor changes * added natural gradient for BetaMC * added docstring to TrustRegion * Added mutation * Bug fix in localization input * Fix indent error --------- Co-authored-by: Mathias Methlie Nilsen <[email protected]> Co-authored-by: Mathias Methlie Nilsen <[email protected]> Co-authored-by: Rolf Johan Lorentzen <[email protected]>
1 parent 49bb2e7 commit 586b98e

File tree

10 files changed

+1971
-218
lines changed

10 files changed

+1971
-218
lines changed

‎popt/cost_functions/ren_npv_co2.py

+196
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,196 @@
1+
''' Net Present Value with Renewable Power and co2 emissions '''
2+
3+
import pandas as pd
4+
import numpy as np
5+
import os
6+
import yaml
7+
8+
from pathlib import Path
9+
from pqdm.processes import pqdm
10+
11+
__all__ = ['ren_npv_co2']
12+
13+
HERE = Path().cwd() # fallback for ipynb's
14+
HERE = HERE.resolve()
15+
16+
def ren_npv_co2(pred_data, keys_opt, report, save_emissions=False):
17+
'''
18+
Net Present Value with Renewable Power and co2 emissions (with eCalc)
19+
20+
Parameters
21+
----------
22+
pred_data : array_like
23+
Ensemble of predicted data.
24+
25+
keys_opt : list
26+
Keys with economic data.
27+
28+
report : list
29+
Report dates.
30+
31+
Returns
32+
-------
33+
objective_values : array_like
34+
Objective function values (NPV) for all ensemble members.
35+
'''
36+
37+
# some globals, for pqdm
38+
global const
39+
global kwargs
40+
global report_dates
41+
global sim_data
42+
43+
# define a data getter
44+
get_data = lambda i, key: pred_data[i+1][key].squeeze() - pred_data[i][key].squeeze()
45+
46+
# ensemble size (ne), number of report-dates (nt)
47+
nt = len(pred_data)
48+
try:
49+
ne = len(get_data(1,'fopt'))
50+
except:
51+
ne = 1
52+
53+
np.save('co2_emissions', np.zeros((ne, nt-1)))
54+
55+
# Economic and other constatns
56+
const = dict(keys_opt['npv_const'])
57+
kwargs = dict(keys_opt['npv_kwargs'])
58+
report_dates = report[1]
59+
60+
# Load energy arrays. These arrays contain the excess windpower used for gas compression,
61+
# and the energy from gas which is used in the water intection.
62+
power_arrays = np.load(kwargs['power']+'.npz')
63+
64+
sim_data = {'fopt': np.zeros((ne, nt-1)),
65+
'fgpt': np.zeros((ne, nt-1)),
66+
'fwpt': np.zeros((ne, nt-1)),
67+
'fwit': np.zeros((ne, nt-1)),
68+
'thp' : np.zeros((ne, nt-1)),
69+
'days': np.zeros(nt-1),
70+
'wind': power_arrays['wind'][:,:-1]}
71+
72+
# loop over pred_data
73+
for t in range(nt-1):
74+
75+
for datatype in ['fopt', 'fgpt', 'fwpt', 'fwit']:
76+
sim_data[datatype][:,t] = get_data(t, datatype)
77+
78+
# days in time-step
79+
sim_data['days'][t] = (report_dates[t+1] - report_dates[t]).days
80+
81+
# get maximum well head pressure (for each ensemble member)
82+
thp_keys = [k for k in keys_opt['datatype'] if 'wthp' in k] # assume only injection wells
83+
thp_vals = []
84+
for key in thp_keys:
85+
thp_vals.append(pred_data[t][key].squeeze())
86+
87+
sim_data['thp'][:,t] = np.max(np.array(thp_vals), axis=0)
88+
89+
# calculate NPV values
90+
npv_values = pqdm(array=range(ne), function=npv, n_jobs=keys_opt['parallel'], disable=True)
91+
92+
if not save_emissions:
93+
os.remove('co2_emissions.npy')
94+
95+
# clear energy arrays
96+
np.savez(kwargs['power']+'.npz', wind=np.zeros((ne, nt)), ren=np.zeros((ne,nt)), gas=np.zeros((ne,nt)))
97+
98+
scaling = 1.0
99+
if 'obj_scaling' in const:
100+
scaling = const['obj_scaling']
101+
102+
return np.asarray(npv_values)/scaling
103+
104+
105+
def emissions(yaml_file="ecalc_config.yaml"):
106+
107+
from libecalc.application.energy_calculator import EnergyCalculator
108+
from libecalc.common.time_utils import Frequency
109+
from libecalc.presentation.yaml.model import YamlModel
110+
111+
# Config
112+
model_path = HERE / yaml_file
113+
yaml_model = YamlModel(path=model_path, output_frequency=Frequency.NONE)
114+
115+
# Compute energy, emissions
116+
model = EnergyCalculator(graph=yaml_model.graph)
117+
consumer_results = model.evaluate_energy_usage(yaml_model.variables)
118+
emission_results = model.evaluate_emissions(yaml_model.variables, consumer_results)
119+
120+
# print power from pump
121+
co2 = []
122+
for identity, component in yaml_model.graph.nodes.items():
123+
if identity in emission_results:
124+
co2.append(emission_results[identity]['co2_fuel_gas'].rate.values)
125+
126+
co2 = np.sum(np.asarray(co2), axis=0)
127+
return co2
128+
129+
130+
def npv(n):
131+
132+
days = sim_data['days']
133+
134+
# config eCalc
135+
pd.DataFrame( {'dd-mm-yyyy' : report_dates[1:],
136+
'OIL_PROD' : sim_data['fopt'][n]/days,
137+
'GAS_PROD' : sim_data['fgpt'][n]/days,
138+
'WATER_INJ' : sim_data['fwit'][n]/days,
139+
'THP_MAX' : sim_data['thp'][n],
140+
'WIND_POWER' : sim_data['wind'][n]*(-1)
141+
} ).to_csv(f'ecalc_input_{n}.csv', index=False)
142+
143+
ecalc_yaml_file = kwargs['yamlfile']+'.yaml'
144+
new_yaml = duplicate_yaml_file(filename=ecalc_yaml_file, member=n)
145+
146+
#calc emissions
147+
co2 = emissions(new_yaml)*days
148+
149+
# save emissions
150+
try:
151+
en_co2 = np.load('co2_emissions.npy')
152+
en_co2[n] = co2
153+
np.save('co2_emissions', en_co2)
154+
except:
155+
import time
156+
time.sleep(1)
157+
158+
#calc npv
159+
gain = const['wop']*sim_data['fopt'][n] + const['wgp']*sim_data['fgpt'][n]
160+
loss = const['wwp']*sim_data['fwpt'][n] + const['wwi']*sim_data['fwit'][n] + const['wem']*co2
161+
disc = (1+const['disc'])**(days/365)
162+
163+
npv_value = np.sum( (gain-loss)/disc )
164+
165+
# delete dummy files
166+
os.remove(new_yaml)
167+
os.remove(f'ecalc_input_{n}.csv')
168+
169+
return npv_value
170+
171+
172+
def duplicate_yaml_file(filename, member):
173+
174+
try:
175+
# Load the YAML file
176+
with open(filename, 'r') as yaml_file:
177+
data = yaml.safe_load(yaml_file)
178+
179+
input_name = data['TIME_SERIES'][0]['FILE']
180+
data['TIME_SERIES'][0]['FILE'] = input_name.replace('.csv', f'_{member}.csv')
181+
182+
# Write the updated content to a new file
183+
new_filename = filename.replace(".yaml", f"_{member}.yaml")
184+
with open(new_filename, 'w') as new_yaml_file:
185+
yaml.dump(data, new_yaml_file, default_flow_style=False)
186+
187+
except FileNotFoundError:
188+
print(f"File '{filename}' not found.")
189+
190+
return new_filename
191+
192+
193+
194+
195+
196+

‎popt/loop/base.py

+163
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
1+
# External imports
2+
import numpy as np
3+
import sys
4+
import warnings
5+
6+
from copy import deepcopy
7+
8+
# Internal imports
9+
from popt.misc_tools import optim_tools as ot
10+
from pipt.misc_tools import analysis_tools as at
11+
from ensemble.ensemble import Ensemble as PETEnsemble
12+
13+
class EnsembleOptimizationBase(PETEnsemble):
14+
'''
15+
Base class for the popt ensemble
16+
'''
17+
def __init__(self, kwargs_ens, sim, obj_func):
18+
'''
19+
Parameters
20+
----------
21+
kwargs_ens : dict
22+
Options for the ensemble class
23+
24+
sim : callable
25+
The forward simulator (e.g. flow)
26+
27+
obj_func : callable
28+
The objective function (e.g. npv)
29+
'''
30+
31+
# Initialize PETEnsemble
32+
super().__init__(kwargs_ens, sim)
33+
34+
self.save_prediction = kwargs_ens.get('save_prediction', None)
35+
self.num_models = kwargs_ens.get('num_models', 1)
36+
self.transform = kwargs_ens.get('transform', False)
37+
self.num_samples = self.ne
38+
39+
# Get bounds and varaince
40+
self.upper_bound = []
41+
self.lower_bound = []
42+
self.bounds = []
43+
self.cov = np.array([])
44+
for name in self.prior_info.keys():
45+
self.state[name] = np.asarray(self.prior_info[name]['mean'])
46+
num_state_var = len(self.state[name])
47+
value_cov = self.prior_info[name]['variance'] * np.ones((num_state_var,))
48+
if 'limits' in self.prior_info[name].keys():
49+
lb = self.prior_info[name]['limits'][0]
50+
ub = self.prior_info[name]['limits'][1]
51+
self.lower_bound.append(lb)
52+
self.upper_bound.append(ub)
53+
if self.transform:
54+
value_cov = value_cov / (ub - lb)**2
55+
np.clip(value_cov, 0, 1, out=value_cov)
56+
self.bounds += num_state_var*[(0, 1)]
57+
else:
58+
self.bounds += num_state_var*[(lb, ub)]
59+
self.cov = np.append(self.cov, value_cov)
60+
else:
61+
self.bounds += num_state_var*[(None, None)]
62+
63+
64+
self._scale_state()
65+
self.cov = np.diag(self.cov)
66+
67+
# Set objective function (callable)
68+
self.obj_func = obj_func
69+
70+
# Objective function values
71+
self.state_func_values = None
72+
self.ens_func_values = None
73+
74+
def get_state(self):
75+
"""
76+
Returns
77+
-------
78+
x : numpy.ndarray
79+
Control vector as ndarray, shape (number of controls, number of perturbations)
80+
"""
81+
x = ot.aug_optim_state(self.state, list(self.state.keys()))
82+
return x
83+
84+
def get_bounds(self):
85+
"""
86+
Returns
87+
-------
88+
bounds : list
89+
(min, max) pairs for each element in x. None is used to specify no bound.
90+
"""
91+
92+
return self.bounds
93+
94+
def function(self, x, *args):
95+
"""
96+
This is the main function called during optimization.
97+
98+
Parameters
99+
----------
100+
x : ndarray
101+
Control vector, shape (number of controls, number of perturbations)
102+
103+
Returns
104+
-------
105+
obj_func_values : numpy.ndarray
106+
Objective function values, shape (number of perturbations, )
107+
"""
108+
self._aux_input()
109+
110+
if len(x.shape) == 1:
111+
self.ne = self.num_models
112+
else:
113+
self.ne = x.shape[1]
114+
115+
self.state = ot.update_optim_state(x, self.state, list(self.state.keys())) # go from nparray to dict
116+
self._invert_scale_state() # ensure that state is in [lb,ub]
117+
run_success = self.calc_prediction(save_prediction=self.save_prediction) # calculate flow data
118+
self._scale_state() # scale back to [0, 1]
119+
if run_success:
120+
func_values = self.obj_func(self.pred_data, self.sim.input_dict, self.sim.true_order)
121+
else:
122+
func_values = np.inf # the simulations have crashed
123+
124+
if len(x.shape) == 1:
125+
self.state_func_values = func_values
126+
else:
127+
self.ens_func_values = func_values
128+
129+
return func_values
130+
131+
def _aux_input(self):
132+
"""
133+
Set the auxiliary input used for multiple geological realizations
134+
"""
135+
136+
nr = 1 # nr is the ratio of samples over models
137+
if self.num_models > 1:
138+
if np.remainder(self.num_samples, self.num_models) == 0:
139+
nr = int(self.num_samples / self.num_models)
140+
self.aux_input = list(np.repeat(np.arange(self.num_models), nr))
141+
else:
142+
print('num_samples must be a multiplum of num_models!')
143+
sys.exit(0)
144+
return nr
145+
146+
def _scale_state(self):
147+
"""
148+
Transform the internal state from [lb, ub] to [0, 1]
149+
"""
150+
if self.transform and (self.upper_bound and self.lower_bound):
151+
for i, key in enumerate(self.state):
152+
self.state[key] = (self.state[key] - self.lower_bound[i])/(self.upper_bound[i] - self.lower_bound[i])
153+
np.clip(self.state[key], 0, 1, out=self.state[key])
154+
155+
def _invert_scale_state(self):
156+
"""
157+
Transform the internal state from [0, 1] to [lb, ub]
158+
"""
159+
if self.transform and (self.upper_bound and self.lower_bound):
160+
for i, key in enumerate(self.state):
161+
if self.transform:
162+
self.state[key] = self.lower_bound[i] + self.state[key]*(self.upper_bound[i] - self.lower_bound[i])
163+
np.clip(self.state[key], self.lower_bound[i], self.upper_bound[i], out=self.state[key])

0 commit comments

Comments
 (0)