diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index ae12d910cab85da37eae5cddcb387c9ac246fecd..5cef5c686c373e9b4240852b8b79f45bf45f48fc 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -32,7 +32,7 @@
SCALAR fitToPreviousAreas controlling parameter 0 means dont fit and 1 is to calibrate the cropAdj in the calibration year
-*$gdxin "/Users/peteralexander/Documents/R_Workspace/temp/GamsTest/_gams_java_189513040.gdx"
+*$gdxin "/Users/peteralexander/Documents/R_Workspace/temp/GamsTest/_gams_java_gdb1.gdx"
$gdxin %gdxincname%
$load location, suitableLandArea, previousArea, demand, landChangeEnergy
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
@@ -75,7 +75,7 @@ $gdxin
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
+ agriLandExpansion(location) addition agricultural land needed as it must be positive it deliberately does not account for abandonment
energy total input energy - energy;
POSITIVE VARIABLE area, fertI, irrigI, feedAmount, agriLandExpansion;
@@ -94,11 +94,11 @@ $gdxin
MAX_FERT_INTENSITY_CONSTRAINT(crop, location) constraint on maximum fertilizer intensity
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity
TOTAL_LAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
- CROPLAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
- PASTURE_MAX_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
- PASTURE_MIN_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
- AGRI_LAND_INCREASE_CONSTRAINT(location) constraint expansion of agricultural land
- AGRI_LAND_DECREASE_CONSTRAINT(location) constraint expansion of agricultural land
+* CROPLAND_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
+* PASTURE_MAX_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
+* PASTURE_MIN_CHANGE_CONSTRAINT(location) constraint on the rate of land use change
+* 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_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
@@ -142,13 +142,11 @@ $gdxin
TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(crop, area(crop, location));
- CROPLAND_CHANGE_CONSTRAINT(location) .. sum(crop_less_pasture, area(crop_less_pasture, location)) =L= (1 + maxLandUseChange) * sum(crop_less_pasture, previousArea(crop_less_pasture, location));
-
- PASTURE_MAX_CHANGE_CONSTRAINT(location) .. area('pasture', location) =L= (1 + maxLandUseChange) * previousArea('pasture', location);
- PASTURE_MIN_CHANGE_CONSTRAINT(location) .. area('pasture', location) =G= (1 - maxLandUseChange) * previousArea('pasture', location);
-
- AGRI_LAND_INCREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) =L= (1 + maxLandUseChange) * sum(crop, previousArea(crop, location));
- AGRI_LAND_DECREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) =G= (1 - maxLandUseChange) * sum(crop, previousArea(crop, location));
+* CROPLAND_CHANGE_CONSTRAINT(location) .. sum(crop_less_pasture, area(crop_less_pasture, location)) =L= (1 + maxLandUseChange) * sum(crop_less_pasture, previousArea(crop_less_pasture, location));
+* PASTURE_MAX_CHANGE_CONSTRAINT(location) .. area('pasture', location) =L= (1 + maxLandUseChange) * previousArea('pasture', location);
+* PASTURE_MIN_CHANGE_CONSTRAINT(location) .. area('pasture', location) =G= (1 - maxLandUseChange) * previousArea('pasture', location);
+* AGRI_LAND_INCREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) =L= (1 + maxLandUseChange) * sum(crop, previousArea(crop, location));
+* AGRI_LAND_DECREASE_CONSTRAINT(location) .. sum(crop, area(crop, location)) =G= (1 - maxLandUseChange) * sum(crop, previousArea(crop, location));
NON_FEED_CROP_CONSTRAINT(not_feed_crop) .. feedAmount(not_feed_crop) =E= 0;
@@ -185,7 +183,7 @@ $gdxin
basePasture = sum(location, previousArea("pasture", location));
baseCropland = sum(location, sum(crop_less_pasture, previousArea(crop_less_pasture, location)));
- while((count le 10 and not (abs(basePasture-newPasture) le 0.2 and abs(baseCropland-newCropland) le 0.2)),
+ while((count le 10 and not (abs(basePasture-newPasture) le 0.1 and abs(baseCropland-newCropland) le 0.1)),
SOLVE LAND_USE USING NLP MINIMIZING energy;
newPasture = sum(location, area.l("pasture", location));
@@ -213,11 +211,11 @@ $gdxin
totalImportEnergy(import_crop) = netImportAmount.l(import_crop) * worldInputEnergy(import_crop);
feedEnergy(feed_crop)$(netImportAmount.l(feed_crop) gt 0) = (totalProdCost(feed_crop) + totalImportEnergy(feed_crop)) * feedAmount.l(feed_crop) / (totalProd(feed_crop) + netImportAmount.l(feed_crop));
- feedEnergy(feed_crop)$(netImportAmount.l(feed_crop) le 0) = totalProdCost(feed_crop) * feedAmount.l(feed_crop) / totalProd(feed_crop);
+ feedEnergy(feed_crop)$(netImportAmount.l(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 totalProd, totalImportEnergy, feedEnergy, totalProdCost
+ display totalProdCost, totalProd, totalImportEnergy, feedEnergy, totalProdCost
Scalar ms 'model status', ss 'solve status';
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index f94e9f19a0b84095be4e3ab6804a8e753136befe..711672f0ad517b3af43b912b1dcb5dc0e864951f 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -2,6 +2,7 @@ package ac.ed.lurg;
import java.io.File;
import java.util.ArrayList;
+import java.util.Arrays;
import java.util.Collection;
import java.util.HashMap;
import java.util.HashSet;
@@ -92,7 +93,7 @@ public class ModelMain {
LogWriter.println("Country " + ca.getCountry());
Collection<RasterKey> countryKeys = countryToKeysMap.get(ca.getCountry());
YieldRaster countryYieldSurfaces = yieldSurfaces.getSubsetRasterForKeys(countryKeys);
-
+
// do the optimization
GamsRasterOutput result = ca.determineProduction(timestep, countryYieldSurfaces, prevWorldInputCost);
@@ -149,11 +150,7 @@ public class ModelMain {
@Override
public int convertToPixelValue(double value) {
- if (value == 4.0)
- return 10;
- else return 1;
-
- //return (int) (value * 3 + 1);
+ return (int) (value * 3 + 1);
}
}.writeOutput(ModelConfig.WRITE_JPEG_IMAGES);
@@ -172,26 +169,31 @@ public class ModelMain {
@Override
public int convertToPixelValue(double value) {
- return (int) (value * 20);
+ return (int) (value * 20) + 1;
}
}.writeOutput(ModelConfig.WRITE_JPEG_IMAGES);
- new RasterOutputer<AreasItem>(cropAreaRaster, "wheatArea" + timestep) {
+ outputAreas(timestep, cropAreaRaster, CropType.WHEAT);
+ outputAreas(timestep, cropAreaRaster, CropType.MAIZE);
+ outputAreas(timestep, cropAreaRaster, CropType.PASTURE);
+ }
+
+ private void outputAreas(int timestep, RasterSet<AreasItem> cropAreaRaster, final CropType crop) {
+ new RasterOutputer<AreasItem>(cropAreaRaster, crop.getGamsName() + "Area" + timestep) {
@Override
public Double getValue(RasterKey location) {
AreasItem area = results.get(location);
if (area == null)
return null;
- return area.getCropArea(CropType.WHEAT);
+ return area.getCropArea(crop);
}
@Override
public int convertToPixelValue(double value) {
- return (int) (value * 100);
+ return value > 0 ? 10 : 1;
}
}.writeOutput(ModelConfig.WRITE_JPEG_IMAGES);
-
}
public Map<Country, List<RasterKey>> createCountryToKeysMap() {
@@ -231,7 +233,7 @@ public class ModelMain {
RasterSet<IrrigationCostItem> allIrrigationCosts = getIrrigationCosts();
Map<Country, Map<CropType, CropUsageData>> cropUsageDataMap = CropUsageData.readCommodityData();
- // HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("French Polynesia", "Cabo Verde", "Samoa", "Saint Vincent and the Grenadines"));
+ HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("Sierra Leone", "Togo", "United Arab Emirates")); //"French Polynesia", "Cabo Verde", "Samoa", "Saint Vincent and the Grenadines"));
for (Map.Entry<Country, List<RasterKey>> entry : countryToKeysMap.entrySet()) {
Country country = entry.getKey();
@@ -239,11 +241,11 @@ public class ModelMain {
// DEBUG code
- if (!(country.getCountryName().equals("United States of America") )) { //|| country.getCountryName().equals("China")
+ if (!(country.getCountryName().equals("United States of America") || country.getCountryName().equals("China"))) { //|| country.getCountryName().equals("China")
continue;
}
- if (demandManager.getPopulation(country, 2010) < 5) { // countryExclusionList.contains(country.getCountryName()
+ if (demandManager.getPopulation(country, 2010) < 5 || countryExclusionList.contains(country.getCountryName())) {
LogWriter.printlnError("Skipping " + country);
continue;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 2699ed0f2d87795196457bb905fd0496a7b3ab74..c04453ef0b01528509d50dd426067e4c674df151 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -36,9 +36,9 @@ public class GamsLocationOptimiser {
// oilcrops 0.90
// soy 0.89
// potatoes 0.21
-
+
private static final boolean DEBUG = true;
-
+
private GamsLocationInput inputData;
public GamsLocationOptimiser(GamsLocationInput inputData) {
@@ -61,14 +61,14 @@ public class GamsLocationOptimiser {
GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.GAMS_MODEL);
GAMSOptions opt = ws.addOptions();
-
+
opt.setAllModelTypes("conopt");
-
+
opt.defines("gdxincname", inDB.getName());
long startTime = System.currentTimeMillis();
gamsJob.run(opt, inDB);
- // cleanup(ws.workingDirectory());
+ // cleanup(ws.workingDirectory());
LogWriter.println("Took " + (System.currentTimeMillis() - startTime) + " ms to run");
return handleResults(gamsJob.OutDB());
@@ -87,7 +87,7 @@ public class GamsLocationOptimiser {
GAMSParameter landP = inDB.addParameter("suitableLandArea", 1);
double totalAgriLand = 0;
-
+
for (Map.Entry<Integer, ? extends AreasItem> entry : inputData.getPreviousAreas().entrySet()) {
Integer locationId = entry.getKey();
String locString = Integer.toString(locationId);
@@ -95,7 +95,8 @@ public class GamsLocationOptimiser {
for (CropType cropType : CropType.values()) {
double d = areasItem.getCropArea(cropType);
- if (DEBUG) if (d != 0) LogWriter.println(String.format(" %d %15s,\t %.1f", locationId, cropType.getGamsName(), d));
+ if (DEBUG && d != 0) LogWriter.println(String.format(" %d %15s,\t %.1f", locationId, cropType.getGamsName(), d));
+
Vector<String> v = new Vector<String>();
v.add(cropType.getGamsName());
v.add(locString);
@@ -164,16 +165,16 @@ public class GamsLocationOptimiser {
addScalar(inDB.addParameter("meatEfficency", 0), countryInput.getMeatEfficiency());
if (DEBUG) LogWriter.println("\nMeatEfficiency: " + countryInput.getMeatEfficiency());
-
+
addScalar(inDB.addParameter("maxLandUseChange", 0), countryInput.getMaxLandUseChange());
if (DEBUG) LogWriter.println("\nMaxLandUseChange: " + countryInput.getMaxLandUseChange());
-
+
addScalar(inDB.addParameter("landChangeEnergy", 0), countryInput.getLandChangeEnergy());
if (DEBUG) LogWriter.println("\nLandChangeEnergy: " + countryInput.getLandChangeEnergy());
-
+
addScalar(inDB.addParameter("minFeedRate", 0), countryInput.getMinFeedRate());
if (DEBUG) LogWriter.println("\nMinFeedRate: " + countryInput.getMinFeedRate());
-
+
double fitToPreviousAreas = countryInput.isCalibrateToObserved() ? 1 : 0;
if (DEBUG) LogWriter.println("\nfitToPreviousAreas: " + fitToPreviousAreas);
addScalar(inDB.addParameter("fitToPreviousAreas", 0), fitToPreviousAreas);
@@ -191,7 +192,7 @@ public class GamsLocationOptimiser {
GAMSVariable varAreas = outDB.getVariable("area");
GAMSVariable varFertIntensities = outDB.getVariable("fertI");
GAMSVariable varIrrigIntensities = outDB.getVariable("irrigI");
-// GAMSVariable varOtherIntensities = outDB.getVariable("otherIntensity");
+ // GAMSVariable varOtherIntensities = outDB.getVariable("otherIntensity");
GAMSVariable varFeedAmount = outDB.getVariable("feedAmount");
GAMSVariable varNetImports = outDB.getVariable("netImportAmount");
GAMSVariable varYields = outDB.getVariable("yield");
@@ -199,7 +200,7 @@ public class GamsLocationOptimiser {
GAMSParameter parmCropAdj = outDB.getParameter("cropAdj");
GAMSParameter parmProd = outDB.getParameter("totalProd");
GAMSParameter parmProdCost = outDB.getParameter("totalProdCost");
-
+
double totalArea = 0;
double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, feedAmount, netImport, yield, unitEnergy, cropAdj, prod, prodCost;
@@ -215,7 +216,7 @@ public class GamsLocationOptimiser {
area = rec.getLevel();
fertIntensity = varFertIntensities.findRecord(itemName, locationName).getLevel();
irrigIntensity = varIrrigIntensities.findRecord(itemName, locationName).getLevel();
-// otherIntensity = varOtherIntensities.findRecord(itemName, locationName).getLevel();
+ // otherIntensity = varOtherIntensities.findRecord(itemName, locationName).getLevel();
yield = varYields.findRecord(itemName, locationName).getLevel();
unitEnergy = varUnitEnergies.findRecord(itemName, locationName).getLevel();
@@ -228,34 +229,35 @@ public class GamsLocationOptimiser {
cropAdj = getParmValue(parmCropAdj, itemName);
prod = getParmValue(parmProd, itemName);
prodCost = getParmValue(parmProdCost, itemName);
-
+
cropUsageData.put(cropType, new CropUsageData(feedAmount, netImport, prod, prodCost));
cropAdjs.put(cropType, cropAdj);
if (DEBUG) LogWriter.println(String.format("\n%s:\tfeedAmount= %.1f,\tnetImports= %.3f,\tprod= %.3f,\tprodCost= %.3f,\tcropAdj= %.3f", itemName, feedAmount, netImport, prod, prodCost, cropAdj));
}
- if (area > 0) {
+ if (area > 0) {
if (DEBUG) LogWriter.println(String.format("\t location %s, %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f",
- locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity));
-
+ locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity));
+
IntensitiesItem intensityItem = intensities.get(locId);
if (intensityItem == null) {
intensityItem = new IntensitiesItem();
intensities.put(locId, intensityItem);
}
intensityItem.setIntensity(cropType, new Intensity(fertIntensity, irrigIntensity, otherIntensity, yield, unitEnergy));
+ }
- AreasItem areaItem = cropAreas.get(locId);
- if (areaItem == null) {
- areaItem = new AreasItem();
- cropAreas.put(locId, areaItem);
- }
- areaItem.setCropArea(cropType, area);
- totalArea += area;
+ AreasItem areaItem = cropAreas.get(locId);
+ if (areaItem == null) {
+ areaItem = new AreasItem();
+ cropAreas.put(locId, areaItem);
}
+ areaItem.setCropArea(cropType, area);
+ totalArea += area;
+
prevCropType = cropType;
}
-
+
netImport = varNetImports.findRecord(CropType.MEAT.getGamsName()).getLevel();
prod = parmProd.findRecord(CropType.MEAT.getGamsName()).getValue();
prodCost = parmProdCost.findRecord(CropType.MEAT.getGamsName()).getValue();
@@ -264,11 +266,11 @@ public class GamsLocationOptimiser {
LogWriter.println(String.format("\n%s:\t\t\t\t\tnetImports= %.3f,\tprod= %.3f,\tprodCost= %.3f", CropType.MEAT.getGamsName(), netImport, prod, prodCost));
LogWriter.println(String.format("\nTotal area= %.1f", totalArea));
}
-
+
GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, cropUsageData, cropAdjs);
return results;
}
-
+
private double getParmValue(GAMSParameter aParm, String itemName) {
try {
double d = aParm.findRecord(itemName).getValue();
@@ -292,7 +294,7 @@ public class GamsLocationOptimiser {
parm.addRecord(entry.getKey().getGamsName()).setValue(d);
}
}
-
+
private void addCommodityMapParm(GAMSParameter parm, Map<CommodityType, Double> itemMap) {
for (Map.Entry<CommodityType, Double> entry : itemMap.entrySet()) {
double d = entry.getValue();