Skip to content

Instantly share code, notes, and snippets.

@jdherman
Last active August 29, 2015 14:07
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jdherman/d519b3ef81840fc27d96 to your computer and use it in GitHub Desktop.
Save jdherman/d519b3ef81840fc27d96 to your computer and use it in GitHub Desktop.
Lake problem (multiobjective) in python with Borg wrapper
from __future__ import division
import numpy as np
from scipy.optimize import brentq as root
from borg import *
nvars = 100 # decision variables: pollution for each timestep
nobjs = 4
nconstr = 1
nsamples = 100 # 100 scenarios of natural inflows
b = 0.42 # decay rate for P in lake (0.42 = irreversible)
q = 2.0 # recycling exponent
alpha = 0.4 # utility from pollution
delta = 0.98 # (1-r) discount rate
Pcrit = root(lambda x: x**q/(1+x**q) - b*x, 0.01, 1.5)
def lake(*decisions):
objs = [0.0]*nobjs
constr = [0.0]
X = np.zeros((nvars,))
average_daily_P = np.zeros((nvars,))
decisions = np.array(decisions)
for s in xrange(nsamples):
X[0] = 0
natural_inflows = np.random.lognormal(-3.52, 0.105, size=nvars)
for t in xrange(1,nvars):
X[t] = (1-b)*X[t-1] + X[t-1]**q/(1+X[t-1]**q) + decisions[t-1] + natural_inflows[t-1]
average_daily_P[t] += X[t]/nsamples
objs[3] -= np.sum(X < Pcrit)/(nsamples*nvars) # Maximize timesteps satisfying X < Pcrit
objs[0] = np.max(average_daily_P) # Minimize the maximum daily P in lake
objs[1] = -1*np.sum(alpha*decisions*np.power(delta,np.arange(nvars))) # Maximize the average sum of discounted benefits
objs[2] = -1*np.sum(np.diff(decisions) > -0.02)/(nvars-1) # Maximize timesteps meeting inertia threshold
constr[0] = Constraint.greaterThan(-1*objs[3], 0.95)
return (objs, constr)
borg = Borg(nvars, nobjs, nconstr, lake)
borg.setBounds(*[[0.0, 0.1]]*nvars)
borg.setEpsilons(*[0.01]*nobjs)
result = borg.solve({"maxEvaluations":10000})
for solution in result:
print solution.getObjectives() + solution.getConstraints()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment