Code updated with possibility to run constrained optimization using e…

…xternal penalty functions

popt/cost_functions/ecalc_npv.py

@@ -9,7 +9,7 @@
 HERE = HERE.resolve()
 
 
-def ecalc_npv(pred_data, keys_opt, report):
+def ecalc_npv(pred_data, **kwargs):
     """
     Net present value cost function using eCalc to calculate emmisions
 
@@ -18,11 +18,14 @@ def ecalc_npv(pred_data, keys_opt, report):
     pred_data : array_like
         Ensemble of predicted data.
 
-    keys_opt : list
-        Keys with economic data.
+    **kwargs : dict
+        Other arguments sent to the npv function
 
-    report : list
-        Report dates.
+        keys_opt : list
+            Keys with economic data.
+
+        report : list
+            Report dates.
 
     Returns
     -------
@@ -34,6 +37,10 @@ def ecalc_npv(pred_data, keys_opt, report):
     from libecalc.common.time_utils import Frequency
     from libecalc.presentation.yaml.model import YamlModel
 
+    # Get the necessary input
+    keys_opt = kwargs.get('input_dict', {})
+    report = kwargs.get('true_order', [])
+
     # Economic values
     npv_const = {}
     for name, value in keys_opt['npv_const']:

popt/cost_functions/ecalc_pareto_npv.py

@@ -11,7 +11,7 @@
 HERE = HERE.resolve()
 
 
-def ecalc_pareto_npv(pred_data, keys_opt, report):
+def ecalc_pareto_npv(pred_data, kwargs):
     """
     Net present value cost function using eCalc to calculate emmisions
 
@@ -20,11 +20,14 @@ def ecalc_pareto_npv(pred_data, keys_opt, report):
     pred_data : array_like
         Ensemble of predicted data.
 
-    keys_opt : list
-        Keys with economic data.
+    **kwargs : dict
+        Other arguments sent to the npv function
 
-    report : list
-        Report dates.
+        keys_opt : list
+            Keys with economic data.
+
+        report : list
+            Report dates.
 
     Returns
     -------
@@ -36,6 +39,10 @@ def ecalc_pareto_npv(pred_data, keys_opt, report):
     from libecalc.common.time_utils import Frequency
     from libecalc.presentation.yaml.model import YamlModel
 
+    # Get the necessary input
+    keys_opt = kwargs.get('input_dict', {})
+    report = kwargs.get('true_order', [])
+
     # Economic values
     npv_const = {}
     for name, value in keys_opt['npv_const']:

popt/cost_functions/npv.py

@@ -2,7 +2,7 @@
 import numpy as np
 
 
-def npv(pred_data, keys_opt, report):
+def npv(pred_data, **kwargs):
     """
     Net present value cost function
 
@@ -11,25 +11,32 @@ def npv(pred_data, keys_opt, report):
     pred_data : array_like
         Ensemble of predicted data.
 
-    keys_opt : list
-        Keys with economic data.
+    **kwargs : dict
+        Other arguments sent to the npv function
 
-            - wop: oil price
-            - wgp: gas price
-            - wwp: water production cost
-            - wwi: water injection cost
-            - disc: discount factor
-            - obj_scaling: used to scale the objective function (negative since all methods are minimizers)
+        keys_opt : list
+            Keys with economic data.
 
-    report : list
-        Report dates.
+                - wop: oil price
+                - wgp: gas price
+                - wwp: water production cost
+                - wwi: water injection cost
+                - disc: discount factor
+                - obj_scaling: used to scale the objective function (negative since all methods are minimizers)
+
+        report : list
+            Report dates.
 
     Returns
     -------
     objective_values : numpy.ndarray
         Objective function values (NPV) for all ensemble members.
     """
 
+    # Get the necessary input
+    keys_opt = kwargs.get('input_dict',{})
+    report = kwargs.get('true_order', [])
+
     # Economic values
     npv_const = {}
     for name, value in keys_opt['npv_const']:

popt/cost_functions/quadratic.py

@@ -1,13 +1,15 @@
 import numpy as np
 
 
-def quadratic(state, *args):
+def quadratic(state, *args, **kwargs):
 	r"""
 	Quadratic objective function
 	.. math::
 		f(x) = ||x - b||^2_A
 	"""
 
+	r = kwargs.get('r', -1)
+
 	x = state[0]['vector']
 	dim, ne = x.shape
 	A = 0.5*np.diag(np.ones(dim))
@@ -17,4 +19,19 @@ def quadratic(state, *args):
 		u = x[:, i] - b
 		f[i] = u.T@A@u
 
+		# check for contraints
+		if r >= 0:
+			c_eq = g(x[:, i])
+			c_iq = h(x[:, i])
+			f[i] += r*0.5*( np.sum(c_eq**2) + np.sum(np.maximum(-c_iq,0)**2) )
+
 	return f
+
+# Equality constraint saying that sum of x should be equal to dimention + 1
+def g(x):
+	return sum(x) - (x.size + 1)
+
+# Inequality constrint saying that x_1 should be equal or less than 0
+def h(x):
+	return -x[0]
+

