diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms index 587cae08aa08bcf09befc977e4684bc24b7f8d7f..2da56a774f6274ddb0d45075e6e4136b06a865f8 100644 --- a/GAMS/IntExtOpt.gms +++ b/GAMS/IntExtOpt.gms @@ -3,13 +3,15 @@ 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*1 /; + SET location / 1*5 /; PARAMETER previous_area(crop, location) areas for previous timestep in ha; 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) yield response to fertilizer parameter; + PARAMETER irrigParam(crop) yield response to irrigation parameter; 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; @@ -20,8 +22,6 @@ 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; - SCALAR fertParam yield response to fertilizer parameter; - SCALAR irrigParam yield response to irrigation parameter; $gdxin %gdxincname% $load previous_area, demand, yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, fertParam, irrigParam, world_input_energy, maxNetImport, minNetImport, meatEfficency, maxLandUseChange, tradeBarrier, landChangeEnergy, minFeedRate @@ -56,7 +56,7 @@ $gdxin 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 + MEAT_DEMAND_CONSTRAINT satisfy demand for crop 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 @@ -69,14 +69,14 @@ $gdxin MIN_FEED_CONSTRAINT constraint on min feed rate ENERGY_EQ total energy objective function; - 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; + UNIT_ENERGY_EQ(crop, location) .. unit_energy(crop, location) =E= power((1+ fertI(crop, location)+irrigI(crop, location)+otherIntensity(crop, location)), 2); YIELD_EQ(crop, location) .. yield(crop, location) =E= ( yieldNone(crop, location) + - (yieldFertOnly(crop, location) - yieldNone(crop, location)) * ((fertI(crop, location) + delta) ** fertParam) + - (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam) + - (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) * ((fertI(crop, location) + delta) ** fertParam) * ((irrigI(crop, location) + delta) ** irrigParam) + (yieldFertOnly(crop, location) - yieldNone(crop, location)) * ((fertI(crop, location) + delta) ** fertParam(crop)) + + (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop)) + + (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) * + ((fertI(crop, location) + delta) ** fertParam(crop)) * ((irrigI(crop, location) + delta) ** irrigParam(crop)) ) * otherIntensity(crop, location); * YIELD_EQ(crop, location) .. yield(crop, location) =E= ( @@ -91,7 +91,7 @@ $gdxin 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, location) .. fertI(crop, location) =L= 1; MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1; @@ -112,9 +112,9 @@ $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, location), area(crop, location)*unit_energy(crop, location)) + 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)); + + sum(crop, (netImportAmount(crop) * tradeBarrier) * world_input_energy(crop))) / 1000000; MODEL LAND_USE /ALL/ ; SOLVE LAND_USE USING NLP MINIMIZING energy; diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java index b86a7590a838aca2333a3b38ea70cd173b5e8086..af3df404cfe5566ee35e1599b4846192364bd14b 100644 --- a/src/ac/ed/lurg/country/CountryAgent.java +++ b/src/ac/ed/lurg/country/CountryAgent.java @@ -89,13 +89,23 @@ public class CountryAgent implements GamsInputData { } @Override - public double getFertParam() { - return 0.4; + public Map<CropType, Double> getFertParam() { + Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>(); + + for (CropType crop : CropType.getAllItems()) + dummyMaxExports.put(crop, 0.4); + + return dummyMaxExports; } @Override - public double getIrrigParam() { - return 0.5; + public Map<CropType, Double> getIrrigParam() { + Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>(); + + for (CropType crop : CropType.getAllItems()) + dummyMaxExports.put(crop, 0.4); + + return dummyMaxExports; } diff --git a/src/ac/ed/lurg/country/gams/GamsInputData.java b/src/ac/ed/lurg/country/gams/GamsInputData.java index 80ac380aa40375c343bf18852a8f280fb81c7a0b..8bfb24e5df20119ac48b310eec197905b4a51da7 100644 --- a/src/ac/ed/lurg/country/gams/GamsInputData.java +++ b/src/ac/ed/lurg/country/gams/GamsInputData.java @@ -21,7 +21,7 @@ public interface GamsInputData { double getTradeBarrier(); double getLandChangeEnergy(); double getMinFeedRate(); - double getFertParam(); - double getIrrigParam(); + Map<CropType, Double> getFertParam(); + Map<CropType, Double> getIrrigParam(); } diff --git a/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java index ab4ae9ef3e67e5875a5ae6d5077c1ded3ba47143..0c6e04845a283a7490dcc33f8f7815b4eac30394 100644 --- a/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsLandUseOptimiser.java @@ -39,7 +39,7 @@ public class GamsLandUseOptimiser { GAMSWorkspace ws = new GAMSWorkspace(wsInfo); GAMSDatabase db = ws.addDatabase(); - int numLocations = 1; + int numLocations = 5; LogWriter.println("\nPrevious areas"); GAMSParameter parm = db.addParameter("previous_area", 2, "the previous area for each land use"); @@ -69,6 +69,14 @@ public class GamsLandUseOptimiser { parm = db.addParameter("yieldBoth", 2, "ref yield t per ha"); for (int i= 1; i<=numLocations; i++) addLocationMapParm(parm, inputData.getYieldBoth(), i); + + LogWriter.println("\nfertParam"); + parm = db.addParameter("fertParam", 1); + addItemMapParm(parm, inputData.getFertParam()); + + LogWriter.println("\nirrigParam"); + parm = db.addParameter("irrigParam", 1); + addItemMapParm(parm, inputData.getIrrigParam()); 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"); @@ -87,12 +95,13 @@ public class GamsLandUseOptimiser { 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("minFeedRate", 0, "minimum rate of feed for producing animal products"), inputData.getMinFeedRate()); - addScalar(db.addParameter("fertParam", 0), inputData.getFertParam()); - addScalar(db.addParameter("irrigParam", 0), inputData.getIrrigParam()); - + GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.GAMS_MODEL); GAMSOptions opt = ws.addOptions(); opt.defines("gdxincname", db.getName()); + + long startTime = System.currentTimeMillis(); + gamsJob.run(opt, db); System.out.println( @@ -121,14 +130,16 @@ public class GamsLandUseOptimiser { netImports = varNetImports.findRecord(itemName).getLevel(); totalArea += rec.getLevel(); - 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)); + if (area > 0) + 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); } LogWriter.println(String.format(" Total area= %.1f", totalArea)); //cleanup(ws.workingDirectory()); + LogWriter.println("Took " + (System.currentTimeMillis() - startTime) + " ms to run"); return output; } @@ -149,9 +160,9 @@ public class GamsLandUseOptimiser { Vector<String> v = new Vector<String>(); v.add(entry.getKey().getGamsName()); v.add(Integer.toString(locationId)); - if (locationId== 2 & entry.getKey().getGamsName().equals("cereals")) - parm.addRecord(v).setValue(0.99); - else + // if (locationId== 2 & entry.getKey().getGamsName().equals("cereals")) + // parm.addRecord(v).setValue(2.5); + // else parm.addRecord(v).setValue(entry.getValue()); } } diff --git a/src/ac/ed/lurg/country/gams/GamsTest.java b/src/ac/ed/lurg/country/gams/GamsTest.java index a28d2411fb772fb20157a2aa4dfc8a18b26b0bfb..1121e726c495d1bfe88c73e7c92ebbae98dc3be8 100644 --- a/src/ac/ed/lurg/country/gams/GamsTest.java +++ b/src/ac/ed/lurg/country/gams/GamsTest.java @@ -23,14 +23,14 @@ public class GamsTest implements GamsInputData { @Override public Map<CropType, Double> getProjectedDemand() { Map<CropType, Double> dummyMap = new HashMap<CropType, Double>(); - dummyMap.put(CropType.CEREALS, 50.0); - dummyMap.put(CropType.FRUIT, 250.0); - dummyMap.put(CropType.OILCROPS, 100.0); - dummyMap.put(CropType.STARCHY_ROOTS, 250.0); - dummyMap.put(CropType.PULSES, 260.0); + dummyMap.put(CropType.CEREALS, 100.0); + dummyMap.put(CropType.FRUIT, 50.0); + dummyMap.put(CropType.OILCROPS, 50.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.TREENUTS, 50.0); + dummyMap.put(CropType.MEAT_OR_PASTURE, 80.0); + dummyMap.put(CropType.TREENUTS, 10.0); return dummyMap; } @@ -56,7 +56,7 @@ public class GamsTest implements GamsInputData { 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.MEAT_OR_PASTURE, 2.0); dummyMap.put(CropType.TREENUTS, 3.0); return dummyMap; } @@ -174,12 +174,22 @@ public class GamsTest implements GamsInputData { } @Override - public double getFertParam() { - return 0.4; + public Map<CropType, Double> getFertParam() { + Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>(); + + for (CropType crop : CropType.getAllItems()) + dummyMaxExports.put(crop, 0.6); + + return dummyMaxExports; } @Override - public double getIrrigParam() { - return 0.3; + public Map<CropType, Double> getIrrigParam() { + Map<CropType, Double> dummyMaxExports = new HashMap<CropType, Double>(); + + for (CropType crop : CropType.getAllItems()) + dummyMaxExports.put(crop, 0.5); + + return dummyMaxExports; } }