From 2bc68eff6df62780d026f6f9f27243edf0c3fffb Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Tue, 15 Feb 2022 14:29:07 +0000
Subject: [PATCH] Annualised wood harvest.
---
GAMS/IntExtOpt.gms | 42 ++++++++--------------
src/ac/ed/lurg/forestry/WoodYieldItem.java | 6 ++--
2 files changed, 16 insertions(+), 32 deletions(-)
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 39ad479a..f0da17b3 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -154,7 +154,7 @@ $gdxin
landCoverArea(land_cover, location) land cover area in Mha
landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha
totalConversionCost(location) land cover conversion cost - $1000 per ha
- newForestPotentialHarvest(location)
+ forestRotaHarvest(location)
naturalDeforestedArea(location)
woodHarvestNat(location)
woodHarvestLuc(location)
@@ -164,7 +164,6 @@ $gdxin
woodSold
carbonFlux(location) total carbon flux - Mt C
supply(all_types)
- npv
* 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;
@@ -202,7 +201,7 @@ $gdxin
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
CONVERSION_COST(location) cost of land cover conversion
- NEW_FOREST_POTENTIAL_HARVEST(location)
+ FOREST_ROTA_HARVEST_CALC(location)
WOOD_HARVEST_NATURAL_CALC(location)
WOOD_HARVEST_NATURAL_CONSTRAINT(location)
WOOD_HARVEST_LUC_CALC(location)
@@ -214,7 +213,7 @@ $gdxin
CARBON_FLUX_CALC(location) calc carbon flux
- NPV_CALC Net Present Value function;
+ COST_CALC Net Present Value function;
**************** Crop cost **************************
@@ -279,7 +278,7 @@ $gdxin
************* Forestry ***********************************
- NEW_FOREST_POTENTIAL_HARVEST(location) .. newForestPotentialHarvest(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location));
+ FOREST_ROTA_HARVEST_CALC(location) .. forestRotaHarvest(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location)) / forestRotationPeriod(location);
WOOD_HARVEST_NATURAL_CALC(location) .. woodHarvestNat(location) =E= naturalDeforestedArea(location) * woodYieldLuc('natural', 'natural', location);
@@ -287,7 +286,7 @@ $gdxin
WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((forested, exc_natural), landCoverChange(forested, exc_natural, location) * woodYieldLuc(forested, exc_natural, location));
- WOOD_SUPPLY_CALC .. woodSupply =E= woodFromRotation + sum(location, woodHarvestNat(location) + woodHarvestLuc(location));
+ WOOD_SUPPLY_CALC .. woodSupply =E= sum(location, forestRotaHarvest(location) + woodHarvestNat(location) + woodHarvestLuc(location));
WOOD_MARKET_CONSTRAINT .. woodSold =L= woodDemand + woodExported - woodImported;
@@ -305,26 +304,13 @@ $gdxin
************ Net Present Value ******************************
- NPV_CALC .. npv =E= {
- [
- sum(import_crop, exportPrices(import_crop) * supply(import_crop)) +
- sum(location, landCoverArea('natural', location)) * naturalAreaValue
- ] * exp(-discountRate) -
- sum((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop)))
- } / (1 - exp(-discountRate)) +
-
- sum(location,
- [
- (woodExportPrice - vegClearingCost) * newForestPotentialHarvest(location) * exp(-discountRate * forestRotationPeriod(location)) -
- landCoverArea('timberForest', location) * forestEstablishmentCost
- ] / (1 - exp(-discountRate * forestRotationPeriod(location)))
- ) +
- (woodExportPrice - vegClearingCost) * woodSold +
- sum(location, carbonFlux(location) * carbonPrice) -
- sum(import_crop, importPrices(import_crop) * importAmount(import_crop)) -
- sum(location, totalConversionCost(location)) -
- woodImportPrice * woodImported -
- sum(location, woodHarvestLuc(location) * vegClearingCost);
+ COST_CALC .. total_cost =E= sum((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) * domesticPriceMarkup +
+ sum(location, totalConversionCost(location)) +
+ sum(import_crop, importPrices(import_crop) * importAmount(import_crop)) +
+ woodImported * woodImportPrice -
+ sum(import_crop, exportPrices(import_crop) * exportAmount(import_crop)) -
+ (woodExportPrice - vegClearingCost) * woodSold;
+
@@ -344,8 +330,8 @@ $gdxin
* display landCoverChange.L
-* SOLVE LAND_USE USING NLP MINIMIZING total_cost;
- SOLVE LAND_USE USING NLP MAXIMIZING npv;
+ SOLVE LAND_USE USING NLP MINIMIZING total_cost;
+* SOLVE LAND_USE USING NLP MAXIMIZING npv;
display previousCropArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, cropArea.L;
display supply.l, demand, ruminantFeed.l, monogastricFeed.l, exportAmount.l;
diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java
index 5dcc34eb..e3da5d72 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldItem.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java
@@ -29,10 +29,6 @@ public class WoodYieldItem implements RasterItem {
LccKey key = entry.getKey();
Double[] yields = entry.getValue();
- if (!landUseTiles.containsKey(key.getFromLc())) {
- int foo = 1;
- }
-
LandCoverTile tiles = landUseTiles.get(key.getFromLc());
double totalArea = tiles.getTotalConvertibleArea(); // Assuming no harvest from protected areas
@@ -79,6 +75,7 @@ public class WoodYieldItem implements RasterItem {
currentRotationHarvest = 0.0;
} else {
currentRotationHarvest += timberTiles.getConvertibleArea(optimalRotation) * yieldAtRotation;
+ /*
timberTiles.resetAreaForAge(optimalRotation); // WARNING: side effect
for (int age=optimalRotation+1; age<=LandCoverTile.getMaxAgeBin(); age++) { // TODO better way of dealing with stands over rotation age?
@@ -86,6 +83,7 @@ public class WoodYieldItem implements RasterItem {
currentRotationHarvest += areaHarvested * yields[age];
timberTiles.resetAreaForAge(age, areaHarvested);
}
+ */
}
}
--
GitLab