forked from ayman510/WASH
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGAMS code
406 lines (282 loc) · 15.8 KB
/
GAMS code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
*Last update: Oct 8, 2015
*Commented by David Rosenberg - Oct 15, 2015. Comments are prefaced by *DER - ...
$Title Watershed Area of Suitable Habitat (WASH) model and an application on the Lower Bear River Watershed
$OnText
WASH is a system optimization model that finds the balance between environmental needs and human needs of available water in the watershed.
WASH receommends allocation of water to improve habitat quality by incorporating watershed-scale habitat quality indexes as objectives
to maximize in a systems optimization model. WASH measures physically-available suitable habitat in three main components of most watersheds, namely
riverine, floodplain and wetland habitats while maintaining water for human beneficial uses. WASH also highlights promising restoration and
conservation sites. In addition, WASH considers, quantifies, and measures multiple sources of uncertainty in management decisions and proposes and rank
alternatives for managers to restore and protect natural habitat. WASH formulation is generic and adaptable to other regulated river systems and for species
of concern. I apply the model to the Lower Bear River, Utah to guide existing habitat conservation efforts, recommend water allocation and show tradeoffs between
human and habitat use of the available flow.
Citation:
??
Documentation on GitHub:
???
####################
Programmed by Ayman H. Alafifi
Dept. of Civil & Env. Engineering and Utah Water Research Lab
Utah State University
Last Updated: May 22, 2015
####################
$OffText
*Declare sets
SETS s sub-indicators
j river netowrk nodes
wt(j) impounded wetlands
dem(j) demand sites
v(j) reservoirs
;
Alias (j,k);
SETS
* link(j,k) river links from node j to k /#j.#k/
dv(j,k) diversions
rf(j,k) return flow
f type of fish species /trout/
n type of vegetation species /cottonwood/
t timesteps in months /1*12/
R_par_indx indexes for the RSI relationship /R_par1, R_par2, R_par3, R_par4/
W_par_indx indexes for the WSI monthly realtionship /W_par1, W_par2/
EleV_par_indx indexes for the reservoir volume elevation curve /EleV_par1, EleV_par2, EleV_par3/
sf_par_indx index for stage-flow relationship parameters /sf1_par, sf2_par/
wf_par_indx index for width-flow relationship parameters /wf1_par, wf2_par/
wsi_par_indx index for wetlands water availbility relationship /wsi_par1, wsi_par2/
;
*Define Parameters and Scalars
PARAMETERS
linkexist(j,k) link from j to k exists (1=yes and 0=no)
envSiteExist(j,k) environmental site exists (1=yes and 0=no) where sensitive habitat is located and data were collected
wght(s,j,t) spatial and temporal weights set by stakeholders [0: no important - 1: important]
theta(j,k) topological slope
aw(wt,t) Impounded wetlands area [m2]
fmin(j,k,n) minimum flood depth required for species n
fmax(j,k,n) maximum flood depth required for species n
bankdepth(j,k,t) Bank depth at site i
lss(j,k,t) net losses on link j entering k [%] - applicable to all nodes
evap(v,t) evaporation losses [m per month]
ar(v,t) reservoir area [m2]
cons(dem,t) consupmtive use fraction expressed as [%] of inflow received at site d
lng(j,k) length of river segment
minstor(v) inactive reservoir storage [m3]
maxstor(v) storage capacity [m3]
dReq(dem,t) demand requirements [m3 per month]
dCap(j,dem,t) diversion link capacity [m3 per month]
dmin(j,k,t) minimum instream flow requirement [m3 per month]
cst(f) unit cost of implementing management objective
vegD(j,k,t,n) Distance between floodplain vegetation and river bank
EleV_par(v, EleV_par_indx) Reservoir elevation volume curve parameters
R_par(f,t,R_par_indx) Riverine Suitability Index equation parameters
W_par(wt,t,W_par_indx) Wetlands Subindicator equation parameters
sf_par(j,k,sf_par_indx) Stage-flow relationship parameters
wf_par(j,k,wf_par_indx) Width-flow relationship parameters
wsi_par(wt,t,wsi_par_indx) Wetlands suitability index monthly relationship parameters
;
SCALARS
b budget [in $] ;
* Define a piecewise polynomial function and its lower boundaries
*This function has two segments, each row defines: FuncInd.SegInd, leftBound, Coef0, Coef1, Coef2 ...
* where CoefX is the Xth degree coefficient of the polynomial corresponding to this segment.
* For example: 0.3X**2 -2.7*X +2.4 will have coef in the order of 2.4, -2.7, 0.3
*In pwdata table:
*equation 1: Riverine Suitability Index in Jan-Feb: RSI01 for Spawn/eggs
*equation 2: RSI in Mar-Aug RSI02 for Fry, equation3: RSI in Aug-Dec RSI03 for Juvenile
*equation4: h index
*equation5: m index
Table pwpdata(*,*,*) '1st index: function number, 2nd index: segment number, 3rd index: degree'
leftBound 0 1 2
1.1 0 0 0 0
1.2 0.25 0.25 0 0
1.3 0.4 0.5 0 0
1.4 0.5 1 0 0
2.1 0 0 0 0
2.2 0.25 0.25 0 0
2.3 0.5 0.5 0 0
2.4 1.0 1 0 0
3.1 0 0 0 0
3.2 0.25 0.75 0 0
3.3 0.4 1.25 0 0
3.4 0.5 1.5 0 0
4.1 0 0.2 0 0
4.2 150 1.0 0 0
5.1 -1.0 0.1 0 0
5.2 0.0 0.2 0 0
5.3 0.2 0.8 0 0
5.4 0.8 0.2 0 0
* Write pwp data to gdx file read by external library
$gdxout pwp.gdx
$unload pwpdata
$gdxout
* Load extrinsic function
$funclibin pwplib pwpcclib
function pwp /pwplib.pwpfunc/ ;
* Read sets and parameter input values from Excel
$CALL GDXXRW.EXE input=WASHData-sample.xlsx output=WASH.gdx Set=s rng=SubInd!A1 Rdim=1 Set=j rng=Nodes!A1 Rdim=1 Set=wt rng=Wetlands!A1 Rdim=1 Set=dem rng=Demand!A1 Rdim=1 Set=v rng=Reservoirs!A1 Rdim=1 par=envSiteExist rng=EnvSite!A1 Rdim=1 Cdim=1 Set=dv rng=Diversions!A1 Rdim=2 Set=rf rng=ReturnFlow!A1 Rdim=2 Set=f rng=FishSpp!A1 Rdim=1 Set=n rng=VegSpp!A1 Rdim=1 Set=t rng=Month!A1 Rdim=1 par=bankdepth rng=Bankfull!A1 Rdim= 2 Cdim=1 par=wght rng=weights!A1 Rdim=3 Cdim=1 par=theta rng=theta!A1 Rdim=1 Cdim=1 par=aw rng= AW Rdim=2 par=fmin rng=Fmin Rdim=2 Cdim=1 par=fmax rng=Fmax Rdim=2 Cdim=1 par=lss rng=lss Rdim=2 Cdim=1 par=evap rng=evap Rdim=1 Cdim=1 par=ar rng=ResA Rdim=1 Cdim=1 par=EleV_par rng=ResElevVol Rdim=1 Cdim=1 par=Cons rng=Cons Rdim=1 Cdim=1 par=lng rng=Length Rdim=1 Cdim=1 par=minstor rng=inactive Rdim=1 par=maxstor rng=capacity Rdim=1 par=dReq rng=demandReq Rdim=1 Cdim=1 par=dCap rng=divCap Rdim=1 Cdim=1 par=dmin rng=Instream Rdim=2 Cdim=1 par=cst rng=UnitCost Rdim=1 par=VegD rng=VegD Rdim=3 Cdim=1 par=Linkexist rng=Connect!A1 Rdim=1 Cdim=1 par=R_par rng=R_par!A1 Rdim=2 Cdim=1 par=W_par rng=W_par!A1 Rdim=2 Cdim=1 par=sf_par rng=StageFlow Rdim=2 Cdim=1 par=wf_par rng=WidthFlow Rdim=2 Cdim=1 par=wsi_par rng=wp Rdim=2 Cdim=1 par=b rng=Budget
$GDXIN WASH.gdx
$LOAD s
$LOAD j
$LOAD wt
$LOAD dem
$LOAD v
$LOAD envSiteExist
$LOAD linkexist
$LOAD dv
$LOAD rf
$LOAD f
$LOAD n
$LOAD t
$LOAD R_par_indx
$LOAD W_par_indx
$LOAD sf_par_indx
$LOAD wf_par_indx
$LOAD wsi_par_indx
$LOAD wght
$LOAD theta
$LOAD aw
$LOAD fmin
$LOAD fmax
$LOAD lss
$LOAD evap
$LOAD EleV_par
$LOAD ar
$LOAD cons
$LOAD lng
$LOAD minstor
$LOAD maxstor
$LOAD dReq
$LOAD dCap
$LOAD dmin
$LOAD cst
$LOAD vegD
$LOAD R_par
$LOAD W_par
$LOAD sf_par
$LOAD wf_par
$LOAD wsi_par
$LOAD b
$GDXIN
*Define Variables
VARIABLES
Z Objective Function value of WASH
*DER - be a bit more specific in description of objective function
Ind(s,i,t) sub-indicator
Q(j,k,t) flow in link at environmental sites [m3 per month]
D(i,t) water depth[m]
A(i,t) Channel surface area [m2]
WD(i,t) Channel width [m]
RV(i,t,n) Area to re-vegetate [m2]
C(i,t,n) vegetation cover [m2]
O(i,t) overbank water level [m]
E(i,t) distance of flood [m]
K(i,t,n) recessiona rate
WA(wt,t) water available to wetlands [ha-m per month]
RR(v,t) Reservoir Releases [m3 per month]
STOR(v,t) Reservoir Storage [m3]
ELEV(v) Reservoir elevation [m]
h(i,t,n) index of flood distance [0.2-1]
m(i,t,n) index of inundation depth [0.2 - 1]
g(i,t,n) index of stream recession rate [0.2 - 1]
WSI(wt,t) Wetlands suitability index
RSI(i,t,f) Rivering Suitability Index
FCI(i,t,n) Floodplain Suitability Index
R(i,t) Riverine habitat sub-indicator
F(i,t) Floodplain connectivity sub-indicator
W(i,t) Wetlands habitat sub-indicator
;
*Define equations
EQUATIONS
EQ1 Z Objective function: Watershed Area of Suitable Habitat [m2]
EQ1a(i,t) ind Summation of the three perfomance indicators [m2]
EQ2(i,t) R Riverine Habitat [m2]
EQ2a(i,1,f) RSI Riversine Suitability Index [0-1] for Jan
EQ2b(i,2,f) RSI for Feb
EQ2c(i,3,f) RSI for Mar
EQ2d(i,4,f) RSI for Apr
EQ2e(i,5,f) RSI for May
EQ2f(i,6,f) RSI for Jun
EQ2g(i,7,f) RSI for Jul
EQ2h(i,8,f) RSI for Aug
EQ2i(i,9,f) RSI for Sep
EQ2g(i,10,f) RSI for Oct
EQ2k(i,11,f) RSI for Nov
EQ2l(i,12,f) RSI for Dec
EQ3(i,t) F Floodplain habitat [m2]
EQ3a(i,t,n) FCI floodplains suiability index [0-1]
EQ3b(i,t,n) h binary index of flood distance
EQ4a(i,t) O overbank level [m]
EQ4(i,t) E flood distance [m]
EQ5(n) m index for floodplain depth [0.2-1]
EQ6(n) g index for flood recession [0.2- 1]
EQ7(i,t,n) K value for g index
EQ7a(wt,t) WSI wetlands suitability index
EQ8(wt,t) Impounded Wetlands [m2]
EQ9(v,t) mass balance at reservoirs
EQ9a(v,t) reservoir elevation - storage relationship
EQ10(j,k,t) mass balance at each node
EQ11(j,d,t) mass balance at demand sites
EQ12(wt,j,t) mass balance at wetlands sites
EQ13(i,t) stage flow relationship
EQ14(i,t) width flow relationship
EQ15(i,t) Channel surface area [m2]
EQ16(i,t,n) Vegetation cover mass balance
EQ17(v,t) Reservoir storage cannot go below inactive zone
EQ17a(v,t) Reservoir storage cannot exceed storage capacirt
EQ18(j,d,t) Diversions to demand sites should meet requirements
EQ19(j,d,t) Diversions to demand sites should not exceed diversion capacity
EQ20(i,t) minimum instream flow requirements
EQ21(s,i,t) Limitation on Budget
;
EQ1.. Z =e= sum((s,i,t), wght(s,i,t) * Ind(s,i,t)) ;
*EQ1a(i,t).. Ind(s,i,t) =e= sum((i,t), R(i,t), F(i,t), W(wt,t)) OLD EQUATION ;
*DER - EQ1a will give several errors. The index s is not controlled (only appears on the left-
* hand side. R and F are for links but W is for a node (wetland units).
* Intead let's rewrite so the set s is controlled and we set R(i,t) equal to the
* first element of the set, F(i,t) equal to the second, etc... I think it will also be easier
* to change the logic so wetland units are defined as links not a node. In that case:
EQ1a(s,i,t) Ind(s,i,t) =e= R(i,t)$(ord(s) eq 1) + F(i,t)$(ord(s) eq 2) + W(i,t)$(ord(s) eq 3);
EQ2(i,t).. R(i,t) =e= prod((f), RSI(i,t,f) * A(i,t)) ;
*EQ2a(i,t,f).. RSI(i,t,f) =e= R_par(f,t,R_par_indx) * D(i,t) OLD EQUATION;
EQ2a(i,1,f).. RSI(i,1,f) =e= PWP(1,D(i,1));
EQ2b(i,2,f).. RSI(i,2,f) =e= PWP(1,D(i,2));
EQ2c(i,3,f).. RSI(i,3,f) =e= PWP(2,D(i,3));
EQ2d(i,4,f).. RSI(i,4,f) =e= PWP(2,D(i,4));
EQ2e(i,5,f).. RSI(i,5,f) =e= PWP(2,D(i,5));
EQ2f(i,6,f).. RSI(i,6,f) =e= PWP(2,D(i,6));
EQ2g(i,7,f).. RSI(i,7,f) =e= PWP(2,D(i,7));
EQ2h(i,8,f).. RSI(i,8,f) =e= PWP(2,D(i,8));
EQ2i(i,9,f).. RSI(i,9,f) =e= PWP(3,D(i,9));
EQ2j(i,10,f).. RSI(i,10,f) =e= PWP(3,D(i,10));
EQ2k(i,11,f).. RSI(i,11,f) =e= PWP(3,D(i,11));
EQ2l(i,12,f).. RSI(i,12,f) =e= PWP(3,D(i,12));
EQ3(i,t).. F(i,t) =e= sum((n),FCI(i,t,n) * C(i,t,n)) ;
EQ3a(i,t,n).. FCI(i,t,n) =e= h(i,t,n) * m(i,t,n) * g(i,t,n) ;
*EQ3b(i,t,n).. h(i,t,n) =e= 0.2$(E(i,t) gt vegD(i,t)) OLD EQUATION ;
EQ3b(i,t,n).. h(i,t,n) =e= PWP(4,E(i,t,n));
EQ4(i,t).. E(i,t) = O(i,t)/theta(i,t) ;
EQ4a(i,t).. O(i,t) = D(i,t) - bankdepth(i,t) ;
EQ5(n).. m(i,t,n) =e= PWP(5,O(i,t,n)) ;
EQ6(n).. g(i,t,n)=e= 0.94* ( exp(-0.5*(ln(K(i,t,n)/1.28))/0.99)**2) ;
EQ7(i,t,n)..
K(i,t,n) = Abs(100*(ln(D(t) - ln(D(t))))/30) ;
EQ8(wt,t).. W(wt,t) =e= WSI(wt,t) * aw(wt,t) ;
EQ8a(wt,t).. WST(wt,t) =e= wsi_par(w,t,wsi_par_indx) * Q(wt,t) ;
EQ9(v,t).. STOR(v,t) =e= STOR(v,t-1)+sum (j, Q(j,v,t)*(1-lss(j,v,t))) - RR(v,t) - (evap(v,t)*ar(v,t)) ;
* modify the equation to specify the initial storage.
EQ9a(v,t).. ELEV(v,t) =e= EleV_par1*sqrt(STOR(v,t)) + EleV_par2 * STOR(v,t) + EleV_par3 ;
EQ10(j,k,t).. sum(k, Q(j,k,t) * (1-lss(j,k,t))) =g= sum (k, Q(j,k,t)) ;
EQ11(j,d,t).. sum (j, Q(j,d,t) * (1-lss(j,d,t)) -sum (j,Q(j,k,t) * (1 - Cons(d,t))) =g= sum(j, Q(d,j,t)) ;
EQ12(wt,j,t).. WA(wt,t) =l= sum (j,Q(j,wt,t) * (1-lss(j,w,t))) ;
EQ13(i,t).. D(i,t)=e= sf_par(i,sf_par_indx)*Q(i,t) ;
EQ14(i,t).. WD(i,t) =e= wf_par(i,wf_par_indx)*Q(i,t) ;
EQ15(i,t).. A(i,t)=e= WD(i,t)* lng(i) ;
EQ16(i,t,n).. C(i,t,n) =e= C(i,t,n) + RV(i,t,n) ;
EQ17(v,t).. STOR(v,t) =g= minstor(v) ;
EQ17a(v,t).. STOR(v,t) =l= maxstor(v) ;
EQ18(j,d,t).. sum(j, Q(j,d,t) * lss(j,d,t)) =g= dReq(d,t) ;
EQ19(j,d,t).. Q(j,d,t)$linkexist(j,k) =l= dcap(j,d,t) ;
EQ20(i,t) Q(i,t) =l= dmin(i,t) ;
EQ21(s,i,t).. sum((i,t,n), cst(f) * RV(i,t,n)) =l= b ;
Model WASH /all/;
Solve WASH using maximizing Z using DNLP;
Display Z.l, Z.m;