Thanks Michael. It’s working now!
On Thursday, December 3, 2015 at 10:49:54 PM UTC-7, Michael Bussieck wrote:
The pwplib is very inflexible when it comes to file names and symbol name. Your files with the piecewise data needs to be called pwp.gdx (not pwp-Nov14.gdx) and the table in the GDX file needs to be called pwpdata. You were picking up some pwp.gdx that was sitting there from some old run. So if you change the unload lines to
- Write pwp data to gdx file read by external library
$gdxout pwp.gdx
$unload pwpdata2=pwpdata
$gdxout
Your model runs fine.
Michael Bussieck - GAMS World coordinator
On Thursday, December 3, 2015 at 5:17:30 PM UTC-5, aaf…@aggiemail.usu.edu wrote:
Here’s my code (also attached):
$Title Watershed Area of Suitable Habitat (WASH) model and an application on the Lower Bear River Watershed
*Declare sets
SETS s sub-indicators
j river netowrk nodes
Alias (j,k);
SETS wt(j) impounded wetlands
dem(j) demand sites
v(j) reservoirs
;
SETS
y type of fish species /trout/
n type of vegetation species /cottonwood/
t timesteps in months /t1*t12/
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/
RA_par_indx indexes for the reservoir volume elevation curve /RA_par1, RA_par2/
sf_par_indx index for stage-flow relationship parameters /sf_par1, sf_par2/
wf_par_indx index for width-flow relationship parameters /wf_par1, wf_par2/
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
returnFlowExist(j,k) define if a return flow exists on a link (1=yes and 0=no)
DiversionExist(dem,j) define if a diversion exists on a link (1=yes and 0=no)
WetlandsExist(wt,j) define if a diversion exists on a link (1=yes and 0=no)
wght(s,j,k,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]
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) diversion link capacity [m3 per month]
dmin(j,k,t) minimum instream flow requirement [m3 per month]
cst(y) unit cost of implementing management objective
vegD(j,k,t,n) Distance between floodplain vegetation and river bank
RA_par(v, RA_par_indx) Reservoir area volume curve parameters
R_par(y,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
RSIequationNum(t) /t1 1,t2 1,t3 2,t4 2,t5 2,t6 2,t7 2,t8 2,t9 3,t10 3,t11 3,t12 3/
;
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 pwpdata 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 pwpdata2(,,) ‘1st index: function number, 2nd index: segment number, 3rd index: degree’
leftBound 0
1.1 0.001 0.001
1.2 0.25 0.25
1.3 0.4 0.5
1.4 0.5 1
2.1 0.001 0.001
2.2 0.25 0.25
2.3 0.5 0.5
2.4 1.0 1
3.1 0.001 0.001
3.2 0.25 0.75
3.3 0.4 1.25
3.4 0.5 1.5
4.1 0.001 0.2
4.2 150 1.0
5.1 -1.0 0.1
5.2 0.0 0.2
5.3 0.2 0.8
5.4 0.8 0.2
;
\
-
Write pwp data to gdx file read by external library
$gdxout pwp-Nov14.gdx
$unload pwpdata2
$gdxout
-
Load extrinsic function
$FuncLibIn pwplib pwpcclib
Function pwp /pwplib.pwpfunc/;
\
-
Read sets and parameter input values from Excel
$CALL GDXXRW.EXE input=WASH-SampleData.xlsx output=WASH-Nov14.gdx Set=s rng=SubInd!A1:A100 Rdim=1 Set=j rng=Nodes!A1:A100 Rdim=1 Set=wt rng=Wetlands!A1:A100 Rdim=1 Set=dem rng=Demand!A1:A100 Rdim=1 Set=v rng=Reservoirs!A1:A100 Rdim=1 par=envSiteExist rng=EnvSite!A1 Rdim=1 Cdim=1 par=DiversionExist rng=Diversions!A1 Rdim=1 Cdim=1 par=returnFlowExist rng=ReturnFlow!A1 Rdim=1 Cdim=1 par=WetlandsExist rng=WetlandsSites!A1 Rdim=1 Cdim=1 Set=f rng=FishSpp!A1:A10 Rdim=1 Set=n rng=VegSpp!A1:A10 Rdim=1 Set=t rng=Month!A1:A100 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!A1 Rdim=2 par=fmin rng=Fmin!A1 Rdim=2 Cdim=1 par=fmax rng=Fmax!A1 Rdim=2 Cdim=1 par=lss rng=lss!A1 Rdim=2 Cdim=1 par=evap rng=evap!A1 Rdim=1 Cdim=1 par=RA_par rng=ResElevVol!A1 Rdim=1 Cdim=1 par=Cons rng=Cons!A1 Rdim=1 Cdim=1 par=lng rng=Length!A1 Rdim=1 Cdim=1 par=minstor rng=inactive!A1 Rdim=1 par=maxstor rng=capacity!A1 Rdim=1 par=dReq rng=demandReq!A1 Rdim=1 Cdim=1 par=dCap rng=divCap!A1 Rdim=1 Cdim=1 par=dmin rng=Instream!A1 Rdim=2 Cdim=1 par=cst rng=UnitCost!A1 Rdim=1 par=VegD rng=VegD!A1 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!A1 Rdim=2 Cdim=1 par=wf_par rng=WidthFlow!A1 Rdim=2 Cdim=1 par=wsi_par rng=wp!A1 Rdim=2 Cdim=1 par=b rng=Budget!A1 Rdim=0
*Load GDX file and values
$GDXIN WASH-Nov14.gdx
$LOAD s
$LOAD j
$LOAD wt
$LOAD dem
$LOAD v
*$LOAD t
$LOAD envSiteExist
$LOAD linkexist
$LOAD DiversionExist
$LOAD returnFlowExist
$LOAD WetlandsExist
$LOAD bankdepth
$LOAD wght
$LOAD theta
$LOAD aw
$LOAD fmin
$LOAD fmax
$LOAD lss
$LOAD evap
$LOAD RA_par
$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,j,k,t) sub-indicator
Q(j,k,t) flow in link at environmental sites [m3 per month]
D(j,k,t) water depth[m]
A(j,k,t) Channel surface area [m2]
WD(j,k,t) Channel width [m]
RV(j,k,t,n) Area to re-vegetate [m2]
C(j,k,t,n) vegetation cover [m2]
O(j,k,t) overbank water level [m]
E(j,k,t) distance of flood [m]
Rec(j,k,t,n) recession 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]
RA(v,t) Reservoir surface area as a func. of storage [m2]
h(j,k,t,n) index of flood distance [0.2-1]
m(j,k,t,n) index of inundation depth [0.2 - 1]
g(j,k,t,n) index of stream recession rate [0.2 - 1]
WSI(wt,t) Wetlands suitability index
RSI(j,k,t,y) Rivering Suitability Index
FCI(j,k,t,n) Floodplain Suitability Index
R(j,k,t) Riverine habitat sub-indicator
F(j,k,t) Floodplain connectivity sub-indicator
W(j,k,t) Wetlands habitat sub-indicator
;
*Define initial values for the decision varaibles
D.L(j,k,t) = 1.60; D.UP(j,k,t) = 10; D.LO(j,k,t) = 0;
Q.L(j,k,t) = 4734.61;
A.L(j,k,t) = 1172000 ;
WD.L(j,k,t) = 10 ;
RV.L(j,k,t,n) = 1000;
C.L(j,k,t,n) =1500;
O.L(j,k,t) = 1.2 ;
E.L(j,k,t) = 10; E.UP(j,k,t) = 250; E.LO(j,k,t) = 0;
Rec.L(j,k,t,n) = 0.08 ;
WA.L(wt,t) = 7000;
RR.L(v,t) = 1000;
STOR.L(v,t) = 800529 ;
RA.L(v,t) = 1405;
h.L(j,k,t,n)=0.4; h.LO(j,k,t,n)=0.2; h.UP(j,k,t,n) =1;
m.L(j,k,t,n)=0.3; m.LO(j,k,t,n)=0.2; m.UP(j,k,t,n) =1;
*Define equations
EQUATIONS
EQ1 Z Objective function: Watershed Area of Suitable Habitat [m2]
EQ1a(s,j,k,t) ind Summation of the three perfomance indicators [m2]
EQ2(j,k,t) R Riverine Habitat [m2]
EQ2a(j,k,t,y) RSI Riversine Suitability Index [0-1]
EQ3(j,k,t) F Floodplain habitat [m2]
EQ3a(j,k,t,n) FCI floodplains suiability index [0-1]
EQ3b(j,k,t,n) h binary index of flood distance
EQ4a(j,k,t) O overbank level [m]
EQ4(j,k,t) E flood distance [m]
EQ5(j,k,t,n) m index for floodplain depth [0.2-1]
EQ6(j,k,t,n) g index for flood recession [0.2- 1]
EQ7(j,k,t,n) Rec recession rate value for g index
EQ7a(j,wt,t) WSI wetlands suitability index
EQ8(wt,j,t) Impounded Wetlands [m2]
EQ9(v,t) mass balance at reservoirs
EQ9a(v,t) reservoir area storage relationship
EQ10(j,t) mass balance at each node
EQ11(j,dem,t) mass balance at demand sites
EQ12(wt,j,t) mass balance at wetlands sites
EQ13(j,k,t) stage flow relationship
EQ14(j,k,t) width flow relationship
EQ15(j,k,t) Channel surface area [m2]
EQ16(j,k,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,dem,t) Diversions to demand sites should meet requirements
EQ19(j,k,dem,t) Diversions to demand sites should not exceed diversion capacity
EQ20(j, k,t) minimum instream flow requirements
EQ21(s,y) Limitation on Budget
;
EQ1… Z =e= sum((s,j,k,t), wght(s,j,k,t) * Ind(s,j,k,t)) ;
EQ1a(s,j,k,t)… Ind(s,j,k,t) =e= R(j,k,t)(ord(s) eq 1) + F(j,k,t)(ord(s) eq 2) + W(j,k,t)$(ord(s) eq 3);
EQ2(j,k,t)$envSiteExist(j,k)… R(j,k,t) =e= prod((y), RSI(j,k,t,y) * A(j,k,t)) ;
EQ2a(j,k,t,y)$envSiteExist(j,k)… RSI(j,k,t,y) =e= PWP(RSIequationNum(t),D(j,k,t)) ;
EQ3(j,k,t)$envSiteExist(j,k)… F(j,k,t) =e= sum((n),FCI(j,k,t,n) * C(j,k,t,n)) ;
EQ3a(j,k,t,n)$envSiteExist(j,k)… FCI(j,k,t,n) =e= h(j,k,t,n) * m(j,k,t,n) * g(j,k,t,n) ;
EQ3b(j,k,t,n)$envSiteExist(j,k)… h(j,k,t,n) =e= PWP(4,E(j,k,t)) ;
EQ4(j,k,t)$envSiteExist(j,k)… E(j,k,t)=e= O(j,k,t)/theta(j,k) ;
EQ4a(j,k,t)$envSiteExist(j,k)… O(j,k,t) =e= D(j,k,t) - bankdepth(j,k,t) ;
EQ5(j,k,t,n)$envSiteExist(j,k)… m(j,k,t,n) =e= PWP(5,O(j,k,t)) ;
EQ6(j,k,t,n)$envSiteExist(j,k)… g(j,k,t,n)=e= 0.94* (exp(-0.5*(log(Rec(j,k,t,n)/1.28))/0.99)**2) ;
EQ7(j,k,t,n)$envSiteExist(j,k)…
Rec(j,k,t,n) =e= Abs(100*(log(D(j,k,t) - log(D(j,k,t))))/30) ;
EQ8(wt,j,t)$WetlandsExist(wt,j)… W(j,wt,t) =e= WSI(wt,t) * aw(wt,t) ;
EQ7a(j,wt,t)$WetlandsExist(wt,j)… WSI(wt,t) =e= wsi_par(wt,t,“wsi_par1”) * Q(j,wt,t) +wsi_par(wt,t,“wsi_par2”) ;
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)*RA(v,t)) ;
EQ9a(v,t)… RA(v,t) =e= RA_par(v, “RA_par1”)*STOR(v,t) + RA_par(v, “RA_par2”) ;
EQ10(j,t)… sum(k, (Q(j,k,t)$linkexist(j,k) * (1-lss(j,k,t)))) =g= sum (k, Q(j,k,t)$linkexist(j,k)) ;
EQ11(j,dem,t)$DiversionExist(dem,j)… sum (k, Q(j,dem,t) * (1-lss(j,dem,t))) - sum (k,Q(j,k,t) * (1 - Cons(dem,t))) =g= sum(k, Q(dem,k,t)) ;
EQ12(wt,j,t)$WetlandsExist(wt,j)… WA(wt,t) =l= sum (k,Q(k,wt,t) * (1-lss(j,wt,t))) ;
EQ13(j,k,t)$envSiteExist(j,k)… D(j,k,t)=e= sf_par(j,k,“sf_par1”)*Q(j,k,t) + sf_par(j,k,“sf_par2”) ;
EQ14(j,k,t)$envSiteExist(j,k)… WD(j,k,t) =e= wf_par(j,k,“wf_par1”)Q(j,k,t) +wf_par(j,k,“wf_par2”) ;
EQ15(j,k,t)$envSiteExist(j,k)… A(j,k,t)=e= WD(j,k,t) lng(j,k) ;
EQ16(j,k,t,n)$envSiteExist(j,k)… C(j,k,t,n) =e= C(j,k,t,n) + RV(j,k,t,n) ;
EQ17(v,t)… STOR(v,t) =g= minstor(v) ;
EQ17a(v,t)… STOR(v,t) =l= maxstor(v) ;
EQ18(j,dem,t)… sum(k,Q(k,dem,t)$DiversionExist(dem,j) * lss(k,dem,t)) =g= dReq(dem,t) ;
EQ19(j,k,dem,t)… Q(j,dem,t)$linkexist(j,k) =l= dcap(j,dem) ;
EQ20(j,k,t)$linkexist(j,k)… Q(j,k,t) =l= dmin(j,k,t) ;
EQ21(s,y)… sum((j,k,t,n), cst(y) * RV(j,k,t,n)) =l= b ;
option limrow = 10000 ;
Model WASH /all/;
Solve WASH maximizing Z using DNLP;
On Thursday, December 3, 2015 at 1:02:27 AM UTC-7, Michael Bussieck wrote:
Your table looks fine. You need to send your entire model to see where the real problem is.
Michael Bussieck - GAMSWorld Coordinator
On Thursday, December 3, 2015 at 2:29:28 AM UTC-5, aaf…@aggiemail.usu.edu wrote:
Hey,
I’m trying to use the Piecewise polynomial library following to the user guide and the pwplib01 test library. I have 5 piecewise functions that the model is reading from the following table entry:
Table pwpdata2(,,*) ‘1st index: function number, 2nd index: segment number, 3rd index: degree’
leftBound 0
1.1 0.001 0.001
1.2 0.25 0.25
1.3 0.4 0.5
1.4 0.5 1
2.1 0.001 0.001
2.2 0.25 0.25
2.3 0.5 0.5
2.4 1.0 1
3.1 0.001 0.001
3.2 0.25 0.75
3.3 0.4 1.25
3.4 0.5 1.5
4.1 0.001 0.2
4.2 150 1.0
5.1 -1.0 0.1
5.2 0.0 0.2
5.3 0.2 0.8
5.4 0.8 0.2
;
However, when I run it, the code can only read functions 1 and 2 but it cannot read 3, 4 or 5 and I’m getting the following error:
*** Error at line 350: unkown function 3. Available 1 to 2.
and the same for functions 4 and 5.
Any reason why the code cannot read those functions from the same table?
Thanks,
Ayman
–
To unsubscribe from this group and stop receiving emails from it, send an email to gamsworld+unsubscribe@googlegroups.com.
To post to this group, send email to gamsworld@googlegroups.com.
Visit this group at http://groups.google.com/group/gamsworld.
For more options, visit https://groups.google.com/d/optout.