From 03105ab6e3aa37d5056fe6e2561ae197d524ad00 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Mon, 17 May 2021 16:00:22 +0100
Subject: [PATCH] Reverted GAMS code to earlier implementation.

Added NATURAL_REGROWTH_TIME parameter.
---
 GAMS/IntExtOpt.gms                            | 86 ++++++++-----------
 src/ac/ed/lurg/InternationalMarket.java       | 28 ++++--
 src/ac/ed/lurg/ModelConfig.java               |  4 +-
 src/ac/ed/lurg/country/ForestManager.java     | 20 +++--
 .../country/gams/GamsLocationOptimiser.java   | 10 ++-
 5 files changed, 82 insertions(+), 66 deletions(-)

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index ecc99819..8938e8fe 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -118,14 +118,15 @@ $gdxin
  baseCost('pasture') = 0.04;
  otherIntCost('pasture') = 0.8 + otherICost;
 
- SCALAR forestExpansionMaxRate / 0.05 /;
- SCALAR restrictedDeforestationCost might fail to optimise if set too high / 100 /;
- unrestrictedArea(land_cover, location) = previousLandCoverArea(land_cover, location) - minimumLandCoverArea(land_cover, location);
+ SCALAR forestExpansionMaxRate / 0.2 /;
  woodYield('carbonForest', 'naturalNew', location) = woodYield('carbonForest', 'timberForest', location);
  woodYield('natural', 'naturalNew', location) = woodYield('natural', 'pasture', location);
  conversionCost(land_cover, 'naturalNew') = conversionCost(land_cover, 'natural');
  conversionCost('natural', 'naturalNew') = 0.05;
 
+ SCALAR loggingCost / 0.05 /;
+
+
  VARIABLES
        cropArea(crop, location)             total area for crops - Mha
        fertI(crop, location)              fertilizer intensity for each crop - factor between 0 and 1
@@ -140,11 +141,6 @@ $gdxin
        net_supply(crop)                   supply after exports and feed
        landCoverArea(land_cover, location)  land cover area in Mha
        landCoverChange(land_cover_before, land_cover_after, location) land cover change
-       restrictedLandCoverArea(land_cover, location)
-       restrictedLandCoverChange(land_cover, land_cover_after, location)
-       prematureDeforestationCost(location)
-       totalLandCoverArea(land_cover, location)
-       totalLandCoverChange(land_cover, land_cover_after, location)
        woodHarvest(location)              total wood harvested
        carbonFlux(location)               total carbon flux
        timberHarvest(location)
@@ -156,9 +152,7 @@ $gdxin
 
  POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
                    totalFeedDM,
-                   prematureDeforestationCost,
-                   landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost,
-                   restrictedLandCoverArea, restrictedLandCoverChange, totalLandCoverArea, totalLandCoverChange;
+                   landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost;
 
 * POSITIVE VARIABLE A;
 
@@ -187,14 +181,9 @@ $gdxin
        PASTURE_LAND_COVER_CALC(location)  pasture area (land cover) equals pasture area (land use)
        LAND_COVER_CHANGE_CALC(land_cover, location)
        LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
-       RESTRICTED_LAND_COVER_CHANGE_CALC(land_cover, location)
-       RESTRICTED_LAND_COVER_CHANGE_CONSTRAINT(land_cover, location)
+       MINIMUM_LAND_COVER_CONSTRAINT(land_cover, location)
        FOREST_EXPANSION_CONSTRAINT(forest, location)
        FOREST_CONTRACTION_CONSTRAINT(forest, location)
-       PREMATURE_DEFORESTATION_COST_CALC(location)
-       TOTAL_LAND_COVER_CALC(land_cover, location)
-       TOTAL_LAND_COVER_CHANGE_CALC(land_cover, land_cover_after, location)
-*       TOTAL_LAND_COVER_CONSTRAINT(location)
        WOOD_HARVEST_CALC(location)                          calc wood harvested
        CARBON_FLUX_CALC(location)                           calc carbon flux
        TIMBER_HARVEST_CALC(location)
@@ -238,10 +227,6 @@ $gdxin
 
  SETASIDE_AREA_CALC(location) .. cropArea('setaside', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) * setAsideRate;
 
- TOTAL_LAND_COVER_CALC(land_cover, location) .. totalLandCoverArea(land_cover, location) =E= landCoverArea(land_cover, location) + restrictedLandCoverArea(land_cover, location);
-
-* TOTAL_LAND_COVER_CONSTRAINT(location) .. sum(land_cover, totalLandCoverArea(land_cover, location)) =E= sum(land_cover, previousLandCoverArea(land_cover, location));
-
 * 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));
 
@@ -252,57 +237,54 @@ $gdxin
 
  IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) * (1.0 - unhandledCropRate) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * cropArea(crop, location));
 
