17 Lake Problem

.title[Simulation-Optimization and the Shallow Lake Problem] <br> .subtitle[BEE 4750/5750] <br> .subtitle[Environmental Systems Analysis, Fall 2022] <hr> .author[Vivek Srikrishnan] <br> .date[October 31, 2022]

---

# Outline

<hr>

1. Project Due Dates
2. HW 4 Due Thursday
3. Questions?
4. Challenges to Mathematical Programming
5. Shallow Lake Model
6. Simulation-Optimization

---

# Poll

<hr>

URL: [https://pollev.com/vsrikrish](https://pollev.com/vsrikrish)

Text: **VSRIKRISH** to 22333, then message]

---

***Any questions?***

---

# Last Class

<hr>

* Multiple Point Sources/Receptors
  * Box Modeling for Airsheds

---

# Mathematical Programming

<hr>

**Previous weeks**: "Mathematical Programming" with formal decision models: objectives and constraints as functions (linear or otherwise) of decision variables.

Mathematically programming is nice when possible because you can guarantee that you can find optima for certain classes of problems.

---

# What If We Can't Write Down a Mathematical Program?

<hr>

But what if our constraints and objectives cannot be written down as a mathematical program?

This often occurs with systems models, due to complex, nonlinear dynamics.

---

# Systems Dynamics: Feedback Loops

<hr>

Dampening feedback loops are associated with *stable* equilibria, while reinforcing feedback loops are associated with *instability*.]

---

# Systems Dynamics: Bifurcations

<hr>

These thresholds are sometimes referred to as "tipping points" and are often associated with a *stable* state equilibrium transitioning to an *unstable* one.]

---

# Systems Dynamics: Emergence

<hr>

**Emergence** refers to when simple rules governing interactions result in macro-scale "ordered" patterns and possibly unexpected behaviors.

* Schelling segregation model (left)
  * Macroeconomic systems
  * Snowflake formation

]

