diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index d0201ff9dd2f850ad8bf9ac090ae5d2772e05647..665ea1ad598139ad53c4373c95a7d63c1f1e1321 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -16,7 +16,7 @@
ALIAS (land_cover, land_cover_after);
SET location;
- PARAMETER suitableLandArea(location) areas of land in Mha;
+ PARAMETER suitableLandArea(location) areas of land in Mha;
PARAMETER previousCropArea(crop, location) areas for previous timestep in Mha;
PARAMETER previousFertIntensity(crop, location);
PARAMETER previousIrrigIntensity(crop, location);
@@ -45,12 +45,6 @@
PARAMETER seedAndWasteRate(all_types) rate of use for seed and waste combined;
PARAMETER subsidyRate(crop) rates of subsidy compared to costs;
- PARAMETER agriExpansionCost(location) price for increase agricultural area varies based on amount of managed or unmanaged forest - 1000$ per ha;
- SCALAR pastureIncCost price for increasing pasture area - 1000$ per ha;
- SCALAR cropIncCost price for increasing crop - 1000$ per ha;
- SCALAR cropDecCost price for decreasing crop - 1000$ per ha;
- SCALAR pastureDecCost price for decreasing pasture - 1000$ per ha;
-
SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
SCALAR fertiliserUnitCost fert cost at max fert rate;
SCALAR otherIParam yield response to other intensity;
@@ -61,17 +55,17 @@
SCALAR maxLandExpansionRate max rate of country land expansion;
PARAMETER previousLandCoverArea(land_cover, location) land cover area in Mha;
- PARAMETER minimumLandCoverArea(land_cover, location) minimum land cover area to constrain conversion;
+ PARAMETER minimumLandCoverArea(land_cover, location) minimum land cover area to constrain conversion in Mha;
PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location) carbon flux - MtC-eq per Mha;
PARAMETER woodYield(land_cover_before, land_cover_after, location) wood yield - MtC-eq per Mha;
- SCALAR carbonPrice price of carbon - $ per tonne or million$ per Mt;
- SCALAR woodPrice price of wood and timber - $ per tonne or million$ per Mt;
+ SCALAR carbonPrice price of carbon - $1000 per tonne;
+ SCALAR woodPrice price of wood and timber - $1000 per tC-eq;
PARAMETER conversionCost(land_cover, land_cover) cost of converting from one land cover to another;
*$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
$gdxin %gdxincname%
-$load location, suitableLandArea, demand, agriExpansionCost, cropIncCost, pastureIncCost, cropDecCost, pastureDecCost
+$load location, suitableLandArea, demand
$load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
@@ -80,10 +74,6 @@ $load previousLandCoverArea, carbonFluxRate, woodYield, carbonPrice, woodPrice
$load minimumLandCoverArea, conversionCost
$gdxin
-*convert units from $/t to billion$/Mt
- carbonPrice = carbonPrice * 0.001;
- woodPrice = woodPrice * 0.001;
-
SCALAR delta "use to smooth power function see 7.5 www.gams.com/dd/docs/solversconopt.pdf" / 0.00000000001 /;
demand(cereal_crop) = demand('cereals') * minDemandPerCereal(cereal_crop);
@@ -120,7 +110,7 @@ $gdxin
baseCost('pasture') = 0.04;
otherIntCost('pasture') = 0.8 + otherICost;
- SCALAR forestExpansionMaxRate / 0.01 /;
+ SCALAR forestExpansionMaxRate / 1 /;
VARIABLES
cropArea(crop, location) total area for crops - Mha
@@ -134,11 +124,6 @@ $gdxin
yield(crop, location) yield per area for each crop - t per ha
unitCost(crop, location) cost per area for each crop - cost
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)
landCoverArea(land_cover, location) land cover area in Mha
landCoverChange(land_cover_before, land_cover_after, location) land cover change
deforestedArea(forest, location) area of forest cut down
@@ -179,12 +164,6 @@ $gdxin
PASTURE_EXPORT_CONSTRAINT
IRRIGATION_CONSTRAINT(location) constraint on water usage
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)
- PASTURE_TOTAL_CHANGE_CONSTRAINT(location)
CROPLAND_LAND_COVER_CALC(location) cropland area equals sum of all crop areas
PASTURE_LAND_COVER_CALC(location) pasture area (land cover) equals pasture area (land use)
* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location) constraint on land cover conversion
@@ -248,14 +227,6 @@ $gdxin
IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) * (1.0 - unhandledCropRate) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * cropArea(crop, location));
- AGRI_LAND_EXPANSION_CALC(location) .. agriLandExpansion(location) =G= sum(crop, cropArea(crop, location) - previousCropArea(crop, location));
-
- CROP_INCREASE_CALC(location) .. cropIncrease(location) =G= sum(crop_less_pasture, cropArea(crop_less_pasture, location) - previousCropArea(crop_less_pasture, location));
- CROP_DECREASE_CALC(location) .. cropDecrease(location) =G= -sum(crop_less_pasture, cropArea(crop_less_pasture, location) - previousCropArea(crop_less_pasture, location));
- PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= cropArea('pasture', location) - previousCropArea('pasture', location);
- PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(cropArea('pasture', location) - previousCropArea('pasture', location));
- PASTURE_TOTAL_CHANGE_CONSTRAINT(location) .. pastureIncrease(location) -pastureDecrease(location) =G= cropArea('pasture', location) - previousCropArea('pasture', location);
-
CROPLAND_LAND_COVER_CALC(location) .. landCoverArea('cropland', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate);
PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
@@ -291,23 +262,16 @@ $gdxin
CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
- CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after));
+* Subtracting cost of not fully grown timberForest reserved under minimumLandCoverArea
+ CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after)) -
+ minimumLandCoverArea('timberForest', location) * conversionCost('timberForest', 'timberForest');
* CROPLAND_INCREASE_CONSTRAINT(location) .. cropIncrease(location)/(1.0 - unhandledCropRate) =L= maxLandChangeRate * ( suitableLandArea(location) -sum(crop_less_pasture, previousCropArea(crop_less_pasture, location))/*(1.0 - unhandledCropRate) - previousCropArea('pasture', location) );
* PASTURE_DECREASE_CONSTRAINT(location) .. pastureDecrease(location) =L= maxLandChangeRate * previousCropArea('pasture', location);
COST_EQ .. total_cost =E=
(
- ( SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) +
-
- SUM(location,
- agriExpansionCost(location) * agriLandExpansion(location) +
- cropIncCost * cropIncrease(location) +
- pastureIncCost * pastureIncrease(location) +
- cropDecCost * cropDecrease(location) +
- pastureDecCost * pastureDecrease(location)
- )
- ) * domesticPriceMarkup +
+ ( SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) ) * domesticPriceMarkup +
SUM(location, carbonFlux(location) * carbonPrice) - SUM(location, woodHarvest(location) * woodPrice) +
@@ -329,7 +293,7 @@ $gdxin
SOLVE LAND_USE USING NLP MINIMIZING total_cost;
- display agriLandExpansion.L, previousCropArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, cropArea.L, cropIncrease.L, cropDecrease.L, pastureIncrease.L, pastureDecrease.L;
+ display previousCropArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, cropArea.L;
display net_supply.l, demand, ruminantFeed.l, monogastricFeed.l, importAmount.l, exportAmount.l;
* Calculate summary information used in Java process
diff --git a/debug_config.properties b/debug_config.properties
index 7e7ff02c108d732e79d99300a7d63993e8f5b1af..680b2e71c94c13468b976d9a018889df5f6305dc 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -4,22 +4,21 @@ YIELD_DIR=C:/Users/Bart/Documents/PhD/LURG/rcp60
OUTPUT_DIR = C:/Users/Bart/git/plumv2/output
-END_TIMESTEP=10
-
-END_FIRST_STAGE_CALIBRATION=10
+BASE_YEAR=2020
+START_TIMESTEP=0
TIMESTEP_SIZE=1
+END_TIMESTEP=10
-IS_CALIBRATION_RUN=true
+IS_CALIBRATION_RUN=false
+END_FIRST_STAGE_CALIBRATION=0
GENERATE_NEW_YIELD_CLUSTERS=false
-DEMAND_PRICE_IMPORT_AND_PROD_COST=false
-
YIELD_FILENAME=yield.out
DEBUG_LIMIT_COUNTRIES=true
-DEBUG_COUNTRY_NAME=Central America
-GAMS_COUNTRY_TO_SAVE=Central America
+DEBUG_COUNTRY_NAME=Brazil
+GAMS_COUNTRY_TO_SAVE=Brazil
-WOOD_PRICE=25
-CARBON_PRICE=10
\ No newline at end of file
+WOOD_PRICE=0.15
+CARBON_PRICE=0.05
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index ebeae63ea52cdeda56da39d2cdf30436ee079f58..5c172641d20c53a6ffe1b343e38df3f2414caa6f 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -424,8 +424,8 @@ public class ModelConfig {
// Forestry parameters
public static final String CONVERSION_COST_FILE = DATA_DIR + File.separator + "conversion_costs.csv";
- public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 10.0); // $/tC-eq
- public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 50.0); // $/tC-eq
+ public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 0.02); // $1000/tC-eq
+ public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.05); // $1000/tC-eq
public static final int FOREST_LOCKIN_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 30); // cannot convert forest after planting for this many years
}
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 4bc0ff980da86e33a03025d54151d3c084fc2100..e069b666428a79b0daf58aef6a11a39e40074235 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -664,16 +664,16 @@ public class ModelMain {
/** Get carbon flux data */
private CarbonFluxRasterSet getCarbonFluxData(Timestep timestep) {
CarbonFluxRasterSet carbonFluxData = new CarbonFluxRasterSet(desiredProjection);
+ LogWriter.println("Reading carbon flux data");
new CarbonFluxReader(carbonFluxData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + "cflux_net.out");
- LogWriter.println("Carbon flux data read");
return carbonFluxData;
}
/** Get carbon wood yield data */
private WoodYieldRasterSet getWoodYieldData(Timestep timestep) {
WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(desiredProjection);
+ LogWriter.println("Reading wood yield data");
new WoodYieldReader(woodYieldData).getRasterDataFromFile(timestep.getYearSubDir(ModelConfig.FOREST_DIR) + File.separator + "wood_yield.out");
- LogWriter.println("Wood yield data read");
return woodYieldData;
}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 6837681363c9402acb676eac7c8c5710c9016e66..55a3a546b7a36babeb9f5b4a571532df08d897c8 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -200,9 +200,10 @@ public class CountryAgent extends AbstractCountryAgent {
changeDown = changeUp = allowedImportChange * maxOfProdOrSupply;
}
+ // avoid getting stuck at 0
if (CropType.ENERGY_CROPS.equals(crop)) {
double croplandArea = LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.CROPLAND);
- double ecMaxExport = gen2EcIncrease * croplandArea / 1500.0 * 2;
+ double ecMaxExport = gen2EcIncrease * croplandArea / 1500.0 * 2; // twice the share
if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
baseTrade = 0.0;
ecMaxExport = 0.0;
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 623413cd800453ee24c48c65df824e98e9423d25..4d03e0d910c8bbf05d45e8733c1c5cd67b6a5964 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -110,7 +110,6 @@ public class GamsLocationOptimiser {
GAMSParameter previousExportAmountP = inDB.addParameter("previousExportAmount", 1);
GAMSParameter landP = inDB.addParameter("suitableLandArea", 1);
- GAMSParameter agriExpansionCostP = inDB.addParameter("agriExpansionCost", 1);
GAMSParameter seedAndWasteRateP = inDB.addParameter("seedAndWasteRate", 1);
GAMSParameter prevLandCoverP = inDB.addParameter("previousLandCoverArea", 2);
@@ -126,12 +125,7 @@ public class GamsLocationOptimiser {
double suitableLand = landUseItem.getSuitableArea();
totalSuitable += suitableLand;
if (DEBUG) LogWriter.println(String.format(" %d %15s,\t %.3f", locationId, "suitableLand", suitableLand));
- setGamsParamValueTruncate(landP.addRecord(locString), suitableLand, 5);
-
- double agriExpansionCost = ModelConfig.AGRI_EXPANSION_COST_BASE + landUseItem.getForestManagedFraction() * ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST;
- if (DEBUG) LogWriter.println(String.format(" %d %15s,\t %.3f", locationId, "agriExpansionCost", agriExpansionCost));
- setGamsParamValue(agriExpansionCostP.addRecord(Integer.toString(locationId)), agriExpansionCost, 3);
-
+ setGamsParamValueTruncate(landP.addRecord(locString), suitableLand, 5);
for (CropType cropType : CropType.getNonMeatTypes()) {
Vector<String> v = new Vector<String>();
@@ -327,11 +321,7 @@ public class GamsLocationOptimiser {
double otherIntCost = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.OTHER_INTENSITY_COST : ModelConfig.getParameter(Parameter.OTHER_INTENSITY_COST, inputData.getTimestep().getYear());
- LogWriter.print("\n");
- addScalar(inDB, "cropIncCost", ModelConfig.CROP_INCREASE_COST, 3);
- addScalar(inDB, "cropDecCost", ModelConfig.CROP_DECREASE_COST, 3);
- addScalar(inDB, "pastureDecCost", ModelConfig.PASTURE_DECREASE_COST, 3);
- addScalar(inDB, "pastureIncCost", ModelConfig.PASTURE_INCREASE_COST, 3);
+ LogWriter.print("\n");
addScalar(inDB, "meatEfficency", meatEff, 3);
addScalar(inDB, "fertiliserUnitCost", fertCost, 3);
addScalar(inDB, "otherICost",otherIntCost, 3);
@@ -364,8 +354,7 @@ public class GamsLocationOptimiser {
v.add(locString);
setGamsParamValue(cFluxRateP.addRecord(v), cFlux.getCarbonFlux(prevLC, newLC), 3);
}
- }
-
+ }
}
// Wood yields
@@ -384,8 +373,7 @@ public class GamsLocationOptimiser {
v.add(locString);
setGamsParamValue(woodYieldP.addRecord(v), wYield.getWoodYield(prevLC, newLC), 3);
}
- }
-
+ }
}
// Minimum land covers
diff --git a/src/ac/ed/lurg/landuse/WoodYieldReader.java b/src/ac/ed/lurg/landuse/WoodYieldReader.java
index 15f77c00caffaa78397d8a739da72cd2c91c39a8..f775fc144dae152e4173aab55431f1a37343167f 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldReader.java
+++ b/src/ac/ed/lurg/landuse/WoodYieldReader.java
@@ -10,7 +10,7 @@ import ac.sac.raster.RasterSet;
public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem> {
private static final int MIN_COLS = 27;
- private static final double conversionFactor = 10.0; // convert kgC/m2 to MtC/Mha
+ private static final double conversionFactor = 10.0; // convert kgC/m2 to tC/ha
public WoodYieldReader(RasterSet<WoodYieldItem> woodYield) {
super("[ |\t]+", MIN_COLS, woodYield);