- CROPLAND_LAND_COVER_CALC(location) .. totalLandCoverArea('cropland', location) =E= sum(crop_less_pasture, cropArea(crop_less_pasture, location)) / (1.0 - unhandledCropRate);
+ 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) .. totalLandCoverArea('pasture', location) =E= cropArea('pasture', location);
+ PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
 
 
- LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= unrestrictedArea(land_cover, location) +
+ LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= previousLandCoverArea(land_cover, location) +
                                                                                          sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) -
                                                                                          sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
 
- RESTRICTED_LAND_COVER_CHANGE_CALC(land_cover, location) .. restrictedLandCoverArea(land_cover, location) =E= minimumLandCoverArea(land_cover, location) +
-                                                                                         sum(land_cover_before, restrictedLandCoverChange(land_cover_before, land_cover, location)) -
-                                                                                         sum(land_cover_after, restrictedLandCoverChange(land_cover, land_cover_after, location));
-
- TOTAL_LAND_COVER_CHANGE_CALC(land_cover, land_cover_after, location) .. totalLandCoverChange(land_cover, land_cover_after, location) =E= landCoverChange(land_cover, land_cover_after, location) +
-                                                                                                                                          restrictedLandCoverChange(land_cover, land_cover_after, location);
-
 * 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= unrestrictedArea(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);
+ LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= previousLandCoverArea(land_cover, location);
 
- NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) .. totalLandCoverChange(exc_natural, 'natural', location) =E= 0;
+* 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);
 
- FOREST_EXPANSION_CONSTRAINT(forest, location) .. totalLandCoverArea(forest, location) =L= previousLandCoverArea(forest, location) * (1 + forestExpansionMaxRate);
- FOREST_CONTRACTION_CONSTRAINT(forest, location) .. totalLandCoverArea(forest, location) =G= previousLandCoverArea(forest, location) * (1 - forestExpansionMaxRate);
+ NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) .. landCoverChange(exc_natural, 'natural', location) =E= 0;
 
- PREMATURE_DEFORESTATION_COST_CALC(location) .. prematureDeforestationCost(location) =E= sum((forest, land_cover)$[not sameAs(forest, land_cover)], restrictedLandCoverChange(forest, land_cover, location)) * restrictedDeforestationCost;
+ 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;
 
 * 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));
 
 * Timber from timberForest rotation
- TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E=  (sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location))) * timberForestYield(location);
+ TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E=  [sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location)) - previousLandCoverArea('timberForest', location)] * timberForestYield(location);
 
 * Future timber from newly planted timberForest
- TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * timberYield(land_cover, 'timberForest', location));
+ TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(land_cover$[not sameAs(land_cover, 'timberForest')], landCoverChange(land_cover, 'timberForest', location) * timberYield(land_cover, 'timberForest', location)) +
+                                                                                 [landCoverChange('timberForest', 'timberForest', location) - minimumLandCoverArea('timberForest', location)] * timberYield('timberForest', 'timberForest', location);
 
- CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(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));
 
 * Subtracting cost of not fully grown timberForest reserved under minimumLandCoverArea
- CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), totalLandCoverChange(land_cover_before, land_cover_after, location) * (conversionCost(land_cover_before, land_cover_after) + infrastructureCost(location)));
+ CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after)$[not (sameAs(land_cover_before, 'timberForest') and sameAs(land_cover_after, 'timberForest'))], landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after)) +
+                                                                 (landCoverChange('timberForest', 'timberForest', location) - minimumLandCoverArea('timberForest', location)) * conversionCost('timberForest', 'timberForest') +
+                                                                 sum(exc_natural, landCoverChange('natural', exc_natural, location) * infrastructureCost(location));
 
  COST_EQ .. total_cost =E=
          (
               (   SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) ) * domesticPriceMarkup +
 
-               SUM(location, carbonFlux(location)) * carbonPrice - SUM(location, woodHarvest(location) + timberForwardContract(location)) * woodPrice +
+               SUM(location, carbonFlux(location)) * carbonPrice - SUM(location, woodHarvest(location)) * (woodPrice - loggingCost) -
+
+               SUM(location, timberForwardContract(location)) * (woodPrice - loggingCost) +
 
                SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
 
-               SUM(location, prematureDeforestationCost(location)) + SUM(location, totalConversionCost(location))
+               SUM(location, totalConversionCost(location))
          );
 
  MODEL LAND_USE /ALL/ ;
@@ -310,13 +292,21 @@ $gdxin
  irrigI.L(crop, location) = previousIrrigIntensity(crop, location);
  otherIntensity.L(crop, location) = previousOtherIntensity(crop, location);
  cropArea.L(crop, location) = previousCropArea(crop, location);
