$EOLCOM // $TITLE RIO GRANDE BASIN HYDROECONOMIC PROTOTYPE $OFFSYMXREF OFFSYMLIST OFFLISTING OFFUPPER OPTION LIMROW=0, LIMCOL = 0; $ONTEXT Endogenous capacity expansion Allocates scarce infrastructure capacity expansions among competing reservoirs 2 reservoirs x 2 agricultural use nodes quadratic total benefits function at each node reservoir capacity expansion is allocated to minimize evaporation for other things equal Prototype example is Rio Grande Basin of Colorado, New Mexico, and Texas (USA) and Mexico * --------------------------------------------------------------------------------------- Frank A. Ward: Dept of Agr Economics/Agr Business New Mexico State University, Las Cruces, NM USA e-mail: fward@nmsu.edu February 17 2008 Rio Grande Basin Model: Expandable Prototype Contains essential elements of full Upper Rio Grande Basin Model. * --------------------------------------------------------------------------------------- Sponsored by Water Resources Research Institutes and Agricultural Experiment Stations of Colorado, New Mexico,and Texas * --------------------------------------------------------------------------------------- Model has these flow nodes: 6 river gauge nodes 2 inflow nodes 2 diversion nodes 2 consumptive use nodes 2 surface water return flow nodes 2 reservoir release node and these stock nodes: 2 reservoir nodes * --------------------------------------------------------------------------------------- FLOWS: Spatial unit for FLOWS is set (index) i. Each element in the set i is assigned to one water use subset (category) Subset categories include: 1. Inflow nodes to the system, inflow(i); 2. Nodes on a river or tributary river(i); 3 Diversion nodes divert(i); 4. Consumptive uses use(i); 5. Return flow nodes directly to the river, return(i); 6. NET reservoir releases from storage, outflow - inflow rel(i); STOCKS: Spatial unit for STOCKS is the set index u. Each element of the set u is assigned to one water use subset (category). Subset categories are: 1. Reservoir nodes, res(u). * --------------------------------------------------------------------------------------- TABLE OF CONTENTS Section 1. Sets Section 2. Data Section 3. Variables Section 4. Equations Section 5. Models Section 6. Solves Section 7. Displays * --------------------------------------------------------------------------------------- **************** Section 1 ************************************************************** * The following sets are specified as indices * * for parameters (data), variables, and equations * ***************************************************************************************** $OFFTEXT SETS ***************************************************************************************** i Flows -- location of important nodes in RG Basin -- CO to MX ***************************************************************************************** / RG_h_f Headwater flow nodes inflow(i) Chama_h_f Lobatos_v_f River gage measurement nodes river(i) Embudo_v_f Chamita_v_f Otowi_v_f Acacia_v_f MX_v_f EBID_d_f Diversion nodes divert(i) TX_d_f EBID_u_f Consumptive use flow nodes use(i) TX_u_f EBID_r_f Surface water return flow nodes return(i) TX_r_f CO_rel_f Reservoir-to-river release flow nodes rel(i) EB_rel_f CO_evp_f Cochiti evap(i) EB_evp_f Elephant Butte / ***************************************************************************************** * Subsets of all Flow nodes above by class (function) ***************************************************************************************** inflow(i) Headwater flow nodes inflow(i) / RG_h_f Rio Grande headwaters CO Chama_h_f Rio Chama headwaters near CO-NM state line / river(i) River gage measurement nodes river(i) / Lobatos_v_f Lobatos gauge on RG at CO-NM state line Embudo_v_f Embudo gauge on RG northern NM Chamita_v_f Chamita gauge on Rio Chama northern NM Otowi_v_f Otowi gauge on RG downstream of Chama RG confluence Acacia_v_f San Acacia gauge central NM downstream of MRGCD ag MX_v_f RG at Mexico near El Paso TX / divert(i) Diversion nodes divert(i) / EBID_d_f Elephant Butte Irr Dist ag near Las Cruces NM TX_d_f West TX irrigated Ag near El Paso TX / use(i) Consumptive use flow nodes = div nodes use(i) / EBID_u_f same as divert(i) TX_u_f / return(i) Surf wat return flow nodes = div nodes return(i) / EBID_r_f same as divert(i) TX_r_f / rel(i) Reservoir to river release flow nodes rel(i) / CO_rel_f Cochiti Reservoir EB_rel_f Elephant Butte reservoir releases to RG / evap(i) Reservoir evaporation flow nodes evap(i) / CO_evp_f Cochiti EB_evp_f Elephant Butte / ***************************************************************************************** u Stocks--location of important nodes on Rio Grande CO to MX ***************************************************************************************** / CO_res_s Cochiti reservoir stock node res(u) EB_res_s Elephant Butte stock node / ***************************************************************************************** * Stock subsets ***************************************************************************************** res(u) Reservoir stock nodes res(u) / CO_res_s Cochiti reservoir storage volume EB_res_s Elephant Butte reservoir storage volume / ***************************************************************************************** t time ***************************************************************************************** / 1*10 # of time periods / tlast(t) terminal period among all periods above tearly(t) early periods tlate(t) late periods ; tlast(t) = yes $ (ord(t) eq card(t)); // GAMS language -- picks last pd tearly(t) = yes $ (ord(t) le card(t)/2); // early half years tlate(t) = yes $ (ord(t) gt card(t)/2); // late half years ALIAS (i,ip); ALIAS (river, riverp); ALIAS (divert, divertp); * lets some tables' nodes be rows or columns **************** Section 2 ************************************************************** * This section defines all data in 3 formats * * 1. Scalars (single numbers), * * 2. Parameters (columns of numbers) or * * 3. Tables (data in rows and columns) * ***************************************************************************************** * Below are several maps summarizing a basin's geometry * By geometry we mean location of mainstems, tributaries, confluence, * source nodes, use nodes, return flow nodes, reservoir nodes, etc. * Basin geometry is summarized through judicious use of numbers 1, -1, and 0 (blank) ***************************************************************************************** * Map #1: * Each column below is a streamgage. Each row is a source or use of water. * Flow at ea gage (column) is directly influenced by at least 1 upstream row. * SOURCE adds to columns flow (+1) * USE deplete from col flow (-1) * BLANK has no effect on col flow ( ) * Geometry accounts for all sources (supplies) and uses (demands) in basin * Map is used to produce coefficients in equations below to define: * X(river) = Bhv * X(inflow) + Bvv * X(river) + Bdv * X(divert) * + Brv * X(return) + Bgv * X(gwflows) + BLv * X(rel) * * These B coeff vectors are stacked below as a single matrix, Bv ***************************************************************************************** TABLE Bv(i,river) Hydrologic Balance Table **************************** Column Heads are River Gauges ************************* Lobatos_v_f Embudo_v_f Chamita_v_f Otowi_v_f Acacia_v_f MX_v_f * ---------------- headwater inflow node rows (+) --------------------------------------- RG_h_f 1 Chama_h_f 1 * ---------------- river gage node rows (+) --------------------------------------------- Lobatos_v_f 1 Embudo_v_f 1 Chamita_v_f 1 Otowi_v_f 1 Acacia_v_f 1 MX_v_f * --------------- diversion nodes (-) ------------------------------------------------- EBID_d_f -1 TX_d_f -1 * -------------- return flow node rows (+) ---------------------------------------------- EBID_r_f 1 TX_r_f 1 * ------------- reservoir release (outflow) to river -- stock-to-flow rows (+) ---------- CO_rel_f 1 EB_rel_f 1 ; ***************************************************************************************** * Map #2: * Enforces nonnegative flows at each use node (wet river) * water sources are rows. Diversion nodes are columns. * For any column, diversion < summed flows from upstream sources (rows) * e.g. SLV Colorado ag use < flows from RG and Conejos headwater sources * X(divert) < Bhd * X(inflow) + Brd * X(river) + Bdd * X(divert) + * Brd * X(return) + Bgd * X(gwflow) + BLd * X(rel) * These B coeff vectors are stacked below as the matrix, Bd ***************************************************************************************** TABLE Bd(i, divert) Wet river table * ------- Col Heads are Diversion nodes ----------------------------------------------- EBID_d_f TX_d_f * ----------------------- headwater inflow nodes (+)------------------------------------- RG_h_f Chama_h_f * ----------------------- river gage nodes ---------------------------------------------- Lobatos_v_f Embudo_v_f Chamita_v_f Otowi_v_f Acacia_v_f 1 1 MX_v_f * ------------------------ diversion nodes (-) ---------------------------------------- EBID_d_f -1 TX_d_f * ---------------------- return flow nodes (+) ------------------------------------------ EBID_r_f 1 TX_r_f * ---------------- reservoir outflow stock-to-flow node row (+)------------------------- CO_rel_f EB_rel_f 1 1 * --------------------------------------------------------------------------------------- ***************************************************************************************** * Map #3: * Defines use (simplistically) as a percentage of diversion * X(use) = Bdu * X(divert) * These B coeffs are shown below as the matrix, Bu ***************************************************************************************** TABLE Bu(divert, use) Defines consumptive use as a percent of diversion * ------------------------- Use nodes -------------------------------------------------- EBID_u_f TX_u_f * ------------------------ divert nodes (+) --------------------------------------------- EBID_d_f 0.5 TX_d_f 0.5 *---------------------------------------------------------------------------------------- ; ***************************************************************************************** * Map #4: * Table defines relation between diversions and return flow nodes * * Tabled entries = proportion return flow by diversion column nodes * (+) means the row diversion contributes to column's ret flow * ( ) means the column diversion makes no cont to row's ret flow *X(return) = Br * X(divert) ***************************************************************************************** TABLE Br(divert, return) Defines surface return flow as a percent of diversion ********************* Column Heads are Return Flow Nodes ******************************* EBID_r_f TX_r_f EBID_d_f 0.5 TX_d_f 0.5 ; ***************************************************************************************** * Map #5: * Table relates reservoir stocks in a period to its prev periods' stocks minus releases. * For any reservoir stock node at the column head * (+1) :added water at flow node -- thru releases -- takes from column's res stock (-) * (-1) :added water at flow node adds to column's reservoir stock * ( ) :added water at flow node has no effect on column's reservoir stock * Z(res(t)) = Z(res(t-1)) + BLv * X(rel(t)) ***************************************************************************************** TABLE BLv(rel, u) Links reservoir releases to downstream flows ********* Column Heads are Reservoir Stocks -- rows are release flows ****************** ********* Table = diagonal matrix for > 1 reservoir--only 1 for now ****************** CO_res_s EB_res_s CO_rel_f 1 EB_rel_f 1 ; * Map #6 TABLE Ber(evap, res) Links reservoir evaporation to volume loss ********** Column Heads are reservoir stocks -- rows are evaporation loss flows ******** ********** Table = diagonal matrix for > 1 reservoir *********************************** CO_res_s EB_res_s CO_evp_f 1 EB_evp_f 1 *---------------------------------------------------------------------------------------- * Map #7 * Defines pan evaporation in inches - translate to af lost per per exposed acre per year by reservoir TABLE Bev(evap, res) CO_res_s EB_res_s CO_evp_f 93.0 EB_evp_f 111.2 ; Bev(evap,res) = Bev(evap,res)/12; // inches to feet ***************************************************************************************** * END OF BASIN GEOMETRY MAPS * ***************************************************************************************** ***************************************************************************************** * NEXT APPEAR BASIN INFLOWS, OTHER FLOWS, FLOW RELATIONSHIPS, AND * * RESERVOIR STARTING VOLUMES, SIMPLE ECONOMIC VALUES PER AC FT WATER USE * ***************************************************************************************** * all water flows are measured in 1000s acre feet per yer * all water stocks are measured in 1000s acre feet instantaneous volume TABLE sourc(inflow,t) annual basin inflows at headwaters -- snowpack or rain ***** Data are from historical or forecast headwater node flows ****************** 1 RG_h_f 500 Chama_h_f 250 ; parameter source(inflow,t) annual basin inflows over all periods ; source(inflow, tearly) = 4.00 * sourc(inflow,'1'); // high flows in early years source(inflow, tlate ) = 0.10 * sourc(inflow,'1'); // severe drought in later years PARAMETERS z0(res) initial reservoir levels at stock nodes ZMAX(res) maximum reservoir capacity ; z0('CO_res_s' ) = 75; // Cochiti Reservoir starting volume z0('EB_res_s' ) = 400; // Elephant Butte Reservoir starting vol *co 75 *eb 400 ZMAX('CO_res_s') = 719; // Cochiti Reservoir max current capacity ZMAX('EB_res_s') = 2065; // Elephant Butte Reservoir max current capacity ; *co 719 *eb 2065 PARAMETERS B0ar(res) Area - Capacity intercept: Intcpt for reservoir area as linear fn of volume = 0 B1ar(res) Area - Capacity slope : (1st order) Slope for res area = linear fn of vol = d(area)\d(vol) ; B0ar('CO_res_s') = 0; //Cochiti intercept max vol = 0.719 kk af B1ar('CO_res_s') = .014798; //Cochiti slope max area = 10.636 k ac B0ar('EB_res_s') = 0; //Elephant Butte inter max vol = 2.065 kk af B1ar('EB_res_s') = .01743; //Elephant Butte slope max area = 36.000 k ac *0 *.014798 *0 *.01743 parameters BC0(res) Total construction cost intercept: fixed cost of expanding reservoir capacity BC1(res) Total construction cost linear: constant marginal cost of expanding reservoir capacity BC2(res) Total construction cost quadratic: rising marginal cost of expanding reservoir capacity ; BC0('CO_res_s') = 0; // Cochiti BC1('CO_res_s') = 1.00; BC2('CO_res_s') = 0.05; BC0('EB_res_s') = 0; // Elephant Butte BC1('EB_res_s') = 1.00; BC2('EB_res_s') = 0.05; ***************************************************************************************** * ECONOMICS MAPS BEGIN * ***************************************************************************************** parameter scales(use) 1000s of acres by ag node Data source USDA 2002 Ag Census / EBID_u_f 100.000 TX_u_f 50.000 / parameter growth(use) annual forecast growth rate of acreage or population by node / EBID_u_f 0.0 TX_u_f 0.0 / parameter scale(use,t) 1000s forecast acres or households by node and time ; scale(use,t) = ((1 + growth(use ))**(ord(t)-1)) * scales(use); TABLE Ben_u_p(use, *) Per acre or per hhold benefit fn params for all uses * ($) ($/ac-af) ($/ac-af2) (af/ac-yr1) (af-yr1) intercept linear quadratic u_mb_0_a u_mb_0 EBID_u_f 0.00 200.00 -20.00 5.00 500 TX_u_f 0.00 200.00 -20.00 5.00 250 ; *----------------------------------------------------------------------------------- **************** Section 3 ************************************************************** * These endogenous (unknown) variables are defined * * Their numerical values are not known til GAMS finds optimal soln * ***************************************************************************************** FREE VARIABLES X(i,t) water flows -- diversion-use-return - etc. Ben_v(use, t) total benefits by use and time Tben_v total benefits over uses and periods dvol_v(res,t) change in reservoir volume TNben_v Total net benefits over reservoirs net of cost of expanding capacity ; POSITIVE VARIABLES Z (u, t) Stocks: reservoir volumes Za(res,t) Stocks: reservoir areas DZ(res) Stock: added reservoir capacity Zadd_cap Stock: added reservoir capacity for region X_capac(res) Stock: defines capacity Z_gain_v Stock: regional capacity gain Cost_v (res) Total cost of expanding a reservoirs capacity Tcost_v Total cost over reservoirs of expanding capacity ; **************** Section 4 ************************************************************** * The following equations state relationships among a basin's * * hydrology, institutions, and economics * ***************************************************************************************** EQUATIONS ***************************************************************************************** * Equations named ***************************************************************************************** * Hydrology Block Inflows (inflow,t) Flows: set source nodes Rivers (i, t) Flows: hydrologic mass balance for each flow node: sources = uses Divs (divert,t) Flows: set divert nodes Returns (return,t) Flows: set return flows Uses (use, t) Flows: define use = diversions - return flows Evaps (evap, t) flows: defies evaporation reservoirs(res, t) Stock: reservoir mass balance accounting area (res, t) Stock: reservoir area dvol_e (res, t) stock: change in reservoir volume def_capacity(res,t) Stock: defines capacity bound_storage(res,t) Stock: bounds storage Z_gain_e Stock: regional storage addition * Institutions Block (empty) * Economics Block Ben_e(use,t) Total benefits by use and time Tben_e Total benefits over use and time Cost_e(res) Total expansion cost of expanding a reservoirs capacity Tcost_e Total cost over reservoirs of expanding TNben_e Total net benefits over reservoirs including cost of expanding capacity ; ***************************************************************************************** * Equations defined algebraiclly using equation names ***************************************************************************************** * Hydrology Block Inflows(inflow,t).. X(inflow,t) =E= source(inflow,t); Rivers(river,t).. X(river,t) =E= sum(inflow, Bv(inflow, river) * X(inflow, t)) + sum(riverp, Bv(riverp, river) * X(riverp, t)) + sum(divert, Bv(divert, river) * X(divert, t)) + sum(return, Bv(return, river) * X(return, t)) + sum(rel, Bv(rel, river) * X(rel, t)) ; Divs(divert,t).. X(divert,t) =L= sum(inflow, Bd(inflow, divert) * X(inflow, t)) + sum(river, Bd(river, divert) * X(river, t)) + sum(divertp, Bd(divertp,divert) * X(divertp,t)) + sum(return, Bd(return, divert) * X(return, t)) + sum(rel, Bd(rel, divert) * X(rel, t)) ; Evaps(evap,t).. X(evap,t) =E= sum(res, Bev(evap, res) * Za(res, t)) ; Uses (use, t).. X(use, t) =E= sum(divert, Bu(divert, use ) * X(divert, t)) ; Returns(return,t).. X(return,t) =E= sum(divert, Br(divert, return) * X(divert, t)) ; reservoirs(res,t).. Z(res,t) =E= z0(res)$(ORD(t) EQ 1) + Z(res,t-1) - SUM(rel, BLv(rel, res) * X(rel, t)) - sum(evap, Ber(evap, res) * X(evap, t)); area(res,t).. Za(res,t) =E= B0ar(res) + B1ar(res) * Z(res, t); dvol_e(res,t).. dvol_v(res,t) =e= Z(res,t) $ (ord(t) gt 1) - z(res,t-1); * this part of hydrology block treats capacity additions as endogenous -- looking for best reservoir to expand def_capacity(res,t).. X_capac(res) =e= ZMAX(res) + DZ(res) ; // capacity = original zmax + endogenous additions bound_storage(res,t).. Z(res,t) =l= X_capac(res) ; // a reservoir cannot exceed its endogenous capacity Z_gain_e.. Z_gain_v =e= sum(res, DZ(res)); ; // regional gain in capacity = total additions (dz) *DZ_equal(res).. DZ('CO_res_s') =e= 0; // both expansions are equal *Z_equal(res,t).. Z('CO_res_s',t) =e= z('EB_res_s',t); // reservoirs are equal *Z_eq_1 (res,t).. Z(res,'1') =e= z(res,'2'); *z_eq_2 (res,t).. z(res,'2') =e= z(res,'3'); *z_eq_3 (res,t).. z(res,'3') =e= z(res,'4'); Cost_e(res).. Cost_v(res) =e= BC0(res) + BC1(res) * DZ(res) + BC2(res) * DZ(res) * DZ(res); * Institutions Block --water laws, compacts, treaties, etc constrains use (Empty for now) * Economics Block Ben_e(use,t).. Ben_v(use, t) =E= scale(use,t) * Ben_u_p(use, 'intercept') + Ben_u_p(use, 'linear') * X(use,t) + (1/scale(use,t)) * Ben_u_p(use, 'quadratic') * X(use,t) * X(use,t) ; Tben_e.. Tben_v =e= sum((t,use), Ben_v(use,t)); Tcost_e.. Tcost_v =e= sum((t, res), Cost_v(res)); TNben_e.. TNben_v =e= Tben_v - Tcost_v ; // total net benefits includes cost of capacity expansion *************************** End of equations ******************************************* **************** Section 5 ************************************************************** * The following section defines models. * * Each model is defined by a set of equations used * * for which one single variable is optimized (min or max) * ***************************************************************************************** * This simple prototype model uses ALL equations defined above. But larger models * may exclude some equations. For example, each of several institution could be defined * by one equation. And each of several model might conduct a single policy experiment * in which that model tries out a single institution. This would require deleting all * institutional equations except the one analyzed. * If you need to EXclude some equations, list INcluded equations where ALL appears below MODEL RIO_PROTOTYPE /ALL/; **************** Section 6 ************************************************************** * The following section defines all solves requested, * Each solve states a single model for which an optimum is requested. * * Upper, lower and fixed bounds on certain variables can also be included here * Bounding variables here gives that variable a non-zero shadow price where the optimal * solution appears at that boundary. If the bound doesn't constrain the model * the variable's shadow price is zero (complementary slackness) ***************************************************************************************** * Non-negative flows at nodes below X.lo(inflow,t) = 0; X.lo(river, t) = 0; X.lo(divert,t) = 0; X.lo(use, t) = 0; X.lo(return,t) = 0; * Sustainability terminal condition -- each water stock (reservoir, aquifer) ends with * terminal volume > starting volume. * It avoids depleting stocks in last period -- saves water for future generations Z.lo('CO_res_s', t) = 1.00 * z0('CO_res_s'); Z.lo('EB_res_s', t) = 1.00 * z0('EB_res_s'); dvol_v.up(res, t) = 50.0; // upper bound on changes in reservoir water volume by pd Z.lo(res, tlast) = Z0(res); // reservoir terminal period vol > starting value Z_gain_v.lo = 0; // added regional capcaity SOLVE RIO_PROTOTYPE USING NLP MAXIMIZING TNben_v; * the end