diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 248f57454e04b65fd96fea66158df11135086337..ca75e43420839214608ae621eec6e717c8c50274 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -2,13 +2,13 @@
SET crop(all_types) / wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots, pasture /;
SET crop_less_pasture(crop) / wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots /;
- SET cereal_crop(crop_less_pasture) / wheat, maize, rice, tropicalCereals /;
- SET non_cereal_crop(crop_less_pasture) / oilcrops, soybean, pulses, starchyRoots /;
- SET feed_crop(crop) / wheat, maize, oilcrops, soybean /;
- SET not_feed_crop(crop) / rice, tropicalCereals, pulses, starchyRoots, pasture /;
+ SET cereal_crop(crop) / wheat, maize, rice, tropicalCereals /;
+ SET non_cereal_crop(crop) / oilcrops, soybean, pulses, starchyRoots, pasture /;
+ SET feed_crop(crop) / wheat, maize, oilcrops, soybean, pasture/;
+ SET feed_crop_less_pasture(feed_crop) / wheat, maize, oilcrops, soybean /;
+ SET not_feed_crop(crop) / rice, tropicalCereals, pulses, starchyRoots /;
SET import_crop(all_types) / meat, wheat, maize, rice, tropicalCereals, oilcrops, soybean, pulses, starchyRoots /;
-
SET location;
PARAMETER suitableLandArea(location) areas of land in Mha;
PARAMETER previousArea(crop, location) areas for previous timestep in Mha;
@@ -23,23 +23,28 @@
PARAMETER maxNetImport(import_crop) maximum net import for each crop based on world market;
PARAMETER minNetImport(import_crop) minimum net import for each crop based on world market;
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?);
-$gdxin %gdxincname%
+ 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_731978163.gdx"
+$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 cropAdj, fitToPreviousAreas
$gdxin
SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /
SCALAR minDemandPerCereal / 0.1 /
-
+
PARAMETER feedDM(crop) kg DM per kg of feed the conversion from feed to meat is done in the R animal product index
/ wheat 0.89
maize 0.7
@@ -47,7 +52,7 @@ $gdxin
soybean 0.7
pasture 1 / ;
* pasture yield is done in DM terms
-
+
PARAMETER baseCost(crop) cost per ha before intensity values just made up at the moment
/ wheat 0.5
maize 0.5
@@ -57,30 +62,25 @@ $gdxin
soybean 0.5
pulses 0.4
starchyRoots 1.0
- pasture 0 / ;
+ pasture 0.2 / ;
SCALAR fertiliserUnitCost / 0.4 /
-
- PARAMETER zeroIntensityYieldFrac(crop) fraction of no fertiliser and irrigiation yield with zero intensity. Most or all arable crops this is 0 but not for pasture
- / pasture 0.3 / ;
-
+
VARIABLES
area(crop, location) total area for each crop - Mha
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 - Mt
netImportAmount(all_types) net imports of crops and meat - Mt
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_less_pasture) supply after exports and feed
+ net_supply(crop) supply after exports and feed
energy total input energy - energy;
- POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount;
+ POSITIVE VARIABLE area, fertI, irrigI, feedAmount;
fertI.L(crop, location) = 0.5;
irrigI.L(crop, location) = 0.5;
- otherIntensity.L(crop, location) = 1.0;
area.L(crop, location) = previousArea(crop, location)
EQUATIONS
@@ -92,31 +92,31 @@ $gdxin
MEAT_DEMAND_CONSTRAINT satisfy demand for meat
MAX_FERT_INTENSITY_CONSTRAINT(crop, location) constraint on maximum fertilizer intensity
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity
- MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) constraint on maximum other 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_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
MIN_FEED_CONSTRAINT constraint on min feed rate
- NET_SUPPLY_EQ(crop_less_pasture) calc net supply for crops
+ NET_SUPPLY_EQ(crop) calc net supply for crops
ENERGY_EQ total energy objective function;
- UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= baseCost(crop) +
- fertiliserUnitCost* fertI(crop, location) +
- (exp(irrigI(crop, location)*2) - 1) * irrigCost(location) +
- exp(otherIntensity(crop, location)*2) - 1;
+ UNIT_ENERGY_EQ(crop, location) .. unitEnergy(crop, location) =E= baseCost(crop) +
+ (fertiliserUnitCost * fertI(crop, location)) +
+ (irrigCost(location) * irrigI(crop, location));
YIELD_EQ(crop, location) .. yield(crop, location) =E= (
yieldNone(crop, location) +
(yieldFertOnly(crop, location) - yieldNone(crop, location)) * ((fertI(crop, location) + delta) ** fertParam(crop, location)) +
(yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * ((irrigI(crop, location) + delta) ** irrigParam(crop, location)) +
(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))
- ) * (zeroIntensityYieldFrac(crop) + otherIntensity(crop, location) * (1-zeroIntensityYieldFrac(crop)));
+ ((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) +
@@ -125,48 +125,77 @@ $gdxin
* (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_less_pasture) .. net_supply(crop_less_pasture) =E= sum(location, area(crop_less_pasture, location) * yield(crop_less_pasture, location)) - feedAmount(crop_less_pasture) + netImportAmount(crop_less_pasture);
+ 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);
+ NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop) + feedAmount(non_cereal_crop);
CEREAL_DEMAND_CONSTRAINT(cereal_crop) .. net_supply(cereal_crop) =G= demand('cereals') * minDemandPerCereal;
TOTAL_CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, net_supply(cereal_crop)) =G= demand('cereals');
- MEAT_DEMAND_CONSTRAINT .. meatEfficency*(sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) + feedDM('pasture') * sum(location, area('pasture', location) * yield('pasture', location))) =G=
- demand('meat') - netImportAmount('meat');
+ MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= demand('meat') - netImportAmount('meat');
MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
- MAX_OTHER_INTENSITY_CONSTRAINT(crop, location) .. otherIntensity(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_CHANGE_CONSTRAINT(location) .. area('pasture', location) =L= (1 + maxLandUseChange) * previousArea('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);
-
MIN_NET_IMPORT_CONSTRAINT(import_crop) .. netImportAmount(import_crop) =G= minNetImport(import_crop);
- MIN_FEED_CONSTRAINT .. sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= minFeedRate * (demand('meat') - netImportAmount('meat'));
+ PASTURE_IMPORT_CONSTRAINT .. netImportAmount('pasture') =E= 0;
+
+ MIN_FEED_CONSTRAINT .. sum(feed_crop_less_pasture, feedDM(feed_crop_less_pasture) * feedAmount(feed_crop_less_pasture)) =G= minFeedRate * (demand('meat') - netImportAmount('meat'));
- ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location))
- + sum(location, sum(crop, area(crop, location)) - sum(crop, previousArea(crop, location)) * landChangeEnergy)
+ ENERGY_EQ .. energy =E= (SUM((crop, location), area(crop, location) * unitEnergy(crop, location))
+ + (sum(location, sum(crop, area(crop, location) - previousArea(crop, location))) * landChangeEnergy)
+ sum(import_crop, (netImportAmount(import_crop)) * worldInputEnergy(import_crop))) / 1000000;
MODEL LAND_USE /ALL/ ;
- SOLVE LAND_USE USING NLP MINIMIZING energy;
-
-display meatEfficency, feedDM, feedAmount.l, area.l, yield.l, demand, netImportAmount.l;
-
+
+
+ scalar basePasture;
+ scalar baseCropland;
+ scalar count / 1 /;
+ scalar newPasture / 0 /;
+ scalar newCropland / 0 /;
+ scalar pastureAdjRate / 0.5 /;
+
+* 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)));
+
+ while((count le 10 and not (abs(basePasture-newPasture) le 0.2 and abs(baseCropland-newCropland) le 0.2)),
+ SOLVE LAND_USE USING NLP MINIMIZING energy;
+
+ newPasture = sum(location, area.l("pasture", location));
+ newCropland = sum(location, sum(crop_less_pasture, area.l(crop_less_pasture, location)));
+
+ cropAdj(crop_less_pasture) = cropAdj(crop_less_pasture) * newCropland / baseCropland;
+ cropAdj("pasture") = cropAdj("pasture") * (pastureAdjRate * newPasture / basePasture + (1-pastureAdjRate));
+
+ display area.l, count, basePasture, baseCropland, newPasture, newCropland, cropAdj;
+ count = count+1 ;
+ );
+ );
+
+ display net_supply.l, demand, feedAmount.l, area.l, yield.l, netImportAmount.l;
+
Scalar ms 'model status', ss 'solve status';
ms=LAND_USE.modelstat;
ss=LAND_USE.solvestat;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 63a07041cb5fd47729e07d277ad154f0e8ffbe67..456700abb09926271a7117b543c731f6d9f7e27a 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -106,7 +106,7 @@ public class ModelConfig {
public static final int START_TIMESTEP = getIntProperty("START_TIMESTEP", 0);
- public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 0);
+ public static final int END_TIMESTEP = getIntProperty("END_TIMESTEP", 1);
public static final int BASE_YEAR = getIntProperty("BASE_YEAR", 2010);
public static final double MAX_IMPORT_CHANGE = getDoubleProperty("MAX_IMPORT_CHANGE", 0.05);
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 50b8487bbe2ed898ba4a5ab4d5ed0a8270342a8f..c10cae0fd8b4988096855286f3185ca032396611 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -2,7 +2,6 @@ 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;
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 4045f40d3415b31474cad1a75428551c44c04aa2..7ab67717c7a382c4f30499c5c778af97b5183ace 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -100,27 +100,27 @@ public class CountryAgent {
private GamsRasterInput getGamsRasterInput(Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy) {
- int previousTimestep = (ModelConfig.START_TIMESTEP == currentTimestep) ? currentTimestep : currentTimestep - 1;
- GamsRasterOutput prevOutput = resultsTimeseries.get(previousTimestep);
-
- Map<CropType, CropUsageData> prevCropUsageData = prevOutput.getCropUsageData();
-
- double allowedImportChange;
- if (currentTimestep == ModelConfig.START_TIMESTEP)
- allowedImportChange = 0;
- else
- allowedImportChange = ModelConfig.MAX_IMPORT_CHANGE;
+ GamsRasterOutput prevOutput;
+ Map<CropType, Double> cropAdjs;
+ boolean calibrate;
+
+ if (ModelConfig.START_TIMESTEP == currentTimestep) { // initialisation time-step
+ calibrate = true;
+ prevOutput = resultsTimeseries.get(currentTimestep);
+ cropAdjs = null;
+ }
+ else { // normal (not the initial) time-step
+ calibrate = false;
+ prevOutput = resultsTimeseries.get(currentTimestep - 1);
+ cropAdjs = prevOutput.getCropAdjustments();
+ }
- Map<CropType, Double> maxNetImport = new HashMap<CropType, Double>();
- Map<CropType, Double> minNetImport = new HashMap<CropType, Double>();
- for (Map.Entry<CropType, CropUsageData> entry : prevCropUsageData.entrySet()) {
- double largerOfProductionOrSupply; // this is what might be better to use than the netImport amount
-
- maxNetImport.put(entry.getKey(), entry.getValue().getNetImports() * (1 + allowedImportChange));
- minNetImport.put(entry.getKey(), entry.getValue().getNetImports() * (1 - allowedImportChange));
+ Map<CropType, Double> baseNetImport = new HashMap<CropType, Double>();
+ for (Map.Entry<CropType, CropUsageData> entry : prevOutput.getCropUsageData().entrySet()) {
+ baseNetImport.put(entry.getKey(), entry.getValue().getNetImports());
}
- GamsCountryInput countryLevelInputs = new GamsCountryInput(country, projectedDemand, worldInputEnergy, maxNetImport, maxNetImport);
+ GamsCountryInput countryLevelInputs = GamsCountryInput.createInput(country, projectedDemand, worldInputEnergy, baseNetImport, cropAdjs, calibrate);
GamsRasterInput input = new GamsRasterInput(countryYieldSurfaces, prevOutput.getCropAreaRaster(), irrigationCostRaster, countryLevelInputs);
return input;
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 0529ba9f587d1dd21962ac97b69a09e3d440414d..1404e490b671395b45074a5db8b5ca827697984e 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -1,7 +1,9 @@
package ac.ed.lurg.country.gams;
+import java.util.HashMap;
import java.util.Map;
+import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.country.Country;
import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
@@ -13,6 +15,8 @@ public class GamsCountryInput {
private Map<CropType, Double> worldInputEnergy;
private Map<CropType, Double> maxNetImport;
private Map<CropType, Double> minNetImport;
+ private Map<CropType, Double> cropAdjustments;
+ private boolean calibrateToObserved;
// limits to areas for each location and crop
@@ -23,15 +27,48 @@ public class GamsCountryInput {
private double landChangeEnergy;*/
- public GamsCountryInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy,
- Map<CropType, Double> maxNetImport, Map<CropType, Double> minNetImport) {
+ private GamsCountryInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy,
+ Map<CropType, Double> maxNetImport, Map<CropType, Double> minNetImport, Map<CropType, Double> cropAdjustments, boolean calibrateToObserved) {
super();
this.country = country;
this.projectedDemand = projectedDemand;
this.worldInputEnergy = worldInputEnergy;
this.maxNetImport = maxNetImport;
this.minNetImport = minNetImport;
+ this.cropAdjustments = cropAdjustments;
+ this.calibrateToObserved = calibrateToObserved;
}
+
+ public static GamsCountryInput createInput(Country country, Map<CommodityType, Double> projectedDemand, Map<CropType, Double> worldInputEnergy,
+ Map<CropType, Double> baseNetImport, Map<CropType, Double> cropAdjustments, boolean calibrateToObserved) {
+
+ Map<CropType, Double> maxNetImport;
+ Map<CropType, Double> minNetImport;
+
+ if (calibrateToObserved) {
+ maxNetImport = baseNetImport;
+ minNetImport = baseNetImport;
+
+ cropAdjustments = new HashMap<CropType, Double>(); // recreate, but might have been null is this case anyway
+ for (CropType c : CropType.getNonMeatTypes()) {
+ cropAdjustments.put(c, 1.0);
+ }
+
+ }
+ else {
+ double allowedImportChange = ModelConfig.MAX_IMPORT_CHANGE;
+
+ maxNetImport = new HashMap<CropType, Double>();
+ minNetImport = new HashMap<CropType, Double>();
+ for (Map.Entry<CropType, Double> entry : baseNetImport.entrySet()) {
+ maxNetImport.put(entry.getKey(), entry.getValue() * (1 + allowedImportChange));
+ minNetImport.put(entry.getKey(), entry.getValue() * (1 - allowedImportChange));
+ }
+ }
+
+ return new GamsCountryInput(country, projectedDemand, worldInputEnergy, maxNetImport, minNetImport, cropAdjustments, calibrateToObserved);
+ }
+
public Country getCountry() {
return country;
@@ -62,7 +99,7 @@ public class GamsCountryInput {
}
public double getMaxLandUseChange() {
- return 0;
+ return 0.95;
}
public double getLandChangeEnergy() {
@@ -72,4 +109,12 @@ public class GamsCountryInput {
public double getMinFeedRate() {
return 0.15;
}
+
+ public Map<CropType, Double> getCropAdjustments() {
+ return cropAdjustments;
+ }
+
+ public boolean isCalibrateToObserved() {
+ return calibrateToObserved;
+ }
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index c03ca9d4bf9b17598d4c0ef94cb980d9515dcb52..96f56b54f5bad88fc6679c97ebd576d777beeb3d 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -12,7 +12,7 @@ public class GamsLocationInput {
private Map<Integer, ? extends AreasItem> previousAreas;
private Map<Integer, ? extends IrrigationCostItem> irrigationCosts;
private GamsCountryInput countryInput;
-
+
public GamsLocationInput(Map<Integer, ? extends YieldResponsesItem> yields, Map<Integer, ? extends AreasItem> previousAreas,
Map<Integer, ? extends IrrigationCostItem> irrigationCosts, GamsCountryInput countryInput) {
super();
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index f42f6cc88cdf067b0ef350c4e6518be70b05728b..84626b7eb0ecda147f7e1fa6c9ce37c998726d7a 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -171,7 +171,14 @@ public class GamsLocationOptimiser {
if (DEBUG) LogWriter.println("\nLandChangeEnergy: " + countryInput.getLandChangeEnergy());
addScalar(inDB.addParameter("minFeedRate", 0), countryInput.getMinFeedRate());
- if (DEBUG) LogWriter.println("\nMinFeedRate:" + 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);
+
+ if (DEBUG) LogWriter.println("\ncropAdj");
+ addItemMapParm(inDB.addParameter("cropAdj", 1), countryInput.getCropAdjustments());
}
private GamsLocationOutput handleResults(GAMSDatabase outDB) {
@@ -188,13 +195,15 @@ public class GamsLocationOptimiser {
GAMSVariable varNetImports = outDB.getVariable("netImportAmount");
GAMSVariable varYields = outDB.getVariable("yield");
GAMSVariable varUnitEnergies = outDB.getVariable("unitEnergy");
+ GAMSParameter parmCropAdj = outDB.getParameter("cropAdj");
double totalArea = 0;
- double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, feedAmount, netImport, yield, unitEnergy;
+ double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, feedAmount, netImport, yield, unitEnergy, cropAdj;
Map<Integer, IntensitiesItem> intensities = new HashMap<Integer, IntensitiesItem>();
Map<Integer, AreasItem> cropAreas = new HashMap<Integer, AreasItem>();
Map<CropType, CropUsageData> cropUsageData = new HashMap<CropType, CropUsageData>();
+ Map<CropType, Double> cropAdjs = new HashMap<CropType, Double>();
CropType prevCropType = null;
for (GAMSVariableRecord rec : varAreas) {
@@ -213,9 +222,11 @@ public class GamsLocationOptimiser {
if (!cropType.equals(prevCropType)) {
feedAmount = varFeedAmount.findRecord(itemName).getLevel();
netImport = cropType.isImportedCrop() ? varNetImports.findRecord(itemName).getLevel() : 0;
-
+ cropAdj = parmCropAdj.findRecord(itemName).getValue();
+
cropUsageData.put(cropType, new CropUsageData(feedAmount, netImport));
- if (DEBUG) LogWriter.println(String.format("\n%s:\tfeedAmount= %.1f,\tnetImports= %.3f", itemName, feedAmount, netImport));
+ cropAdjs.put(cropType, cropAdj);
+ if (DEBUG) LogWriter.println(String.format("\n%s:\tfeedAmount= %.1f,\tnetImports= %.3f,\tcropAdj= %.3f", itemName, feedAmount, netImport,cropAdj));
}
if (area > 0) {
@@ -247,7 +258,7 @@ public class GamsLocationOptimiser {
LogWriter.println(String.format("\nTotal area= %.1f", totalArea));
}
- GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, cropUsageData);
+ GamsLocationOutput results = new GamsLocationOutput(modelStatus, intensities, cropAreas, cropUsageData, cropAdjs);
return results;
}
@@ -258,7 +269,7 @@ public class GamsLocationOptimiser {
private void addItemMapParm(GAMSParameter parm, Map<CropType, Double> itemMap) {
for (Map.Entry<CropType, Double> entry : itemMap.entrySet()) {
double d = entry.getValue();
- if (DEBUG) LogWriter.println(String.format(" %15s,\t %.1f", entry.getKey().getGamsName(), d));
+ if (DEBUG) LogWriter.println(String.format(" %15s,\t %.3f", entry.getKey().getGamsName(), d));
if (!Double.isNaN(d))
parm.addRecord(entry.getKey().getGamsName()).setValue(d);
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 1a91e9151eb51fb7e2c21b6556df6b5f4883c90b..327662bb21ffa2770852e3151613063251bf92c4 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -13,16 +13,19 @@ public class GamsLocationOutput {
Map<Integer, IntensitiesItem> intensities; // data mapped from id (not raster)
Map<Integer, AreasItem> cropAreas;
private Map<CropType, CropUsageData> cropUsageData;
+ Map<CropType, Double> cropAdjustments;
public GamsLocationOutput(int status,
Map<Integer, IntensitiesItem> intensities,
Map<Integer, AreasItem> cropAreas,
- Map<CropType, CropUsageData> cropUsageData) {
+ Map<CropType, CropUsageData> cropUsageData,
+ Map<CropType, Double> cropAdjustments) {
super();
this.status = status;
this.intensities = intensities;
this.cropAreas = cropAreas;
this.cropUsageData = cropUsageData;
+ this.cropAdjustments = cropAdjustments;
}
public int getStatus() {
@@ -37,4 +40,9 @@ public class GamsLocationOutput {
public Map<CropType, CropUsageData> getCommoditiesData() {
return cropUsageData;
}
+
+ public Map<CropType, Double> getCropAdjustments() {
+ return cropAdjustments;
+ }
+
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationTest.java b/src/ac/ed/lurg/country/gams/GamsLocationTest.java
index d9112d7670061e846302e89796c307949f386cb7..3954a5c4c378ce4c1ca7a634c83e94d7d6cf7647 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationTest.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationTest.java
@@ -21,7 +21,7 @@ public class GamsLocationTest {
}
private void run() {
- GamsCountryInput countryLevelInputs = new GamsCountryInput(new Country("Test", "TES", 123), getProjectedDemand(), getWorldInputEnergy(), getMaxNetImport(), getMinNetImport());
+ GamsCountryInput countryLevelInputs = GamsCountryInput.createInput(new Country("Test", "TES", 123), getProjectedDemand(), getWorldInputEnergy(), getBaseNetImport(), null, true);
GamsLocationInput gamsInput = new GamsLocationInput(getYields(), getPreviousArea(), getIrrigationCosts(), countryLevelInputs);
GamsLocationOptimiser opti = new GamsLocationOptimiser(gamsInput);
@@ -123,7 +123,7 @@ public class GamsLocationTest {
return dummyMap;
}
- Map<CropType, Double> getMaxNetImport() {
+ Map<CropType, Double> getBaseNetImport() {
Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
dummyMap.put(CropType.WHEAT, 24.0);
dummyMap.put(CropType.MAIZE, 5.0);
@@ -136,18 +136,4 @@ public class GamsLocationTest {
dummyMap.put(CropType.PASTURE, 5.0);
return dummyMap;
}
-
- Map<CropType, Double> getMinNetImport() {
- Map<CropType, Double> dummyMap = new HashMap<CropType, Double>();
- dummyMap.put(CropType.WHEAT, -24.0);
- dummyMap.put(CropType.MAIZE, -5.0);
- dummyMap.put(CropType.OILCROPS, -5.0);
- dummyMap.put(CropType.TROPICAL_CEREALS, 0.0);
- dummyMap.put(CropType.SOYBEAN, -0.6);
- dummyMap.put(CropType.OILCROPS, -10.0);
- dummyMap.put(CropType.PULSES, -20.0);
- dummyMap.put(CropType.STARCHY_ROOTS, -5.0);
- dummyMap.put(CropType.PASTURE, -5.0);
- return dummyMap;
- }
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index fa0f2c0ddfd2100aed5d1086b9ccd8b5dd15003e..cbe5400d4ab4622c4f91fd94c298bf663f4ce48c 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -11,7 +11,7 @@ public class GamsRasterInput {
private RasterSet<AreasItem> previousAreas;
private RasterSet<IrrigationCostItem> irrigationCost;
private GamsCountryInput countryInput;
-
+
public GamsRasterInput(YieldRaster yields, RasterSet<AreasItem> previousAreas,
RasterSet<IrrigationCostItem> irrigationCost, GamsCountryInput countryInput) {
super();
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index be7e66bc60e337e35f559898c06c47e8e8a5a273..5c6c4c77c3f72dd32f59339894f105b89a04ec2f 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -50,7 +50,7 @@ public class GamsRasterOptimiser {
RasterSet<AreasItem> newAreaRaster = allocAreas(gamsInput.getPreviousAreas(), gamsOutput);
RasterSet<IntensitiesItem> newIntensityRaster = allocIntensities(gamsOutput);
- return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData());
+ return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, newAreaRaster, gamsOutput.getCommoditiesData(), gamsOutput.getCropAdjustments());
}
private RasterSet<AreasItem> allocAreas(Map<Integer, ? extends AreasItem> prevAreasAgg, GamsLocationOutput gamsOutput) {
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
index 78f2cd1252969589a6c01c71a83884bec29bb607..fe6d166851d4327724147f572af58ef1cd8ea84a 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
@@ -15,17 +15,19 @@ public class GamsRasterOutput {
private RasterSet<IntensitiesItem> intensityRaster;
private RasterSet<AreasItem> cropAreaRaster;
private Map<CropType, CropUsageData> cropUsageData;
-
+ private Map<CropType, Double> cropAdjustments;
+
public GamsRasterOutput(RasterSet<AreasItem> cropAreaRaster, Map<CropType, CropUsageData> cropUsageData) {
super();
this.cropAreaRaster = cropAreaRaster;
this.cropUsageData = cropUsageData;
}
- public GamsRasterOutput(int status, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, Map<CropType, CropUsageData> cropUsageData) {
+ public GamsRasterOutput(int status, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, Map<CropType, CropUsageData> cropUsageData, Map<CropType, Double> cropAdjustments) {
this(cropAreaRaster, cropUsageData);
this.status = status;
this.intensityRaster = intensityRaster;
+ this.cropAdjustments = cropAdjustments;
}
public int getStatus() {
@@ -40,4 +42,8 @@ public class GamsRasterOutput {
public Map<CropType, CropUsageData> getCropUsageData() {
return cropUsageData;
}
+
+ public Map<CropType, Double> getCropAdjustments() {
+ return cropAdjustments;
+ }
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterTest.java b/src/ac/ed/lurg/country/gams/GamsRasterTest.java
index 6ce825fd9f7760f2af83d5e1575b10812afcdde5..902b1bedad74c849de874f11842e3f658eb4e1e1 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterTest.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterTest.java
@@ -17,7 +17,7 @@ public class GamsRasterTest extends GamsLocationTest {
}
private void run() {
- GamsCountryInput countryLevelInputs = new GamsCountryInput(new Country("Test", "TES", 123), getProjectedDemand(), getWorldInputEnergy(), getMaxNetImport(), getMinNetImport());
+ GamsCountryInput countryLevelInputs = GamsCountryInput.createInput(new Country("Test", "TES", 123), getProjectedDemand(), getWorldInputEnergy(), getBaseNetImport(), null, true);
GamsRasterInput input = new GamsRasterInput(getYieldRaster(), getPreviousAreaRaster(), getIrrigationCost(), countryLevelInputs);
GamsRasterOptimiser opti = new GamsRasterOptimiser(input);
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index ad89e96c0a732771533f71413f334b47a73ca862..19554b58a22385510e3a3842d7a0e9ec4589f69a 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/src/ac/ed/lurg/yield/YieldResponse.java
@@ -29,14 +29,14 @@ public class YieldResponse {
public double getFertParam() {
if (fertParm == 0)
- 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));
+ 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));
return fertParm;
}
public double getIrrigParam() {
if (irrigParm == 0)
- irrigParm = calcParam (0, 0.7, 1); // we don't have a mid irrigation figure, so lets assume 60% at mid point
+ irrigParm = calcParam (0, 0.8, 1); // we don't have a mid irrigation figure, so lets assume 60% at mid point
return irrigParm;
}