From 320026059644aee8f2cf64ded14a62974f02e7f6 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Mon, 8 Feb 2021 14:47:31 +0000
Subject: [PATCH] Initial changes to include wood yields and carbon fluxes.
---
GAMS/IntExtOpt.gms | 275 ++++++++++--------
.../lurg/country/gams/GamsLocationInput.java | 17 +-
.../ed/lurg/country/gams/GamsRasterInput.java | 18 +-
src/ac/ed/lurg/landuse/CarbonFluxItem.java | 22 ++
.../ed/lurg/landuse/CarbonFluxRasterSet.java | 28 ++
src/ac/ed/lurg/landuse/WoodYieldItem.java | 21 ++
src/ac/ed/lurg/types/NonCropType.java | 53 ++++
7 files changed, 316 insertions(+), 118 deletions(-)
create mode 100644 src/ac/ed/lurg/landuse/CarbonFluxItem.java
create mode 100644 src/ac/ed/lurg/landuse/CarbonFluxRasterSet.java
create mode 100644 src/ac/ed/lurg/landuse/WoodYieldItem.java
create mode 100644 src/ac/ed/lurg/types/NonCropType.java
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 6b558131..33495b6e 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -1,18 +1,24 @@
- SET all_types / monogastrics, ruminants, cereals, oilcropspulses, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside /;
-
- SET crop(all_types) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside /;
- SET crop_less_pasture(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, setaside/;
- SET cereal_crop(crop) / wheat, maize, rice /;
- SET oilpulse_crop(crop) / oilcrops, pulses /;
- SET feed_crop(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, pasture/;
- SET feed_crop_less_pasture(feed_crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg /;
- SET import_crop(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops /;
-
+
+ SET all_types / monogastrics, ruminants, cereals, oilcropspulses, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside, carbonForest, timberForest, natural /;
+
+ SET crop(all_types) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside /;
+ SET crop_less_pasture(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, setaside/;
+ SET cereal_crop(crop) / wheat, maize, rice /;
+ SET oilpulse_crop(crop) / oilcrops, pulses /;
+ SET feed_crop(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, pasture/;
+ SET feed_crop_less_pasture(feed_crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg /;
+ SET import_crop(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops/;
+ SET wood_producing(all_types) / timberForest, carbonForest, natural/;
+* SET forest(all_types) / timberForest, carbonForest /;
+
+ ALIAS (all_types, lu_type_previous);
+
SET location;
PARAMETER suitableLandArea(location) areas of land in Mha;
- PARAMETER previousArea(crop, location) areas for previous timestep in Mha;
+ PARAMETER previousCropArea(crop, location) areas for previous timestep in Mha;
+ PARAMETER previousNonCropArea(wood_producing, location);
PARAMETER previousFertIntensity(crop, location);
- PARAMETER previousIrrigIntensity(crop, location);
+ PARAMETER previousIrrigIntensity(crop, location);
PARAMETER previousOtherIntensity(crop, location);
PARAMETER previousRuminantFeed(crop);
PARAMETER previousMonogastricFeed(crop);
@@ -25,7 +31,7 @@
PARAMETER yieldShock(crop, location) rate of yield shock;
PARAMETER fertParam(crop, location) yield response to fertilizer parameter;
PARAMETER irrigParam(crop, location) yield response to irrigation parameter;
- PARAMETER demand(all_types) in t;
+ PARAMETER demand(all_types) in t;
PARAMETER exportPrices(import_crop) prices for exports;
PARAMETER importPrices(import_crop) prices for imports;
PARAMETER maxNetImport(import_crop) maximum net import for each crop based on world market;
@@ -37,13 +43,13 @@
PARAMETER minDemandPerOilcrop(oilpulse_crop) min demand for oilcrop pulses as factor of total;
PARAMETER seedAndWasteRate(all_types) rate of use for seed and waste combined;
PARAMETER subsidyRate(crop) rates of subsidy compared to costs;
-
+
PARAMETER agriExpansionCost(location) price for increase agricultural area varies based on amount of managed or unmanaged forest;
SCALAR pastureIncCost price for increasing pasture area;
SCALAR cropIncCost price for increasing crop;
SCALAR cropDecCost price for decreasing crop;
SCALAR pastureDecCost price for decreasing pasture;
-
+
SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
SCALAR fertiliserUnitCost fert cost at max fert rate;
SCALAR otherIParam yield response to other intensity;
@@ -51,44 +57,50 @@
SCALAR unhandledCropRate rate of fruit veg and other crops not modelled;
SCALAR setAsideRate rate of set aside and failed crop;
SCALAR domesticPriceMarkup factor price increased from cost of production;
- SCALAR maxLandExpansionRate max rate of country land expansion;
+ SCALAR maxLandExpansionRate max rate of country land expansion;
+
+ PARAMETER carbonFluxRate(all_types, location) carbon flux;
+ PARAMETER woodYield(wood_producing, location) wood yield;
+ SCALAR carbonPrice price of carbon;
+ SCALAR woodPrice price of wood and timber;
*$gdxin "/Users/peteralexander/Documents/R_Workspace/UNPLUM/temp/GamsTmp/t1.gdx"
-$gdxin %gdxincname%
+$gdxin %gdxincname%
$load location, suitableLandArea, demand, agriExpansionCost, cropIncCost, pastureIncCost, cropDecCost, pastureDecCost
-$load previousArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
+$load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
-$gdxin
-
+$load carbonFluxRate, woodYield, carbonPrice, woodPrice
+$gdxin
+
SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /;
-
+
demand(cereal_crop) = demand('cereals') * minDemandPerCereal(cereal_crop);
demand(oilpulse_crop) = demand('oilcropspulses') * minDemandPerOilcrop(oilpulse_crop);
-
- previousArea(crop_less_pasture, location) = previousArea(crop_less_pasture, location) * (1.0 - unhandledCropRate);
-
+
+ previousCropArea(crop_less_pasture, location) = previousCropArea(crop_less_pasture, location) * (1.0 - unhandledCropRate);
+
PARAMETER cropDM(crop) kg DM per kg of feed the conversion from feed to meat is done in the R animal product index. Pasture is in DM terms already
- / wheat 0.87
- maize 0.86
+ / wheat 0.87
+ maize 0.86
rice 0.89
- oilcrops 0.96
+ oilcrops 0.96
pulses 0.31
starchyRoots 0.25
fruitveg 0.1
sugar 1
- pasture 1 / ;
-
+ pasture 1 / ;
+
PARAMETER prodCost(crop) cost per ha before intensity values not including fertiliser or irrigation i.e. seed spray and other. In 1000 $ per ha
- / wheat 0.32
- maize 0.31
+ / wheat 0.32
+ maize 0.31
rice 0.36
oilcrops 0.2
pulses 0.31
- starchyRoots 3.14
+ starchyRoots 3.14
fruitveg 4.0
- sugar 3.0
+ sugar 3.0
energyCrops 0.34 / ;
PARAMETER baseCost(crop);
@@ -97,9 +109,13 @@ $gdxin
otherIntCost(crop) = baseCost(crop)*0.7 + otherICost;
baseCost('pasture') = 0.04;
otherIntCost('pasture') = 0.8 + otherICost;
-
+
+*PARAMETER lcTransitionCost(lc_type_previous, lc_type, location) cost of land transitions including carbon cost and conversion cost;
+*lcTransitionCost(lc_type_previous, lc_type, location) = carbonFluxRate(lc_type_previous, lc_type, location) * carbonPrice + agriExpansionCost(location);
+
VARIABLES
- area(crop, location) total area for each crop - Mha
+ cropArea(crop, location) total area for crops - Mha
+ nonCropArea(wood_producing, location) forest and natural area
fertI(crop, location) fertilizer intensity for each crop - factor between 0 and 1
irrigI(crop, location) irrigation intensity for each crop - factor between 0 and 1
otherIntensity(crop, location)
@@ -111,18 +127,21 @@ $gdxin
unitCost(crop, location) cost per area for each crop - cost
net_supply(crop) supply after exports and feed
agriLandExpansion(location) addition agricultural land needed as it must be positive it deliberately does not account for abandonment
- cropIncrease(location)
- cropDecrease(location)
- pastureIncrease(location)
- pastureDecrease(location)
- totalFeedDM
+ cropIncrease(location)
+ cropDecrease(location)
+ pastureIncrease(location)
+ pastureDecrease(location)
+ woodHarvest(location) total wood harvested
+ carbonFlux(location) total carbon flux
+ totalFeedDM
total_cost total cost of domestic supply including net imports;
-
- POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
- agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease, totalFeedDM;
-
+
+ POSITIVE VARIABLE cropArea, nonCropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
+ agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease, totalFeedDM,
+ woodHarvest, carbonFlux;
+
EQUATIONS
- UNIT_COST_EQ(crop, location) cost per area
+ UNIT_COST_EQ(crop, location) cost per area
YIELD_EQ(crop, location) yield given chosen intensity
CROP_DEMAND_CONSTRAINT(crop) satisfy demand for individual crops
TOTAL_CEREAL_DEMAND_CONSTRAINT satisfy demand for combined cereals
@@ -135,108 +154,134 @@ $gdxin
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity
MAX_OTHER_INTENSITY_CONSTRAINT(crop, location)
SETASIDE_AREA_CALC(location)
- TOTAL_LAND_CHANGE_CONSTRAINT(location) constraint on suitable land use
+ TOTAL_LAND_CHANGE_CONSTRAINT(location) constraint on suitable land use
MAX_NET_IMPORT_CONSTRAINT(import_crop) constraint on max net imports
MIN_NET_IMPORT_CONSTRAINT(import_crop) constraint on min net imports
PASTURE_IMPORT_CONSTRAINT constraint to not import pasture
PASTURE_EXPORT_CONSTRAINT
+ TIMBER_IMPORT_CONSTRAINT constraint to not import timber (for now)
+ TIMBER_EXPORT_CONSTRAINT
+ CARBON_IMPORT_CONSTRAINT constraint to not import carbon
+ CARBON_EXPORT_CONSTRAINT
+ NATURAL_IMPORT_CONSTRAINT constraint to not import natural
+ NATURAL_EXPORT_CONSTRAINT
IRRIGATION_CONSTRAINT(location) constraint no water usage
NET_SUPPLY_EQ(crop) calc net supply for crops
AGRI_LAND_EXPANSION_CALC(location) calc agriLandExpansion
- CROP_INCREASE_CALC(location)
- CROP_DECREASE_CALC(location)
- PASTURE_INCREASE_CONV_CALC(location)
- PASTURE_DECREASE_CONV_CALC(location)
- PASTURE_TOTAL_CHANGE_CONSTRAINT(location)
+ CROP_INCREASE_CALC(location)
+ CROP_DECREASE_CALC(location)
+ PASTURE_INCREASE_CONV_CALC(location)
+ PASTURE_DECREASE_CONV_CALC(location)
+ PASTURE_TOTAL_CHANGE_CONSTRAINT(location)
+ WOOD_HARVEST_CALC(location) calc wood harvested
+ CARBON_FLUX_CALC(location) calc carbon flux
COST_EQ total cost objective function;
-
- UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E= ( baseCost(crop) +
- fertiliserUnitCost * fertI(crop, location) +
+
+ UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E= ( baseCost(crop) +
+ fertiliserUnitCost * fertI(crop, location) +
irrigCost(location) * irrigMaxRate(crop, location) * irrigI(crop, location) +
otherIntCost(crop) * otherIntensity(crop, location)
) ;
-
+
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
- yieldNone(crop, location) +
+ yieldNone(crop, location) +
(yieldFertOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-fertParam(crop, location)*fertI(crop, location))) +
(yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location))) +
(yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
(1 - exp(-fertParam(crop, location)*fertI(crop, location))) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location)))
) * (1 - exp(-otherIntensity(crop, location)*otherIParam));
-
- NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) * (1 - seedAndWasteRate(crop)) - ruminantFeed(crop) - monogastricFeed(crop) + importAmount(crop) - exportAmount(crop);
-
+
+ NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, cropArea(crop, location) * yield(crop, location)) * (1 - seedAndWasteRate(crop)) - ruminantFeed(crop) - monogastricFeed(crop) + importAmount(crop) - exportAmount(crop);
+
CROP_DEMAND_CONSTRAINT(crop) .. net_supply(crop) =G= demand(crop);
-
+
TOTAL_CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, net_supply(cereal_crop)) =G= demand('cereals');
-
+
TOTAL_OIL_PULSE_DEMAND_CONSTRAINT .. sum(oilpulse_crop, net_supply(oilpulse_crop)) =G= demand('oilcropspulses');
-
+
RUMINANT_DEMAND_CONSTRAINT .. meatEfficency * sum(feed_crop, ruminantFeed(feed_crop) * cropDM(feed_crop)) * (1 - seedAndWasteRate('ruminants')) =G= (demand('ruminants') - importAmount('ruminants') + exportAmount('ruminants'));
-
+
MONOGASTRICS_DEMAND_CONSTRAINT .. meatEfficency * sum(feed_crop_less_pasture, monogastricFeed(feed_crop_less_pasture) * cropDM(feed_crop_less_pasture)) * (1 - seedAndWasteRate('monogastrics')) =G= (demand('monogastrics') - importAmount('monogastrics') + exportAmount('monogastrics'));
-
+
TOTAL_NON_PASTURE_FEED_DM_CALC .. totalFeedDM =E= sum(feed_crop_less_pasture, (ruminantFeed(feed_crop_less_pasture) + monogastricFeed(feed_crop_less_pasture)) * cropDM(feed_crop_less_pasture));
FEED_MIX_CONSTRAINT(feed_crop_less_pasture) .. (ruminantFeed(feed_crop_less_pasture) + monogastricFeed(feed_crop_less_pasture)) * cropDM(feed_crop_less_pasture) =L= totalFeedDM * 0.7;
-
+
MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) .. otherIntensity(crop, location) =L= 1;
-
- SETASIDE_AREA_CALC(location) .. area('setaside', location) =E= sum(crop_less_pasture, area(crop_less_pasture, location)) * setAsideRate;
-
- TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, area(crop_less_pasture, location)) / (1.0 - unhandledCropRate) + area('pasture', location);
-
+
+ SETASIDE_AREA_CALC(location) .. cropArea('setaside', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) * setAsideRate;
+
+ TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate) + cropArea('pasture', location) +
+ sum(wood_producing, nonCropArea(wood_producing, location));
+
MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
PASTURE_IMPORT_CONSTRAINT .. importAmount('pasture') =E= 0;
PASTURE_EXPORT_CONSTRAINT .. exportAmount('pasture') =E= 0;
-
- IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) * (1.0 - unhandledCropRate) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * area(crop, location));
-
- AGRI_LAND_EXPANSION_CALC(location) .. agriLandExpansion(location) =G= sum(crop, area(crop, location) - previousArea(crop, location));
-
- CROP_INCREASE_CALC(location) .. cropIncrease(location) =G= sum(crop_less_pasture, area(crop_less_pasture, location) - previousArea(crop_less_pasture, location));
- CROP_DECREASE_CALC(location) .. cropDecrease(location) =G= -sum(crop_less_pasture, area(crop_less_pasture, location) - previousArea(crop_less_pasture, location));
- PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= area('pasture', location) - previousArea('pasture', location);
- PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(area('pasture', location) - previousArea('pasture', location));
- PASTURE_TOTAL_CHANGE_CONSTRAINT(location) .. pastureIncrease(location) -pastureDecrease(location) =G= area('pasture', location) - previousArea('pasture', location);
-
-* CROPLAND_INCREASE_CONSTRAINT(location) .. cropIncrease(location)/(1.0 - unhandledCropRate) =L= maxLandChangeRate * ( suitableLandArea(location) -sum(crop_less_pasture, previousArea(crop_less_pasture, location))/*(1.0 - unhandledCropRate) - previousArea('pasture', location) );
-* PASTURE_DECREASE_CONSTRAINT(location) .. pastureDecrease(location) =L= maxLandChangeRate * previousArea('pasture', location);
-
- COST_EQ .. total_cost =E=
- (
- ( SUM((crop, location), area(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) +
-
- SUM(location,
+
+ MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
+ MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
+
+ TIMBER_IMPORT_CONSTRAINT .. importAmount('timberForest') =E= 0;
+ TIMBER_EXPORT_CONSTRAINT .. exportAmount('timberForest') =E= 0;
+ CARBON_IMPORT_CONSTRAINT .. importAmount('carbonForest') =E= 0;
+ CARBON_EXPORT_CONSTRAINT .. exportAmount('carbonForest') =E= 0;
+ NATURAL_IMPORT_CONSTRAINT .. importAmount('natural') =E= 0;
+ NATURAL_EXPORT_CONSTRAINT .. exportAmount('natural') =E= 0;
+
+ IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) * (1.0 - unhandledCropRate) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * cropArea(crop, location));
+
+ AGRI_LAND_EXPANSION_CALC(location) .. agriLandExpansion(location) =G= sum(crop, cropArea(crop, location) - previousCropArea(crop, location));
+
+ CROP_INCREASE_CALC(location) .. cropIncrease(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location) - previousCropArea(crop_less_pasture, location));
+ CROP_DECREASE_CALC(location) .. cropDecrease(location) =G= -sum(crop_less_pasture, cropArea(crop_less_pasture, location) - previousCropArea(crop_less_pasture, location));
+ PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= cropArea('pasture', location) - previousCropArea('pasture', location);
+ PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(cropArea('pasture', location) - previousCropArea('pasture', location));
+ PASTURE_TOTAL_CHANGE_CONSTRAINT(location) .. pastureIncrease(location) -pastureDecrease(location) =G= cropArea('pasture', location) - previousCropArea('pasture', location);
+
+ WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum(wood_producing, nonCropArea(wood_producing, location) * woodYield(wood_producing, location));
+
+ CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum(wood_producing, nonCropArea(wood_producing, location) * carbonFluxRate(wood_producing, location));
+
+* CROPLAND_INCREASE_CONSTRAINT(location) .. cropIncrease(location)/(1.0 - unhandledCropRate) =L= maxLandChangeRate * ( suitableLandArea(location) -sum(crop_less_pasture, previousCropArea(crop_less_pasture, location))/*(1.0 - unhandledCropRate) - previousCropArea('pasture', location) );
+* PASTURE_DECREASE_CONSTRAINT(location) .. pastureDecrease(location) =L= maxLandChangeRate * previousCropArea('pasture', location);
+
+ COST_EQ .. total_cost =E=
+ (
+ ( SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) +
+
+ SUM(location,
agriExpansionCost(location) * agriLandExpansion(location) +
- cropIncCost * cropIncrease(location) +
- pastureIncCost * pastureIncrease(location) +
- cropDecCost * cropDecrease(location) +
+ cropIncCost * cropIncrease(location) +
+ pastureIncCost * pastureIncrease(location) +
+ cropDecCost * cropDecrease(location) +
pastureDecCost * pastureDecrease(location)
)
- ) * domesticPriceMarkup +
+ ) * domesticPriceMarkup +
+
+ SUM(location, carbonFlux(location) * carbonPrice) - SUM(location, woodHarvest(location) * woodPrice) +
- SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop))
+ SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop))
);
-
+
MODEL LAND_USE /ALL/ ;
fertI.L(crop, location) = previousFertIntensity(crop, location);
irrigI.L(crop, location) = previousIrrigIntensity(crop, location);
otherIntensity.L(crop, location) = previousOtherIntensity(crop, location);
- area.L(crop, location) = previousArea(crop, location);
+ cropArea.L(crop, location) = previousCropArea(crop, location);
+ nonCropArea.L(wood_producing, location) = previousNonCropArea(wood_producing, location);
ruminantFeed.L(feed_crop) = previousRuminantFeed(feed_crop);
monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop);
importAmount.L(all_types) = previousImportAmount(all_types);
exportAmount.L(all_types) = previousExportAmount(all_types);
-
+
SOLVE LAND_USE USING NLP MINIMIZING total_cost;
-
- display agriLandExpansion.L, previousArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, area.L, cropIncrease.L, cropDecrease.L, pastureIncrease.L, pastureDecrease.L;
+
+ display agriLandExpansion.L, previousCropArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, cropArea.L, cropIncrease.L, cropDecrease.L, pastureIncrease.L, pastureDecrease.L;
display net_supply.l, demand, ruminantFeed.l, monogastricFeed.l, importAmount.l, exportAmount.l;
-* Calculate summary information used in Java process
+* Calculate summary information used in Java process
parameter totalProd(all_types);
parameter totalProdCost(all_types);
parameter totalArea(crop);
@@ -245,32 +290,32 @@ $gdxin
parameter netImportCost(all_types);
parameter feedCostRate(feed_crop);
parameter productionShock(all_types);
-
+
* Production quantities based on smaller area (before unhandledCropArea adjustment applied)
- totalProd(crop) = sum(location, area.l(crop, location) * yield.l(crop, location));
- productionShock(crop) = sum(location, area.l(crop, location) * yield.l(crop, location) * yieldShock(crop, location));
+ totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location));
+ productionShock(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location) * yieldShock(crop, location));
totalProd('ruminants') = meatEfficency*(sum(feed_crop, ruminantFeed.l(feed_crop) * cropDM(feed_crop)));
totalProd('monogastrics') = meatEfficency*(sum(feed_crop, monogastricFeed.l(feed_crop) * cropDM(feed_crop)));
-
+
* Cost based on adjusted area
- area.l(crop_less_pasture, location) = area.l(crop_less_pasture, location) / (1.0 - unhandledCropRate);
- totalProdCost(crop) = sum(location, unitCost.l(crop, location) * area.l(crop, location));
- totalCropland(location) = sum(crop_less_pasture, area.l(crop_less_pasture, location));
- totalArea(crop) = sum(location, area.l(crop, location));
-
+ cropArea.l(crop_less_pasture, location) = cropArea.l(crop_less_pasture, location) / (1.0 - unhandledCropRate);
+ totalProdCost(crop) = sum(location, unitCost.l(crop, location) * cropArea.l(crop, location));
+ totalCropland(location) = sum(crop_less_pasture, cropArea.l(crop_less_pasture, location));
+ totalArea(crop) = sum(location, cropArea.l(crop, location));
+
feedCostRate(feed_crop)$[totalProd(feed_crop) <> 0] = totalProdCost(feed_crop) / totalProd(feed_crop);
totalProdCost('ruminants') = sum(feed_crop, ruminantFeed.l(feed_crop) * feedCostRate(feed_crop));
totalProdCost('monogastrics') = sum(feed_crop, monogastricFeed.l(feed_crop) * feedCostRate(feed_crop));
-
+
netImportCost(import_crop) = importAmount.l(import_crop) * importPrices(import_crop) - exportAmount.l(import_crop) * exportPrices(import_crop);
netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
-
+
Scalar totalCostsLessLU;
-
+
totalCostsLessLU = sum(crop,totalProdCost(crop))+ sum(import_crop,netImportCost(import_crop));
-
- Scalar ms 'model status', ss 'solve status';
+
+ Scalar ms 'model status', ss 'solve status';
ms=LAND_USE.modelstat;
ss=LAND_USE.solvestat;
-
\ No newline at end of file
+
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index 218d6086..d00a90bf 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -3,8 +3,10 @@ package ac.ed.lurg.country.gams;
import java.util.Map;
import ac.ed.lurg.Timestep;
+import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.WoodYieldItem;
import ac.ed.lurg.yield.YieldResponsesItem;
public class GamsLocationInput {
@@ -13,15 +15,20 @@ public class GamsLocationInput {
private Map<Integer, ? extends YieldResponsesItem> yields;
private Map<Integer, ? extends LandUseItem> previousLandUse;
private Map<Integer, ? extends IrrigationItem> irrigationCosts;
+ private Map<Integer, ? extends CarbonFluxItem> carbonFluxes;
+ private Map<Integer, ? extends WoodYieldItem> woodYields;
private GamsCountryInput countryInput;
public GamsLocationInput(Timestep timestep, Map<Integer, ? extends YieldResponsesItem> yields, Map<Integer, ? extends LandUseItem> previousLandUse,
- Map<Integer, ? extends IrrigationItem> irrigationCosts, GamsCountryInput countryInput) {
+ Map<Integer, ? extends IrrigationItem> irrigationCosts, Map<Integer, ? extends CarbonFluxItem> carbonFluxes,
+ Map<Integer, ? extends WoodYieldItem> woodYields, GamsCountryInput countryInput) {
super();
this.timestep = timestep;
this.yields = yields;
this.previousLandUse = previousLandUse;
this.irrigationCosts = irrigationCosts;
+ this.carbonFluxes = carbonFluxes;
+ this.woodYields = woodYields;
this.countryInput = countryInput;
}
@@ -37,6 +44,14 @@ public class GamsLocationInput {
return irrigationCosts;
}
+ public Map<Integer, ? extends CarbonFluxItem> getCarbonFluxes() {
+ return carbonFluxes;
+ }
+
+ public Map<Integer, ? extends WoodYieldItem> getWoodYields() {
+ return woodYields;
+ }
+
public GamsCountryInput getCountryInput() {
return countryInput;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index 072d0754..f6ef7f85 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -2,6 +2,8 @@ package ac.ed.lurg.country.gams;
import ac.ed.lurg.Timestep;
import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.WoodYieldItem;
+import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.yield.YieldRaster;
import ac.sac.raster.RasterSet;
@@ -12,15 +14,19 @@ public class GamsRasterInput {
private YieldRaster yields;
private RasterSet<LandUseItem> previousLandsUses;
private RasterSet<IrrigationItem> irrigationCost;
+ private RasterSet<CarbonFluxItem> carbonFluxes;
+ private RasterSet<WoodYieldItem> woodYields;
private GamsCountryInput countryInput;
- public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses,
- RasterSet<IrrigationItem> irrigationCost, GamsCountryInput countryInput) {
+ public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses, RasterSet<IrrigationItem> irrigationCost,
+ RasterSet<CarbonFluxItem> carbonFluxes, RasterSet<WoodYieldItem> woodYields, GamsCountryInput countryInput) {
super();
this.timestep = timestep;
this.yields = yields;
this.previousLandsUses = previousLandsUses;
this.irrigationCost = irrigationCost;
+ this.carbonFluxes = carbonFluxes;
+ this.woodYields = woodYields;
this.countryInput = countryInput;
}
@@ -36,6 +42,14 @@ public class GamsRasterInput {
return irrigationCost;
}
+ public RasterSet<CarbonFluxItem> getCarbonFluxes() {
+ return carbonFluxes;
+ }
+
+ public RasterSet<WoodYieldItem> getWoodYields() {
+ return woodYields;
+ }
+
public GamsCountryInput getCountryInput() {
return countryInput;
}
diff --git a/src/ac/ed/lurg/landuse/CarbonFluxItem.java b/src/ac/ed/lurg/landuse/CarbonFluxItem.java
new file mode 100644
index 00000000..e83ebf77
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/CarbonFluxItem.java
@@ -0,0 +1,22 @@
+package ac.ed.lurg.landuse;
+
+import ac.sac.raster.RasterItem;
+import java.util.HashMap;
+import java.util.Map;
+
+import ac.ed.lurg.types.LandCoverType;
+
+public class CarbonFluxItem implements RasterItem {
+ Map<LandCoverType, Map<LandCoverType, Double>> carbonFluxes = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
+
+ public void setCarbonFlux(LandCoverType previousLandCover, LandCoverType newLandCover, Double carbonFlux) {
+ Map<LandCoverType, Double> cfMap = new HashMap<LandCoverType, Double>();
+ cfMap.put(newLandCover, carbonFlux);
+ carbonFluxes.put(previousLandCover, cfMap);
+ }
+
+ public Double getCarbonFlux(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ return carbonFluxes.get(previousLandCover).get(newLandCover);
+ }
+
+}
diff --git a/src/ac/ed/lurg/landuse/CarbonFluxRasterSet.java b/src/ac/ed/lurg/landuse/CarbonFluxRasterSet.java
new file mode 100644
index 00000000..059bfac4
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/CarbonFluxRasterSet.java
@@ -0,0 +1,28 @@
+package ac.ed.lurg.landuse;
+
+import java.util.Collection;
+
+import ac.sac.raster.RasterHeaderDetails;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+
+public class CarbonFluxRasterSet extends RasterSet<CarbonFluxItem> {
+
+ private static final long serialVersionUID = 4986837799421955684L;
+
+ public CarbonFluxRasterSet(RasterHeaderDetails header) {
+ super(header);
+ }
+
+ protected CarbonFluxItem createRasterData() {
+ return new CarbonFluxItem();
+ }
+
+ // not very efficient, we could keep the mapping of country to area somewhere.
+ @Override
+ public CarbonFluxRasterSet createSubsetForKeys(Collection<RasterKey> keys) {
+ CarbonFluxRasterSet subsetCarbonFluxRaster = new CarbonFluxRasterSet(getHeaderDetails());
+ popSubsetForKeys(subsetCarbonFluxRaster, keys);
+ return subsetCarbonFluxRaster;
+ }
+}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldItem.java b/src/ac/ed/lurg/landuse/WoodYieldItem.java
new file mode 100644
index 00000000..ff7c23ae
--- /dev/null
+++ b/src/ac/ed/lurg/landuse/WoodYieldItem.java
@@ -0,0 +1,21 @@
+package ac.ed.lurg.landuse;
+
+import java.util.HashMap;
+import java.util.Map;
+
+import ac.ed.lurg.types.LandCoverType;
+import ac.sac.raster.RasterItem;
+
+public class WoodYieldItem implements RasterItem {
+ Map<LandCoverType, Map<LandCoverType, Double>> woodYields = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
+
+ public void setWoodYield(LandCoverType previousLandCover, LandCoverType newLandCover, Double woodYield) {
+ Map<LandCoverType, Double> cfMap = new HashMap<LandCoverType, Double>();
+ cfMap.put(newLandCover, woodYield);
+ woodYields.put(previousLandCover, cfMap);
+ }
+
+ public Double getWoodYield(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ return woodYields.get(previousLandCover).get(newLandCover);
+ }
+}
diff --git a/src/ac/ed/lurg/types/NonCropType.java b/src/ac/ed/lurg/types/NonCropType.java
new file mode 100644
index 00000000..ad5f9a6d
--- /dev/null
+++ b/src/ac/ed/lurg/types/NonCropType.java
@@ -0,0 +1,53 @@
+package ac.ed.lurg.types;
+
+import java.util.Collection;
+import java.util.HashSet;
+
+public enum NonCropType {
+ TIMBER_FOREST("timberForest", true, true),
+ CARBON_FOREST("carbonForest", true, true),
+ NATURAL("natural", true, false);
+
+ private String gamsName;
+ private boolean isWoodProducing;
+ private boolean isManagedForest;
+
+ private NonCropType (String gamsName, boolean isWoodProducing, boolean isManagedForest) {
+ this.gamsName = gamsName;
+ this.isWoodProducing = isWoodProducing;
+ this.isManagedForest = isManagedForest;
+ }
+
+ public String getGamsName() {
+ return gamsName;
+ }
+
+ public boolean isWoodProducing() {
+ return isWoodProducing;
+ }
+
+ public boolean isManagedForest() {
+ return isManagedForest;
+ }
+
+ public static Collection<NonCropType> getWoodProducingTypes() {
+ Collection<NonCropType> lu = new HashSet<NonCropType>();
+
+ for (NonCropType c : values())
+ if (c.isWoodProducing)
+ lu.add(c);
+
+ return lu;
+ }
+
+ public static Collection<NonCropType> getManagedForestTypes() {
+ Collection<NonCropType> lu = new HashSet<NonCropType>();
+
+ for (NonCropType c : values())
+ if (c.isManagedForest)
+ lu.add(c);
+
+ return lu;
+ }
+
+}
--
GitLab