-
Peter Alexander authoredPeter Alexander authored
IntExtOpt.gms 14.92 KiB
SET all_types / meat, cereals, wheat, maize, rice, oilcrops, pulses, starchyRoots, pasture /;
SET crop(all_types) / wheat, maize, rice, oilcrops, pulses, starchyRoots, pasture /;
SET crop_less_pasture(crop) / wheat, maize, rice, oilcrops, pulses, starchyRoots /;
SET cereal_crop(crop) / wheat, maize, rice /;
SET non_cereal_crop(crop) / oilcrops, pulses, starchyRoots, pasture /;
SET feed_crop(crop) / wheat, maize, oilcrops, pasture/;
SET feed_crop_less_pasture(feed_crop) / wheat, maize, oilcrops /;
SET not_feed_crop(crop) / rice, pulses, starchyRoots /;
SET import_crop(all_types) / meat, wheat, maize, rice, oilcrops, pulses, starchyRoots /;
SET location;
PARAMETER suitableLandArea(location) areas of land in Mha;
PARAMETER previousArea(crop, location) areas for previous timestep in Mha;
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 fertParam(crop, location) yield response to fertilizer parameter;
PARAMETER irrigParam(crop, location) yield response to irrigation parameter;
PARAMETER demand(all_types) in t;
PARAMETER worldExportPrices(import_crop) prices for exports;
PARAMETER worldImportPrices(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 energy 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 cropAdj(crop) this is the fudge factor that allows the model to fit the observed data
SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
SCALAR landChangeEnergy energy required to add ha of agricultural land;
SCALAR minFeedRate minimum rate of feed for producing animal products (why is this needed?);
SCALAR fertiliserUnitCost fert cost at max fert rate;
*$gdxin "/Users/peteralexander/Documents/R_Workspace/temp/GamsTest/_gams_java_1091455539.gdx"
$gdxin %gdxincname%
$load location, suitableLandArea, previousArea, demand, landChangeEnergy
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
$load fertParam, irrigParam, worldExportPrices, worldImportPrices, maxNetImport, minNetImport
$load meatEfficency, minFeedRate, irrigCost, irrigMaxRate, irrigConstraint, cropAdj, fertiliserUnitCost
$gdxin
SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /
SCALAR minDemandPerCereal / 0.1 /;
SCALAR unhandledCropArea includes fruit veg forage crops set aside and failed crop / 0.4 /;
previousArea(crop_less_pasture, location) = previousArea(crop_less_pasture, location) * (1.0 - unhandledCropArea);
PARAMETER feedDM(crop) kg DM per kg of feed the conversion from feed to meat is done in the R animal product index
/ wheat 0.89
maize 0.7
oilcrops 0.9
pasture 1 / ;
* pasture yield is done in DM terms
PARAMETER baseCost(crop) cost per ha before intensity values just made up at the moment
/ wheat 1.1
maize 0.9
rice 1.1
oilcrops 0.7
pulses 0.5
starchyRoots 4.0
pasture 0.0 / ;
VARIABLES
area(crop, location) total area for each crop - 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 otherIntensity(crop, location)
feedAmount(crop) amount of feed for each crop - Mt
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
unitEnergy(crop, location) energy per area for each crop - energy
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)
energy total input energy - energy;
POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount, importAmount, exportAmount,
agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease;
EQUATIONS
UNIT_ENERGY_EQ(crop, location) energy per area
YIELD_EQ(crop, location) yield given chosen intensity
NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) satisfy demand for non-cereal crop
CEREAL_DEMAND_CONSTRAINT(cereal_crop) satisfy demand for cereal so that exports can at least be met. Could also allow min. proporation of cereal consumption of each type
TOTAL_CEREAL_DEMAND_CONSTRAINT satisfy demand for cereal
MEAT_DEMAND_CONSTRAINT satisfy demand for meat
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)
TOTAL_LAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
NON_FEED_CROP_CONSTRAINT(not_feed_crop) constraint to set non feed crop feed usage to zero
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
MIN_FEED_CONSTRAINT constraint on min feed rate
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)
ENERGY_EQ total energy objective function;
UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= ( baseCost(crop) +
fertiliserUnitCost * fertI(crop, location) +
irrigCost(location) * irrigMaxRate(crop, location) * irrigI(crop, location) +
SQR(otherIntensity(crop, location)) * 1
) ;
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yieldNone(crop, location) +
* (yieldFertOnly(crop, location) - yieldNone(crop, location)) * ((fertI(crop, location) + delta) ** fertParam(crop, location)) +
* (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) +
* (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
* ((fertI(crop, location) + delta) ** fertParam(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(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)))
) * cropAdj(crop) * otherIntensity(crop, location);
NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + importAmount(crop) - exportAmount(crop);
NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop) ;
CEREAL_DEMAND_CONSTRAINT(cereal_crop) .. net_supply(cereal_crop) =G= demand('cereals') * minDemandPerCereal;
TOTAL_CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, net_supply(cereal_crop)) =G= demand('cereals');
MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= demand('meat') - importAmount('meat') - exportAmount('meat');
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;
TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop_less_pasture, area(crop_less_pasture, location)) / (1.0 - unhandledCropArea) + area('pasture', location);
NON_FEED_CROP_CONSTRAINT(not_feed_crop) .. feedAmount(not_feed_crop) =E= 0;
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') - exportAmount('pasture') =E= 0;
MIN_FEED_CONSTRAINT .. sum(feed_crop_less_pasture, feedDM(feed_crop_less_pasture) * feedAmount(feed_crop_less_pasture)) =G= minFeedRate * (demand('meat') - importAmount('meat') + exportAmount('meat'));
IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * area(crop, location)) / suitableLandArea(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));
ENERGY_EQ .. energy =E=
(
SUM((crop, location), area(crop, location) * unitEnergy(crop, location)) +
sum(location,
agriLandExpansion(location) +
0.45 * SQR(cropIncrease(location) + cropDecrease(location)) / (suitableLandArea(location) * 0.01 + sum(crop_less_pasture, previousArea(crop_less_pasture, location))) +
0.05 * (cropIncrease(location) + cropDecrease(location)) +
0.45 * SQR(pastureIncrease(location) + pastureDecrease(location)) / (suitableLandArea(location) * 0.01 + previousArea('pasture', location)) +
0.05 * (pastureIncrease(location) + pastureDecrease(location))
) * landChangeEnergy +
sum(import_crop, importAmount(import_crop) * worldImportPrices(import_crop) - exportAmount(import_crop) * worldExportPrices(import_crop))
) / 1000000;
MODEL LAND_USE /ALL/ ;
fertI.L(crop, location) = 0.5;
irrigI.L(crop, location) = 0.5;
otherIntensity.L(crop, location) = 0.1;
importAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) gt 0) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
exportAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) lt 0) = -(maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
area.L(crop, location) = previousArea(crop, location)
SOLVE LAND_USE USING NLP MINIMIZING energy;
* display agriLandExpansion.l, previousArea, area.l, net_supply.l, demand, feedAmount.l, yield.l, netImportAmount.l;
area.l(crop_less_pasture, location) = area.l(crop_less_pasture, location) / (1.0 - unhandledCropArea);
parameter totalProd(all_types);
parameter totalProdCost(all_types);
parameter totalCropland(location);
parameter netImportAmount(all_types);
parameter netImportEnergy(all_types);
parameter feedEnergy(all_types);
totalProd(crop) = sum(location, area.l(crop, location) * yield.l(crop, location));
totalProd('meat') = meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount.l(feed_crop));
totalProdCost(crop) = sum(location, unitEnergy.l(crop, location) * area.l(crop, location));
totalCropland(location) = sum(crop_less_pasture, area.l(crop_less_pasture, location));
netImportEnergy(import_crop) = importAmount.l(import_crop) * worldImportPrices(import_crop) - exportAmount.l(import_crop) * worldExportPrices(import_crop);
netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
feedEnergy(feed_crop)$(netImportAmount(feed_crop) gt 0) = (totalProdCost(feed_crop) + netImportEnergy(feed_crop)) * feedAmount.l(feed_crop) / (totalProd(feed_crop) + netImportAmount(feed_crop));
feedEnergy(feed_crop)$(netImportAmount(feed_crop) le 0 and totalProd(feed_crop) gt 0) = totalProdCost(feed_crop) * feedAmount.l(feed_crop) / totalProd(feed_crop);
totalProdCost('meat') = sum(feed_crop, feedEnergy(feed_crop));
* display totalProdCost, totalProd, totalImportEnergy, feedEnergy, totalProdCost
Scalar ms 'model status', ss 'solve status';
ms=LAND_USE.modelstat;
ss=LAND_USE.solvestat;