Skip to content
Snippets Groups Projects

Bringing Bart's recent changes into the inequality branch

Closed Bringing Bart's recent changes into the inequality branch
master into subnational_inequality
Closed s1713176 requested to merge master into subnational_inequality
Viewing commit 12de71fd
Prev Next
Show latest version
10 files
+ 54
− 371
Compare changes
  • Side-by-side
  • Inline
Files
10
GAMS/IntExtOpt.gms
+ 14
− 22
@@ -6,7 +6,7 @@
SET crop_less_pasture_setaside(crop) / wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops/;
SET feed_crop(crop) / wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, pasture/;
SET feed_crop_less_pasture(feed_crop) / wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg /;
SET import_types(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, wood, carbonCredits/;
SET import_types(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, wood, carbonCredits, energyPV, energyAV/;
SET animal(all_types) / monogastrics, ruminants /;
SET agri_types(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture /;
SET not_feed(all_types) / monogastrics, ruminants, sugar, energycrops /;
@@ -63,7 +63,7 @@
PARAMETER maxLandCoverArea(land_cover, protection, location);
PARAMETER solarEnergyDensity(solar_land, location) MWh per ha;
PARAMETER solarNpvRate(solar_land, location) 1000$ per ha;
PARAMETER solarCostRate(solar_land, location) 1000$ per ha;
PARAMETER agrivoltaicsYieldFactor(crop, farming_type, location);
PARAMETER agrivoltaicsCropCostFactor(farming_type) crop unit cost is lower since crop area is smaller due to space taken by PV panels;
@@ -93,7 +93,7 @@
SCALAR discountRate;
SCALAR energyPrice;
*$gdxin "/Users/bart/Documents/PLUM_Output/calib/GamsTmp/Brazil2020landuse.gdx"
*$gdxin "/Users/bart/Documents/PLUM_Output/test/GamsTmp/China2020landuse.gdx"
$gdxin %gdxincname%
$load location, suitableLandArea, demand
$load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
@@ -103,7 +103,7 @@ $load meatEfficency, otherICost, irrigCost, irrigConstraint, irrigMaxRate, ferti
$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMaxParam, woodYieldSlopeParam, woodYieldShapeParam, tradeAdjustmentCostRate
$load forestManagementCost, carbonForestMaxProportion, previousRotationIntensity, maxLandCoverChange, maxLandCoverArea
$load fertLimRestricted, irrigLimRestricted, otherLimRestricted, fertLimConventional, irrigLimConventional, otherLimConventional
$load solarEnergyDensity, agrivoltaicsYieldFactor, prodCost, animalProdCost, solarNpvRate
$load solarEnergyDensity, agrivoltaicsYieldFactor, prodCost, animalProdCost, solarCostRate
$load discountRate, energyPrice, monoculturePenalty, agrivoltaicsCropCostFactor, previousProduction, pastureIntensityPenalty
$gdxin
@@ -156,7 +156,6 @@ $gdxin
netImportIncrease(import_types)
netImportDecrease(import_types)
excessImports(import_types)
importFraction(all_types)
landCoverArea(land_cover, protection, location) land cover area in Mha
landCoverChange(land_cover_before, land_cover_after, protection, location) land cover change in Mha
@@ -165,7 +164,6 @@ $gdxin
woodYieldRota(location) wood yield m3 per ha
rotationIntensity(location) equivalent to 1 divided by rotation period
woodSupply total wood supply - million m3
forestryNPV(location) NPV of wood production
carbonCredits(location) carbon credits generated - Mt C
@@ -174,7 +172,7 @@ $gdxin
;
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, yield, feedAmount, importAmount, exportAmount, cropProduction, excessImports, importFraction,
POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, yield, feedAmount, importAmount, exportAmount, cropProduction, excessImports,
landCoverArea, landCoverChange, woodYieldRota, rotationIntensity, woodSupply, totalProduction, netImportIncrease, netImportDecrease;
@@ -195,7 +193,6 @@ $gdxin
MIN_NET_IMPORT_CONSTRAINT(import_types)
MAX_NET_IMPORT_CONSTRAINT(import_types)
NET_IMPORTS_CHANGE_CALC(import_types)
IMPORT_FRACTION_CALC(all_types) fraction of supply that is imported
SETASIDE_AREA_CALC(farming_type, location) setaside area constraint
CROPLAND_LAND_COVER_CALC(protection, farming_type, location) cropland area equals sum of all crop areas
@@ -209,7 +206,6 @@ $gdxin
WOOD_YIELD_CALC(location) wood yield at rotation - tC per ha
WOOD_SUPPLY_CALC wood supply calculation
WOOD_DEMAND_CONSTRAINT wood demand constraint
FORESTRY_NPV_CALC(location) net present value of wood production (Faustmann's formula)
PV_ENERGY_DEMAND_CONSTRAINT energy produced from photovoltaics
AV_ENERGY_DEMAND_CONSTRAINT energy produced from agrivoltaics
@@ -322,8 +318,6 @@ $gdxin
WOOD_DEMAND_CONSTRAINT .. woodSupply =E= demand('wood') + exportAmount('wood') - importAmount('wood');
FORESTRY_NPV_CALC(location) .. forestryNPV(location) =E= (-forestManagementCost) * rotationIntensity(location);
*********** Carbon ***********************************
CARBON_CREDIT_CALC(location) .. carbonCredits(location) =E= sum((land_cover_before, land_cover_after, protection), landCoverChange(land_cover_before, land_cover_after, protection, location) * carbonCreditRate(land_cover_before, land_cover_after, location));
@@ -348,8 +342,6 @@ $gdxin
NET_IMPORTS_CHANGE_CALC(import_types) .. netImportIncrease(import_types) - netImportDecrease(import_types) =E= netImportAmount(import_types) - previousNetImport(import_types);
IMPORT_FRACTION_CALC(all_types) .. importFraction(all_types) =E= importAmount(all_types) / (importAmount(all_types) + totalProduction(all_types) + delta);
* No imports or exports of pasture allowed
importAmount.FX('pasture') = 0;
exportAmount.FX('pasture') = 0;
@@ -371,10 +363,10 @@ $gdxin
- sum(animal, animalProdCost(animal) * totalProduction(animal) * (1 - subsidyRate(animal)))
* Cost of feed imports
- sum(import_types, importPrices(import_types) * importAmount(import_types))
* NPV of wood production
+ sum(location, forestryNPV(location) * landCoverArea('timberForest', 'unprotected', location))
* NPV of solar energy
+ sum((solar_land, location), solarNpvRate(solar_land, location) * landCoverArea(solar_land, 'unprotected', location))
* Cost of wood production
- sum(location, forestManagementCost * rotationIntensity(location) * landCoverArea('timberForest', 'unprotected', location))
* Cost of solar energy
- sum((solar_land, location), solarCostRate(solar_land, location) * landCoverArea(solar_land, 'unprotected', location))
* Land cover conversion costs
- sum(location, totalConversionCost(location))
* Value of carbon credits
@@ -387,7 +379,7 @@ $gdxin
- sum(location, pastureIntensityPenalty * sqr(cropProduction('pasture', location))) / (sum(location, cropProduction('pasture', location)) + delta)
* Small reward for increasing natural area to allow land abandonment
+ sum((land_cover, protection, location)$(not sameAs(land_cover, 'natural')), 0.00001 * landCoverChange(land_cover, 'natural', protection, location))
* Feed diversity penalty
* Feed diversity penalty
- sum(animal, feedDiversityCost(animal) * sum(feed_crop, sqr(feedAmount(animal, feed_crop) * cropDM(feed_crop))) / (sum(feed_crop, feedAmount(animal, feed_crop) * cropDM(feed_crop)) + delta))
* Net import change cost
- sum(import_types, tradeAdjustmentCostRate(import_types) * (netImportIncrease(import_types) + netImportDecrease(import_types)))
@@ -432,7 +424,7 @@ $gdxin
parameter totalProdCost(all_types);
parameter totalCropArea(crop);
parameter netImportCost(all_types);
* parameter importFraction(all_types);
parameter importFraction(all_types);
parameter prodCostRate(crop);
parameter productionShock(all_types);
parameter woodProdCost;
@@ -445,12 +437,12 @@ $gdxin
netImportCost(import_types) = importAmount.L(import_types) * importPrices(import_types) - exportAmount.l(import_types) * exportPrices(import_types);
* importFraction(import_types) $ (importAmount.l(import_types) > 0) = importAmount.l(import_types) / (totalProduction.l(import_types) + importAmount.l(import_types));
importFraction(import_types) $ (importAmount.l(import_types) > 0) = importAmount.l(import_types) / (totalProduction.l(import_types) + importAmount.l(import_types));
prodCostRate(crop)$[totalProduction.l(crop) <> 0] = totalProdCost(crop) / totalProduction.l(crop);
totalProdCost(animal) = sum(feed_crop, importFraction.l(feed_crop) * feedAmount.l(animal, feed_crop) * importPrices(feed_crop)
+ (1 - importFraction.l(feed_crop)) * prodCostRate(feed_crop)) + totalProduction.l(animal) * animalProdCost(animal);
totalProdCost(animal) = sum(feed_crop, importFraction(feed_crop) * feedAmount.l(animal, feed_crop) * importPrices(feed_crop)
+ (1 - importFraction(feed_crop)) * prodCostRate(feed_crop)) + totalProduction.l(animal) * animalProdCost(animal);
If(woodSupply.L > 0,
woodProdCost = sum(location, forestManagementCost * rotationIntensity.L(location) * landCoverArea.L('timberForest', 'unprotected', location));