From ac5d6df9f7abd1e09cf0e01071faf03525415728 Mon Sep 17 00:00:00 2001 From: s1924442 <b.arendarczyk@sms.ed.ac.uk> Date: Fri, 19 Mar 2021 18:48:52 +0000 Subject: [PATCH] Added land cover conversion cost. --- GAMS/IntExtOpt.gms | 21 +++++--- debug_config.properties | 11 ++-- src/ac/ed/lurg/ModelConfig.java | 1 + src/ac/ed/lurg/country/CountryAgent.java | 50 ++++++++--------- .../country/gams/GamsLocationOptimiser.java | 23 ++++++-- .../ed/lurg/country/gams/GamsRasterInput.java | 8 +-- .../country/gams/GamsRasterOptimiser.java | 21 -------- .../ed/lurg/landuse/ConversionCostReader.java | 53 +++++++++++++++++++ src/ac/ed/lurg/landuse/WoodYieldReader.java | 12 ++--- 9 files changed, 124 insertions(+), 76 deletions(-) create mode 100644 src/ac/ed/lurg/landuse/ConversionCostReader.java diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms index 4fa6d238..d0201ff9 100644 --- a/GAMS/IntExtOpt.gms +++ b/GAMS/IntExtOpt.gms @@ -67,6 +67,8 @@ SCALAR carbonPrice price of carbon - $ per tonne or million$ per Mt; SCALAR woodPrice price of wood and timber - $ per tonne or million$ per Mt; + 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 @@ -75,7 +77,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, woodYield, carbonPrice, woodPrice -$load minimumLandCoverArea +$load minimumLandCoverArea, conversionCost $gdxin *convert units from $/t to billion$/Mt @@ -118,9 +120,6 @@ $gdxin baseCost('pasture') = 0.04; otherIntCost('pasture') = 0.8 + otherICost; -*PARAMETER lcTransitionCost(lc_type_previous, lc_type, location) cost of land transitions including carbon cost and conversion cost; -*lcTransitionCost(lc_type_previous, lc_type, location) = carbonFluxRate(lc_type_previous, lc_type, location) * carbonPrice + agriExpansionCost(location); - SCALAR forestExpansionMaxRate / 0.01 /; VARIABLES @@ -147,6 +146,7 @@ $gdxin prematureDeforestationCost woodHarvest(location) total wood harvested carbonFlux(location) total carbon flux + totalConversionCost(location) totalFeedDM * A "artificial variable https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/" total_cost total cost of domestic supply including net imports; @@ -195,8 +195,10 @@ $gdxin DEFORESTED_AREA_CALC(forest, location) RESERVED_AREA_DEFORESTED_CALC(forest, location) PREMATURE_DEFORESTATION_COST_CALC + UNMANAGED_FOREST_CONSTRAINT(location) unmanagedForest cannot be created WOOD_HARVEST_CALC(location) calc wood harvested CARBON_FLUX_CALC(location) calc carbon flux + CONVERSION_COST(location) COST_EQ total cost objective function; UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E= ( baseCost(crop) + @@ -258,11 +260,12 @@ $gdxin PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location); -* landCoverChange(A, A, location) defined as unchanged land cover + 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)); +* 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); @@ -280,12 +283,16 @@ $gdxin sqrt(sqr(deforestedArea(forest, location) - [previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)] + delta)) + sqr(delta)] / 2; - PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, location), reservedAreaDeforested(forest, location) * 100); + PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, location), reservedAreaDeforested(forest, location) * woodYield('pasture', forest, location) * woodPrice); + + UNMANAGED_FOREST_CONSTRAINT(location) .. sum((land_cover)$[not sameAs(land_cover, 'unmanagedForest')], landCoverChange(land_cover, 'unmanagedForest', location)) =E= 0; 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)); 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)); + * 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); @@ -306,7 +313,7 @@ $gdxin SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) + - prematureDeforestationCost + prematureDeforestationCost + SUM(location, totalConversionCost(location)) ); MODEL LAND_USE /ALL/ ; diff --git a/debug_config.properties b/debug_config.properties index c161ff3a..7e7ff02c 100644 --- a/debug_config.properties +++ b/debug_config.properties @@ -11,12 +11,15 @@ TIMESTEP_SIZE=1 IS_CALIBRATION_RUN=true -GENERATE_NEW_YIELD_CLUSTERS=true +GENERATE_NEW_YIELD_CLUSTERS=false DEMAND_PRICE_IMPORT_AND_PROD_COST=false YIELD_FILENAME=yield.out -DEBUG_LIMIT_COUNTRIES=false -DEBUG_COUNTRY_NAME=Peru & Ecuador -GAMS_COUNTRY_TO_SAVE=Peru & Ecuador \ No newline at end of file +DEBUG_LIMIT_COUNTRIES=true +DEBUG_COUNTRY_NAME=Central America +GAMS_COUNTRY_TO_SAVE=Central America + +WOOD_PRICE=25 +CARBON_PRICE=10 \ No newline at end of file diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java index ef6786e0..ebeae63e 100755 --- a/src/ac/ed/lurg/ModelConfig.java +++ b/src/ac/ed/lurg/ModelConfig.java @@ -423,6 +423,7 @@ public class ModelConfig { public static final String TARGET_DIET_FILE = getProperty("TARGET_DIET_FILE", DATA_DIR + File.separator + "TargetDiet.txt"); // 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 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/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java index 421a98a3..68376813 100644 --- a/src/ac/ed/lurg/country/CountryAgent.java +++ b/src/ac/ed/lurg/country/CountryAgent.java @@ -228,28 +228,11 @@ public class CountryAgent extends AbstractCountryAgent { } } } - - // Logged timberForest area - to be converted to otherNatural - DoubleMap<Integer, LandCoverType, Double> loggedForest = new DoubleMap<Integer, LandCoverType, Double>(); //<Location, LandCoverType, Area> - - if (currentTimestep.getTimestep() != ModelConfig.START_TIMESTEP) { - for (Entry<Integer, DoubleMap<Integer, LandCoverType, Double>> yearMap : minimumLandCover.getMap().entrySet()) { - if (yearMap.getKey() <= currentTimestep.getYear() && yearMap.getKey() > currentTimestep.getPreviousYear()) { - for (Entry<Integer, Map<LandCoverType, Double>> locMap : yearMap.getValue().getMap().entrySet()) { - for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) { - loggedForest.put(locMap.getKey(), coverMap.getKey(), coverMap.getValue()); - } - } - } - } - } - - GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates); GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, - carbonFluxData, woodYieldData, minimumLandCoverForTimestep, loggedForest, countryLevelInputs); + carbonFluxData, woodYieldData, minimumLandCoverForTimestep, countryLevelInputs); return input; } @@ -374,23 +357,36 @@ public class CountryAgent extends AbstractCountryAgent { if (!reservedDeforested.getMap().isEmpty()) { for (Entry<Integer, Map<LandCoverType, Double>> locMap : reservedDeforested.getMap().entrySet()) { int locId = locMap.getKey(); + for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) { LandCoverType forestType = coverMap.getKey(); double deforestedArea = coverMap.getValue(); List<Integer> expiryYearsFound = new ArrayList<Integer>(minimumLandCover.getMap().keySet()) ; Collections.sort(expiryYearsFound); // ascending order, remove area from most mature forest - while (deforestedArea > 0) { - for (int year : expiryYearsFound) { - if (year > currentTimestep.getYear()) { - double currentMinimumArea = minimumLandCover.get(year, locId, forestType); - if (currentMinimumArea <= 0) - LogWriter.printlnError("currentMinimumArea is <=0"); - double areaToSubtract = Math.min(currentMinimumArea, deforestedArea); - minimumLandCover.put(year, locId, forestType, currentMinimumArea - areaToSubtract); - deforestedArea -= areaToSubtract; + + for (int year : expiryYearsFound) { + if (year > currentTimestep.getYear()) { + double currentMinimumArea = minimumLandCover.get(year, locId, forestType); + if (currentMinimumArea < 0.0) + LogWriter.printlnError("CountryAgent.updateMinimumLandCover(): currentMinimumArea is <=0"); + else if (currentMinimumArea == 0.0) { + continue; + } + + double areaToSubtract = Math.min(currentMinimumArea, deforestedArea); + if (currentMinimumArea - areaToSubtract < 0.0) { // floating point errors + areaToSubtract = currentMinimumArea; } + + minimumLandCover.put(year, locId, forestType, currentMinimumArea - areaToSubtract); + deforestedArea -= areaToSubtract; + if (deforestedArea < 1E-6) + break; } } + + if (deforestedArea >= 1E-6) + LogWriter.printlnError("CountryAgent.updateMinimumLandCover(): unable to subtract all deforested area"); } } } diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java index 815eb097..623413cd 100644 --- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java @@ -24,6 +24,7 @@ import ac.ed.lurg.ModelConfig; import ac.ed.lurg.country.CountryPrice; import ac.ed.lurg.country.TradeConstraint; import ac.ed.lurg.landuse.CarbonFluxItem; +import ac.ed.lurg.landuse.ConversionCostReader; import ac.ed.lurg.landuse.CropUsageData; import ac.ed.lurg.landuse.Intensity; import ac.ed.lurg.landuse.IrrigationItem; @@ -42,7 +43,7 @@ import ac.ed.lurg.yield.YieldResponsesItem; public class GamsLocationOptimiser { - private static final boolean DEBUG = true; + private static final boolean DEBUG = false; private GamsLocationInput inputData; @@ -388,7 +389,7 @@ public class GamsLocationOptimiser { } // Minimum land covers - GAMSParameter minimumLandCover = inDB.addParameter("minimumLandCoverArea", 2); + GAMSParameter minimumLandCoverP = inDB.addParameter("minimumLandCoverArea", 2); for (Entry<Integer, Map<LandCoverType, Double>> locMap : inputData.getMinimumLandCover().getMap().entrySet()) { Integer locationId = locMap.getKey(); @@ -403,11 +404,25 @@ public class GamsLocationOptimiser { double previousCover = inputData.getPreviousLandUse().get(locationId).getUnprotectedLandCoverArea(landCover); // Correction for floating point errors double minCoverCorrected = Math.min(minCover, previousCover); - setGamsParamValueTruncate(minimumLandCover.addRecord(v), minCoverCorrected, 5); + setGamsParamValueTruncate(minimumLandCoverP.addRecord(v), minCoverCorrected, 5); } } - + // Land cover conversion cost + GAMSParameter conversionCostP = inDB.addParameter("conversionCost", 2); + DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = new ConversionCostReader().read(); + + for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap : conversionCosts.getMap().entrySet()) { + String fromName = fromMap.getKey().getName(); + for (Entry<LandCoverType, Double> toMap : fromMap.getValue().entrySet()) { + String toName = toMap.getKey().getName(); + double cost = toMap.getValue(); + Vector<String> v = new Vector<String>(); + v.add(fromName); + v.add(toName); + setGamsParamValue(conversionCostP.addRecord(v), cost, 5); + } + } } private void addScalar(GAMSDatabase gamsDb, String recordName, double val, int places) { diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java index e4e1171b..759107b3 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java @@ -19,12 +19,11 @@ public class GamsRasterInput { private RasterSet<CarbonFluxItem> carbonFluxes; private RasterSet<WoodYieldItem> woodYields; private DoubleMap<Integer, LandCoverType, Double> minimumLandCover; - private DoubleMap<Integer, LandCoverType, Double> loggedForest; private GamsCountryInput countryInput; public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses, RasterSet<IrrigationItem> irrigationCost, RasterSet<CarbonFluxItem> carbonFluxes, RasterSet<WoodYieldItem> woodYields, DoubleMap<Integer, LandCoverType, Double> minimumLandCover, - DoubleMap<Integer, LandCoverType, Double> loggedForest, GamsCountryInput countryInput) { + GamsCountryInput countryInput) { super(); this.timestep = timestep; this.yields = yields; @@ -33,7 +32,6 @@ public class GamsRasterInput { this.carbonFluxes = carbonFluxes; this.woodYields = woodYields; this.minimumLandCover = minimumLandCover; - this.loggedForest = loggedForest; this.countryInput = countryInput; } @@ -61,10 +59,6 @@ public class GamsRasterInput { return minimumLandCover; } - public DoubleMap<Integer, LandCoverType, Double> getLoggedForest() { - return loggedForest; - } - public GamsCountryInput getCountryInput() { return countryInput; } diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java index b7dcecda..ca040d51 100644 --- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java +++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java @@ -384,27 +384,6 @@ public class GamsRasterOptimiser { checkedTotalAreas(landUseRaster, "before"); checkedTotalAreas(aggregatedAreas, "after"); - - // Auto convert timberForest to otherNatural once matured and cut down - if (!rasterInputData.getLoggedForest().getMap().isEmpty()) { - for (Entry<Integer, Map<LandCoverType, Double>> locMap : rasterInputData.getLoggedForest().getMap().entrySet()) { - int locId = locMap.getKey(); - - Set<RasterKey> keys = new HashSet<RasterKey>(); - for (Entry<RasterKey, IntegerRasterItem> mapEntry : mapping.entrySet()) { - IntegerRasterItem iri = mapEntry.getValue(); - if (iri != null && locId == iri.getInt()) - keys.add(mapEntry.getKey()); - } - - for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) { - // move land for GAMS input - aggregatedAreas.get(locId).moveUnprotectedArea(LandCoverType.OTHER_NATURAL, coverMap.getKey(), coverMap.getValue()); - // move land for raster - allocAllLandCropsTop(rasterInputData.getPreviousLandUses(), keys, LandCoverType.OTHER_NATURAL, coverMap.getKey(), coverMap.getValue(), locId); - } - } - } return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getMinimumLandCover(), rasterInputData.getCountryInput()); diff --git a/src/ac/ed/lurg/landuse/ConversionCostReader.java b/src/ac/ed/lurg/landuse/ConversionCostReader.java new file mode 100644 index 00000000..05aec14e --- /dev/null +++ b/src/ac/ed/lurg/landuse/ConversionCostReader.java @@ -0,0 +1,53 @@ +package ac.ed.lurg.landuse; + +import java.io.BufferedReader; +import java.io.FileReader; +import java.io.IOException; + +import ac.ed.lurg.ModelConfig; +import ac.ed.lurg.types.LandCoverType; +import ac.ed.lurg.utils.DoubleMap; +import ac.ed.lurg.utils.LogWriter; + +public class ConversionCostReader { + + private static final int FROM_COL = 0; + private static final int TO_COL = 1; + private static final int COST_COL = 2; + + public DoubleMap<LandCoverType, LandCoverType, Double> read() { + DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = new DoubleMap<LandCoverType, LandCoverType, Double>(); + String filename = ModelConfig.CONVERSION_COST_FILE; + try { + BufferedReader reader = new BufferedReader(new FileReader(filename)); + String line, fromName, toName; + Double cost; + reader.readLine(); // read header + + while ((line=reader.readLine()) != null) { + String[] tokens = line.split(","); + + if (tokens.length < 3) + LogWriter.printlnError("Too few columns in " + filename + ", " + line); + + fromName = tokens[FROM_COL]; + toName = tokens[TO_COL]; + cost = Double.parseDouble(tokens[COST_COL]); + + LandCoverType fromLC = LandCoverType.getForName(fromName); + LandCoverType toLC = LandCoverType.getForName(toName); + + conversionCosts.put(fromLC, toLC, cost); + } + reader.close(); + + } catch (IOException e) { + LogWriter.printlnError("Failed in reading conversion costs file"); + LogWriter.print(e); + } + LogWriter.println("Processed " + filename); + + return conversionCosts; + } + +} diff --git a/src/ac/ed/lurg/landuse/WoodYieldReader.java b/src/ac/ed/lurg/landuse/WoodYieldReader.java index 76f92aaf..15f77c00 100644 --- a/src/ac/ed/lurg/landuse/WoodYieldReader.java +++ b/src/ac/ed/lurg/landuse/WoodYieldReader.java @@ -31,12 +31,12 @@ public class WoodYieldReader extends AbstractTabularRasterReader<WoodYieldItem> item.setWoodYield(LandCoverType.PASTURE, LandCoverType.CROPLAND, getValueForCol(rowValues, "past_to_crop") * conversionFactor); item.setWoodYield(LandCoverType.PASTURE, LandCoverType.PASTURE, getValueForCol(rowValues, "past_to_past") * conversionFactor); - item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.PASTURE, getValueForCol(rowValues, "ntrl_to_past") * conversionFactor); - item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.TIMBER_FOREST, getValueForCol(rowValues, "ntrl_to_forT") * conversionFactor); - item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.CARBON_FOREST, getValueForCol(rowValues, "ntrl_to_forC") * conversionFactor); - item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.UNMANAGED_FOREST, getValueForCol(rowValues, "ntrl_to_ntrl") * conversionFactor); - item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.CROPLAND, getValueForCol(rowValues, "ntrl_to_crop") * conversionFactor); - item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "ntrl_to_ntrl") * conversionFactor); + item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.PASTURE, 0.0); + item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.TIMBER_FOREST, 0.0); + item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.CARBON_FOREST, 0.0); + item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.UNMANAGED_FOREST, 0.0); + item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.CROPLAND, 0.0); + item.setWoodYield(LandCoverType.OTHER_NATURAL, LandCoverType.OTHER_NATURAL, 0.0); item.setWoodYield(LandCoverType.TIMBER_FOREST, LandCoverType.PASTURE, getValueForCol(rowValues, "forT_to_past") * conversionFactor); item.setWoodYield(LandCoverType.TIMBER_FOREST, LandCoverType.OTHER_NATURAL, getValueForCol(rowValues, "forT_to_ntrl") * conversionFactor); -- GitLab