.right-column[ .center[![Schelling Model Animation](https://upload.wikimedia.org/wikipedia/commons/e/e8/Schellings.gif)] .cite[Source: Blaqdolphin - Own work, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=91228415] ]

---

# Example: Shallow Lake Problem

<hr>

.left-column[ Suppose we have a town which emits phosphorous into a lake as a by-product of economic activity.

The lake also receives phosphorous from non-point source runoff. There is some nutrient cycling in the lake which naturally reduces phosphorous concentrations. ]

---

# Shallow Lake Model

<hr>

Suppose the lake dynamics are given by:

$$
X_{t+1} = X_t + a_t + y_t + \frac{X_t^q}{1 + X_t^q} - bX_t, \quad y_t \sim LogNormal(\mu, \sigma^2)
$$

---

# Impacts of Parameters on Equilibria

<hr>

.center[![Impact On Flux](figures/lake-nox.svg)] .cite[Adapted from [Quinn et al (2017)](https://dx.doi.org/10.1016/j.envsoft.2017.02.017)] ]

---

# Can We Write Down a Mathematical Program For This Problem?

<hr>

Our objective might be to maximize $\sum_{t=1}^T a_t$ (as a proxy for economic activity), while keeping a low probability of eutrophication.

Straightforward to compute probability of eutrophication for a given input sequence $a_t$ using Monte Carlo simulation for $y_t$, but highly nontrivial to write this out in analytic form.

Even if we did, it'd be *very* nonlinear and definitely nonconvex, which rules out most forms of mathematical programming.

---

# Simulation-Optimization

<hr>

*So how do we make decisions?*

Recall our previous comment about Monte Carlo!

We can use a search algorithm to explore the exceedance probabilities using our simulation model, and find the "optimum" sequence of releases that keeps the probability of eutrophification at some specified level.

---

# Why "Optimum" in Quotes?

<hr>

Unlike mathematical programming, **there is usually no guarantee that a solution found this way is truly optimal**.

Simulation-optimization methods typically rely on *heuristics* to decide that a solution is good enough. These can include

* number of evaluations/iterations; or
  * lack of improvement of solution.

---

# Challenges to Simulation-Optimization In General

<hr>

* **Monte Carlo Error**: If constraints or objective is probabilistic, how many Monte Carlo runs are needed to ensure difference in function values is "real" and not stochastic noise.
  * **Computational**: Can be expensive depending on the simulation model.
  * **Local vs. Global Optima**: Depending on type of search algorithm, may not be able to guarantee more than a local optimum.

But, like Monte Carlo, often other methods are impossible or worse.

---

# Gradient Estimation

<hr>

.left-column[Find estimate of gradient near current point and step in positive/negative direction (depending on max/min).

$$
x_{n+1} = x_n \pm \alpha_n \nabla f(x_n)
$$

But: may not find global optimum; stepsize plays a big role in convergence. ]

---

# Random Search

<hr>

* Evolutionary Algorithms fall into this category.
  * No guarantees! Need to evaluate lots and convince yourself you've done enough.]

---

# Random Search with Constraints

<hr>

---

# Random Search In Julia

<hr>

Two main packages:

* [`Metaheuristics.jl`](https://docs.juliahub.com/Metaheuristics/aJ70z/3.2.12/)
  * [`BlackBoxOptim.jl`](https://github.com/robertfeldt/BlackBoxOptim.jl)

---

# Quick Example (Adapted From `Metaheuristics.jl` Tutorial)

<hr>

```julia
using Metaheuristics
using Distributions
using Random
Random.seed!(1)
# define function to optimize
function f(x)
    lnorm = LogNormal(0.25, 1)
    y = rand(lnorm)
    return sum(x .- y)^2
end
fobj(x) = mean([f(x) for i in 1:1000])
```

---

# Quick Example (Adapted From `Metaheuristics.jl` Tutorial)

<hr>

```julia
# set bounds
bounds = [0.0 1000.0]'
# optimize
results = optimize(fobj, bounds, DE())
results.best_sol
```

```
(f = 3.6109e+00, x = [2.368314566931006])
```

---

# `Metaheuristics.jl` Algorithms

<hr>

`Metaheuristics.jl` contains [a number of algorithms](https://docs.juliahub.com/Metaheuristics/aJ70z/3.2.12/algorithms/), covering a number of single-objective and multi-objective algorithms.

We won't go into details here, and will just stick with `DE()` (differential evolution) in our examples.

---

# Lake Example (Including Constraints)

<hr>

Lake Problem: Let's try to keep the probability of exceeding the critical value (leading to eutrophication) less than 80%.

```julia
using Metaheuristics
using Distributions
using Plots
using Random
using Roots
Random.seed!(1)

# set parameters
q = 2.5
b = 0.4
T = 100
nsamples = 1000
```

```
1000
```

---

# Lake Example (Including Constraints)

<hr>

```julia
# get NPS inflows
lnorm = LogNormal(log(0.03), 0.1)
y = rand(lnorm, (T, nsamples))
# find critical value
crit(x) = (x^q/(1+x^q)) - b*x
Xcrit = find_zero(crit, (0.1, 1.5))
```

```
0.6678778690448219
```

]

---

# Lake Example (Including Constraints)

<hr>

```julia
# define lake model
function lake(a, y, q, b, T)
    X = zeros(T+1, size(y, 2))
    # calculate states
    for t = 1:T
        X[t+1, :] = X[t, :] .+ a[t] .+ y[t, :] .+ (X[t, :].^q./(1 .+ X[t, :].^q)) .- b.*X[t, :]
    end
    return X
end
function lake_opt(a, y, q, b, T, Xcrit)
    X = lake(a, y, q, b, T)
    # calculate exceedance of critical value
    Pexceed = sum(X[101, :] .> Xcrit) / size(X, 2)
    failconst = [Pexceed - 0.2]
    return mean(a), failconst, [0.0]
end
```

---

# Lake Example (Including Constraints)

<hr>

```julia
# set bounds
bounds = [0.0ones(T) 0.1ones(T)]'
# optimize
obj(a) = lake_opt(a, y, q, b, T, Xcrit)
options = Options(f_calls_limit=1000)
results = optimize(obj, bounds, DE(options=options))
```

<pre>
+=========== RESULT ==========+
  iteration: 1
    minimum: 0.0442926
  minimizer: [0.039240054578424256, 0.09338124781042632, 0.08835418983139154, 0.011427041624799706, 0.048924893262125936, 0.037265780904226796, 0.09599174108992913, 0.016868255218518348, 0.06501532126834804, 0.012504347730588805, 0.000746213129172535, 0.08767941825273518, 0.06444447819690417, 0.054203305489866606, 0.04515010867664222, 0.03718958851984027, 0.025678307010852875, 0.002614556880381591, 0.031664324511728625, 0.07500180739034422, 0.043465369178345536, 0.02618137451907561, 0.012125732365595233, 0.0007835638382379329, 0.050504057095811676, 0.01973042349246107, 0.033136951386881955, 0.01677370068553733, 0.03257252223625352, 0.03630833605591843, 0.09431288747936745, 0.09808233362986095, 0.00500268791860089, 0.05647886662206815, 0.03155775927695148, 0.07254766454510762, 0.004762093575137772, 0.027549599629570133, 0.07984938749477144, 0.05387033975708925, 0.0830215466837763, 0.03715250339291026, 0.06260338936336694, 0.05665233610796105, 0.020867171309697265, 0.003031407257133101, 0.07038442196539017, 0.081973209766261, 0.049986184715054574, 0.0024662374832810086, 0.006323147080792702, 0.009657881913163725, 0.02315890085132412, 0.09193974477638128, 0.03760366519862346, 0.07564759057190282, 0.030848913630144338, 0.01669100231588613, 0.015936053575238196, 0.07402535946544608, 0.05861471118712115, 0.024710497229536024, 0.04584411954893829, 0.006433931745784916, 0.02525147898785858, 0.0927269006025585, 0.007840973215041902, 0.03206206338343415, 0.09203388652557211, 0.013111902500402828, 0.004544827274546015, 0.08026903850057003, 0.020570018234940746, 0.011876941026341604, 0.004183774888220338, 0.09852820698013408, 0.057271685561092404, 0.029796571502513525, 0.015045512458381129, 0.08260028598701494, 0.0036321355165284453, 0.059813074523417824, 0.0862832791944941, 0.075865373638978, 0.04496204337475345, 0.04088387618008386, 0.013879043219367927, 0.05961114061643075, 0.039901062880284735, 0.08933411431006494, 0.004006065449934948, 0.012758921045441619, 0.048643482508816366, 0.0924091515586886, 0.07363173417161059, 0.08219603133985784, 0.008796614369613464, 0.036447873444808855, 0.09373048347884377, 0.08024559983252538]
    f calls: 1000
  feasibles: 1 / 1000 in final population
 total time: 27.1178 s
stop reason: Maximum objective function calls exceeded.
+============================+
</pre>

---

# Lake Example (Including Constraints)

<hr>

```julia
using Plots

a = results.best_sol.x
X = lake(a, y, q, b, T)
plot(X, alpha=0.1,
    guidefontsize=14, tickfontsize=12,
    legendfontsize=12, label=:false,
    legend=:topleft)
hline!([Xcrit], color=:red, linestyle=:dot,
    label="Critical Value")
plot!(size=(400, 400))
```

]

---

# Lake Example (Including Constraints)

<hr>

```julia
histogram(X[101, :], xlabel="End P Concentration", ylabel="Count",
    guidefontsize=14, tickfontsize=12,
    legendfontsize=12, label=:false)
vline!([Xcrit], color=:red, linestyle=:dot,
    label="Critical Value")
plot!(size=(400, 400))
```

]

---

# Key Takeaways

<hr>

* Many challenges to using mathematical programming for general systems analysis.
  * Can use simulation-optimization approaches.
  * But be careful of computational expense and convergence!

---

<hr>

# Next Class

<hr>

* In-Class Coding Exercise with the Lake Problem