From d85106f28a00584d54dec46bcc66df875e31ca98 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Mon, 8 Mar 2021 16:23:04 +0000
Subject: [PATCH] Change to how landCoverChange is calculated in GAMS.

---
 GAMS/IntExtOpt.gms | 31 +++++++++++--------------------
 1 file changed, 11 insertions(+), 20 deletions(-)

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index c915b7de..6e72b652 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -77,23 +77,10 @@ $load previousLandCoverArea, carbonFluxRate, woodYield, carbonPrice, woodPrice
 $load minimumLandCoverArea
 $gdxin
 
-* convert units from $/t to billion$/Mt
+*convert units from $/t to billion$/Mt
  carbonPrice = carbonPrice * 0.001;
  woodPrice = woodPrice * 0.001;
 
-* convert units from 1000$/ha to million$/Mha
-* agriExpansionCost(location) = agriExpansionCost(location) * 1000;
-* pastureIncCost = pastureIncCost * 1000;
-* cropIncCost = cropIncCost * 1000;
-* cropDecCost = cropDecCost * 1000;
-* pastureDecCost = pastureDecCost * 1000;
-
-* covert units from 1000$/t to million$/Mt
-* exportPrices(import_crop) = exportPrices(import_crop) * 1000;
-* importPrices(import_crop) = importPrices(import_crop) * 1000;
-
-
-
  SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /;
 
  demand(cereal_crop) = demand('cereals') * minDemandPerCereal(cereal_crop);
@@ -162,7 +149,7 @@ $gdxin
                    landCoverArea, landCoverChange, woodHarvest;
 
  EQUATIONS
-       UNIT_COST_EQ(crop, location)                     cost per area - million$ per Mha
+       UNIT_COST_EQ(crop, location)                     cost per area - $1000 per ha or $billion per Mha
        YIELD_EQ(crop, location)                         yield given chosen intensity - tonnes per hectare
        CROP_DEMAND_CONSTRAINT(crop)                     satisfy demand for individual crops
        TOTAL_CEREAL_DEMAND_CONSTRAINT                   satisfy demand for combined cereals
@@ -192,8 +179,8 @@ $gdxin
        PASTURE_LAND_COVER_CALC(location)  pasture area (land cover) equals pasture area (land use)
        MINIMUM_LAND_COVER_CONSTRAINT(location) constraint on land cover conversion
        LAND_COVER_CHANGE_CALC(land_cover, location)
-       LAND_COVER_CHANGE_CONSTRAINT(land_cover, location)
-       LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location)
+       LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) can't convert more land than was previously available
+       LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) calculate no change area (e.g. cropland to cropland)
        WOOD_HARVEST_CALC(location)                          calc wood harvested
        CARBON_FLUX_CALC(location)                           calc carbon flux
        COST_EQ                                        total cost objective function;
@@ -254,14 +241,19 @@ $gdxin
 
  PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
 
- MINIMUM_LAND_COVER_CONSTRAINT(location, land_cover) .. landCoverArea(land_cover, location)) =G= minimumLandCover(land_cover_after, location);
+ MINIMUM_LAND_COVER_CONSTRAINT(location, land_cover) .. landCoverArea(land_cover, location) =G= minimumLandCover(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));
 
+* LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) -
+*                                                                                         sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
+
+
+
  LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =L= previousLandCoverArea(land_cover, location);
 
- LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =E= 0;
+ LAND_COVER_SELF_CHANGE_CONSTRAINT(land_cover, location) .. landCoverChange(land_cover, land_cover, location) =E= previousLandCoverArea(land_cover, location) - sum(land_cover_after$[not sameAs(land_cover, land_cover_after)], landCoverChange(land_cover, land_cover_after, location));
 
  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));
 
@@ -341,4 +333,3 @@ $gdxin
  Scalar ms 'model status', ss 'solve status';
  ms=LAND_USE.modelstat;
  ss=LAND_USE.solvestat;
-
-- 
GitLab