# Lake model

see also [this general introduction to the workbench](https://waterprogramming.wordpress.com/2017/11/01/using-the-exploratory-modelling-workbench/) as a source of inspiration for completing the assignment below

## the lake model
The exploratory modeling workbench includes an example folder. This folder contains a variety of examples that demonstrate the functionality of the workbench. Many of these examples have been drawn from published cases. Here, we use the Lake Problem as an example for demonstrating some of the key functionality of the workbench. 

We demonstrate some of the key capabilities of the exploratory modeling workbench using the Lake problem. The lake problem is a stylized and hypothetical decision problem where the population of a city has to decide on the amount of annual pollution it will put into a lake. It the pollution in the lake passes a threshold, it will suffer irreversible eutrophication. 

\begin{equation}
    X_{(t+1)}=X_t+a_t+\frac{(X_t^q)}{(1+X_t^q )}- bX_t+\epsilon_t
\end{equation}

where $X_t$ is the pollution at time $t$, $a_t$ is the rate of anthropogenic pollution at time $t$, $b$ is the lake’s natural removal rate, $q$ is the lake's natural recycling rate, $\epsilon_t$ is the rate of natural pollution at time $t$. The rate of anthropogenic pollution $a_t$ is the decision variable and is somewhere between 0, and 0.1. So $a_t \in [0,0.1]$. The natural pollution $\epsilon_t$ is modeled, following Singh et al. (2015), as a log normal distribution with mean $\mu$ and standard deviation $\sigma$. 
There are four outcomes of interest. The first is the average concentration of phosphor in the lake. 

\begin{equation}
    f_{phosphorus}=  \frac{1}{\left\vert{T}\right\vert} \sum\limits_{t\in{T}} X_t 
\end{equation}

where $\left\vert{T}\right\vert$ is the cardinality of the set of points in time. 
The second objective is the economic benefit derived from polluting the lake. Following Singh et al. (2015), this is defined as the discounted benefit of pollution mines the costs of having a polluted lake

\begin{equation}
    f_{economic} = \sum\limits_{t \in {T}}\alpha a_t \delta^t 
\end{equation}

where $\alpha$ is the utility derived from polluting and $\delta$ is the discount rate. By default, $\alpha$ is 0.04.
The third objective is related to the year over year change in the anthropogenic pollution rate. 

\begin{equation}
    f_{inertia} =\frac{1}{\left\vert{T}\right\vert-1} \sum\limits_{t=1}^{\left\vert{T}\right\vert} I(|a_{t}-a_{t-1} |>\tau)   
\end{equation}

where $I$ is an indicator function that is 0 if the statement is false, and 1 if the statement is true, $\tau$ is the threshold that is deemed undesirable, and is for illustrative purposes et to 0.2. Effectively, f_{inertia} is the fraction of years where the absolute value of the change in anthropogenic pollution is larger then $\tau$.
The fourth objective is the fraction of years where the pollution in the lake is below the critical threshold.

\begin{equation}
    f_{reliability} =  \frac{1}{\left\vert{T}\right\vert} \sum\limits_{t \in T}I(X_{t}<X_{crit} ) 
\end{equation}

where $I$ is an indicator function that is 0 if the statement is false, and 1 if the statement is true, $X_{crit}$ is the critical threshold of pollution and is a function of both $b$ and $q$.

The lake problem is characterized by both stochastic uncertainty and deep uncertainty. The stochastic uncertainty arises from the natural inflow. To reduce this stochastic uncertainty, multiple replications are performed and the average over the replication is taken. Deep uncertainty is presented by uncertainty about the mean $\mu$ and standard deviation $sigma$ of the lognormal distribution characterizing the natural inflow, the natural removal rate of the lake $\beta$, the natural recycling rate of the lake $q$, and the discount rate $\delta$. The table below specifies the ranges for the deeply uncertain factors, as well as their best estimate or default values. 


## Assignment
1. Given the Python implementation of the lake problem in lakemodel_function.py, adapt this code and connect it to the workbench

for the uncertainties, use the following table

|Parameter	|Range	        |Default value|
|-----------|--------------:|------------:|
|$\mu$    	|0.01 – 0.05	|0.02         |
|$\sigma$	|0.001 – 0.005 	|0.0017       |
|$b$      	|0.1 – 0.45	    |0.42         |
|$q$	    |2 – 4.5	    |2            |
|$\delta$	|0.93 – 0.99	|0.98         |

For now, assume that for each year a release decision is made. The release is between 0 and 0.1. Carefully look at line 23 in lake_model.py to identify the name to use for each lever.




2. Explore the behavior of the system in the absence of any release using 1000 scenarios, and the default sampling approach.
    * visualize the outcomes of interest, are there any apparent trade-offs?
    * can you visually identify the uncertainties that drive system behavior?



```python
import numpy as np
from lakemodel_function import lake_problem
from ema_workbench import (RealParameter, ScalarOutcome, Constant, Model)
import pandas as pd
import plotly_express as px
import seaborn as sns
from ema_workbench.analysis import pairs_plotting
import matplotlib.pyplot as plt
```

### 1. Defining a Model Without Levers


```python
model = Model('lakeproblem', function=lake_problem)
 
#specify uncertainties
model.uncertainties = [RealParameter('b', 0.1, 0.45),
                       RealParameter('q', 2.0, 4.5),
                       RealParameter('mean', 0.01, 0.05),
                       RealParameter('stdev', 0.001, 0.005),
                       RealParameter('delta', 0.93, 0.99)]

#specify outcomes
model.outcomes = [ScalarOutcome('max_P'),
                  ScalarOutcome('utility'),
                  ScalarOutcome('inertia'),
                  ScalarOutcome('reliability')]
```

* Run the model with a sequential evaluator.


```python
%%time
from ema_workbench import (SequentialEvaluator, ema_logging, perform_experiments)

ema_logging.log_to_stderr(ema_logging.INFO)
 
with SequentialEvaluator(model) as evaluator:
    experiments, outcomes = evaluator.perform_experiments(scenarios=1000, policies=0)
```

    [MainProcess/INFO] performing 1000 scenarios * 1 policies * 1 model(s) = 1000 experiments
    [MainProcess/INFO] performing experiments sequentially
    [MainProcess/INFO] 100 cases completed
    [MainProcess/INFO] 200 cases completed
    [MainProcess/INFO] 300 cases completed
    [MainProcess/INFO] 400 cases completed
    [MainProcess/INFO] 500 cases completed
    [MainProcess/INFO] 600 cases completed
    [MainProcess/INFO] 700 cases completed
    [MainProcess/INFO] 800 cases completed
    [MainProcess/INFO] 900 cases completed
    [MainProcess/INFO] 1000 cases completed
    [MainProcess/INFO] experiments finished


    Wall time: 42.3 s


### 1.1 Visualizing results with pair plots


```python
fig, axes = pairs_plotting.pairs_scatter(experiments, outcomes, group_by="policy")

fig.set_size_inches(14,14)

plt.show()
```

### 1.2 Evaluating the Effects of Uncertainties
* We visualized the relation between parameters 'b' and 'q' in terms of the impact on outcome max_P.
* It appears that with low 'b' values the phosphorus levels are high.


```python
P = outcomes['max_P']

fig, ax = plt.subplots(figsize=(6,6))

sc = ax.scatter(experiments.b, experiments.q, c=P)
plt.colorbar(sc)
ax.set_xlabel('b')
ax.set_ylabel('q')
plt.show()
```

* No apparently discernible relation between 'delta' and 'mean' variables in terms of impact on outcome max_P.


```python
P = outcomes['max_P']

fig, ax = plt.subplots(figsize=(6,6))

sc = ax.scatter(experiments.delta, experiments["mean"], c=P)
plt.colorbar(sc)
ax.set_xlabel('delta')
ax.set_ylabel('mean')
plt.show()
```

### 2. Defining a Model with 100 Policy Levers (bounded between 0 and 0.1)


```python
model = Model('lakeproblem', function=lake_problem)

model.levers = [RealParameter('l' + str(i), 0, 0.1) for i in range(4)]
 
#specify uncertainties
model.uncertainties = [RealParameter('b', 0.1, 0.45),
                       RealParameter('q', 2.0, 4.5),
                       RealParameter('mean', 0.01, 0.05),
                       RealParameter('stdev', 0.001, 0.005),
                       RealParameter('delta', 0.93, 0.99)]

#specify outcomes
model.outcomes = [ScalarOutcome('max_P'),
                  ScalarOutcome('utility'),
                  ScalarOutcome('inertia'),
                  ScalarOutcome('reliability')]
```

* Run the model with a sequential evaluator.


```python
%%time
from ema_workbench import (SequentialEvaluator, ema_logging, perform_experiments)

ema_logging.log_to_stderr(ema_logging.INFO)
 
with SequentialEvaluator(model) as evaluator:
    experiments, outcomes = evaluator.perform_experiments(scenarios=1000, policies=4)
```

    [MainProcess/INFO] performing 1000 scenarios * 4 policies * 1 model(s) = 4000 experiments
    [MainProcess/INFO] performing experiments sequentially
    [MainProcess/INFO] 400 cases completed
    [MainProcess/INFO] 800 cases completed
    [MainProcess/INFO] 1200 cases completed
    [MainProcess/INFO] 1600 cases completed
    [MainProcess/INFO] 2000 cases completed
    [MainProcess/INFO] 2400 cases completed
    [MainProcess/INFO] 2800 cases completed
    [MainProcess/INFO] 3200 cases completed
    [MainProcess/INFO] 3600 cases completed
    [MainProcess/INFO] 4000 cases completed
    [MainProcess/INFO] experiments finished


    Wall time: 2min 45s


3. Explore the behavior of the system over 1000 scenarios for 4 randomly sampled candidate strategies.
    * visualize the outcomes of interest
    * what can you say about how the release decision influences the system?

### 2.1 Visualize Results


```python
outcomes_df = pd.DataFrame.from_dict(outcomes)
# melted = pd.melt(outcomes_df, value_vars = outcomes_df.columns)
```


```python
sns.jointplot(data=outcomes_df, y="max_P", x="utility", kind='kde')
```

* Visualize results with pair plots.


```python
fig, axes = pairs_plotting.pairs_scatter(experiments, outcomes, group_by="policy")

fig.set_size_inches(12,12)

plt.show()
```

### 3. Defining a Model with 100 Policy Levers (bounded between 0 and 0.5)


```python
model = Model('lakeproblem', function=lake_problem)

model.levers = [RealParameter('l' + str(i), 0, 0.5) for i in range(100)]
 
#specify uncertainties
model.uncertainties = [RealParameter('b', 0.1, 0.45),
                       RealParameter('q', 2.0, 4.5),
                       RealParameter('mean', 0.01, 0.05),
                       RealParameter('stdev', 0.001, 0.005),
                       RealParameter('delta', 0.93, 0.99)]
#specify outcomes
model.outcomes = [ScalarOutcome('max_P'),
                  ScalarOutcome('utility'),
                  ScalarOutcome('inertia'),
                  ScalarOutcome('reliability')]
```

* Running the model with a MultiprocessingEvaluator


```python
%%time
from ema_workbench import (MultiprocessingEvaluator, ema_logging,
                           perform_experiments)
ema_logging.log_to_stderr(ema_logging.INFO)
 
with MultiprocessingEvaluator(model) as evaluator:
    experiments, outcomes = evaluator.perform_experiments(scenarios=1000, policies=4)
```

    [MainProcess/INFO] pool started
    [MainProcess/INFO] performing 1000 scenarios * 4 policies * 1 model(s) = 4000 experiments
    [MainProcess/INFO] 400 cases completed
    [MainProcess/INFO] 800 cases completed
    [MainProcess/INFO] 1200 cases completed
    [MainProcess/INFO] 1600 cases completed
    [MainProcess/INFO] 2000 cases completed
    [MainProcess/INFO] 2400 cases completed
    [MainProcess/INFO] 2800 cases completed
    [MainProcess/INFO] 3200 cases completed
    [MainProcess/INFO] 3600 cases completed
    [MainProcess/INFO] 4000 cases completed
    [MainProcess/INFO] experiments finished
    [MainProcess/INFO] terminating pool


    Wall time: 1min 14s


#### The multiprocessing evaluator was much faster than the sequential evaluator.


```python
fig, axes = pairs_plotting.pairs_scatter(experiments, outcomes, group_by="policy")

fig.set_size_inches(12,12)

plt.show()
```

*Interpretation of Results:*
* Inertia is obviously different per policy because the amount of release is different per policy; inertia is directly affected by release.
* Inertia does not affect Max_Phosphorus, reliability and utility when release is between 0 and 0.1. However, when release >0.1, it appears to affect utility and max phosphorus (see graph below for policy between 0 and 0.5).

4. If you have not used parallelization in the foregoing, try to adapt your code to use parallelization. The workbench comes with two evaluators for parallelization. The `MultiProcessingingEvaluator` and the `IpyparallelEvaluator`. When can you use each? Adapt your code from above and sue the `MultiProcessingingEvaluator`. Use the time library to check how much faster the computation for 1000 scenarios completes. 

Using multiprocessing within a jupyter notebook is tricky. On linux it will work in general just fine. On mac it depends on the version of Mac OS and the version of Python. If you are on the latest version of Mac OS in combination with Python 3.8, it might work but no guarantees. On older versions of Python it should work fine. On Windows it is always a problem. 

The underlying explanation is quite technical. It has to do with how your operating system creates the additional python processes. On windows, and the latest version of Mac OS in combination with Python 3.8. A completely new Python process is spawned. This new process does **not** inherit what is defined in memory of the parent process. The new child process will try to replicate what is in memory of the parent process by executing many of the import statements that have also been executed within the python process. Thus, if you define a model in the main process, it is not guaranteed to be known in the child processes. This is in particular true if you define the model within a jupyter notebook. Then the child processes will **never** know this function. Within jupyter notebooks, therefore, the best practice is to define your model within a .py file and import this .py file into the notebook. Now, each of the child processes will also execute this import statement and thus know the function. 


83eb85a453f8282dc978c7c8a3629b81825ac55c
ipynb
Jupyter Notebook
epa1361_open/Week 1-2 - general intro to exploratory modelling/assignment 2 - lakemodel.ipynb
lnicoletti/EPA1361-G21
e8dce7169101fd620a3e36f7c38df7f4db598a69
MIT
epa1361_open/Week 1-2 - general intro to exploratory modelling/assignment 2 - lakemodel.ipynb
lnicoletti/EPA1361-G21
e8dce7169101fd620a3e36f7c38df7f4db598a69
MIT
2020-04-28T16:00:57.000Z
2020-04-28T16:17:25.000Z
epa1361_open/Week 1-2 - general intro to exploratory modelling/assignment 2 - lakemodel.ipynb
lnicoletti/EPA1361-G21
e8dce7169101fd620a3e36f7c38df7f4db598a69
MIT
2021-04-19T11:14:49.000Z
2021-04-19T11:14:49.000Z
Qwen/Qwen-72B
1. YES
2. YES