From dd03d90fa2f10b025fe86efc2a7a4a7c80f4a6fe Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Mon, 31 May 2021 09:20:50 +0100
Subject: [PATCH] Fixed issue with timber forest yields not updating.
---
GAMS/IntExtOpt.gms | 159 +++++++++---------
data/timber_demand.csv | 42 ++---
src/ac/ed/lurg/InternationalMarket.java | 2 +-
src/ac/ed/lurg/ModelConfig.java | 8 +-
src/ac/ed/lurg/country/CountryAgent.java | 9 +-
src/ac/ed/lurg/country/ForestManager.java | 16 +-
.../country/gams/GamsLocationOptimiser.java | 26 +--
.../country/gams/GamsRasterOptimiser.java | 22 +--
src/ac/ed/lurg/landuse/WoodYieldItem.java | 40 ++---
src/ac/ed/lurg/landuse/WoodYieldReader.java | 24 +--
10 files changed, 176 insertions(+), 172 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 4e144983..e73089b8 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -10,12 +10,10 @@
SET import_crop(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops/;
SET land_cover / cropland, pasture, timberForest, carbonForest, natural, naturalNew /;
- SET forest(land_cover) / timberForest, carbonForest, natural /;
+ SET forested(land_cover) / timberForest, carbonForest, natural /;
SET natural_forest(land_cover) / carbonForest, natural /;
- SET non_forest(land_cover) / cropland, pasture /;
- SET non_timber_forest_before(land_cover) / cropland, pasture, carbonForest, natural, naturalNew /;
- SET exc_natural(land_cover) / cropland, pasture, timberForest, carbonForest, naturalNew /;
- ALIAS(non_timber_forest_before, non_timber_forest_after);
+ SET non_timber_forest(land_cover) / cropland, pasture, carbonForest, natural, naturalNew /;
+ SET exc_natural(land_cover) / cropland, pasture, timberForest, naturalNew /;
ALIAS (land_cover, land_cover_before);
ALIAS (land_cover, land_cover_after);
@@ -59,11 +57,10 @@
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 in Mha;
-* PARAMETER unrestrictedArea(land_cover, location)
+
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;
- PARAMETER timberYield(land_cover_before, land_cover_after, location)
+ PARAMETER woodYieldLuc(land_cover_before, land_cover_after, location) wood yield - MtC-eq per Mha;
+ PARAMETER woodYieldRota(land_cover_before, land_cover_after, location)
SCALAR carbonPrice price of carbon - $1000 per tonne;
SCALAR woodPrice price of wood and timber - $1000 per tC-eq;
PARAMETER timberForestYield(location);
@@ -84,8 +81,8 @@ $load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousO
$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, woodYield, carbonPrice, woodPrice
-$load conversionCost, timberForestYield, timberYield, infrastructureCost, netPresentValueFactor, forestRotationPeriod, timberMaxNetImport, timberMinNetImport, loggingCost
+$load previousLandCoverArea, carbonFluxRate, woodYieldLuc, carbonPrice, woodPrice
+$load conversionCost, timberForestYield, woodYieldRota, infrastructureCost, netPresentValueFactor, forestRotationPeriod, timberMaxNetImport, timberMinNetImport, loggingCost
$gdxin
SCALAR delta "use to smooth power function see 7.5 www.gams.com/dd/docs/solversconopt.pdf" / 0.00000000001 /;
@@ -124,14 +121,21 @@ $gdxin
baseCost('pasture') = 0.04;
otherIntCost('pasture') = 0.8 + otherICost;
- SCALAR forestExpansionMaxRate / 1 /;
- woodYield('carbonForest', 'naturalNew', location) = woodYield('carbonForest', 'timberForest', location);
- woodYield('natural', 'naturalNew', location) = woodYield('natural', 'pasture', location);
+ SCALAR forestExpansionMaxRate / 0.0175 /;
+ woodYieldLuc('carbonForest', 'naturalNew', location) = woodYieldLuc('carbonForest', 'pasture', location);
+ woodYieldLuc('natural', 'naturalNew', location) = woodYieldLuc('natural', 'pasture', location);
conversionCost(land_cover, 'naturalNew') = conversionCost(land_cover, 'natural');
conversionCost('natural', 'naturalNew') = conversionCost('carbonForest', 'naturalNew');
+ woodYieldLuc('carbonForest', land_cover, location) = 0;
+
+ carbonFluxRate('natural', 'naturalNew', location) = carbonFluxRate('natural', 'pasture', location);
+ carbonFluxRate('carbonForest', 'naturalNew', location) = carbonFluxRate('natural', 'pasture', location);
+ carbonFluxRate('pasture', 'naturalNew', location) = carbonFluxRate('pasture', 'natural', location);
+ carbonFluxRate('cropland', 'naturalNew', location) = carbonFluxRate('cropland', 'natural', location);
+
VARIABLES
- cropArea(crop, location) total area for crops - Mha
+ cropArea(crop, location) total area for crops - 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)
@@ -142,22 +146,21 @@ $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
+ totalFeedDM total feed dry matter
landCoverArea(land_cover, location) land cover area in Mha
- landCoverChange(land_cover_before, land_cover_after, location) land cover change
- woodHarvest(location) total wood harvested
- carbonFlux(location) total carbon flux
- timberHarvest(location)
- timberForwardContract(location)
- totalConversionCost(location)
- totalFeedDM
- totalTimberHarvest
- timberSold
+ landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha
+ 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
+ totalWoodHarvest total wood harvested from all activities - Mt C-eq
+ woodExported total wood sold - Mt C-eq
+ carbonFlux(location) total carbon flux - Mt C
* A "artificial variable for debugging https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
total_cost total cost of domestic supply including net imports;
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
- totalFeedDM,
- landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost, totalTimberHarvest, timberSold;
+ totalFeedDM, landCoverArea, landCoverChange, totalConversionCost, woodHarvestLuc, woodHarvestRota,
+ totalWoodHarvest, woodExported
* POSITIVE VARIABLE A;
@@ -175,32 +178,33 @@ $gdxin
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) constraint on maximum irrigation intensity
MAX_OTHER_INTENSITY_CONSTRAINT(crop, location)
SETASIDE_AREA_CALC(location)
-* TOTAL_LAND_CHANGE_CONSTRAINT(location) constraint on suitable land use
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 constraint to not import pasture
- PASTURE_EXPORT_CONSTRAINT
+ PASTURE_EXPORT_CONSTRAINT constraint to not export pasture
IRRIGATION_CONSTRAINT(location) constraint on water usage
NET_SUPPLY_EQ(crop) calc net supply for crops
- 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)
- LAND_COVER_CHANGE_CALC(land_cover, 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)
+ LAND_COVER_CHANGE_CALC(land_cover, location) calc land cover change
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
-* MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location)
- FOREST_EXPANSION_CONSTRAINT(forest, location)
- FOREST_CONTRACTION_CONSTRAINT(forest, location)
- WOOD_HARVEST_CALC(location) calc wood harvested
- CARBON_FLUX_CALC(location) calc carbon flux
-* TIMBER_HARVEST_CALC(location)
- TIMBER_FORWARD_CONTRACT_CALC(location)
- CONVERSION_COST(location)
- NATURAL_CONVERSION_CONSTRAINT(exc_natural, location)
-* TIMBER_MAX_EXPORT_CONSTRAINT
- TIMBER_MIN_EXPORT_CONSTRAINT
- TOTAL_TIMBER_HARVEST_CALC
- TIMBER_SOLD_CONSTRAINT_A
- TIMBER_SOLD_CONSTRAINT_B
- TIMBER_SOLD_CONSTRAINT_C
+ NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) can't create primary natural area
+* FOREST_EXPANSION_CONSTRAINT(forested, location) constraint on forest expansion
+* FOREST_CONTRACTION_CONSTRAINT(forested, location) constraint on forest contraction
+ NATURAL_CONTRACTION_CONSTRAINT(location)
+ CONVERSION_COST(location) cost of land cover conversion
+
+ WOOD_HARVEST_LUC_CALC(location) calc wood harvested from LUC
+ WOOD_HARVEST_ROTA_CALC(location) calc wood harvested from timber forest rotation
+ TOTAL_WOOD_HARVEST_CALC total wood harvested
+* WOOD_MIN_EXPORT_CONSTRAINT constraint on minimum wood export
+ WOOD_MAX_EXPORT_CONSTRAINT constraint on maximum wood export
+* WOOD_MAX_TRADE_CONSTRAINT
+ WOOD_ROTA_EXPORT_CONSTRAINT must export all wood produced from timber forest rotation
+* WOOD_ROTA_TO_LUC_CONSTRAINT
+
+ CARBON_FLUX_CALC(location) calc carbon flux
COST_EQ total cost objective function;
UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E= ( baseCost(crop) +
@@ -238,9 +242,6 @@ $gdxin
SETASIDE_AREA_CALC(location) .. cropArea('setaside', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) * setAsideRate;
-* Needs inequality due to floating point errors but in theory should be equal
-* TOTAL_LAND_CHANGE_CONSTRAINT(location) .. suitableLandArea(location) =G= sum(land_cover, totalLandCoverArea(land_cover, location));
-
MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
PASTURE_IMPORT_CONSTRAINT .. importAmount('pasture') =E= 0;
@@ -260,41 +261,39 @@ $gdxin
* landCoverChange(A, A, location) defined as unchanged land cover
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= previousLandCoverArea(land_cover, location);
-* RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, restrictedLandCoverChange(land_cover, land_cover_after, location)) =E= minimumLandCoverArea(land_cover, location);
-* MINIMUM_LAND_COVER_CONSTRAINT(forest, location) .. landCoverChange(forest, forest, location) =G= minimumLandCoverArea(forest, location);
-
NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) .. landCoverChange(exc_natural, 'natural', location) =E= 0;
- FOREST_EXPANSION_CONSTRAINT(forest, location) .. landCoverArea(forest, location) =L= previousLandCoverArea(forest, location) + suitableLandArea(location) * forestExpansionMaxRate;
- FOREST_CONTRACTION_CONSTRAINT(forest, location) .. landCoverArea(forest, location) =G= previousLandCoverArea(forest, location) - suitableLandArea(location) * forestExpansionMaxRate;
+* FOREST_EXPANSION_CONSTRAINT(forested, location) .. landCoverArea(forested, location) =L= previousLandCoverArea(forested, location) + suitableLandArea(location) * forestExpansionMaxRate;
+* FOREST_CONTRACTION_CONSTRAINT(forested, location) .. landCoverArea(forested, location) =G= previousLandCoverArea(forested, location) - suitableLandArea(location) * forestExpansionMaxRate;
-* Timber from land cover change
- WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYield(land_cover_before, land_cover_after, location));
+ NATURAL_CONTRACTION_CONSTRAINT(location) .. landCoverArea('natural', location) =G= previousLandCoverArea('natural', location) * (1 - forestExpansionMaxRate);
-* Timber from timberForest rotation
-* TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E= [sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location))] * timberForestYield(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)) +
+ sum(exc_natural, landCoverChange('natural', exc_natural, location) * infrastructureCost(location));
-* Potential future timber from newly planted timberForest plus from existing timberForest
- TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(non_timber_forest_before, landCoverChange(non_timber_forest_before, 'timberForest', location) * timberYield(non_timber_forest_before, 'timberForest', location)) +
- landCoverChange('timberForest', 'timberForest', location) * timberForestYield(location);
+* Wood from land cover change
+ WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYieldLuc(land_cover_before, land_cover_after, location));
- 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));
+* Wood from timberForest rotation = harvest from new forest + existing forest
+ WOOD_HARVEST_ROTA_CALC(location) .. woodHarvestRota(location) =E= [ sum(non_timber_forest, landCoverChange(non_timber_forest, 'timberForest', location) * woodYieldRota(non_timber_forest, 'timberForest', location)) +
+ landCoverChange('timberForest', 'timberForest', location) * timberForestYield(location)
+ ] / forestRotationPeriod;
-* 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)) +
- sum(exc_natural, landCoverChange('natural', exc_natural, location) * infrastructureCost(location));
+ TOTAL_WOOD_HARVEST_CALC .. totalWoodHarvest =E= sum(location, woodHarvestLuc(location) + woodHarvestRota(location));
-* TIMBER_MAX_EXPORT_CONSTRAINT .. 0 - sum(location, woodHarvest(location) + timberForwardContract(location) / forestRotationPeriod) =L= timberMaxNetImport;
+* WOOD_MIN_EXPORT_CONSTRAINT .. 0 - woodExported =L= timberMaxNetImport;
- TIMBER_MIN_EXPORT_CONSTRAINT .. 0 - sum(location, woodHarvest(location) + timberForwardContract(location) / forestRotationPeriod) =G= timberMinNetImport;
+ WOOD_MAX_EXPORT_CONSTRAINT .. woodExported =L= totalWoodHarvest;
- TOTAL_TIMBER_HARVEST_CALC .. totalTimberHarvest =E= sum(location, woodHarvest(location) + [timberForwardContract(location) / forestRotationPeriod]);
+* WOOD_MAX_TRADE_CONSTRAINT .. 0 - woodExported =G= timberMinNetImport;
-* Equivalent of timberSold =E= MIN(totalTimberHarvest, -timberMinNetImport)
- TIMBER_SOLD_CONSTRAINT_A .. timberSold =L= totalTimberHarvest;
- TIMBER_SOLD_CONSTRAINT_B .. timberSold =L= -timberMinNetImport;
* Must sell timber producted in timberForest
- TIMBER_SOLD_CONSTRAINT_C .. timberSold =G= sum(location, timberForwardContract(location)) / forestRotationPeriod;
+ WOOD_ROTA_EXPORT_CONSTRAINT .. woodExported =G= sum(location, woodHarvestRota(location));
+
+* WOOD_ROTA_TO_LUC_CONSTRAINT .. sum(location, woodHarvestLuc(location)) =L= sum(location, woodHarvestRota(location)) * 2;
+
+ 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=
(
@@ -304,11 +303,9 @@ $gdxin
SUM(location, totalConversionCost(location)) -
- timberSold * (woodPrice - loggingCost) -
+ woodExported * woodPrice +
SUM(location, carbonFlux(location)) * carbonPrice
-
-
);
MODEL LAND_USE /ALL/ ;
@@ -322,11 +319,9 @@ $gdxin
monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop);
importAmount.L(all_types) = previousImportAmount(all_types);
exportAmount.L(all_types) = previousExportAmount(all_types);
-* timberHarvest.L(location) = [previousLandCoverArea('timberForest', location) - minimumLandCoverArea('timberForest', location)] * timberForestYield(location);
- timberForwardContract.L(location) = previousLandCoverArea('timberForest', location) * timberForestYield(location);
- woodHarvest.L(location) = 0;
- carbonFlux.L(location) = 0;
-* totalConversionCost.L(location) = [previousLandCoverArea('timberForest', location) - minimumLandCoverArea('timberForest', location)] * conversionCost('timberForest', 'timberForest');
+ woodHarvestRota.L(location) = previousLandCoverArea('timberForest', location) * timberForestYield(location) / forestRotationPeriod;
+ totalWoodHarvest.L = sum(location, woodHarvestRota.L(location));
+ totalConversionCost.L(location) = previousLandCoverArea('timberForest', location) * conversionCost('timberForest', 'timberForest');
LAND_USE.OptFile = 1;
@@ -346,9 +341,7 @@ $gdxin
parameter netImportCost(all_types);
parameter feedCostRate(feed_crop);
parameter productionShock(all_types);
- parameter reservedAreaDeforested(forest, location);
scalar netCarbonFlux;
- scalar totalTimberProduction;
* Production quantities based on smaller area (before unhandledCropArea adjustment applied)
totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location));
@@ -369,9 +362,7 @@ $gdxin
netImportCost(import_crop) = importAmount.l(import_crop) * importPrices(import_crop) - exportAmount.l(import_crop) * exportPrices(import_crop);
netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
- reservedAreaDeforested(forest, location) = 0;
netCarbonFlux = SUM(location, carbonFlux.L(location));
- totalTimberProduction = timberSold.L;
Scalar totalCostsLessLU;
diff --git a/data/timber_demand.csv b/data/timber_demand.csv
index 89ab3786..47e52ee5 100644
--- a/data/timber_demand.csv
+++ b/data/timber_demand.csv
@@ -1,22 +1,22 @@
year,value
-2000,70
-2005,70
-2010,70
-2015,70
-2020,70
-2025,70
-2030,70
-2035,70
-2040,70
-2045,70
-2050,70
-2055,70
-2060,70
-2065,70
-2070,70
-2075,70
-2080,70
-2085,70
-2090,70
-2095,70
-2100,70
+2000,1048
+2005,1059
+2010,1070
+2015,1080
+2020,1091
+2025,1102
+2030,1113
+2035,1124
+2040,1135
+2045,1147
+2050,1158
+2055,1170
+2060,1181
+2065,1193
+2070,1205
+2075,1217
+2080,1229
+2085,1242
+2090,1254
+2095,1267
+2100,1279
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index 0d08b814..cbd77efa 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -42,7 +42,7 @@ public class InternationalMarket {
}
carbonPrice = GlobalPrice.createInitial(ModelConfig.CARBON_PRICE, 0);
- timberPrice = GlobalPrice.createInitial(ModelConfig.WOOD_PRICE, 0);
+ timberPrice = GlobalPrice.createInitial(ModelConfig.WOOD_PRICE, 1000);
}
else {
List<Object> deserializedPrices = deserializeGlobalPrice();
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 9aeb957e..97881150 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -430,10 +430,10 @@ public class ModelConfig {
public static final double CARBON_PRICE = getDoubleProperty("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 WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.2); // $1000/tC-eq
- public static final int FOREST_ROTATION_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 30); // cannot convert forest after planting for this many years
- public static final int NATURAL_REGROWTH_TIME = getIntProperty("NATURAL_REGROWTH_TIME", 100); // years
+ public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.05); // $1000/tC-eq
+ public static final int FOREST_ROTATION_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 20); // 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 DISCOUNT_RATE = getDoubleProperty("DISCOUNT_RATE", 0.05);
public static final double TREE_PLANTING_COST = getDoubleProperty("TREE_PLANTING_COST", 1.0); // $1000/ha
- public static final double TREE_LOGGING_COST = getDoubleProperty("TREE_PLANTING_COST", 0.05); // $1000/tC-eq
+ public static final double TREE_LOGGING_COST = getDoubleProperty("TREE_LOGGING_COST", 0.02); // $1000/tC-eq
}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 73d8216a..277042c3 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -233,7 +233,14 @@ public class CountryAgent extends AbstractCountryAgent {
// Timber import/export constraints
TradeConstraint timberTradeConstraint;
{
- double baseTrade = -getTimberProduction(); // net imports negative as all exported
+ double baseTrade;
+ if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.isInitialTimestep()) {
+ double forestedAndNaturalArea = LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.NATURAL) +
+ LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.TIMBER_FOREST);
+ baseTrade = -timberDemandIncrease * forestedAndNaturalArea / 4000;
+ } else {
+ baseTrade = -getTimberProduction(); // net imports negative as all exported
+ }
double countryArea = LandUseItem.getTotalLandArea(previousGamsRasterOutput.getLandUses().values());
double changeUp = 0.0;
double changeDown = 0.0;
diff --git a/src/ac/ed/lurg/country/ForestManager.java b/src/ac/ed/lurg/country/ForestManager.java
index 4b4c02d8..5d8b4bf2 100644
--- a/src/ac/ed/lurg/country/ForestManager.java
+++ b/src/ac/ed/lurg/country/ForestManager.java
@@ -66,7 +66,7 @@ public class ForestManager {
WoodYieldItem wYItem = initWoodYieldData.get(key);
double timberYieldThisTime;
// TODO Assuming forest was also forest before
- timberYieldThisTime = (wYItem != null) ? initWoodYieldData.get(key).getTimberYield(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST) : 0.0;
+ timberYieldThisTime = (wYItem != null) ? initWoodYieldData.get(key).getWoodYieldRota(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST) : 0.0;
double naturalForestFraction = (luItem.getLandCoverArea(LandCoverType.NATURAL) > 0) ? luItem.getInitialForestedNaturalArea() / luItem.getLandCoverArea(LandCoverType.NATURAL) : 0;
unmanagedForestArea.put(clusterId, unmanagedForestAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.NATURAL) * naturalForestFraction);
otherNaturalArea.put(clusterId, otherNaturalAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.NATURAL) * (1 - naturalForestFraction)); // Natural is other natural + unmanaged forest so need to subtract the latter
@@ -143,8 +143,6 @@ public class ForestManager {
public void updateForestHistory(ArrayList<LandCoverChangeItem> gamsLandCoverChanges, Map<Integer, Double> baseTimberYield) {
- LogWriter.println("foo " + naturalArea.get(13));
-
// Update timber forest yield
for (Integer locId : timberForestYield.keySet()) {
double area = timberForestArea.get(locId);
@@ -156,7 +154,7 @@ public class ForestManager {
// Removed timber forest
for (LandCoverChangeItem item : gamsLandCoverChanges) {
- if (!item.getToLandCover().equals(LandCoverType.TIMBER_FOREST.toString()) && item.getFromLandCover().equals(LandCoverType.TIMBER_FOREST.toString())) {
+ if (!item.getToLandCover().equals(LandCoverType.TIMBER_FOREST.getName()) && item.getFromLandCover().equals(LandCoverType.TIMBER_FOREST.getName())) {
int locId = item.getLocation();
double currentArea = timberForestArea.get(locId);
double removedArea = item.getArea();
@@ -165,10 +163,14 @@ public class ForestManager {
}
// New timber forest
for (LandCoverChangeItem item : gamsLandCoverChanges) {
- if (item.getToLandCover().equals(LandCoverType.TIMBER_FOREST.toString()) && !item.getFromLandCover().equals(LandCoverType.TIMBER_FOREST.toString())) {
+ if (item.getToLandCover().equals(LandCoverType.TIMBER_FOREST.getName()) && !item.getFromLandCover().equals(LandCoverType.TIMBER_FOREST.getName())) {
int locId = item.getLocation();
- double currentYield = timberForestYield.get(locId);
- double currentArea = timberForestArea.get(locId);
+ double currentYield = 0;
+ double currentArea = 0;
+ if (timberForestArea.containsKey(locId)) {
+ currentYield = timberForestYield.get(locId);
+ currentArea = timberForestArea.get(locId);
+ }
double newArea = item.getArea();
double newYield = item.getTimberYield();
double updatedYield = (currentYield * currentArea + newYield * newArea) / (currentArea + newArea); // area weighted average
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 56ee2d25..0d44762f 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -379,15 +379,15 @@ public class GamsLocationOptimiser {
double factor = entry.getValue();
for (LandCoverType toLC : LandCoverType.getConvertibleTypes()) {
WoodYieldItem item = woodYieldData.get(locId);
- if (item.checkForWoodKeys(LandCoverType.NATURAL, toLC)) {
- double maxYield = item.getWoodYield(LandCoverType.NATURAL, toLC);
- item.setWoodYield(LandCoverType.NATURAL, toLC, maxYield * factor);
+ if (item.checkForWoodYieldLucKeys(LandCoverType.NATURAL, toLC)) {
+ double maxYield = item.getWoodYieldLuc(LandCoverType.NATURAL, toLC);
+ item.setWoodYieldLuc(LandCoverType.NATURAL, toLC, maxYield * factor);
}
}
}
- GAMSParameter woodYieldP = inDB.addParameter("woodYield", 3);
- GAMSParameter timberYieldP = inDB.addParameter("timberYield", 3);
+ GAMSParameter woodYieldLucP = inDB.addParameter("woodYieldLuc", 3);
+ GAMSParameter woodYieldRotaP = inDB.addParameter("woodYieldRota", 3);
for (Entry<Integer, ? extends WoodYieldItem> entry : woodYieldData.entrySet()) {
Integer locationId = entry.getKey();
@@ -400,10 +400,10 @@ public class GamsLocationOptimiser {
v.add(prevLC.getName());
v.add(newLC.getName());
v.add(locString);
- if (wYield.checkForWoodKeys(prevLC, newLC))
- setGamsParamValue(woodYieldP.addRecord(v), wYield.getWoodYield(prevLC, newLC), 3);
- if (wYield.checkForTimberKeys(prevLC, newLC))
- setGamsParamValue(timberYieldP.addRecord(v), wYield.getTimberYield(prevLC, newLC), 3);
+ if (wYield.checkForWoodYieldLucKeys(prevLC, newLC))
+ setGamsParamValue(woodYieldLucP.addRecord(v), wYield.getWoodYieldLuc(prevLC, newLC), 3);
+ if (wYield.checkForWoodYieldRotaKeys(prevLC, newLC))
+ setGamsParamValue(woodYieldRotaP.addRecord(v), wYield.getWoodYieldRota(prevLC, newLC), 3);
}
}
}
@@ -491,9 +491,13 @@ public class GamsLocationOptimiser {
if (toLc.equals(fromLc))
cost = 0;
+ // Forest rotation cost
if (toLc.isManagedForest() && fromLc.isManagedForest() & toLc.equals(fromLc))
cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD;
+ if (!toLc.isNatural() && fromLc.isNatural())
+ cost += ModelConfig.AGRI_EXPANSION_COST_BASE;
+
Vector<String> v = new Vector<String>();
v.add(fromLc.getName());
v.add(toLc.getName());
@@ -669,7 +673,7 @@ public class GamsLocationOptimiser {
}
// Data for updating natural areas and forest
- double potentialTimberYield = (toLc.equals(LandCoverType.TIMBER_FOREST)) ? inputData.getWoodYields().get(locId).getTimberYield(fromLc, toLc) : 0;
+ double potentialTimberYield = (toLc.equals(LandCoverType.TIMBER_FOREST)) ? inputData.getWoodYields().get(locId).getWoodYieldRota(fromLc, toLc) : 0;
gamsLandCoverChanges.add(new LandCoverChangeItem(locId, fromGamsLc, toGamsLc, change, potentialTimberYield));
}
@@ -677,7 +681,7 @@ public class GamsLocationOptimiser {
double netCarbonFlux = outDB.getParameter("netCarbonFlux").getFirstRecord().getValue();
// Timber harvest
- double totalTimberProduction = outDB.getParameter("totalTimberProduction").getFirstRecord().getValue();
+ double totalTimberProduction = outDB.getVariable("woodExported").getFirstRecord().getLevel();
GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, gamsLandCoverChanges, netCarbonFlux, totalTimberProduction);
return results;
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index bc488761..c5fc95b8 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -55,7 +55,7 @@ public class GamsRasterOptimiser {
Map<Integer, Double> baseTimberYield = new HashMap<Integer, Double>(); // For updating timber yields of existing forest
for (int locId : gamsInput.getWoodYields().keySet()) {
- baseTimberYield.put(locId, gamsInput.getWoodYields().get(locId).getTimberYield(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST));
+ baseTimberYield.put(locId, gamsInput.getWoodYields().get(locId).getWoodYieldRota(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST));
}
return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getGamsLandCoverChanges(),
@@ -334,8 +334,8 @@ public class GamsRasterOptimiser {
for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
double woodYield;
if (woodYieldItem != null) {
- if (woodYieldItem.checkForWoodKeys(prevLC, newLC)) {
- woodYield = woodYieldItem.getWoodYield(prevLC, newLC);
+ if (woodYieldItem.checkForWoodYieldLucKeys(prevLC, newLC)) {
+ woodYield = woodYieldItem.getWoodYieldLuc(prevLC, newLC);
} else {
continue;
}
@@ -343,10 +343,10 @@ public class GamsRasterOptimiser {
woodYield = 0.0; // if missing data assume 0
}
- if (!aggWYield.checkForWoodKeys(prevLC, newLC)) { // if not seen yet
- aggWYield.setWoodYield(prevLC, newLC, woodYield);
+ if (!aggWYield.checkForWoodYieldLucKeys(prevLC, newLC)) { // if not seen yet
+ aggWYield.setWoodYieldLuc(prevLC, newLC, woodYield);
} else {
- aggWYield.setWoodYield(prevLC, newLC, aggregateMean(aggWYield.getWoodYield(prevLC, newLC), suitableAreaSoFar,
+ aggWYield.setWoodYieldLuc(prevLC, newLC, aggregateMean(aggWYield.getWoodYieldLuc(prevLC, newLC), suitableAreaSoFar,
woodYield, suitableAreaThisTime));
}
}
@@ -357,8 +357,8 @@ public class GamsRasterOptimiser {
for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
double timberYield;
if (woodYieldItem != null) {
- if (woodYieldItem.checkForTimberKeys(prevLC, newLC)) {
- timberYield = woodYieldItem.getTimberYield(prevLC, newLC);
+ if (woodYieldItem.checkForWoodYieldRotaKeys(prevLC, newLC)) {
+ timberYield = woodYieldItem.getWoodYieldRota(prevLC, newLC);
} else {
continue;
}
@@ -366,10 +366,10 @@ public class GamsRasterOptimiser {
timberYield = 0.0; // if missing data assume 0
}
- if (!aggWYield.checkForTimberKeys(prevLC, newLC)) { // if not seen yet
- aggWYield.setTimberYield(prevLC, newLC, timberYield);
+ if (!aggWYield.checkForWoodYieldRotaKeys(prevLC, newLC)) { // if not seen yet
+ aggWYield.setWoodYieldRota(prevLC, newLC, timberYield);
} else {
- aggWYield.setTimberYield(prevLC, newLC, aggregateMean(aggWYield.getTimberYield(prevLC, newLC), suitableAreaSoFar,
+ aggWYield.setWoodYieldRota(prevLC, newLC, aggregateMean(aggWYield.getWoodYieldRota(prevLC, newLC), suitableAreaSoFar,
timberYield, suitableAreaThisTime));
}
}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldItem.java b/src/ac/ed/lurg/landuse/WoodYieldItem.java
index f259c170..3ddf14b9 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldItem.java
+++ b/src/ac/ed/lurg/landuse/WoodYieldItem.java
@@ -7,48 +7,48 @@ import ac.ed.lurg.types.LandCoverType;
import ac.sac.raster.RasterItem;
public class WoodYieldItem implements RasterItem {
- Map<LandCoverType, Map<LandCoverType, Double>> woodYields = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
- Map<LandCoverType, Map<LandCoverType, Double>> timberYields = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
+ Map<LandCoverType, Map<LandCoverType, Double>> woodYieldsLuc = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
+ Map<LandCoverType, Map<LandCoverType, Double>> woodYieldsRota = new HashMap<LandCoverType, Map<LandCoverType, Double>>();
- public void setWoodYield(LandCoverType previousLandCover, LandCoverType newLandCover, double woodYield) {
- if (woodYields.containsKey(previousLandCover)) {
- woodYields.get(previousLandCover).put(newLandCover, woodYield);
+ public void setWoodYieldLuc(LandCoverType previousLandCover, LandCoverType newLandCover, double woodYield) {
+ if (woodYieldsLuc.containsKey(previousLandCover)) {
+ woodYieldsLuc.get(previousLandCover).put(newLandCover, woodYield);
} else {
Map<LandCoverType, Double> cfMap = new HashMap<LandCoverType, Double>();
cfMap.put(newLandCover, woodYield);
- woodYields.put(previousLandCover, cfMap);
+ woodYieldsLuc.put(previousLandCover, cfMap);
}
}
- public void setTimberYield(LandCoverType previousLandCover, LandCoverType newLandCover, double woodYield) {
- if (timberYields.containsKey(previousLandCover)) {
- timberYields.get(previousLandCover).put(newLandCover, woodYield);
+ public void setWoodYieldRota(LandCoverType previousLandCover, LandCoverType newLandCover, double woodYield) {
+ if (woodYieldsRota.containsKey(previousLandCover)) {
+ woodYieldsRota.get(previousLandCover).put(newLandCover, woodYield);
} else {
Map<LandCoverType, Double> cfMap = new HashMap<LandCoverType, Double>();
cfMap.put(newLandCover, woodYield);
- timberYields.put(previousLandCover, cfMap);
+ woodYieldsRota.put(previousLandCover, cfMap);
}
}
- public double getWoodYield(LandCoverType previousLandCover, LandCoverType newLandCover) {
- return woodYields.get(previousLandCover).get(newLandCover);
+ public double getWoodYieldLuc(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ return woodYieldsLuc.get(previousLandCover).get(newLandCover);
}
- public double getTimberYield(LandCoverType previousLandCover, LandCoverType newLandCover) {
- return timberYields.get(previousLandCover).get(newLandCover);
+ public double getWoodYieldRota(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ return woodYieldsRota.get(previousLandCover).get(newLandCover);
}
- public boolean checkForWoodKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
- if (woodYields.containsKey(previousLandCover)) {
- return woodYields.get(previousLandCover).containsKey(newLandCover);
+ public boolean checkForWoodYieldLucKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ if (woodYieldsLuc.containsKey(previousLandCover)) {
+ return woodYieldsLuc.get(previousLandCover).containsKey(newLandCover);
} else {
return false;
}
}
- public boolean checkForTimberKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
- if (timberYields.containsKey(previousLandCover)) {
- return timberYields.get(previousLandCover).containsKey(newLandCover);
+ public boolean checkForWoodYieldRotaKeys(LandCoverType previousLandCover, LandCoverType newLandCover) {
+ if (woodYieldsRota.containsKey(previousLandCover)) {
+ return woodYieldsRota.get(previousLandCover).containsKey(newLandCover);
} else {
return false;
}
diff --git a/src/ac/ed/lurg/landuse/WoodYieldReader.java b/src/ac/ed/lurg/landuse/WoodYieldReader.java
index cfbaa2aa..3df4fec3 100644
--- a/src/ac/ed/lurg/landuse/WoodYieldReader.java
+++ b/src/ac/ed/lurg/landuse/WoodYieldReader.java
@@ -17,19 +17,19 @@ public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem>
}
protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
- item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "forC_to_crop") * conversionFactor);
- item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT") * conversionFactor);
- item.setWoodYield(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "forC_to_past") * conversionFactor);
- item.setWoodYield(LandCoverType.NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "ntrl_to_crop") * conversionFactor);
- item.setWoodYield(LandCoverType.NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "ntrl_to_forC") * conversionFactor);
- item.setWoodYield(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT") * conversionFactor);
- item.setWoodYield(LandCoverType.NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "ntrl_to_past") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.CARBON_FOREST, LandCoverType.CROPLAND, getValueForCol(rowValues, "forC_to_crop") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.CARBON_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "forC_to_past") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "ntrl_to_crop") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "ntrl_to_forC") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT") * conversionFactor);
+ item.setWoodYieldLuc(LandCoverType.NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "ntrl_to_past") * conversionFactor);
- item.setTimberYield(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "crop_to_forT_rot") * conversionFactor);
- item.setTimberYield(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "past_to_forT_rot") * conversionFactor);
- item.setTimberYield(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forT_to_forT_rot") * conversionFactor);
- item.setTimberYield(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT_rot") * conversionFactor);
- item.setTimberYield(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT_rot") * conversionFactor);
+ item.setWoodYieldRota(LandCoverType.CROPLAND, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "crop_to_forT_rot") * conversionFactor);
+ item.setWoodYieldRota(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "past_to_forT_rot") * conversionFactor);
+ item.setWoodYieldRota(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forT_to_forT_rot") * conversionFactor);
+ item.setWoodYieldRota(LandCoverType.CARBON_FOREST, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "forC_to_forT_rot") * conversionFactor);
+ item.setWoodYieldRota(LandCoverType.NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT_rot") * conversionFactor);
/*
item.setWoodYield(LandCoverType.CROPLAND, LandCoverType.NATURAL, getValueForCol(rowValues, "crop_to_ntrl") * conversionFactor);
--
GitLab