Last active
August 29, 2015 14:07
-
-
Save jdherman/d519b3ef81840fc27d96 to your computer and use it in GitHub Desktop.
Lake problem (multiobjective) in python with Borg wrapper
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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