- landCoverArea.L(land_cover, location) = unrestrictedArea(land_cover, location);
- restrictedLandCoverArea.L(land_cover, location) = minimumLandCoverArea(land_cover, location);
- totalLandCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, location);
+ landCoverArea.L(land_cover, location) = previousLandCoverArea(land_cover, location);
+ landCoverChange.L(land_cover, land_cover, location) = previousLandCoverArea(land_cover, location);
  ruminantFeed.L(feed_crop) = previousRuminantFeed(feed_crop);
  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) - minimumLandCoverArea('timberForest', location)] * woodYield('timberForest', 'timberForest', location);
+ woodHarvest.L(location) = 0;
+ carbonFlux.L(location) = 0;
+ totalConversionCost.L(location) = [previousLandCoverArea('timberForest', location) - minimumLandCoverArea('timberForest', location)] * conversionCost('timberForest', 'timberForest');
+
+ LAND_USE.OptFile = 1;
+
+* display landCoverChange.L
 
  SOLVE LAND_USE USING NLP MINIMIZING total_cost;
 
@@ -355,7 +345,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) = sum((land_cover)$[not sameAs(forest, land_cover)], restrictedLandCoverChange.L(forest, land_cover, location));
+ reservedAreaDeforested(forest, location) = 0;
  netCarbonFlux = SUM(location, carbonFlux.L(location));
  totalTimberProduction = SUM(location, woodHarvest.L(location) + timberHarvest.L(location));
 
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index 028b0a23..0d08b814 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -159,13 +159,27 @@ public class InternationalMarket {
 			outputFile.newLine();
 		}
 		// Carbon price
-		StringBuffer sbData = new StringBuffer();
-		sbData.append(String.format("%d,%s", timestep.getYear(), "carbon"));
-		sbData.append(String.format(",%.1f,%.1f", carbonPrice.getImportAmount(), carbonPrice.getExportsAfterTransportLosses()));
-		sbData.append(String.format(",%.8f,%.3f", carbonPrice.getExportPrice(), carbonPrice.getStockLevel()));
-
-		outputFile.write(sbData.toString());
-		outputFile.newLine();
+		{
+			StringBuffer sbData = new StringBuffer();
+			sbData.append(String.format("%d,%s", timestep.getYear(), "carbon"));
+			sbData.append(String.format(",%.1f,%.1f", carbonPrice.getImportAmount(), carbonPrice.getExportsAfterTransportLosses()));
+			sbData.append(String.format(",%.8f,%.3f", carbonPrice.getExportPrice(), carbonPrice.getStockLevel()));
+	
+			outputFile.write(sbData.toString());
+			outputFile.newLine();
+		}
+		
+		// Timber price
+		{
+			StringBuffer sbData = new StringBuffer();
+			sbData.append(String.format("%d,%s", timestep.getYear(), "timber"));
+			sbData.append(String.format(",%.1f,%.1f", timberPrice.getImportAmount(), timberPrice.getExportsAfterTransportLosses()));
+			sbData.append(String.format(",%.8f,%.3f", timberPrice.getExportPrice(), timberPrice.getStockLevel()));
+	
+			outputFile.write(sbData.toString());
+			outputFile.newLine();
+		}
+		
 	}
 
 	public void serializeGlobalPrices() {
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 59f788e7..d589d2aa 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -431,6 +431,6 @@ public class ModelConfig {
 	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_LOCKIN_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 30); // cannot convert forest after planting for this many years
-	
+	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);
 }
