From f089362a229bff4c0cb59bbb2e59b579c38a3003 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Wed, 3 Dec 2014 16:03:14 +0000
Subject: [PATCH] Multiple locations per country
---
GAMS/IntExtOpt.gms | 80 ++++++++++---------
.../country/gams/GamsLandUseOptimiser.java | 51 ++++++++----
src/ac/ed/lurg/country/gams/GamsTest.java | 4 +-
3 files changed, 80 insertions(+), 55 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 9be8ae43..7874b8a2 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -3,13 +3,14 @@
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 not_feed_crops(crop) / fruits, starchyRoots, treenuts, vegetables, pulses, pasture_or_meat /;
+ SET location / 1*2 /;
- PARAMETER previous_area(crop) areas for previous timestep in ha;
- PARAMETER demand(crop) yield in t per ha;
- PARAMETER yield00(crop) yield in t per ha;
- PARAMETER yield10(crop) yield in t per ha;
- PARAMETER yield01(crop) yield in t per ha;
- PARAMETER yield11(crop) yield in t per ha;
+ PARAMETER previous_area(crop, location) areas for previous timestep in ha;
+ PARAMETER yield00(crop, location) yield in t per ha;
+ PARAMETER yield10(crop, location) yield in t per ha;
+ PARAMETER yield01(crop, location) yield in t per ha;
+ PARAMETER yield11(crop, location) yield in t per ha;
+ PARAMETER demand(crop) in t;
PARAMETER world_input_energy(crop) average input energy from world exports used to determine if we should import or export;
PARAMETER maxNetImport(crop) maximum net import for each crop based on world market;
PARAMETER minNetImport(crop) minimum net import for each crop based on world market;
@@ -23,9 +24,9 @@
$gdxin %gdxincname%
$load previous_area, demand, yield00, yield10, yield01, yield11, world_input_energy, maxNetImport, minNetImport, meatEfficency, maxLandUseChange, tradeBarrier, landChangeEnergy, minFeedRate
$gdxin
-
- DISPLAY landChangeEnergy
-
+
+DISPLAY yield00, yield10, yield01, yield11
+
PARAMETER feedEnergy(crop) energy from feed in MJ per t
/ cereals 13.7
oilcrops 18.6
@@ -37,58 +38,59 @@ $gdxin
pasture_or_meat 1.0 / ;
VARIABLES
- area(crop) total area for each crop - ha
+ area(crop, location) total area for each crop - ha
+ 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) other intensity for each crop - unit less
feedAmount(crop) amount of feed for each crop - t
- fertI(crop) fertilizer intensity for each crop - factor between 0 and 1
- irrigI(crop) irrigation intensity for each crop - factor between 0 and 1
- otherIntensity(crop) other intensity for each crop - unit less
netImportAmount(crop) export of crop - t
- yield(crop) yield per area for each crop - t per ha
- unit_energy(crop) energy per area for each crop - energy
+ yield(crop, location) yield per area for each crop - t per ha
+ unit_energy(crop, location) energy per area for each crop - energy
energy total input energy - energy;
POSITIVE VARIABLE area, fertI, irrigI, otherI, feedAmount;
EQUATIONS
- UNIT_ENERGY_EQ(crop) energy per area
- YIELD_EQ(crop) yield given chosen intensity
+ UNIT_ENERGY_EQ(crop, location) energy per area
+ YIELD_EQ(crop, location) yield given chosen intensity
CROP_DEMAND_CONSTRAINT(crop_less_pasture) satisfy demand for crop
MEAT_DEMAND_CONSTRAINT satisfy demand for crop
- MAX_FERT_INTENSITY_CONSTRAINT(crop) constraint on maximum fertilizer intensity
- MAX_IRRIG_INTENSITY_CONSTRAINT(crop) constraint on maximum irrigation intensity
- CROPLAND_CHANGE_CONSTRAINT constraint on the rate of land use change
- PASTURE_CHANGE_CONSTRAINT constraint on the rate of land use change
- AGRI_LAND_CHANGE_CONSTRAINT constraint expansion of agricultural land
+ MAX_FERT_INTENSITY_CONSTRAINT(crop, location) constraint on maximum fertilizer intensity
+ MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity
+ CROPLAND_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
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
MIN_NET_IMPORT_CONSTRAINT(crop) constraint on min net imports
MIN_FEED_CONSTRAINT constraint on min feed rate
ENERGY_EQ total energy objective function;
- UNIT_ENERGY_EQ(crop) .. unit_energy(crop) =E= (fertI(crop)+irrigI(crop)+otherIntensity(crop)) + sqr(fertI(crop)+irrigI(crop)+otherIntensity(crop)) / 8;
+ UNIT_ENERGY_EQ(crop, location) .. unit_energy(crop, location) =E= fertI(crop, location)+irrigI(crop, location)+otherIntensity(crop, location) +
+ sqr(fertI(crop, location)+irrigI(crop, location)+otherIntensity(crop, location)) / 2;
- YIELD_EQ(crop) .. yield(crop) =E= (
- yield00(crop)*(1-fertI(crop))*(1-irrigI(crop)) +
- yield10(crop)*fertI(crop)*(1-irrigI(crop)) +
- yield01(crop)*(1-fertI(crop))*irrigI(crop) +
- yield11(crop)*fertI(crop)*irrigI(crop)
- ) * otherIntensity(crop);
+ YIELD_EQ(crop, location) .. yield(crop, location) =E= (
+ yield00(crop, location)*(1-fertI(crop, location))*(1-irrigI(crop, location)) +
+ yield10(crop, location)*fertI(crop, location)*(1-irrigI(crop, location)) +
+ yield01(crop, location)*(1-fertI(crop, location))*irrigI(crop, location) +
+ yield11(crop, location)*fertI(crop, location)*irrigI(crop, location)
+ ) * otherIntensity(crop, location);
- CROP_DEMAND_CONSTRAINT(crop_less_pasture) .. area(crop_less_pasture) * yield(crop_less_pasture) - 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=
demand(crop_less_pasture) - netImportAmount(crop_less_pasture);
- MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) + feedDM('pasture_or_meat') * area('pasture_or_meat') * yield('pasture_or_meat') =G=
+ MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) + feedDM('pasture_or_meat') * sum(location, area('pasture_or_meat', location) * yield('pasture_or_meat', location)) =G=
demand('pasture_or_meat') - netImportAmount('pasture_or_meat');
- MAX_FERT_INTENSITY_CONSTRAINT(crop) .. fertI(crop) =L= 1;
+ MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
- MAX_IRRIG_INTENSITY_CONSTRAINT(crop) .. irrigI(crop) =L= 1;
+ MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
- CROPLAND_CHANGE_CONSTRAINT .. 1 - sum(crop_less_pasture, area(crop_less_pasture))/sum(crop_less_pasture, previous_area(crop_less_pasture)) =L= maxLandUseChange;
+ CROPLAND_CHANGE_CONSTRAINT(location) .. 1 - sum(crop_less_pasture, area(crop_less_pasture, location))/sum(crop_less_pasture, previous_area(crop_less_pasture, location)) =L= maxLandUseChange;
- PASTURE_CHANGE_CONSTRAINT .. 1 - area('pasture_or_meat')/previous_area('pasture_or_meat') =L= maxLandUseChange;
+ PASTURE_CHANGE_CONSTRAINT(location) .. 1 - area('pasture_or_meat', location)/previous_area('pasture_or_meat', location) =L= maxLandUseChange;
- AGRI_LAND_CHANGE_CONSTRAINT .. abs(sum(crop, area(crop)) / sum(crop, previous_area(crop)) - 1) =L= maxLandUseChange;
+ AGRI_LAND_CHANGE_CONSTRAINT(location) .. abs(sum(crop, area(crop, location)) / sum(crop, previous_area(crop, location)) - 1) =L= maxLandUseChange;
NON_FEED_CROP_CONSTRAINT(not_feed_crops) .. feedAmount(not_feed_crops) =E= 0;
@@ -98,13 +100,15 @@ $gdxin
MIN_FEED_CONSTRAINT .. sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= minFeedRate * (demand('pasture_or_meat') - netImportAmount('pasture_or_meat'));
- ENERGY_EQ .. energy =E= SUM(crop, area(crop)*unit_energy(crop))
- + (sum(crop, area(crop)) - sum(crop, previous_area(crop))) * landChangeEnergy
+ ENERGY_EQ .. energy =E= SUM((crop, location), area(crop, location)*unit_energy(crop, location))
+ + sum(location, sum(crop, area(crop, location)) - sum(crop, previous_area(crop, location)) * landChangeEnergy)
+ sum(crop, (netImportAmount(crop) * tradeBarrier) * world_input_energy(crop));
MODEL LAND_USE /ALL/ ;
SOLVE LAND_USE USING DNLP MINIMIZING energy;
+ DISPLAY yield00, yield10, yield01, yield11
+
Scalar ms 'model status', ss 'solve status';
ms=LAND_USE.modelstat;
ss=LAND_USE.solvestat;
diff --git a/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java
index 953336f2..bd83bc77 100644
--- a/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java
@@ -2,6 +2,7 @@ package ac.ed.lurg.country.gams;
import java.io.File;
import java.util.Map;
+import java.util.Vector;
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.country.CropAreas;
@@ -39,28 +40,33 @@ public class GamsLandUseOptimiser {
GAMSDatabase db = ws.addDatabase();
LogWriter.println("\nPrevious areas");
- GAMSParameter parm = db.addParameter("previous_area", 1, "the previous area for each land use");
- addItemMapParm(parm, inputData.getPreviousCropArea());
+ GAMSParameter parm = db.addParameter("previous_area", 2, "the previous area for each land use");
+ addLocationMapParm(parm, inputData.getPreviousCropArea(), "1");
+ addLocationMapParm(parm, inputData.getPreviousCropArea(), "2");
LogWriter.println("\nDemand");
parm = db.addParameter("demand", 1, "demand for crop");
addItemMapParm(parm, inputData.getProjectedDemand());
LogWriter.println("\nYield");
- parm = db.addParameter("yield00", 1, "ref yield t per ha");
- addItemMapParm(parm, inputData.getRefYield());
+ parm = db.addParameter("yield00", 2, "ref yield t per ha");
+ addLocationMapParm(parm, inputData.getRefYield(), "1");
+ addLocationMapParm(parm, inputData.getRefYield(), "2");
LogWriter.println("\nYield");
- parm = db.addParameter("yield10", 1, "ref yield t per ha");
- addItemMapParm(parm, inputData.getRefYield());
+ parm = db.addParameter("yield10", 2, "ref yield t per ha");
+ addLocationMapParm(parm, inputData.getRefYield(), "1");
+ addLocationMapParm(parm, inputData.getRefYield(), "2");
LogWriter.println("\nYield");
- parm = db.addParameter("yield01", 1, "ref yield t per ha");
- addItemMapParm(parm, inputData.getRefYield());
+ parm = db.addParameter("yield01", 2, "ref yield t per ha");
+ addLocationMapParm(parm, inputData.getRefYield(), "1");
+ addLocationMapParm(parm, inputData.getRefYield(), "2");
LogWriter.println("\nYield");
- parm = db.addParameter("yield11", 1, "ref yield t per ha");
- addItemMapParm(parm, inputData.getRefYield());
+ parm = db.addParameter("yield11", 2, "ref yield t per ha");
+ addLocationMapParm(parm, inputData.getRefYield(), "1");
+ addLocationMapParm(parm, inputData.getRefYield(), "2");
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");
@@ -102,16 +108,17 @@ public class GamsLandUseOptimiser {
for (GAMSVariableRecord rec : varAreas) {
String itemName = rec.getKeys()[0];
+ String locationName = rec.getKeys()[1];
area = rec.getLevel();
- fertIntensity = varFertIntensities.findRecord(itemName).getLevel();
- irrigIntensity = varIrrigIntensities.findRecord(itemName).getLevel();
- otherIntensity = varOtherIntensities.findRecord(itemName).getLevel();
+ fertIntensity = varFertIntensities.findRecord(itemName, locationName).getLevel();
+ irrigIntensity = varIrrigIntensities.findRecord(itemName, locationName).getLevel();
+ otherIntensity = varOtherIntensities.findRecord(itemName, locationName).getLevel();
feedAmount = varFeedAmount.findRecord(itemName).getLevel();
netImports = varNetImports.findRecord(itemName).getLevel();
totalArea += rec.getLevel();
- LogWriter.println(String.format("%15s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f,\tfeedAmount= %.0f,\tnetImports= %.0f",
- itemName, area, fertIntensity, irrigIntensity, otherIntensity, feedAmount, netImports));
+ LogWriter.println(String.format("%15s, l%s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f,\tfeedAmount= %.0f,\tnetImports= %.0f",
+ itemName, locationName, area, fertIntensity, irrigIntensity, otherIntensity, feedAmount, netImports));
output.addCropArea(CropType.getForGamsName(itemName), area, otherIntensity, feedAmount, netImports);
}
@@ -132,6 +139,20 @@ public class GamsLandUseOptimiser {
}
}
+ private void addLocationMapParm(GAMSParameter parm, Map<CropType, Double> itemMap, String locationId) {
+ for (Map.Entry<CropType, Double> entry : itemMap.entrySet()) {
+ LogWriter.println(String.format(" %15s,\t %.1f", entry.getKey().getGamsName(), entry.getValue()));
+ Vector<String> v = new Vector<String>();
+ v.add(entry.getKey().getGamsName());
+ v.add(locationId);
+ if (locationId.equals("2") & entry.getKey().getGamsName().equals("cereals"))
+ parm.addRecord(v).setValue(0.99);
+ else
+ parm.addRecord(v).setValue(entry.getValue());
+ }
+ }
+
+
@SuppressWarnings("unused")
private void cleanup(String directory) {
File directoryToDelete = new File(directory);
diff --git a/src/ac/ed/lurg/country/gams/GamsTest.java b/src/ac/ed/lurg/country/gams/GamsTest.java
index b780158c..94b5901a 100644
--- a/src/ac/ed/lurg/country/gams/GamsTest.java
+++ b/src/ac/ed/lurg/country/gams/GamsTest.java
@@ -23,13 +23,13 @@ public class GamsTest implements GamsInputData {
@Override
public Map<CropType, Double> getProjectedDemand() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
- dummyMap.put(CropType.CEREALS, 390.0);
+ dummyMap.put(CropType.CEREALS, 290.0);
dummyMap.put(CropType.FRUIT, 50.0);
dummyMap.put(CropType.OILCROPS, 100.0);
dummyMap.put(CropType.STARCHY_ROOTS, 50.0);
dummyMap.put(CropType.PULSES, 60.0);
dummyMap.put(CropType.VEGETABLES, 20.0);
- dummyMap.put(CropType.MEAT_OR_PASTURE, 800.0);
+ dummyMap.put(CropType.MEAT_OR_PASTURE, 80.0);
dummyMap.put(CropType.TREENUTS, 5.0);
return dummyMap;
}
--
GitLab