diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index d30d5af775d4bfe7c415b1bd34bafb76e15cbeaa..c960e1ed6742b819eaeacd14ab238c18f4da169f 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -70,6 +70,8 @@
SCALAR forestRotationPeriod;
SCALAR woodMaxNetImport;
SCALAR woodMinNetImport;
+ SCALAR loggingCostBase;
+ SCALAR naturalWoodMaxHarvestRate;
*$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
$gdxin %gdxincname%
@@ -79,7 +81,7 @@ $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
$load previousLandCoverArea, carbonFluxRate, woodYieldLuc, carbonPrice, woodPrice
-$load conversionCost, woodYieldRota, infrastructureCost, forestRotationPeriod, woodMaxNetImport, woodMinNetImport
+$load conversionCost, woodYieldRota, infrastructureCost, forestRotationPeriod, woodMaxNetImport, woodMinNetImport, loggingCostBase, naturalWoodMaxHarvestRate
$gdxin
SCALAR delta "use to smooth power function see 7.5 www.gams.com/dd/docs/solversconopt.pdf" / 0.00000000001 /;
@@ -118,12 +120,6 @@ $gdxin
baseCost('pasture') = 0.04;
otherIntCost('pasture') = 0.8 + otherICost;
- SCALAR naturalWoodHarvestMaxRate avaiable wood max harvest proportion / 0.02 /;
-
- woodYieldLuc('carbonForest', land_cover, location) = 0;
-
- SCALAR forestRotationCost / 0.073 /;
-
VARIABLES
cropArea(crop, location) total area for crops - Mha
fertI(crop, location) fertilizer intensity for each crop - factor between 0 and 1
@@ -142,6 +138,7 @@ $gdxin
totalConversionCost(location) land cover conversion cost - $1000 per ha
woodHarvestLuc(location) wood harvested due to land cover change - Mt C-eq
woodHarvestRota(location) wood harvested from timber forest rotation - Mt C-eq
+ forestForwardContract(location) present value of timber forest
totalWoodHarvest(location) total wood harvested from all activities - Mt C-eq
woodExported total wood sold - Mt C-eq
carbonFlux(location) total carbon flux - Mt C
@@ -151,7 +148,7 @@ $gdxin
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
totalFeedDM, landCoverArea, landCoverChange, totalConversionCost, woodHarvestLuc, woodHarvestRota,
- totalWoodHarvest, woodExported;
+ totalWoodHarvest, woodExported, loggingCost;
* POSITIVE VARIABLE A;
@@ -184,6 +181,7 @@ $gdxin
WOOD_HARVEST_LUC_CALC(location) calc wood harvested from LUC
WOOD_HARVEST_ROTA_CALC(location) calc wood harvested from timber forest rotation
+ FOREST_FORWARD_CONTRACT_CALC(location)
TOTAL_WOOD_HARVEST_CALC(location) total wood harvested
WOOD_MAX_EXPORT_CONSTRAINT constraint on maximum wood export
WOOD_MIN_TRADE_CONSTRAINT constraint on minimum wood export
@@ -251,38 +249,38 @@ $gdxin
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));
* Wood from land cover change
- WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =L= woodYieldLuc('natural', 'pasture', location) * landCoverChange('natural', 'natural', location) * naturalWoodHarvestMaxRate;
+ WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =L= woodYieldLuc('natural', 'pasture', location) * landCoverChange('natural', 'natural', location) * naturalWoodMaxHarvestRate;
* Wood from timberForest rotation and new timberForest
WOOD_HARVEST_ROTA_CALC(location) .. woodHarvestRota(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location)) / forestRotationPeriod;
+ FOREST_FORWARD_CONTRACT_CALC(location) .. forestForwardContract(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location));
+
TOTAL_WOOD_HARVEST_CALC(location) .. totalWoodHarvest(location) =E= woodHarvestLuc(location) + woodHarvestRota(location);
WOOD_MAX_EXPORT_CONSTRAINT .. woodExported =L= sum(location, totalWoodHarvest(location));
- WOOD_MIN_TRADE_CONSTRAINT .. 0 - woodExported =L= woodMaxNetImport;
+ WOOD_MIN_TRADE_CONSTRAINT .. 0 - woodExported =L= woodMaxNetImport / 5;
- WOOD_MAX_TRADE_CONSTRAINT .. 0 - woodExported =G= woodMinNetImport;
+ WOOD_MAX_TRADE_CONSTRAINT .. 0 - woodExported =G= woodMinNetImport * 5;
* Must sell timber producted in timberForest
WOOD_ROTA_EXPORT_CONSTRAINT .. woodExported =G= sum(location, woodHarvestRota(location));
- LOGGING_COST_CALC(location) .. loggingCost(location) =E= 0.03 + 0.03 * woodHarvestLuc(location) / [woodYieldLuc('natural', 'pasture', location) * landCoverChange('natural', 'natural', location) * naturalWoodHarvestMaxRate + delta];
+ LOGGING_COST_CALC(location) .. loggingCost(location) =E= loggingCostBase + 0.03 * woodHarvestLuc(location) / [woodYieldLuc('natural', 'pasture', location) * landCoverChange('natural', 'natural', location) * naturalWoodMaxHarvestRate + delta];
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));
COST_EQ .. total_cost =E=
(
- ( SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) ) * domesticPriceMarkup +
+ [ SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) ] * domesticPriceMarkup +
SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
SUM(location, totalConversionCost(location)) +
- [SUM((managed_forest, location), landCoverArea(managed_forest, location) * forestRotationCost) - woodExported * woodPrice] +
-
- SUM(location, loggingCost(location) * woodHarvestLuc(location)) +
+ SUM(location, (loggingCost(location) - woodPrice) * woodHarvestLuc(location) - woodPrice * forestForwardContract(location)) +
SUM(location, carbonFlux(location)) * carbonPrice
);
@@ -298,9 +296,9 @@ $gdxin
monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop);
importAmount.L(all_types) = previousImportAmount(all_types);
exportAmount.L(all_types) = previousExportAmount(all_types);
- totalConversionCost.L(location) = [previousLandCoverArea('timberForest', location) + previousLandCoverArea('carbonForest', location)] * forestRotationCost;
+ totalConversionCost.L(location) = sum(land_cover, landCoverChange.L(land_cover, land_cover, location) * conversionCost(land_cover, land_cover));
woodHarvestRota.L(location) = previousLandCoverArea('timberForest', location) * woodYieldRota('timberForest', 'timberForest', location) / forestRotationPeriod;
- woodHarvestLuc.L(location) = previousLandCoverArea('natural', location) * woodYieldLuc('natural', 'pasture', location) * naturalWoodHarvestMaxRate;
+ woodHarvestLuc.L(location) = previousLandCoverArea('natural', location) * woodYieldLuc('natural', 'pasture', location) * naturalWoodMaxHarvestRate;
totalWoodHarvest.L(location) = woodHarvestRota.L(location) + woodHarvestLuc.L(location);
woodExported.L = sum(location, totalWoodHarvest.L(location));
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index d09e1dc576a431035df6a3952bad7978be756b99..415f57975157102694726479643afac8b90b30e1 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -425,19 +425,18 @@ public class ModelConfig {
// Forestry and carbon market parameters
public static final String CONVERSION_COST_FILE = DATA_DIR + File.separator + "conversion_costs.csv"; // cost of converting from one land cover to another, $1000/ha
- public static final String SERIALIZED_CARBON_MARKET_FILE = CALIB_DIR + File.separator + "CarbonMarket.ser";
public static final String CARBON_DEMAND_FILENAME = getProperty("CARBON_DEMAND_FILENAME", "carbon_demand.csv");
public static final String CARBON_DEMAND_FILE = getProperty("CARBON_DEMAND_FILE", DATA_DIR + File.separator + CARBON_DEMAND_FILENAME);
public static final double INIT_CARBON_PRICE = getDoubleProperty("INIT_CARBON_PRICE", 0.02); // $1000/tC-eq
public static final String TIMBER_DEMAND_FILENAME = getProperty("TIMBER_DEMAND_FILENAME", "timber_demand.csv");
public static final String TIMBER_DEMAND_FILE = getProperty("TIMBER_DEMAND_FILE", DATA_DIR + File.separator + TIMBER_DEMAND_FILENAME);
- public static final double INIT_WOOD_PRICE = getDoubleProperty("INIT_WOOD_PRICE", 0.05); // $1000/tC-eq
- public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 0.0); // $1000/tC-eq
+ public static final double INIT_WOOD_PRICE = getDoubleProperty("INIT_WOOD_PRICE", 0.06); // $1000/tC-eq
+ public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 1000.0); // tC-eq
public static final int FOREST_ROTATION_PERIOD = getIntProperty("FOREST_ROTATION_PERIOD", 30); // cannot convert forest after planting for this many years
- public static final int NATURAL_REGROWTH_TIME = getIntProperty("NATURAL_REGROWTH_TIME", 50); // years
- public static final double NATURAL_REGROWTH_PARAM_P = getDoubleProperty("NATURAL_REGROWTH_PARAM_P", 1.0); // Chapman-Richards paramater. Controls gradient of growth curve
- public static final double NATURAL_REGROWTH_PARAM_K = -Math.log(1 - Math.pow(0.95, (1 / NATURAL_REGROWTH_PARAM_P))) / NATURAL_REGROWTH_TIME; // Calc Chapman-Richards paramater. 95% yield reached at NATURAL_REGROWTH_TIME
- public static final double DISCOUNT_RATE = getDoubleProperty("DISCOUNT_RATE", 0.05);
- public static final double TREE_PLANTING_COST = getDoubleProperty("TREE_PLANTING_COST", 1.2); // $1000/ha
- public static final double TREE_LOGGING_COST = getDoubleProperty("TREE_LOGGING_COST", 0.02); // $1000/tC-eq
+ public static final int NATURAL_REGROWTH_TIME = getIntProperty("NATURAL_REGROWTH_TIME", 60); // years
+ public static final double NATURAL_WOOD_MAX_HARVEST_RATE = getDoubleProperty("NATURAL_WOOD_MAX_HARVEST_RATE", 1.0 / NATURAL_REGROWTH_TIME); // $1000/ha
+ public static final double FOREST_MANAGEMENT_COST = getDoubleProperty("FOREST_MANAGEMENT_COST", 0.025); // $1000/ha
+ public static final double TREE_LOGGING_COST = getDoubleProperty("TREE_LOGGING_COST", 0.04); // $1000/tC-eq
+ public static final double MANAGED_FOREST_INCREASE_COST = getDoubleProperty("MANAGED_FOREST_INCREASE_COST", 10 * LAND_CHANGE_COST); // $1000/ha
+ public static final double MANAGED_FOREST_DECREASE_COST = getDoubleProperty("MANAGED_FOREST_DECREASE_COST", 0.05 * LAND_CHANGE_COST); // $1000/ha
}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 348c0abd53440f9a06d85d33ab500534ce4f0219..1289c2206d267c23017a8a7f76df78f2b19a10d4 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -239,9 +239,9 @@ public class CountryAgent extends AbstractCountryAgent {
double forestedAndNaturalArea = LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.NATURAL) +
LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.TIMBER_FOREST);
changeDown = -timberDemandIncrease * forestedAndNaturalArea / 4000;
- if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {changeDown = -timberDemandIncrease;}
+ if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {changeDown = -timberDemandIncrease;} // for debugging purposes only
} else {
- changeDown = Math.min(baseTrade * allowedImportChange, -(timberDemandIncrease * countryArea / 15000 * 2));
+ changeDown = Math.min(baseTrade * allowedImportChange, -(timberDemandIncrease * countryArea / 15000 * 2)); // allow export in case baseTrade is 0
}
changeUp = baseTrade * allowedImportChange;
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index ce790c677b6b1a747ea135047476ccf4f9079092..e6d63ea8cb23806c26d24a634230870744b4b03f 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -391,9 +391,40 @@ public class GamsLocationOptimiser {
// Land cover conversion cost
GAMSParameter conversionCostP = inDB.addParameter("conversionCost", 2);
-
+ /*
DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = inputData.getConversionCosts();
+ */
+
+ DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = new DoubleMap<LandCoverType, LandCoverType, Double>();
+
+ conversionCosts.put(LandCoverType.CROPLAND, LandCoverType.PASTURE, ModelConfig.CROP_DECREASE_COST + ModelConfig.PASTURE_INCREASE_COST);
+ conversionCosts.put(LandCoverType.CROPLAND, LandCoverType.NATURAL, ModelConfig.CROP_DECREASE_COST);
+ conversionCosts.put(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, ModelConfig.CROP_DECREASE_COST + ModelConfig.MANAGED_FOREST_INCREASE_COST);
+ conversionCosts.put(LandCoverType.CROPLAND, LandCoverType.CARBON_FOREST, ModelConfig.CROP_DECREASE_COST + ModelConfig.MANAGED_FOREST_INCREASE_COST);
+
+ conversionCosts.put(LandCoverType.PASTURE, LandCoverType.CROPLAND, ModelConfig.PASTURE_DECREASE_COST + ModelConfig.CROP_INCREASE_COST);
+ conversionCosts.put(LandCoverType.PASTURE, LandCoverType.NATURAL, ModelConfig.PASTURE_DECREASE_COST);
+ conversionCosts.put(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, ModelConfig.PASTURE_DECREASE_COST + ModelConfig.MANAGED_FOREST_INCREASE_COST);
+ conversionCosts.put(LandCoverType.PASTURE, LandCoverType.CARBON_FOREST, ModelConfig.PASTURE_DECREASE_COST + ModelConfig.MANAGED_FOREST_INCREASE_COST);
+
+ conversionCosts.put(LandCoverType.NATURAL, LandCoverType.CROPLAND, ModelConfig.AGRI_EXPANSION_COST_BASE + ModelConfig.CROP_INCREASE_COST);
+ conversionCosts.put(LandCoverType.NATURAL, LandCoverType.PASTURE, ModelConfig.AGRI_EXPANSION_COST_BASE + ModelConfig.PASTURE_INCREASE_COST);
+ conversionCosts.put(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, ModelConfig.MANAGED_FOREST_INCREASE_COST);
+ conversionCosts.put(LandCoverType.NATURAL, LandCoverType.CARBON_FOREST, ModelConfig.MANAGED_FOREST_INCREASE_COST);
+
+ conversionCosts.put(LandCoverType.TIMBER_FOREST, LandCoverType.CROPLAND, ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST + ModelConfig.CROP_INCREASE_COST);
+ conversionCosts.put(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE, ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST + ModelConfig.PASTURE_INCREASE_COST);
+ conversionCosts.put(LandCoverType.TIMBER_FOREST, LandCoverType.NATURAL, ModelConfig.MANAGED_FOREST_DECREASE_COST);
+ conversionCosts.put(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST, ModelConfig.MANAGED_FOREST_INCREASE_COST);
+ conversionCosts.put(LandCoverType.TIMBER_FOREST, LandCoverType.CARBON_FOREST, ModelConfig.FOREST_MANAGEMENT_COST * 0.05);
+
+ conversionCosts.put(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST + ModelConfig.CROP_INCREASE_COST);
+ conversionCosts.put(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST + ModelConfig.PASTURE_INCREASE_COST);
+ conversionCosts.put(LandCoverType.CARBON_FOREST, LandCoverType.NATURAL, ModelConfig.AGRI_EXPANSION_COST_BASE_MANAGED_FOREST);
+ conversionCosts.put(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, ModelConfig.FOREST_MANAGEMENT_COST);
+ conversionCosts.put(LandCoverType.CARBON_FOREST, LandCoverType.CARBON_FOREST, ModelConfig.FOREST_MANAGEMENT_COST);
+
for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap : conversionCosts.getMap().entrySet()) {
String fromName = fromMap.getKey().getName();
for (Entry<LandCoverType, Double> toMap : fromMap.getValue().entrySet()) {
@@ -405,78 +436,11 @@ public class GamsLocationOptimiser {
setGamsParamValue(conversionCostP.addRecord(v), cost, 5);
}
}
- /*
- for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
- for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
- double cost = 0;
-
- switch (fromLc) {
- case PASTURE:
- cost += ModelConfig.PASTURE_DECREASE_COST;
- break;
- case CROPLAND:
- cost += ModelConfig.CROP_DECREASE_COST;
- break;
- case TIMBER_FOREST:
- cost += ModelConfig.LAND_CHANGE_COST;
- break;
- case CARBON_FOREST:
- cost += ModelConfig.LAND_CHANGE_COST;
- break;
- case NATURAL:
- cost += ModelConfig.LAND_CHANGE_COST;
- default:
- break;
- }
-
- switch (toLc) {
- case PASTURE:
- cost += ModelConfig.PASTURE_INCREASE_COST;
- break;
- case CROPLAND:
- cost += ModelConfig.CROP_INCREASE_COST;
- break;
- case TIMBER_FOREST:
- cost += ModelConfig.TREE_PLANTING_COST;
- break;
- case CARBON_FOREST:
- cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD;
- break;
- case NATURAL:
- cost += 0;
- default:
- break;
- }
-
- if (toLc.equals(fromLc))
- cost = 0;
-
- if (toLc.equals(LandCoverType.NATURAL) && fromLc.isManagedForest())
- cost = 0;
-
- // Forest rotation cost
- if (toLc.isManagedForest() && fromLc.isManagedForest() & toLc.equals(fromLc))
- cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD;
-
- if (toLc.equals(LandCoverType.NATURAL) && fromLc.equals(LandCoverType.TIMBER_FOREST))
- cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD * 0.9;
-
- //if (toLc.equals(LandCoverType.TIMBER_FOREST))
- // cost = cost / ModelConfig.FOREST_ROTATION_PERIOD;
-
-
- cost = 0;
-
- Vector<String> v = new Vector<String>();
- v.add(fromLc.getName());
- v.add(toLc.getName());
- setGamsParamValue(conversionCostP.addRecord(v), cost, 5);
- }
- }
- */
addScalar(inDB, "forestRotationPeriod", ModelConfig.FOREST_ROTATION_PERIOD, 3);
- addScalar(inDB, "loggingCost", ModelConfig.TREE_LOGGING_COST, 5);
+ addScalar(inDB, "loggingCostBase", ModelConfig.TREE_LOGGING_COST, 5);
+ addScalar(inDB, "naturalWoodMaxHarvestRate", ModelConfig.NATURAL_WOOD_MAX_HARVEST_RATE, 5);
+
}
private void addScalar(GAMSDatabase gamsDb, String recordName, double val, int places) {