diff --git a/src/ac/ed/lurg/country/ForestManager.java b/src/ac/ed/lurg/country/ForestManager.java
index 6e0b6e45..c8b3656d 100644
--- a/src/ac/ed/lurg/country/ForestManager.java
+++ b/src/ac/ed/lurg/country/ForestManager.java
@@ -92,8 +92,8 @@ public class ForestManager {
 		
 		// Assume unmanaged forest is fully mature (max wood yield) and other natural is new (no wood yield).
 		int otherNaturalPseudoYear = ModelConfig.BASE_YEAR;
-		int forestYearPlanted = ModelConfig.BASE_YEAR - ModelConfig.FOREST_LOCKIN_PERIOD;
-		int forestAge = ModelConfig.FOREST_LOCKIN_PERIOD;
+		int forestYearPlanted = ModelConfig.BASE_YEAR - ModelConfig.FOREST_ROTATION_PERIOD;
+		int forestAge = ModelConfig.FOREST_ROTATION_PERIOD;
 		int otherNaturalAge = 0;
 		
 		/*
@@ -137,7 +137,7 @@ public class ForestManager {
 	// Area of forest not fully grown
 	public DoubleMap<Integer, LandCoverType, Double> getReservedAreaForTimestep(Timestep timestep) {
 		DoubleMap<Integer, LandCoverType, Double> minimumLandCoverForTimestep = new DoubleMap<Integer, LandCoverType, Double>();
-		double minYear = timestep.getYear() - ModelConfig.FOREST_LOCKIN_PERIOD;
+		double minYear = timestep.getYear() - ModelConfig.FOREST_ROTATION_PERIOD;
 		for (ForestRecord record : forestHistory) {
 			if (record.getYearPlanted() > minYear && record.getForestType().isManagedForest()) {
 				minimumLandCoverForTimestep.addTo(record.getLocation(), record.getForestType(), record.getArea());
@@ -150,7 +150,7 @@ public class ForestManager {
 	public DoubleMap<Integer, LandCoverType, Double> getTimberYieldForTimestep(Timestep timestep) {
 		DoubleMap<Integer, LandCoverType, Double> yieldTimesArea = new DoubleMap<Integer, LandCoverType, Double>();
 		DoubleMap<Integer, LandCoverType, Double> area = new DoubleMap<Integer, LandCoverType, Double>();
-		double yearPlanted = timestep.getYear() - ModelConfig.FOREST_LOCKIN_PERIOD;
+		double yearPlanted = timestep.getYear() - ModelConfig.FOREST_ROTATION_PERIOD;
 		for (ForestRecord record : forestHistory) {
 			if (record.getYearPlanted() == yearPlanted && record.getForestType().isManagedForest()) {
 				yieldTimesArea.addTo(record.getLocation(), record.getForestType(), record.getPotentialYield() * record.getArea());
@@ -174,7 +174,7 @@ public class ForestManager {
 	
 	// Calculates wood yield factor for natural areas based on age
 	private double naturalWoodYieldFunction(double age) {
-		return Math.min(age, ModelConfig.FOREST_LOCKIN_PERIOD) / ModelConfig.FOREST_LOCKIN_PERIOD;
+		return Math.min(age, ModelConfig.NATURAL_REGROWTH_TIME) / ModelConfig.NATURAL_REGROWTH_TIME;
 	}
 	
 	/*
@@ -353,6 +353,14 @@ public class ForestManager {
 			naturalArea.put(locId, currentArea + newArea);
 			naturalMeanAge.put(locId, newAge);
 		}
+		
+		// Increment natural age
+		for (Entry<Integer, Double> entry : naturalMeanAge.entrySet()) {
+			int locId = entry.getKey();
+			double currentAge = entry.getValue();
+			naturalMeanAge.put(locId, currentAge + ModelConfig.TIMESTEP_SIZE);
+		}
+		
 
 		/*
 		for (ForestRecord record: forestHistory) {
@@ -376,7 +384,7 @@ public class ForestManager {
 					yearsFound.add(item.getYearPlanted());
 				}
 				Collections.sort(yearsFound, Collections.reverseOrder()); // descending order, remove area from newest forest first
-				int minYear = timestep.getYear() - ModelConfig.FOREST_LOCKIN_PERIOD;
+				int minYear = timestep.getYear() - ModelConfig.FOREST_ROTATION_PERIOD;
 				
 				for (int year : yearsFound) {
 					if (year > minYear) {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 22c099e9..c1644b84 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -178,7 +178,7 @@ public class GamsLocationOptimiser {
 				Vector<String> v = new Vector<String>();
 				v.add(lc.getName());
 				v.add(locString);
-				setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(lc), 5);
+				setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(lc), 4);
 			}
 			
 			// Infrastructure cost
@@ -448,7 +448,11 @@ public class GamsLocationOptimiser {
 				double previousCover = inputData.getPreviousLandUse().get(locationId).getUnprotectedLandCoverArea(landCover);
 				// Correction for floating point errors
 				double minCoverCorrected = Math.min(minCover, previousCover);
-				setGamsParamValueTruncate(minimumLandCoverP.addRecord(v), minCoverCorrected, 5);
+				
+				if (previousCover - minCoverCorrected <= 1E-5) 
+					minCoverCorrected = previousCover;
+				
+				setGamsParamValueTruncate(minimumLandCoverP.addRecord(v), minCoverCorrected, 4);
 			}
 		}
 		
@@ -596,7 +600,7 @@ public class GamsLocationOptimiser {
 				totalCropArea+totalPastureArea, totalCropArea, totalPastureArea));
 			
 		// Land cover change
-		GAMSVariable varLandCoverChange = outDB.getVariable("totalLandCoverChange");
+		GAMSVariable varLandCoverChange = outDB.getVariable("landCoverChange");
 		
 		TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges = new TripleMap<Integer, LandCoverType, LandCoverType, Double>();
 		DoubleMap<Integer, LandCoverType, Double> totalArea = new DoubleMap<Integer, LandCoverType, Double>();
-- 
GitLab