popt/cost_functions/ren_npv.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 
-def ren_npv(pred_data, keys_opt, report):
+def ren_npv(pred_data, kwargs):
     """
     Net present value cost function with injection from RENewable energy
 
@@ -12,18 +12,25 @@ def ren_npv(pred_data, keys_opt, report):
     pred_data_en : ndarray
         Ensemble of predicted data.
 
-    keys_opt : list
-        Keys with economic data.
+    **kwargs : dict
+        Other arguments sent to the npv function
 
-    report : list
-        Report dates.
+        keys_opt : list
+            Keys with economic data.
+
+        report : list
+            Report dates.
 
     Returns
     -------
     objective_values : ndarray
         Objective function values (NPV) for all ensemble members.
     """
 
+    # Get the necessary input
+    keys_opt = kwargs.get('input_dict', {})
+    report = kwargs.get('true_order', [])
+
     # Economic values
     npv_const = {}
     for name, value in keys_opt['npv_const']:

popt/cost_functions/rosenbrock.py

@@ -1,7 +1,7 @@
 """Rosenbrock objective function."""
 
 
-def rosenbrock(state, *args):
+def rosenbrock(state, *args, **kwargs):
     """
     Rosenbrock: http://en.wikipedia.org/wiki/Rosenbrock_function
     """

popt/loop/ensemble.py

@@ -103,9 +103,9 @@ def __set__variable(var_name=None, defalut=None):
                     self.bounds += num_state_var*[(0, 1)]
                 else:
                     self.bounds += num_state_var*[(lb, ub)]
-                self.cov = np.append(self.cov, value_cov)
             else:
                 self.bounds += num_state_var*[(None, None)]
+            self.cov = np.append(self.cov, value_cov)
 
         self._scale_state()
         self.cov = np.diag(self.cov)
@@ -189,7 +189,7 @@ def get_final_state(self, return_dict=False):
             x = self.get_state()
         return x
 
-    def function(self, x, *args):
+    def function(self, x, *args, **kwargs):
         """
         This is the main function called during optimization.
 
@@ -215,7 +215,8 @@ def function(self, x, *args):
         run_success = self.calc_prediction()  # calculate flow data
         self._scale_state()  # scale back to [0, 1]
         if run_success:
-            func_values = self.obj_func(self.pred_data, self.sim.input_dict, self.sim.true_order)
+            func_values = self.obj_func(self.pred_data, input_dict=self.sim.input_dict,
+                                        true_order=self.sim.true_order, **kwargs)
         else:
             func_values = np.inf  # the simulations have crashed
 
@@ -226,7 +227,7 @@ def function(self, x, *args):
 
         return func_values
 
-    def gradient(self, x, *args):
+    def gradient(self, x, *args, **kwargs):
         r"""
         Calculate the preconditioned gradient associated with ensemble, defined as:
 
@@ -277,7 +278,7 @@ def gradient(self, x, *args):
         self.state = self._gen_state_ensemble()
 
         state_ens = at.aug_state(self.state, list(self.state.keys()))
-        self.function(state_ens)
+        self.function(state_ens, **kwargs)
 
         # If bias correction is used we need to calculate the bias factors, J(u_j,m_j)/J(u_j,m)
         if self.bias_file is not None:  # use bias corrections
@@ -356,7 +357,7 @@ def genopt_mutation_gradient(self, x=None, *args, **kwargs):
         return self.genopt.ensemble_mutation_gradient(return_ensembles=kwargs['return_ensembles'])
     '''
 
-    def calc_ensemble_weights(self, x, *args):
+    def calc_ensemble_weights(self, x, *args, **kwargs):
         r"""
         Calculate weights used in sequential monte carlo optimization.
 
@@ -398,9 +399,12 @@ def calc_ensemble_weights(self, x, *args):
         self._invert_scale_state()  # ensure that state is in [lb,ub]
         self.calc_prediction()  # calculate flow data
         self._scale_state()  # scale back to [0, 1]
-        self.ens_func_values = self.obj_func(self.pred_data, self.sim.input_dict, self.sim.true_order)
-        self.ens_func_values = np.array(self.ens_func_values)
+
+        #self.ens_func_values = self.obj_func(self.pred_data, self.sim.input_dict, self.sim.true_order)
+        #self.ens_func_values = np.array(self.ens_func_values)
         state_ens = at.aug_state(self.state, list(self.state.keys()))
+        self.function(state_ens, **kwargs)
+
         self.particles = np.hstack((self.particles, state_ens))
         self.particle_values = np.hstack((self.particle_values,self.ens_func_values))
 
@@ -421,8 +425,8 @@ def calc_ensemble_weights(self, x, *args):
         index = np.argmin(self.particle_values)
         best_ens = self.particles[:, index]
         best_func = self.particle_values[index]
-        resample_index = np.random.choice(self.num_samples,int(np.round(self.num_samples-self.num_samples*self.survival_factor)),
-                                 replace=True,p=weights)
+        resample_index = np.random.choice(self.num_samples,int(np.round(self.num_samples-
+                                          self.num_samples*self.survival_factor)),replace=True,p=weights)
         self.particles = self.particles[:, resample_index]
         self.particle_values = self.particle_values[resample_index]
         return sens_matrix, best_ens, best_func