Skip to content
Snippets Groups Projects
Commit 35be1dcb authored by Peter Alexander's avatar Peter Alexander
Browse files

Diminishing returns on fertiliser/irrigation

parent f089362a
No related branches found
No related tags found
No related merge requests found
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
SET crop_less_pasture(crop) arable crops types includes treenuts but not pasture / cereals, fruits, oilcrops, starchyRoots, treenuts, vegetables, pulses /; SET crop_less_pasture(crop) arable crops types includes treenuts but not pasture / cereals, fruits, oilcrops, starchyRoots, treenuts, vegetables, pulses /;
SET feed_crop(crop) / cereals, oilcrops /; SET feed_crop(crop) / cereals, oilcrops /;
SET not_feed_crops(crop) / fruits, starchyRoots, treenuts, vegetables, pulses, pasture_or_meat /; SET not_feed_crops(crop) / fruits, starchyRoots, treenuts, vegetables, pulses, pasture_or_meat /;
SET location / 1*2 /; SET location / 1*1 /;
PARAMETER previous_area(crop, location) areas for previous timestep in ha; PARAMETER previous_area(crop, location) areas for previous timestep in ha;
PARAMETER yield00(crop, location) yield in t per ha; PARAMETER yield00(crop, location) yield in t per ha;
...@@ -20,13 +20,17 @@ ...@@ -20,13 +20,17 @@
SCALAR tradeBarrier trade barrier which adjust energy cost of imports; SCALAR tradeBarrier trade barrier which adjust energy cost of imports;
SCALAR landChangeEnergy energy required to add ha of agricultural land; SCALAR landChangeEnergy energy required to add ha of agricultural land;
SCALAR minFeedRate minimum rate of feed for producing animal products; SCALAR minFeedRate minimum rate of feed for producing animal products;
SCALAR alpha yield response to fertilizer parameter;
SCALAR beta yield response to irrigation parameter;
$gdxin %gdxincname% $gdxin %gdxincname%
$load previous_area, demand, yield00, yield10, yield01, yield11, world_input_energy, maxNetImport, minNetImport, meatEfficency, maxLandUseChange, tradeBarrier, landChangeEnergy, minFeedRate $load previous_area, demand, yield00, yield10, yield01, yield11, alpha, beta, world_input_energy, maxNetImport, minNetImport, meatEfficency, maxLandUseChange, tradeBarrier, landChangeEnergy, minFeedRate
$gdxin $gdxin
DISPLAY yield00, yield10, yield01, yield11 DISPLAY yield00, yield10, yield01, yield11
SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.000000001 /
PARAMETER feedEnergy(crop) energy from feed in MJ per t PARAMETER feedEnergy(crop) energy from feed in MJ per t
/ cereals 13.7 / cereals 13.7
oilcrops 18.6 oilcrops 18.6
...@@ -48,7 +52,7 @@ DISPLAY yield00, yield10, yield01, yield11 ...@@ -48,7 +52,7 @@ DISPLAY yield00, yield10, yield01, yield11
unit_energy(crop, location) energy per area for each crop - energy unit_energy(crop, location) energy per area for each crop - energy
energy total input energy - energy; energy total input energy - energy;
POSITIVE VARIABLE area, fertI, irrigI, otherI, feedAmount; POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount;
EQUATIONS EQUATIONS
UNIT_ENERGY_EQ(crop, location) energy per area UNIT_ENERGY_EQ(crop, location) energy per area
...@@ -59,7 +63,8 @@ DISPLAY yield00, yield10, yield01, yield11 ...@@ -59,7 +63,8 @@ DISPLAY yield00, yield10, yield01, yield11
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity
CROPLAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change CROPLAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
PASTURE_CHANGE_CONSTRAINT(location) constraint on the rate of land use change PASTURE_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
AGRI_LAND_CHANGE_CONSTRAINT(location) constraint expansion of agricultural land AGRI_LAND_INCREASE_CONSTRAINT(location) constraint expansion of agricultural land
AGRI_LAND_DECREASE_CONSTRAINT(location) constraint expansion of agricultural land
NON_FEED_CROP_CONSTRAINT(not_feed_crops) constraint to set non feed crop feed usage to zero NON_FEED_CROP_CONSTRAINT(not_feed_crops) constraint to set non feed crop feed usage to zero
MAX_NET_IMPORT_CONSTRAINT(crop) constraint on max net imports MAX_NET_IMPORT_CONSTRAINT(crop) constraint on max net imports
MIN_NET_IMPORT_CONSTRAINT(crop) constraint on min net imports MIN_NET_IMPORT_CONSTRAINT(crop) constraint on min net imports
...@@ -70,10 +75,10 @@ DISPLAY yield00, yield10, yield01, yield11 ...@@ -70,10 +75,10 @@ DISPLAY yield00, yield10, yield01, yield11
sqr(fertI(crop, location)+irrigI(crop, location)+otherIntensity(crop, location)) / 2; sqr(fertI(crop, location)+irrigI(crop, location)+otherIntensity(crop, location)) / 2;
YIELD_EQ(crop, location) .. yield(crop, location) =E= ( YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yield00(crop, location)*(1-fertI(crop, location))*(1-irrigI(crop, location)) + yield00(crop, location) +
yield10(crop, location)*fertI(crop, location)*(1-irrigI(crop, location)) + (yield10(crop, location) - yield00(crop, location)) * ((fertI(crop, location) + delta) ** alpha) +
yield01(crop, location)*(1-fertI(crop, location))*irrigI(crop, location) + (yield10(crop, location) - yield00(crop, location)) * ((irrigI(crop, location) + delta) ** beta) +
yield11(crop, location)*fertI(crop, location)*irrigI(crop, location) (yield11(crop, location) + yield00(crop, location) - yield10(crop, location) - yield01(crop, location)) * ((fertI(crop, location) + delta) ** alpha) * ((irrigI(crop, location) + delta) ** beta)
) * otherIntensity(crop, location); ) * otherIntensity(crop, location);
CROP_DEMAND_CONSTRAINT(crop_less_pasture) .. sum(location, area(crop_less_pasture, location) * yield(crop_less_pasture, location)) - feedAmount(crop_less_pasture) =G= CROP_DEMAND_CONSTRAINT(crop_less_pasture) .. sum(location, area(crop_less_pasture, location) * yield(crop_less_pasture, location)) - feedAmount(crop_less_pasture) =G=
...@@ -90,7 +95,9 @@ DISPLAY yield00, yield10, yield01, yield11 ...@@ -90,7 +95,9 @@ DISPLAY yield00, yield10, yield01, yield11
PASTURE_CHANGE_CONSTRAINT(location) .. 1 - area('pasture_or_meat', location)/previous_area('pasture_or_meat', location) =L= maxLandUseChange; PASTURE_CHANGE_CONSTRAINT(location) .. 1 - area('pasture_or_meat', location)/previous_area('pasture_or_meat', location) =L= maxLandUseChange;
AGRI_LAND_CHANGE_CONSTRAINT(location) .. abs(sum(crop, area(crop, location)) / sum(crop, previous_area(crop, location)) - 1) =L= maxLandUseChange; AGRI_LAND_INCREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) / sum(crop, previous_area(crop, location)) =L= 1 + maxLandUseChange;
AGRI_LAND_DECREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) / sum(crop, previous_area(crop, location)) =G= 1 - maxLandUseChange;
NON_FEED_CROP_CONSTRAINT(not_feed_crops) .. feedAmount(not_feed_crops) =E= 0; NON_FEED_CROP_CONSTRAINT(not_feed_crops) .. feedAmount(not_feed_crops) =E= 0;
...@@ -105,7 +112,7 @@ DISPLAY yield00, yield10, yield01, yield11 ...@@ -105,7 +112,7 @@ DISPLAY yield00, yield10, yield01, yield11
+ sum(crop, (netImportAmount(crop) * tradeBarrier) * world_input_energy(crop)); + sum(crop, (netImportAmount(crop) * tradeBarrier) * world_input_energy(crop));
MODEL LAND_USE /ALL/ ; MODEL LAND_USE /ALL/ ;
SOLVE LAND_USE USING DNLP MINIMIZING energy; SOLVE LAND_USE USING NLP MINIMIZING energy;
DISPLAY yield00, yield10, yield01, yield11 DISPLAY yield00, yield10, yield01, yield11
......
...@@ -69,10 +69,36 @@ public class CountryAgent implements GamsInputData { ...@@ -69,10 +69,36 @@ public class CountryAgent implements GamsInputData {
} }
@Override @Override
public Map<CropType, Double> getRefYield() { public Map<CropType, Double> getYield00() {
return currentRefYield; return currentRefYield;
} }
@Override
public Map<CropType, Double> getYield10() {
return currentRefYield;
}
@Override
public Map<CropType, Double> getYield01() {
return currentRefYield;
}
@Override
public Map<CropType, Double> getYield11() {
return currentRefYield;
}
@Override
public double getAlpha() {
return 0.4;
}
@Override
public double getBeta() {
return 0.5;
}
@Override @Override
public Map<CropType, Double> getPreviousCropArea() { public Map<CropType, Double> getPreviousCropArea() {
Map<CropType, Double> previousCropAreas = new HashMap<CropType, Double>(); Map<CropType, Double> previousCropAreas = new HashMap<CropType, Double>();
......
...@@ -7,7 +7,10 @@ import ac.ed.lurg.types.CropType; ...@@ -7,7 +7,10 @@ import ac.ed.lurg.types.CropType;
public interface GamsInputData { public interface GamsInputData {
Map<CropType, Double> getProjectedDemand(); Map<CropType, Double> getProjectedDemand();
Map<CropType, Double> getRefYield(); Map<CropType, Double> getYield00();
Map<CropType, Double> getYield10();
Map<CropType, Double> getYield01();
Map<CropType, Double> getYield11();
Map<CropType, Double> getPreviousCropArea(); Map<CropType, Double> getPreviousCropArea();
Map<CropType, Double> getWorldInputEnergy(); Map<CropType, Double> getWorldInputEnergy();
Map<CropType, Double> getMaxNetImport(); Map<CropType, Double> getMaxNetImport();
...@@ -18,4 +21,7 @@ public interface GamsInputData { ...@@ -18,4 +21,7 @@ public interface GamsInputData {
double getTradeBarrier(); double getTradeBarrier();
double getLandChangeEnergy(); double getLandChangeEnergy();
double getMinFeedRate(); double getMinFeedRate();
double getAlpha();
double getBeta();
} }
...@@ -39,10 +39,12 @@ public class GamsLandUseOptimiser { ...@@ -39,10 +39,12 @@ public class GamsLandUseOptimiser {
GAMSWorkspace ws = new GAMSWorkspace(wsInfo); GAMSWorkspace ws = new GAMSWorkspace(wsInfo);
GAMSDatabase db = ws.addDatabase(); GAMSDatabase db = ws.addDatabase();
int numLocations = 1;
LogWriter.println("\nPrevious areas"); LogWriter.println("\nPrevious areas");
GAMSParameter parm = db.addParameter("previous_area", 2, "the previous area for each land use"); GAMSParameter parm = db.addParameter("previous_area", 2, "the previous area for each land use");
addLocationMapParm(parm, inputData.getPreviousCropArea(), "1"); for (int i= 1; i<=numLocations; i++)
addLocationMapParm(parm, inputData.getPreviousCropArea(), "2"); addLocationMapParm(parm, inputData.getPreviousCropArea(), i);
LogWriter.println("\nDemand"); LogWriter.println("\nDemand");
parm = db.addParameter("demand", 1, "demand for crop"); parm = db.addParameter("demand", 1, "demand for crop");
...@@ -50,23 +52,23 @@ public class GamsLandUseOptimiser { ...@@ -50,23 +52,23 @@ public class GamsLandUseOptimiser {
LogWriter.println("\nYield"); LogWriter.println("\nYield");
parm = db.addParameter("yield00", 2, "ref yield t per ha"); parm = db.addParameter("yield00", 2, "ref yield t per ha");
addLocationMapParm(parm, inputData.getRefYield(), "1"); for (int i= 1; i<=numLocations; i++)
addLocationMapParm(parm, inputData.getRefYield(), "2"); addLocationMapParm(parm, inputData.getYield00(), i);
LogWriter.println("\nYield"); LogWriter.println("\nYield");
parm = db.addParameter("yield10", 2, "ref yield t per ha"); parm = db.addParameter("yield10", 2, "ref yield t per ha");
addLocationMapParm(parm, inputData.getRefYield(), "1"); for (int i= 1; i<=numLocations; i++)
addLocationMapParm(parm, inputData.getRefYield(), "2"); addLocationMapParm(parm, inputData.getYield10(), i);
LogWriter.println("\nYield"); LogWriter.println("\nYield");
parm = db.addParameter("yield01", 2, "ref yield t per ha"); parm = db.addParameter("yield01", 2, "ref yield t per ha");
addLocationMapParm(parm, inputData.getRefYield(), "1"); for (int i= 1; i<=numLocations; i++)
addLocationMapParm(parm, inputData.getRefYield(), "2"); addLocationMapParm(parm, inputData.getYield01(), i);
LogWriter.println("\nYield"); LogWriter.println("\nYield");
parm = db.addParameter("yield11", 2, "ref yield t per ha"); parm = db.addParameter("yield11", 2, "ref yield t per ha");
addLocationMapParm(parm, inputData.getRefYield(), "1"); for (int i= 1; i<=numLocations; i++)
addLocationMapParm(parm, inputData.getRefYield(), "2"); addLocationMapParm(parm, inputData.getYield11(), i);
LogWriter.println("\nWorld input energy"); LogWriter.println("\nWorld input energy");
parm = db.addParameter("world_input_energy", 1, "average input energy from world exports used to determine if we should import or export energy per t"); parm = db.addParameter("world_input_energy", 1, "average input energy from world exports used to determine if we should import or export energy per t");
...@@ -85,6 +87,8 @@ public class GamsLandUseOptimiser { ...@@ -85,6 +87,8 @@ public class GamsLandUseOptimiser {
addScalar(db.addParameter("tradeBarrier", 0, "trade barrier which adjust energy cost of imports"), inputData.getTradeBarrier()); addScalar(db.addParameter("tradeBarrier", 0, "trade barrier which adjust energy cost of imports"), inputData.getTradeBarrier());
addScalar(db.addParameter("landChangeEnergy", 0, "energy required to add ha of agricultural land"), inputData.getLandChangeEnergy()); addScalar(db.addParameter("landChangeEnergy", 0, "energy required to add ha of agricultural land"), inputData.getLandChangeEnergy());
addScalar(db.addParameter("minFeedRate", 0, "minimum rate of feed for producing animal products"), inputData.getMinFeedRate()); addScalar(db.addParameter("minFeedRate", 0, "minimum rate of feed for producing animal products"), inputData.getMinFeedRate());
addScalar(db.addParameter("alpha", 0), inputData.getAlpha());
addScalar(db.addParameter("beta", 0), inputData.getBeta());
GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.GAMS_MODEL); GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.GAMS_MODEL);
GAMSOptions opt = ws.addOptions(); GAMSOptions opt = ws.addOptions();
...@@ -139,13 +143,13 @@ public class GamsLandUseOptimiser { ...@@ -139,13 +143,13 @@ public class GamsLandUseOptimiser {
} }
} }
private void addLocationMapParm(GAMSParameter parm, Map<CropType, Double> itemMap, String locationId) { private void addLocationMapParm(GAMSParameter parm, Map<CropType, Double> itemMap, int locationId) {
for (Map.Entry<CropType, Double> entry : itemMap.entrySet()) { for (Map.Entry<CropType, Double> entry : itemMap.entrySet()) {
LogWriter.println(String.format(" %15s,\t %.1f", entry.getKey().getGamsName(), entry.getValue())); LogWriter.println(String.format(" %15s,\t %.1f", entry.getKey().getGamsName(), entry.getValue()));
Vector<String> v = new Vector<String>(); Vector<String> v = new Vector<String>();
v.add(entry.getKey().getGamsName()); v.add(entry.getKey().getGamsName());
v.add(locationId); v.add(Integer.toString(locationId));
if (locationId.equals("2") & entry.getKey().getGamsName().equals("cereals")) if (locationId== 2 & entry.getKey().getGamsName().equals("cereals"))
parm.addRecord(v).setValue(0.99); parm.addRecord(v).setValue(0.99);
else else
parm.addRecord(v).setValue(entry.getValue()); parm.addRecord(v).setValue(entry.getValue());
......
...@@ -23,28 +23,67 @@ public class GamsTest implements GamsInputData { ...@@ -23,28 +23,67 @@ public class GamsTest implements GamsInputData {
@Override @Override
public Map<CropType, Double> getProjectedDemand() { public Map<CropType, Double> getProjectedDemand() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>(); Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
dummyMap.put(CropType.CEREALS, 290.0); dummyMap.put(CropType.CEREALS, 50.0);
dummyMap.put(CropType.FRUIT, 50.0); dummyMap.put(CropType.FRUIT, 250.0);
dummyMap.put(CropType.OILCROPS, 100.0); dummyMap.put(CropType.OILCROPS, 100.0);
dummyMap.put(CropType.STARCHY_ROOTS, 50.0); dummyMap.put(CropType.STARCHY_ROOTS, 250.0);
dummyMap.put(CropType.PULSES, 60.0); dummyMap.put(CropType.PULSES, 260.0);
dummyMap.put(CropType.VEGETABLES, 20.0); dummyMap.put(CropType.VEGETABLES, 20.0);
dummyMap.put(CropType.MEAT_OR_PASTURE, 80.0); dummyMap.put(CropType.MEAT_OR_PASTURE, 800.0);
dummyMap.put(CropType.TREENUTS, 5.0); dummyMap.put(CropType.TREENUTS, 50.0);
return dummyMap; return dummyMap;
} }
@Override @Override
public Map<CropType, Double> getRefYield() { public Map<CropType, Double> getYield00() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>(); Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
dummyMap.put(CropType.CEREALS, 1.0); dummyMap.put(CropType.CEREALS, 1.0);
dummyMap.put(CropType.FRUIT, 5.0); dummyMap.put(CropType.FRUIT, 1.0);
dummyMap.put(CropType.OILCROPS, 1.0);
dummyMap.put(CropType.STARCHY_ROOTS, 1.0);
dummyMap.put(CropType.PULSES, 1.0);
dummyMap.put(CropType.VEGETABLES, 1.0);
dummyMap.put(CropType.MEAT_OR_PASTURE, 1.0);
dummyMap.put(CropType.TREENUTS, 1.0);
return dummyMap;
}
public Map<CropType, Double> getYield10() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
dummyMap.put(CropType.CEREALS, 3.0);
dummyMap.put(CropType.FRUIT, 3.0);
dummyMap.put(CropType.OILCROPS, 3.0);
dummyMap.put(CropType.STARCHY_ROOTS, 3.0);
dummyMap.put(CropType.PULSES, 3.0);
dummyMap.put(CropType.VEGETABLES, 3.0);
dummyMap.put(CropType.MEAT_OR_PASTURE, 3.0);
dummyMap.put(CropType.TREENUTS, 3.0);
return dummyMap;
}
public Map<CropType, Double> getYield01() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
dummyMap.put(CropType.CEREALS, 2.0);
dummyMap.put(CropType.FRUIT, 2.0);
dummyMap.put(CropType.OILCROPS, 2.0); dummyMap.put(CropType.OILCROPS, 2.0);
dummyMap.put(CropType.STARCHY_ROOTS, 5.0); dummyMap.put(CropType.STARCHY_ROOTS, 2.0);
dummyMap.put(CropType.PULSES, 20.0); dummyMap.put(CropType.PULSES, 2.0);
dummyMap.put(CropType.VEGETABLES, 20.0); dummyMap.put(CropType.VEGETABLES, 2.0);
dummyMap.put(CropType.MEAT_OR_PASTURE, 8.0); dummyMap.put(CropType.MEAT_OR_PASTURE, 2.0);
dummyMap.put(CropType.TREENUTS, 5.0); dummyMap.put(CropType.TREENUTS, 2.0);
return dummyMap;
}
public Map<CropType, Double> getYield11() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
dummyMap.put(CropType.CEREALS, 4.0);
dummyMap.put(CropType.FRUIT, 4.0);
dummyMap.put(CropType.OILCROPS, 4.0);
dummyMap.put(CropType.STARCHY_ROOTS, 4.0);
dummyMap.put(CropType.PULSES, 4.0);
dummyMap.put(CropType.VEGETABLES, 4.0);
dummyMap.put(CropType.MEAT_OR_PASTURE, 4.0);
dummyMap.put(CropType.TREENUTS, 4.0);
return dummyMap; return dummyMap;
} }
...@@ -133,4 +172,14 @@ public class GamsTest implements GamsInputData { ...@@ -133,4 +172,14 @@ public class GamsTest implements GamsInputData {
public double getMinFeedRate() { public double getMinFeedRate() {
return 0.15; return 0.15;
} }
@Override
public double getAlpha() {
return 0.4;
}
@Override
public double getBeta() {
return 0.5;
}
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment