From a4bfd5551452334449fa2f234eda81a5d61deebe Mon Sep 17 00:00:00 2001 From: Peter Alexander <peter@blackhillock.co.uk> Date: Mon, 15 Dec 2014 12:40:28 +0000 Subject: [PATCH] Optimising from / to raster set seems to be working --- GAMS/IntExtOpt.gms | 10 ++- src/ac/ed/lurg/ModelConfig.java | 6 +- src/ac/ed/lurg/country/CountryAgent.java | 2 +- .../country/gams/GamsLandUseOptimiser.java | 23 ++++-- .../country/gams/GamsRasterOptimiser.java | 49 +++++++++--- .../ed/lurg/country/gams/GamsRasterTest.java | 60 ++++++++++++++ src/ac/ed/lurg/country/gams/GamsTest.java | 79 +++++++++++-------- src/ac/ed/lurg/landuse/AreasItem.java | 4 +- .../ed/lurg/yield/AveragingYieldResponse.java | 10 ++- 9 files changed, 179 insertions(+), 64 deletions(-) create mode 100644 src/ac/ed/lurg/country/gams/GamsRasterTest.java diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms index 64029a6a..fb7c287b 100644 --- a/GAMS/IntExtOpt.gms +++ b/GAMS/IntExtOpt.gms @@ -3,8 +3,8 @@ 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*5 /; - + + SET location; PARAMETER landArea(location) areas of land in ha; PARAMETER previous_area(crop, location) areas for previous timestep in ha; PARAMETER yieldNone(crop, location) yield in t per ha; @@ -23,13 +23,15 @@ SCALAR tradeBarrier trade barrier which adjust energy cost of imports; SCALAR landChangeEnergy energy required to add ha of agricultural land; SCALAR minFeedRate minimum rate of feed for producing animal products; - + $gdxin %gdxincname% -$load landArea, previous_area, demand, yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, fertParam, irrigParam, world_input_energy, maxNetImport, minNetImport, meatEfficency, maxLandUseChange, tradeBarrier, landChangeEnergy, minFeedRate +$load location, landArea, previous_area, demand, yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, fertParam, irrigParam, world_input_energy, maxNetImport, minNetImport, meatEfficency, maxLandUseChange, tradeBarrier, landChangeEnergy, minFeedRate $gdxin SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 / + display previous_area + PARAMETER feedEnergy(crop) energy from feed in MJ per t / cereals 13.7 oilcrops 18.6 diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java index 1894eb05..8995025f 100644 --- a/src/ac/ed/lurg/ModelConfig.java +++ b/src/ac/ed/lurg/ModelConfig.java @@ -78,7 +78,9 @@ public class ModelConfig { String v = configFile.getProperty(prop); return v==null ? defaultBoolean : Boolean.valueOf(v); } - + + public static final boolean DEBUG = true; + public static final boolean SUPPRESS_STD_OUTPUT = getBooleanProperty("SUPPRESS_STD_OUTPUT", Boolean.FALSE); public static final String BASE_DIR = getProperty("OUTPUT_DIR", "/Users/peteralexander/Documents/R_Workspace"); @@ -101,6 +103,4 @@ public class ModelConfig { public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 1); public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2010); - public static final int NUM_LOCATIONS_PER_COUNTRY = 5; - } \ No newline at end of file diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java index c5b97f39..5291b92a 100644 --- a/src/ac/ed/lurg/country/CountryAgent.java +++ b/src/ac/ed/lurg/country/CountryAgent.java @@ -98,7 +98,7 @@ public class CountryAgent { Map<Integer, AreasItem> returnMap = new HashMap<Integer, AreasItem>(); - for (int i= 1; i<=ModelConfig.NUM_LOCATIONS_PER_COUNTRY; i++) + for (int i= 1; i<=5; i++) returnMap.put(i, previousCropAreas); return returnMap; diff --git a/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java index a1759eba..d1801cf2 100644 --- a/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java @@ -20,6 +20,7 @@ import com.gams.api.GAMSGlobals; import com.gams.api.GAMSJob; import com.gams.api.GAMSOptions; import com.gams.api.GAMSParameter; +import com.gams.api.GAMSSet; import com.gams.api.GAMSVariable; import com.gams.api.GAMSVariableRecord; import com.gams.api.GAMSWorkspace; @@ -59,7 +60,13 @@ public class GamsLandUseOptimiser { } private void setupInDB(GAMSDatabase inDB) { - + LogWriter.println("\nLocation set"); + GAMSSet locationSet = inDB.addSet("location", 1); + for (Integer locId : inputData.getPreviousAreas().keySet()) { + LogWriter.println(" " + locId); + locationSet.addRecord(locId.toString()); + } + LogWriter.println("\nPrevious areas"); addAreaParms(inDB, inputData.getPreviousAreas()); @@ -135,7 +142,8 @@ public class GamsLandUseOptimiser { Map<Integer, AreasItem> cropAreas = new HashMap<Integer, AreasItem>(); Map<CropType, Double> feedAmounts = new HashMap<CropType, Double>(); Map<CropType, Double> netImports = new HashMap<CropType, Double>(); - + CropType prevCropType = null; + for (GAMSVariableRecord rec : varAreas) { String itemName = rec.getKeys()[0]; String locationName = rec.getKeys()[1]; @@ -146,7 +154,7 @@ public class GamsLandUseOptimiser { int locId = Integer.parseInt(locationName); CropType cropType = CropType.getForGamsName(itemName); - if (locId == 1) { + if (!cropType.equals(prevCropType)) { feedAmount = varFeedAmount.findRecord(itemName).getLevel(); netImport = varNetImports.findRecord(itemName).getLevel(); feedAmounts.put(cropType, feedAmount); @@ -155,8 +163,8 @@ public class GamsLandUseOptimiser { } if (area > 0) { - LogWriter.println(String.format("\t location %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f", - locationName, area, fertIntensity, irrigIntensity, otherIntensity)); + LogWriter.println(String.format("\t location %s, %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f", + locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity)); IntensitiesItem intensityItem = intensities.get(locId); if (intensityItem == null) { @@ -171,10 +179,9 @@ public class GamsLandUseOptimiser { cropAreas.put(locId, areaItem); } areaItem.setCropArea(cropType, area); - totalArea += area; } - + prevCropType = cropType; } LogWriter.println(String.format("\nTotal area= %.1f", totalArea)); //cleanup(ws.workingDirectory()); @@ -206,7 +213,7 @@ public class GamsLandUseOptimiser { for (CropType cropType : CropType.values()) { double d = areasItem.getCropArea(cropType); - LogWriter.println(String.format(" %15s,\t %.1f", cropType.getGamsName(), d)); + 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); diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java index 53354777..8c5ff657 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java @@ -9,6 +9,7 @@ import java.util.Map; import java.util.Map.Entry; import java.util.Set; +import ac.ed.lurg.ModelConfig; import ac.ed.lurg.landuse.AreasItem; import ac.ed.lurg.landuse.IntensitiesItem; import ac.ed.lurg.types.CropToDouble; @@ -61,7 +62,8 @@ public class GamsRasterOptimiser { Set<RasterKey> keys = mapping.get(locId); CropToDouble cropChangeTotals = newAreaAggItem.getCropChanges(prevAreaAggItem); - + LogWriter.println("Processing location id " + locId); + double shortfall = allocAllAreasForAgg(newAreaRaster, prevAreaRaster, keys, cropChangeTotals); if (shortfall > 0) LogWriter.printlnError("This should never happen, due to GAMS constraint. Not able to incorporate all changes, as not enough forest or natural areas left: " + shortfall); @@ -77,6 +79,7 @@ public class GamsRasterOptimiser { CropToDouble avgCropChange = cropChange.scale(1.0/keys.size()); for (RasterKey key : keys) { + LogWriter.println(" Processing raster key " + key); AreasItem prevAreasRasterItem = prevAreaRaster.get(key); AreasItem newAreasRasterItem = new AreasItem(); @@ -99,6 +102,7 @@ public class GamsRasterOptimiser { return totalShortfall; } + /** returns the net area of crop change that not been able to be allocated */ private double allocRasterArea(AreasItem newAreasRasterItem, AreasItem prevAreasRasterItem, CropToDouble cropChange) { double netAddedCrop = cropChange.getTotal(); @@ -107,19 +111,25 @@ public class GamsRasterOptimiser { double prevForest = prevAreasRasterItem.getLandCoverArea(LandDataType.FOREST); double prevNatural = prevAreasRasterItem.getLandCoverArea(LandDataType.OTHER_NATURAL); - if (prevForest < forestRemoved) { + double ratioOfChanges = 1; + + if (netAddedCrop > prevForest + prevNatural) { + ratioOfChanges = (prevForest + prevNatural) / netAddedCrop; + forestRemoved = prevForest; + naturalRemoved = prevNatural; + LogWriter.println("Not able to incorporate all changes for raster location: " + netAddedCrop); + } + else if (prevForest < forestRemoved) { forestRemoved = prevForest; naturalRemoved = netAddedCrop - forestRemoved; } - - if (prevNatural < naturalRemoved) { + else if (prevNatural < naturalRemoved) { naturalRemoved = prevNatural; + forestRemoved = netAddedCrop - naturalRemoved; + } + else { + // enough space to keep 50/50 allocation } - - double ratioOfChanges = (forestRemoved + naturalRemoved)/netAddedCrop; - double shortfall = netAddedCrop - forestRemoved - naturalRemoved; - if (shortfall > 0) - LogWriter.println("Not able to incorporate all changes for raster location: " + shortfall); newAreasRasterItem.setLandCoverArea(LandDataType.LAND, prevAreasRasterItem.getLandCoverArea(LandDataType.LAND)); newAreasRasterItem.setLandCoverArea(LandDataType.FOREST, prevForest - forestRemoved); @@ -127,11 +137,15 @@ public class GamsRasterOptimiser { for (CropType crop : CropType.values()) { double prevCropArea = prevAreasRasterItem.getCropArea(crop); - double additionalArea = cropChange.get(crop) * ratioOfChanges; + Double d = cropChange.get(crop); + + double additionalArea = (d == null) ? 0 : d * ratioOfChanges; + // LogWriter.println(String.format("%s, %.1f", crop, additionalArea)); + newAreasRasterItem.setCropArea(crop, prevCropArea + additionalArea); } - return shortfall; + return netAddedCrop - prevForest - prevNatural; } private RasterSet<IntensitiesItem> allocIntensities(GamsOutput gamsOutput) { @@ -180,11 +194,22 @@ public class GamsRasterOptimiser { for (YieldType yieldType : YieldType.values()) { avgYResp.setYield(yieldType, crop, yresp.getYield(yieldType, crop)); } - aggAreas.setCropArea(crop, aggAreas.getCropArea(crop) + cropAreas.getCropArea(crop)); + double areaSoFar = aggAreas.getCropArea(crop); + double areaThisTime = cropAreas.getCropArea(crop); + aggAreas.setCropArea(crop, areaSoFar + areaThisTime); } aggAreas.setLandCoverArea(LandDataType.LAND, aggAreas.getLandCoverArea(LandDataType.LAND) + cropAreas.getLandCoverArea(LandDataType.LAND)); } + + if (ModelConfig.DEBUG) { + for (Map.Entry<Integer, Set<RasterKey>> e : mapping.entrySet()) { + LogWriter.println(e.getKey() + " category has " + e.getValue().size() + " raster areas"); + for (RasterKey key : e.getValue()) { + LogWriter.println(String.format("\t%s: yields pasture=%.1f, cereal=%.1f", key, yieldRaster.get(key).getYieldMax(CropType.MEAT_OR_PASTURE), yieldRaster.get(key).getYieldMax(CropType.CEREALS))); + } + } + } return new GamsInput(aggregatedYield, aggregatedAreas, rasterInputData.getCountryInput()); } diff --git a/src/ac/ed/lurg/country/gams/GamsRasterTest.java b/src/ac/ed/lurg/country/gams/GamsRasterTest.java new file mode 100644 index 00000000..136a8be5 --- /dev/null +++ b/src/ac/ed/lurg/country/gams/GamsRasterTest.java @@ -0,0 +1,60 @@ +package ac.ed.lurg.country.gams; + +import ac.ed.lurg.landuse.AreasItem; +import ac.ed.lurg.yield.YieldRaster; +import ac.ed.lurg.yield.YieldResponsesItem; +import ac.sac.raster.RasterKey; +import ac.sac.raster.RasterSet; + +public class GamsRasterTest extends GamsTest { + + public static void main(String[] args) { + GamsRasterTest aGamsRasterTest = new GamsRasterTest(); + aGamsRasterTest.run(); + } + + private void run() { + GamsCountryInput countryLevelInputs = new GamsCountryInput(getProjectedDemand(), getWorldInputEnergy(), getMaxNetImport(), getMinNetImport()); + GamsRasterInput input = new GamsRasterInput(getYieldRaster(), getPreviousAreaRaster(), countryLevelInputs); + + GamsRasterOptimiser opti = new GamsRasterOptimiser(input); + GamsRasterOutput output = opti.run(); + } + + private static final int X_RANGE = 3; + private static final int Y_RANGE = 3; + + public YieldRaster getYieldRaster() { + YieldRaster testYieldRaster= new YieldRaster(null); + + for (int i = 0; i<X_RANGE; i++) { + for (int j = 0; j<Y_RANGE; j++) { + RasterKey key = new RasterKey(i,j); + YieldResponsesItem yresp; + if (i==2 && j==2) + yresp = getYieldRespItem(0.5*i + 0.3*j, 0.7); + else + yresp = getYieldRespItem(2.0/(i+1) + 0.3*j, 1.2); + testYieldRaster.put(key, yresp); + } + } + + return testYieldRaster; + } + + + public RasterSet<AreasItem> getPreviousAreaRaster() { + RasterSet<AreasItem> testAreaRaster= new RasterSet<AreasItem>(); + + for (int i = 0; i<X_RANGE; i++) { + for (int j = 0; j<Y_RANGE; j++) { + RasterKey key = new RasterKey(i,j); + AreasItem areaItem = getAreaItem(i * 20, j * 30); + testAreaRaster.put(key, areaItem); + } + } + + + return testAreaRaster; + } +} diff --git a/src/ac/ed/lurg/country/gams/GamsTest.java b/src/ac/ed/lurg/country/gams/GamsTest.java index 1b491ff0..e0464fa8 100644 --- a/src/ac/ed/lurg/country/gams/GamsTest.java +++ b/src/ac/ed/lurg/country/gams/GamsTest.java @@ -3,7 +3,6 @@ package ac.ed.lurg.country.gams; import java.util.HashMap; import java.util.Map; -import ac.ed.lurg.ModelConfig; import ac.ed.lurg.landuse.AreasItem; import ac.ed.lurg.types.CropType; import ac.ed.lurg.types.LandDataType; @@ -11,7 +10,8 @@ import ac.ed.lurg.types.YieldType; import ac.ed.lurg.yield.YieldResponsesItem; public class GamsTest { - + public static final int NUM_LOCATIONS_PER_COUNTRY = 5; + public static void main(String[] args) { GamsTest aGamsTest = new GamsTest(); aGamsTest.run(); @@ -25,7 +25,7 @@ public class GamsTest { opti.run(); } - public Map<CropType, Double> getProjectedDemand() { + Map<CropType, Double> getProjectedDemand() { Map<CropType, Double> dummyMap = new HashMap<CropType, Double>(); dummyMap.put(CropType.CEREALS, 300.0); dummyMap.put(CropType.FRUIT, 50.0); @@ -39,49 +39,64 @@ public class GamsTest { } - public Map<Integer, YieldResponsesItem> getYields() { - - YieldResponsesItem yresp = new YieldResponsesItem(); - - for (CropType crop : CropType.getAllItems()) { - yresp.setYield(YieldType.NONE, crop, 1); - yresp.setYield(YieldType.FERT_MID_ONLY, crop, 2.4); - yresp.setYield(YieldType.FERT_MAX_ONLY, crop, 3); - yresp.setYield(YieldType.IRRIG_MID_ONLY, crop, 1.8); - yresp.setYield(YieldType.IRRIG_MAX_ONLY, crop, 2); - yresp.setYield(YieldType.FERT_IRRIG_MAX, crop, 4); - } + Map<Integer, YieldResponsesItem> getYields() { + YieldResponsesItem yresp = getYieldRespItem(1, 1); Map<Integer, YieldResponsesItem> returnMap = new HashMap<Integer, YieldResponsesItem>(); - for (int i= 1; i<=ModelConfig.NUM_LOCATIONS_PER_COUNTRY; i++) + for (int i= 1; i<=NUM_LOCATIONS_PER_COUNTRY; i++) returnMap.put(i, yresp); return returnMap; } - public Map<Integer, AreasItem> getPreviousArea() { - AreasItem dummyMap = new AreasItem(); - dummyMap.setCropArea(CropType.CEREALS, 9.0); - dummyMap.setCropArea(CropType.FRUIT, 5.0); - dummyMap.setCropArea(CropType.OILCROPS, 10.0); - dummyMap.setCropArea(CropType.STARCHY_ROOTS, 5.0); - dummyMap.setCropArea(CropType.PULSES, 20.0); - dummyMap.setCropArea(CropType.VEGETABLES, 20.0); - dummyMap.setCropArea(CropType.MEAT_OR_PASTURE, 80.0); - dummyMap.setCropArea(CropType.TREENUTS, 5.0); - - dummyMap.setLandCoverArea(LandDataType.LAND, 300); + YieldResponsesItem getYieldRespItem(double base, double pastureFactor) { + YieldResponsesItem yresp = new YieldResponsesItem(); + for (CropType crop : CropType.getAllItems()) { + double factor = base; + if (CropType.MEAT_OR_PASTURE == crop) + factor = base * pastureFactor; + + yresp.setYield(YieldType.NONE, crop, factor); + yresp.setYield(YieldType.FERT_MID_ONLY, crop, factor*2.2); + yresp.setYield(YieldType.FERT_MAX_ONLY, crop, factor*3); + yresp.setYield(YieldType.IRRIG_MID_ONLY, crop, factor*1.8); + yresp.setYield(YieldType.IRRIG_MAX_ONLY, crop, factor*2); + yresp.setYield(YieldType.FERT_IRRIG_MAX, crop, factor*4); + } + return yresp; + } + + + private Map<Integer, AreasItem> getPreviousArea() { + AreasItem dummyMap = getAreaItem(0,0); Map<Integer, AreasItem> returnMap = new HashMap<Integer, AreasItem>(); - for (int i= 1; i<=ModelConfig.NUM_LOCATIONS_PER_COUNTRY; i++) + for (int i= 1; i<=NUM_LOCATIONS_PER_COUNTRY; i++) returnMap.put(i, dummyMap); return returnMap; } - public Map<CropType, Double> getWorldInputEnergy() { + AreasItem getAreaItem(double cereal, double pasture) { + AreasItem dummyMap = new AreasItem(); + dummyMap.setCropArea(CropType.CEREALS, 5.0 + cereal); + dummyMap.setCropArea(CropType.FRUIT, 5.0); + dummyMap.setCropArea(CropType.OILCROPS, 5.0); + dummyMap.setCropArea(CropType.STARCHY_ROOTS, 5.0); + dummyMap.setCropArea(CropType.PULSES, 5.0); + dummyMap.setCropArea(CropType.VEGETABLES, 5.0); + dummyMap.setCropArea(CropType.MEAT_OR_PASTURE, 5.0 + pasture); + dummyMap.setCropArea(CropType.TREENUTS, 5.0); + + dummyMap.setLandCoverArea(LandDataType.LAND, 400); + dummyMap.setLandCoverArea(LandDataType.FOREST, 200 - pasture - cereal); + dummyMap.setLandCoverArea(LandDataType.OTHER_NATURAL, 30); + return dummyMap; + } + + Map<CropType, Double> getWorldInputEnergy() { Map<CropType, Double> dummyMap = new HashMap<CropType, Double>(); dummyMap.put(CropType.CEREALS, 0.4); dummyMap.put(CropType.FRUIT, 5.0); @@ -94,7 +109,7 @@ public class GamsTest { return dummyMap; } - public Map<CropType, Double> getMaxNetImport() { + Map<CropType, Double> getMaxNetImport() { Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>(); dummyMaxExports.put(CropType.CEREALS, 24.0); dummyMaxExports.put(CropType.FRUIT, 5.0); @@ -107,7 +122,7 @@ public class GamsTest { return dummyMaxExports; } - public Map<CropType, Double> getMinNetImport() { + Map<CropType, Double> getMinNetImport() { Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>(); dummyMaxExports.put(CropType.CEREALS, -24.0); dummyMaxExports.put(CropType.FRUIT, -5.0); diff --git a/src/ac/ed/lurg/landuse/AreasItem.java b/src/ac/ed/lurg/landuse/AreasItem.java index e062e2ad..e2b10c1c 100644 --- a/src/ac/ed/lurg/landuse/AreasItem.java +++ b/src/ac/ed/lurg/landuse/AreasItem.java @@ -19,7 +19,7 @@ public class AreasItem implements RasterItem { public double getCropArea(CropType c) { Double d = cropAreas.get(c); - return d == null ? Double.NaN : d; + return d == null ? 0.0 : d; } public void setCropArea(CropType c, double area) { @@ -29,7 +29,7 @@ public class AreasItem implements RasterItem { public double getLandCoverArea(LandDataType c) { Double d = landCoverAreas.get(c); - return d == null ? Double.NaN : d; + return d == null ? 0.0 : d; } public void setLandCoverArea(LandDataType c, double d) { diff --git a/src/ac/ed/lurg/yield/AveragingYieldResponse.java b/src/ac/ed/lurg/yield/AveragingYieldResponse.java index 7f35276b..502d15ee 100644 --- a/src/ac/ed/lurg/yield/AveragingYieldResponse.java +++ b/src/ac/ed/lurg/yield/AveragingYieldResponse.java @@ -10,8 +10,14 @@ public class AveragingYieldResponse extends YieldResponse { // we are assuming that the areas are all the same size, i.e. countries boundaries follow grid structure public void setYield(YieldType type, double yield) { - double previousYield = yields.get(type); - int prevNumAreas = numAreas.get(type); + + double previousYield = 0; + if (yields.containsKey(type)) + previousYield = yields.get(type); + + int prevNumAreas = 0; + if (numAreas.containsKey(type)) + prevNumAreas = numAreas.get(type); double newAvgYield = (previousYield*prevNumAreas + yield)/(prevNumAreas+1); -- GitLab