Newer
Older
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 land_cover / cropland, pasture, timberForest, carbonForest, unmanagedForest, otherNatural /;
SET forest(land_cover) / timberForest, carbonForest, unmanagedForest /;
SET non_forest(land_cover) / cropland, pasture, otherNatural /;
ALIAS (land_cover, land_cover_before);
ALIAS (land_cover, land_cover_after);
PARAMETER suitableLandArea(location) areas of land in Mha;
PARAMETER previousCropArea(crop, location) areas for previous timestep in Mha;
Peter Alexander
committed
PARAMETER previousFertIntensity(crop, location);
PARAMETER previousIrrigIntensity(crop, location);
Peter Alexander
committed
PARAMETER previousOtherIntensity(crop, location);
PARAMETER previousRuminantFeed(crop);
PARAMETER previousMonogastricFeed(crop);
PARAMETER previousImportAmount(all_types);
PARAMETER previousExportAmount(all_types);
PARAMETER yieldNone(crop, location) yield in t per ha;
PARAMETER yieldFertOnly(crop, location) yield in t per ha;
PARAMETER yieldIrrigOnly(crop, location) yield in t per ha;
PARAMETER yieldBoth(crop, location) yield in t per ha;
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 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;
PARAMETER minNetImport(import_crop) minimum net import for each crop based on world market;
PARAMETER irrigCost(location) irrigation cost in cost per 1000 Mlitre or Mha for each litre per m2;
PARAMETER irrigMaxRate(crop, location) max water application rate irrigation in litre per m2;
PARAMETER irrigConstraint(location) max water available for irrigation in litre per m2;
PARAMETER minDemandPerCereal(cereal_crop) min demand for each cereal crop as factor of all cereals;
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 - 1000$ per ha;
SCALAR pastureIncCost price for increasing pasture area - 1000$ per ha;
SCALAR cropIncCost price for increasing crop - 1000$ per ha;
SCALAR cropDecCost price for decreasing crop - 1000$ per ha;
SCALAR pastureDecCost price for decreasing pasture - 1000$ per ha;
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;
SCALAR otherICost cost of other intensity;
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;
PARAMETER previousLandCoverArea(land_cover, location) land cover area in Mha;
PARAMETER minimumLandCoverArea(land_cover, location) minimum land cover area to constrain conversion;
PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location) carbon flux - MtC-eq per Mha;
PARAMETER woodYield(land_cover_before, land_cover_after, location) wood yield - MtC-eq per Mha;
SCALAR carbonPrice price of carbon - $ per tonne or million$ per Mt;
SCALAR woodPrice price of wood and timber - $ per tonne or million$ per Mt;
*$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
$gdxin %gdxincname%
$load location, suitableLandArea, demand, agriExpansionCost, cropIncCost, pastureIncCost, cropDecCost, pastureDecCost
$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
Peter Alexander
committed
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
$load previousLandCoverArea, carbonFluxRate, woodYield, carbonPrice, woodPrice
$load minimumLandCoverArea
*convert units from $/t to billion$/Mt
carbonPrice = carbonPrice * 0.001;
woodPrice = woodPrice * 0.001;
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);
previousCropArea(crop_less_pasture, location) = previousCropArea(crop_less_pasture, location) * (1.0 - unhandledCropRate);
Peter Alexander
committed
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
oilcrops 0.96
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
starchyRoots 3.14
PARAMETER baseCost(crop);
PARAMETER otherIntCost(crop);
baseCost(crop) = prodCost(crop)*0.3;
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);
SCALAR forestExpansionMaxRate / 0.01 /;
cropArea(crop, location) total area for crops - Mha
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
ruminantFeed(crop) amount of feed for ruminant animals - Mt
monogastricFeed(crop) amount of feed for monogatric animals - Mt
Peter Alexander
committed
importAmount(all_types) imports of crops and meat - Mt
exportAmount(all_types) exports of crops and meat - Mt
yield(crop, location) yield per area for each crop - t per ha
unitCost(crop, location) cost per area for each crop - cost
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)
landCoverArea(land_cover, location) land cover area in Mha
landCoverChange(land_cover_before, land_cover_after, location) land cover change
prematureDeforestationCost(location)
woodHarvest(location) total wood harvested
carbonFlux(location) total carbon flux
totalFeedDM
* A "artificial variable https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
total_cost total cost of domestic supply including net imports;
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease, totalFeedDM,
landCoverArea, landCoverChange, woodHarvest;
UNIT_COST_EQ(crop, location) cost per area - $1000 per ha or $billion per Mha
YIELD_EQ(crop, location) yield given chosen intensity - tonnes per hectare
CROP_DEMAND_CONSTRAINT(crop) satisfy demand for individual crops
TOTAL_CEREAL_DEMAND_CONSTRAINT satisfy demand for combined cereals
TOTAL_OIL_PULSE_DEMAND_CONSTRAINT satisfy demand for combined oilcrops and pulses
RUMINANT_DEMAND_CONSTRAINT satisfy demand for ruminant products
MONOGASTRICS_DEMAND_CONSTRAINT satisfy demand for monogastric products
TOTAL_NON_PASTURE_FEED_DM_CALC calc total feed dry matter not including pasture
FEED_MIX_CONSTRAINT(feed_crop_less_pasture) limit amount of feed for each feed crop
MAX_FERT_INTENSITY_CONSTRAINT(crop, location) constraint on maximum fertilizer intensity
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
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
IRRIGATION_CONSTRAINT(location) constraint on water usage
CROP_INCREASE_CALC(location)
CROP_DECREASE_CALC(location)
PASTURE_INCREASE_CONV_CALC(location)
PASTURE_DECREASE_CONV_CALC(location)
PASTURE_TOTAL_CHANGE_CONSTRAINT(location)
CROPLAND_LAND_COVER_CALC(location) cropland area equals sum of all crop areas
PASTURE_LAND_COVER_CALC(location) pasture area (land cover) equals pasture area (land use)
* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) constraint on land cover conversion
LAND_COVER_CHANGE_CALC(land_cover, location)
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) can't convert more land than was previously available
LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) calculate no change area (e.g. cropland to cropland)
FOREST_EXPANSION_CONSTRAINT
PREMATURE_DEFORESTATION_COST_CALC(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) +
irrigCost(location) * irrigMaxRate(crop, location) * irrigI(crop, location) +
otherIntCost(crop) * otherIntensity(crop, location)
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
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, 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');
Peter Alexander
committed
RUMINANT_DEMAND_CONSTRAINT .. meatEfficency * sum(feed_crop, ruminantFeed(feed_crop) * cropDM(feed_crop)) * (1 - seedAndWasteRate('ruminants')) =G= (demand('ruminants') - importAmount('ruminants') + exportAmount('ruminants'));
Peter Alexander
committed
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) .. 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);
* Needs inequality due to floating point errors
TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(land_cover, landCoverArea(land_cover, location));
Peter Alexander
committed
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) * 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);
CROPLAND_LAND_COVER_CALC(location) .. landCoverArea('cropland', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate);
PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= previousLandCoverArea(land_cover, location) +
sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) - sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =L= previousLandCoverArea(land_cover, location);
LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =E= previousLandCoverArea(land_cover, location) -
sum(land_cover_after$[not sameAs(land_cover, land_cover_after)], landCoverChange(land_cover, land_cover_after, location));
FOREST_EXPANSION_CONSTRAINT .. sum((non_forest, forest, location), landCoverChange(non_forest, forest, location)) =L= sum(location, suitableLandArea(location)) * forestExpansionMaxRate;
* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =G= minimumLandCoverArea(land_cover, location);
PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum(forest, (sum((land_cover_after)$[not sameAs(forest, land_cover_after)], landCoverChange(forest, land_cover_after, location)) -
previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location))) * 1000;
WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYield(land_cover_before, land_cover_after, location));
CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, 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) +
pastureDecCost * pastureDecrease(location)
) * 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(location, prematureDeforestationCost(location))
Peter Alexander
committed
fertI.L(crop, location) = previousFertIntensity(crop, location);
irrigI.L(crop, location) = previousIrrigIntensity(crop, location);
otherIntensity.L(crop, location) = previousOtherIntensity(crop, location);
cropArea.L(crop, location) = previousCropArea(crop, location);
landCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, 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, 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
parameter totalProd(all_types);
parameter totalProdCost(all_types);
Peter Alexander
committed
parameter totalCropland(location);
parameter netImportAmount(all_types);
parameter netImportCost(all_types);
parameter productionShock(all_types);
* Production quantities based on smaller area (before unhandledCropArea adjustment applied)
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));
Peter Alexander
committed
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
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);
Peter Alexander
committed
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';
ms=LAND_USE.modelstat;
ss=LAND_USE.solvestat;