Source code for SALib.test_functions.lake_problem

from typing import Union
import numpy as np
from numpy import log, sqrt
from scipy.optimize import brentq


FLOAT_OR_ARRAY = Union[float, np.array]


[docs] def lake_problem( X: FLOAT_OR_ARRAY, a: FLOAT_OR_ARRAY = 0.1, q: FLOAT_OR_ARRAY = 2.0, b: FLOAT_OR_ARRAY = 0.42, eps: FLOAT_OR_ARRAY = 0.02, ) -> float: """Lake Problem as given in Hadka et al., (2015) and Kwakkel (2017) modified for use as a test function. The `mean` and `stdev` parameters control the log normal distribution of natural inflows (`epsilon` in [1] and [2]). References ---------- .. [1] Hadka, D., Herman, J., Reed, P., Keller, K., (2015). "An open source framework for many-objective robust decision making." Environmental Modelling & Software 74, 114-129. doi:10.1016/j.envsoft.2015.07.014 .. [2] Kwakkel, J.H, (2017). "The Exploratory Modeling Workbench: An open source toolkit for exploratory modeling, scenario discovery, and (multi-objective) robust decision making." Environmental Modelling & Software 96, 239-250. doi:10.1016/j.envsoft.2017.06.054 .. [3] Singh, R., Reed, P., Keller, K., (2015). "Many-objective robust decision making for managing an ecosystem with a deeply uncertain threshold response." Ecology and Society 20. doi:10.5751/ES-07687-200312 Parameters ---------- X : float or np.ndarray normalized concentration of Phosphorus at point in time a : float or np.ndarray rate of anthropogenic pollution (0.0 to 0.1) q : float or np.ndarray exponent controlling recycling rate (2.0 to 4.5). b : float or np.ndarray decay rate for phosphorus (0.1 to 0.45, where default 0.42 is irreversible, as described in [1]) eps : float or np.ndarray natural inflows of phosphorus (pollution), see [3] Returns ------- float, phosphorus pollution for a point in time """ Xq = X**q X_t1 = X + a + (Xq / (1.0 + Xq)) - (b * X) + eps return X_t1
[docs] def evaluate_lake(values: np.ndarray, seed=101) -> np.ndarray: """Evaluate the Lake Problem with an array of parameter values. References ---------- .. [1] Hadka, D., Herman, J., Reed, P., Keller, K., (2015). "An open source framework for many-objective robust decision making." Environmental Modelling & Software 74, 114–129. doi:10.1016/j.envsoft.2015.07.014 .. [2] Singh, R., Reed, P., Keller, K., (2015). "Many-objective robust decision making for managing an ecosystem with a deeply uncertain threshold response." Ecology and Society 20. doi:10.5751/ES-07687-200312 Parameters ---------- values : np.ndarray, model inputs in the (column) order of a, q, b, mean, stdev where * `a` is rate of anthropogenic pollution * `q` is exponent controlling recycling rate * `b` is decay rate for phosphorus * `mean` and * `stdev` set the log normal distribution of `eps`, see [2] Returns ------- np.ndarray, of Phosphorus pollution over time `t` """ rng = np.random.default_rng(seed) nvars = values.shape[0] a, q, b, mean, stdev = values.T sq_mean = mean**2 sq_std = stdev**2 eps = rng.lognormal( log(sq_mean / sqrt(sq_std + sq_mean)), sqrt(log(1.0 + sq_std / sq_mean)), size=nvars, ) Y = np.zeros((nvars, nvars)) for t in range(nvars): # First X value will be last Y value (should be 0 as we are filling # an array of zeros) Y[t] = lake_problem(Y[t - 1], a, q, b, eps) return Y
[docs] def evaluate(values: np.ndarray, nvars: int = 100, seed=101): """Evaluate the Lake Problem with an array of parameter values. Parameters ---------- values : np.ndarray, model inputs in the (column) order of a, q, b, mean, stdev, delta, alpha nvars : int, number of decision variables to simulate (default: 100) Returns ------- np.ndarray : max_P, utility, inertia, reliability """ a, q, b, _, __, delta, alpha = values.T nsamples = len(a) Y = np.empty((nsamples, 4)) for i in range(nsamples): tmp = evaluate_lake(values[i, :5], seed=seed) a_i = a[i] q_i = q[i] Pcrit = brentq(lambda x: x**q_i / (1.0 + x**q_i) - b[i] * x, 0.01, 1.5) reliability = len(tmp[tmp < Pcrit]) / nvars max_P = np.max(tmp) utility = np.sum(alpha[i] * a_i * np.power(delta[i], np.arange(nvars))) # In practice, `a` will be set by a separate decision model # See [2] in docstring for `lake_problem()` # Here, it is a constant for a given scenario. # The value for `tau` (0.02) is taken from [2]. inertia = len(a_i[a_i < 0.02]) / nvars Y[i, :] = max_P, utility, inertia, reliability return Y
if __name__ == "__main__": from SALib.sample import latin from SALib.analyze import delta SEED_VAL = 101 LAKE_SPEC = { "num_vars": 7, "names": ["a", "q", "b", "mean", "stdev", "delta", "alpha"], "bounds": [ [0.0, 0.1], [2.0, 4.5], [0.1, 0.45], [0.01, 0.05], [0.001, 0.005], [0.93, 0.99], [0.2, 0.5], ], "outputs": ["max_P", "Utility", "Inertia", "Reliability"], } latin_samples = latin.sample(LAKE_SPEC, 1000, seed=SEED_VAL) Y = evaluate(latin_samples) for i, name in enumerate(LAKE_SPEC["outputs"]): Si = delta.analyze(LAKE_SPEC, latin_samples, Y[:, i]) print(name) print(Si.to_df())