diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index de004b9d0157751300daab6823c6d10751a7ee44..ecc45a225ce955c5f6c5ce5f18855ed92431aac2 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -25,10 +25,9 @@
PARAMETER irrigCost(location) irrigation cost;
PARAMETER cropAdj(crop) this is the fudge factor that allows the model to fit the observed data
- SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
- SCALAR maxLandUseChange max rate of land use change;
- SCALAR landChangeEnergy energy required to add ha of agricultural land;
- SCALAR minFeedRate minimum rate of feed for producing animal products (why is this needed?);
+ SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
+ SCALAR landChangeEnergy energy required to add ha of agricultural land;
+ SCALAR minFeedRate minimum rate of feed for producing animal products (why is this needed?);
SCALAR fitToPreviousAreas controlling parameter 0 means dont fit and 1 is to calibrate the cropAdj in the calibration year
@@ -37,7 +36,7 @@ $gdxin %gdxincname%
$load location, suitableLandArea, previousArea, demand, landChangeEnergy
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
$load fertParam, irrigParam, worldInputEnergy, maxNetImport, minNetImport
-$load meatEfficency, maxLandUseChange, minFeedRate, irrigCost
+$load meatEfficency, minFeedRate, irrigCost
$load cropAdj, fitToPreviousAreas
$gdxin
@@ -54,17 +53,17 @@ $gdxin
* pasture yield is done in DM terms
PARAMETER baseCost(crop) cost per ha before intensity values just made up at the moment
- / wheat 0.8
- maize 1.0
- rice 1.0
- tropicalCereals 0.5
- oilcrops 0.5
- soybean 0.5
- pulses 0.4
- starchyRoots 1.0
+ / wheat 1.1
+ maize 0.9
+ rice 1.1
+ tropicalCereals 0.9
+ oilcrops 0.7
+ soybean 0.6
+ pulses 0.5
+ starchyRoots 4.0
pasture 0.1 / ;
- SCALAR fertiliserUnitCost / 0.4 /
+ SCALAR fertiliserUnitCost / 1.5 /
VARIABLES
area(crop, location) total area for each crop - Mha
@@ -76,9 +75,13 @@ $gdxin
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
+ cropIncrease(location)
+ cropDecrease(location)
+ pastureIncrease(location)
+ pastureDecrease(location)
energy total input energy - energy;
- POSITIVE VARIABLE area, fertI, irrigI, feedAmount, agriLandExpansion;
+ POSITIVE VARIABLE area, fertI, irrigI, feedAmount, agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease;
fertI.L(crop, location) = 0.5;
irrigI.L(crop, location) = 0.5;
@@ -94,23 +97,22 @@ $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
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
- PASTURE_IMPORT_CONSTRAINT
+ PASTURE_IMPORT_CONSTRAINT constraint to not import pasture
MIN_FEED_CONSTRAINT constraint on min feed rate
NET_SUPPLY_EQ(crop) calc net supply for crops
AGRI_LAND_EXPANSION_CALC(location) calc agriLandExpansion
+ CROP_INCREASE_CALC(location)
+ CROP_DECREASE_CALC(location)
+ PASTURE_INCREASE_CONV_CALC(location)
+ PASTURE_DECREASE_CONV_CALC(location)
ENERGY_EQ total energy objective function;
- UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= baseCost(crop) +
+ UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= ( baseCost(crop) +
(fertiliserUnitCost * fertI(crop, location)) +
- (irrigCost(location) * irrigI(crop, location));
+ (irrigCost(location) * irrigI(crop, location)) ) ;
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yieldNone(crop, location) +
@@ -119,17 +121,10 @@ $gdxin
(yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
((fertI(crop, location) + delta) ** fertParam(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location))
) * cropAdj(crop);
-
-* YIELD_EQ(crop, location) .. yield(crop, location) =E= (
-* yieldNone(crop, location) +
-* (yieldFertOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-10*fertI(crop, location)) +
-* (yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-5*irrigI(crop, location)) +
-* (yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) * (1 - exp(-fertParam(crop, location)*fertI(crop, location))) * (1 - exp(-irrigParam(crop, location)*irrigI(crop, location)))
-* ) * (zeroIntensityYieldFrac(crop) + otherIntensity(crop, location) * (1-zeroIntensityYieldFrac(crop)));
-
+
NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + netImportAmount(crop);
- NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop) + feedAmount(non_cereal_crop);
+ NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop);
CEREAL_DEMAND_CONSTRAINT(cereal_crop) .. net_supply(cereal_crop) =G= demand('cereals') * minDemandPerCereal;
@@ -141,13 +136,7 @@ $gdxin
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
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));
-
+
NON_FEED_CROP_CONSTRAINT(not_feed_crop) .. feedAmount(not_feed_crop) =E= 0;
MAX_NET_IMPORT_CONSTRAINT(import_crop) .. netImportAmount(import_crop) =L= maxNetImport(import_crop);
@@ -158,10 +147,15 @@ $gdxin
AGRI_LAND_EXPANSION_CALC(location) .. agriLandExpansion(location) =G= sum(crop, area(crop, location) - previousArea(crop, location));
-* CROP_TO_PASTURE_CONV_CALC(location) .. cropToPasture(location) =E= sum(crop_less_pasture, area(crop_less_pasture, location)) - area('pasture', location)
+ CROP_INCREASE_CALC(location) .. cropIncrease(location) =G= sum(crop_less_pasture, area(crop_less_pasture, location) - previousArea(crop_less_pasture, location));
+ CROP_DECREASE_CALC(location) .. cropDecrease(location) =G= -sum(crop_less_pasture, area(crop_less_pasture, location) - previousArea(crop_less_pasture, location));
+ PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= area('pasture', location) - previousArea('pasture', location);
+ PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -area('pasture', location) - previousArea('pasture', location);
ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location))
+ (sum(location, agriLandExpansion(location)) * landChangeEnergy)
+ + (sum(location, (((cropIncrease(location) + cropDecrease(location))/sum(crop_less_pasture, previousArea(crop_less_pasture, location)))) * landChangeEnergy/5 * sum(crop_less_pasture, previousArea(crop_less_pasture, location))) )
+ + (sum(location, (((pastureIncrease(location) + pastureDecrease(location))/previousArea('pasture', location))) * landChangeEnergy/5 *previousArea('pasture', location)) )
+ sum(import_crop, (netImportAmount(import_crop)) * worldInputEnergy(import_crop))) / 1000000;
MODEL LAND_USE /ALL/ ;
@@ -172,25 +166,37 @@ $gdxin
scalar count / 1 /;
scalar newPasture / 0 /;
scalar newCropland / 0 /;
- scalar pastureAdjRate / 0.6 /;
+
+ parameter temp(location);
* if fitToPreviousAreas is zero then we just optimise
if (fitToPreviousAreas eq 0,
SOLVE LAND_USE USING NLP MINIMIZING energy;
+
+
* if fitToPreviousAreas is not zero then we try to find a cropAdj factors that give the previous areas ie we attempt to calibrate
else
basePasture = sum(location, previousArea("pasture", location));
- baseCropland = sum(location, sum(crop_less_pasture, previousArea(crop_less_pasture, location)));
+ baseCropland = sum((location, crop_less_pasture), previousArea(crop_less_pasture, location));
- while((count le 10 and not (abs(basePasture-newPasture) le 0.1 and abs(baseCropland-newCropland) le 0.1)),
+ while((count le 1 and not (abs(basePasture-newPasture) le 0.1 and abs(baseCropland-newCropland) le 0.1)),
SOLVE LAND_USE USING NLP MINIMIZING energy;
+ display pastureIncrease.l, pastureDecrease.l, cropDecrease.l, cropIncrease.l;
+
+ temp(location) = SQRT(SQR(cropIncrease.l(location)/sum(crop_less_pasture, previousArea(crop_less_pasture, location)))) * landChangeEnergy/5 * sum(crop_less_pasture, previousArea(crop_less_pasture, location))
+
+ display temp;
+
newPasture = sum(location, area.l("pasture", location));
- newCropland = sum(location, sum(crop_less_pasture, area.l(crop_less_pasture, location)));
+ newCropland = sum((location, crop_less_pasture), area.l(crop_less_pasture, location));
+
+ cropAdj(crop_less_pasture) = cropAdj(crop_less_pasture) * (newCropland / baseCropland);
+ cropAdj("pasture") = cropAdj("pasture") * (newPasture / basePasture);
- cropAdj(crop_less_pasture) = cropAdj(crop_less_pasture) * newCropland / baseCropland;
- cropAdj("pasture") = cropAdj("pasture") * (pastureAdjRate * newPasture / basePasture + (1-pastureAdjRate));
+* cropAdj(crop) = max(cropAdj(crop), 0.2);
+* cropAdj(crop) = min(cropAdj(crop), 5);
display area.l, count, basePasture, baseCropland, newPasture, newCropland, cropAdj;
count = count+1 ;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index c304e5e967ac4b351ca011c292ee2a8f70b294e4..f0b75b426840edd619856041307721a44068475a 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -103,10 +103,11 @@ public class ModelConfig {
public static final String INITAL_LAND_COVER_DIR = SPATIAL_DATA_DIR + File.separator + "initLandUse";
public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries.asc";
public static final String IRRIGATION_COST_FILE = SPATIAL_DATA_DIR + File.separator + "irrigation_cost.asc";;
+ public static final double IRRIG_COST_SCALE_FACTOR = 3.0;
public static final int START_TIMESTEP = getIntProperty("START_TIMESTEP", 0);
- public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 1);
+ public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 0);
public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2010);
public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", 0.2);
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 9500e56e2c05c43d7fe55f5e95e1f791e73e32c6..27bbc0c54a3708356c20e3f691bba2b1baa50227 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -241,7 +241,7 @@ public class ModelMain {
// DEBUG code
- if (!(country.getCountryName().equals("United States of America") || country.getCountryName().equals("China"))) { //|| country.getCountryName().equals("China")
+ if (!(country.getCountryName().equals("United States of AmericaXX") || country.getCountryName().equals("China"))) { //|| country.getCountryName().equals("China")
continue;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 51e5f44c0d6cb733fc2814649fd32ba5fc35656c..33b511768c7afb9e6d941784da69f6e129892a78 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -97,15 +97,11 @@ public class GamsCountryInput {
// }
public double getMeatEfficiency() {
- return 0.5; // this is already handled by the feed conversion efficiency for each animal product
- }
-
- public double getMaxLandUseChange() {
- return 0.95;
+ return 0.6; // this is already handled by the feed conversion efficiency for each animal product
}
public double getLandChangeEnergy() {
- return 0.0;
+ return 1.0;
}
public double getMinFeedRate() {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 8bdf7f14c3a12872aa9990cb3a6a50c97ce96e18..8610ceb3ba6ae9b5880c634525efb750c78181d4 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -169,9 +169,6 @@ 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());
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index fac3693e49376d461408eac346f7338488631091..94d189dc4ee5423b6183abf6edd42ec649edd3ae 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -8,11 +8,9 @@ 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.landuse.IrrigationCostItem;
-import ac.ed.lurg.output.RasterOutputer;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.YieldType;
@@ -251,8 +249,8 @@ public class GamsRasterOptimiser {
}
- int numCerealCats = 2;
- int numPastureCats = 2;
+ int numCerealCats = 5;
+ int numPastureCats = 1;
int thisShouldLookAtCropsOtherThanJustWheat; // need to consider other crops, and perhaps other yieldTypes as well
List<Double> wheatlDivisions = getDivisions(yieldRaster, CropType.WHEAT, numCerealCats);
diff --git a/src/ac/ed/lurg/landuse/IrrigiationCostReader.java b/src/ac/ed/lurg/landuse/IrrigiationCostReader.java
index 74f609eba77081192deb2dfa06d78af29bebba9c..742bb88886f21ec1f73c8df6a50ac27dc2528a57 100644
--- a/src/ac/ed/lurg/landuse/IrrigiationCostReader.java
+++ b/src/ac/ed/lurg/landuse/IrrigiationCostReader.java
@@ -1,5 +1,6 @@
package ac.ed.lurg.landuse;
+import ac.ed.lurg.ModelConfig;
import ac.sac.raster.AbstractRasterReader;
import ac.sac.raster.RasterSet;
@@ -12,6 +13,6 @@ public class IrrigiationCostReader extends AbstractRasterReader<IrrigationCostIt
@Override
public void setData(IrrigationCostItem item, String token) {
double irrigCost = Double.parseDouble(token);
- item.setIrrigCost(irrigCost);
+ item.setIrrigCost(irrigCost * ModelConfig.IRRIG_COST_SCALE_FACTOR);
}
}
diff --git a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
index 158165ec72a53d0f654604d17bd0aab3136d74d5..efb48f754314254356875bc36908fe31ab675355 100644
--- a/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
+++ b/src/ac/ed/lurg/yield/LPJYieldResponseMapReader.java
@@ -123,8 +123,15 @@ public class LPJYieldResponseMapReader {
item.setYield(maxIrrigYieldType, CropType.PULSES, record.teSWirr);
item.setYield(noIrrigYieldType, CropType.STARCHY_ROOTS, record.teSW * 1.1);
item.setYield(maxIrrigYieldType, CropType.STARCHY_ROOTS, record.teSWirr * 1.1);
- item.setYield(noIrrigYieldType, CropType.PASTURE, record.teSW * 2);
- item.setYield(maxIrrigYieldType, CropType.PASTURE, record.teSW * 2.1);
+
+ if (noIrrigYieldType.getFertiliserRate() == FertiliserRate.MAX_FERT) {
+ item.setYield(noIrrigYieldType, CropType.PASTURE, record.teSW*1.2);
+ item.setYield(maxIrrigYieldType, CropType.PASTURE, record.teSWirr*1.2);
+ }
+ else {
+ item.setYield(noIrrigYieldType, CropType.PASTURE, record.teSW*0.6);
+ item.setYield(maxIrrigYieldType, CropType.PASTURE, record.teSWirr*0.7);
+ }
}
class YieldRecord {
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index 19554b58a22385510e3a3842d7a0e9ec4589f69a..9fa8eada5649ab3efc2a4a0f8181d4e9565d9726 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/src/ac/ed/lurg/yield/YieldResponse.java
@@ -29,7 +29,7 @@ public class YieldResponse {
public double getFertParam() {
if (fertParm == 0)
- fertParm = calcParam (0, 0.8, 1); // we do have MID fert data, but looks wrong calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG));
+ fertParm = calcParam (0, 0.7, 1); // we do have MID fert data, but looks wrong calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG));
return fertParm;
}
diff --git a/src/ac/sac/raster/RasterSet.java b/src/ac/sac/raster/RasterSet.java
index 4537e751dfd1c941fc36a2c34f72021c74669662..768c30e7e8427f7cbf50e91972d50bf221d2de14 100755
--- a/src/ac/sac/raster/RasterSet.java
+++ b/src/ac/sac/raster/RasterSet.java
@@ -5,7 +5,6 @@ import java.util.Collection;
import java.util.HashMap;
import java.util.Map;
-import ac.ed.lurg.landuse.AreasItem;
import ac.ed.lurg.utils.LogWriter;
/* Class which holds raster data. The generics used to defines the type of raster data held.