From eb96c77d6b0a10295a7a03c58c183d61942cd667 Mon Sep 17 00:00:00 2001
From: Bart Arendarczyk <s1924442@ed.ac.uk>
Date: Wed, 13 Mar 2024 11:31:41 +0000
Subject: [PATCH] Photovoltaics and agrivoltaics.
---
GAMS/IntExtOpt.gms | 68 +++++++++----
data/agrivoltaic_yield_params.csv | 12 +++
src/ac/ed/lurg/ModelConfig.java | 34 +++++++
src/ac/ed/lurg/ModelMain.java | 27 +++++-
src/ac/ed/lurg/country/CountryAgent.java | 11 ++-
.../ed/lurg/country/CountryAgentManager.java | 7 +-
.../lurg/country/gams/GamsLocationInput.java | 10 +-
.../country/gams/GamsLocationOptimiser.java | 71 ++++++++++++--
.../ed/lurg/country/gams/GamsRasterInput.java | 11 ++-
.../country/gams/GamsRasterOptimiser.java | 29 +++++-
.../ed/lurg/landuse/ConversionCostReader.java | 48 +++++----
src/ac/ed/lurg/landuse/LandUseItem.java | 9 ++
src/ac/ed/lurg/output/LandUseOutputer.java | 16 ++-
.../lurg/solar/AgrivoltaicsYieldManager.java | 42 ++++++++
src/ac/ed/lurg/solar/SolarPotentialItem.java | 97 +++++++++++++++++++
.../ed/lurg/solar/SolarPotentialReader.java | 37 +++++++
src/ac/ed/lurg/types/FarmingType.java | 2 +-
src/ac/ed/lurg/types/LandCoverType.java | 2 +
18 files changed, 468 insertions(+), 65 deletions(-)
create mode 100644 data/agrivoltaic_yield_params.csv
create mode 100644 src/ac/ed/lurg/solar/AgrivoltaicsYieldManager.java
create mode 100644 src/ac/ed/lurg/solar/SolarPotentialItem.java
create mode 100644 src/ac/ed/lurg/solar/SolarPotentialReader.java
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 35cde672..f1e8bab1 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -1,5 +1,5 @@
- SET all_types / monogastrics, ruminants, cerealsStarchyRoots, wheat, maize, rice, oilcrops, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside, roundwood, fuelwood, carbonCredits/;
+ SET all_types / monogastrics, ruminants, cerealsStarchyRoots, wheat, maize, rice, oilcrops, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside, roundwood, fuelwood, carbonCredits, energy /;
SET crop(all_types) / wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside /;
SET crop_less_pasture(crop) / wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops , setaside /;
@@ -9,19 +9,20 @@
SET import_types(all_types) / monogastrics, ruminants, wheat, maize, rice, oilcropsNFix, oilcropsOther, pulses, starchyRoots, fruitveg, sugar, energycrops, roundwood, fuelwood, carbonCredits/;
SET animal(all_types) / monogastrics, ruminants /;
- SET commodity(all_types) / ruminants, monogastrics, cerealsStarchyRoots, oilcrops, pulses, fruitveg, sugar, energycrops, roundwood, fuelwood, carbonCredits /;
+ SET commodity(all_types) / ruminants, monogastrics, cerealsStarchyRoots, oilcrops, pulses, fruitveg, sugar, roundwood, fuelwood, carbonCredits, energycrops, energy /;
SET animal_commodity(commodity) / ruminants, monogastrics /;
SET crop_commodity(commodity) / cerealsStarchyRoots, oilcrops, pulses, fruitveg, sugar, energycrops /;
SET demand_map(commodity, all_types);
- SET land_cover / cropland, pasture, timberForest, carbonForest, natural /;
+ SET land_cover / cropland, pasture, timberForest, carbonForest, natural, agrivoltaics, photovoltaics /;
SET protection / unprotected, protected /;
SET managed_forest(land_cover) / timberForest, carbonForest /;
+ SET solar_land(land_cover) / agrivoltaics, photovoltaics /;
ALIAS (land_cover, land_cover_before);
ALIAS (land_cover, land_cover_after);
- SET farming_type / conventional, restricted /;
+ SET farming_type / conventional, restricted, agrivoltaics /;
SET farming_map(protection, farming_type) / unprotected.conventional, protected.restricted /;
SET wood_type(commodity) / roundwood, fuelwood /;
@@ -66,6 +67,10 @@
PARAMETER forestBaseCost(managed_forest);
PARAMETER previousRotationIntensity(location);
PARAMETER maxLandCoverChange(land_cover_before, land_cover_after, protection, location);
+
+ PARAMETER solarEnergyDensity(solar_land, location) MWh per ha;
+ PARAMETER solarCost(solar_land, location) 1000$ per ha;
+ PARAMETER agrivoltaicsYieldFactor(crop, farming_type, location);
SCALAR meatEfficency efficiency of converting feed and pasture into animal products;
SCALAR fertiliserUnitCost fert cost at max fert rate;
@@ -85,10 +90,12 @@
SCALAR otherLimConventional;
SCALAR forestManagementCost cost $1000 per ha;
-
SCALAR carbonForestMaxProportion maximum proportion of land cover as carbon forest;
+
+ SCALAR energycropsEnergyDensity;
+ SCALAR bioenergyGenerationCost;
-*$gdxin "./_gams_java_gdb1.gdx"
+*$gdxin "/Users/bart/Documents/PLUM_Output/calib_forest/GamsTmp/loc4965753693958365911/_gams_java_gdb1.gdx"
$gdxin %gdxincname%
$load location, suitableLandArea, demand, demand_map
$load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
@@ -98,6 +105,7 @@ $load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, ferti
$load previousLandCoverArea, maxCroplandArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam, forestBaseCost, tradeAdjustmentCostRate
$load forestManagementCost, carbonForestMaxProportion, maxFertChange, maxIrrigChange, previousRotationIntensity, maxLandCoverChange
$load fertLimRestricted, irrigLimRestricted, otherLimRestricted, fertLimConventional, irrigLimConventional, otherLimConventional
+$load solarEnergyDensity, agrivoltaicsYieldFactor, solarCost, energycropsEnergyDensity, bioenergyGenerationCost
$gdxin
SCALAR delta "use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf" / 0.00000000001 /;
@@ -160,6 +168,8 @@ $gdxin
fuelwood 0.03
/;
+demand('energy') = demand('energycrops') * energycropsEnergyDensity;
+demand('energycrops') = 0;
VARIABLES
unitCost(crop, farming_type, location) cost per area for each crop - cost
@@ -188,11 +198,16 @@ $gdxin
forestryCost(managed_forest, location) total cost in 1000$
carbonCredits(location)
+ bioenergyFeedstock(crop)
+
+ energyProdCost
+
+
* 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;
- POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, netSupply, ruminantFeed, monogastricFeed, importAmount, exportAmount, totalFeedDM,
- landCoverArea, landCoverChange, totalConversionCost, animalProd, woodYieldRota, rotationIntensity, woodSupply;
+ POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, netSupply, ruminantFeed, monogastricFeed, importAmount, exportAmount, totalFeedDM, bioenergyFeedstock,
+ landCoverArea, landCoverChange, totalConversionCost, animalProd, woodYieldRota, rotationIntensity, woodSupply, energyProd, energyProdCost;
* POSITIVE VARIABLE A;
@@ -201,8 +216,6 @@ $gdxin
UNIT_COST_EQ(crop, farming_type, location) cost per area - $1000 per ha or $billion per Mha
YIELD_EQ(crop, farming_type, location) yield given chosen intensity - tonnes per hectare
IRRIGATION_CONSTRAINT(location) constraint on water usage
-* FERT_RATE_INCREASE_CONSTRAINT
-* IRRIG_RATE_INCREASE_CONSTRAINT
NET_SUPPLY_EQ(crop) calc net supply for crops
CROP_DEMAND_CONSTRAINT(crop_commodity, crop) satisfy demand for crop products
@@ -220,6 +233,7 @@ $gdxin
SETASIDE_AREA_CALC(farming_type, location)
CROPLAND_LAND_COVER_CALC(protection, farming_type, location) cropland area equals sum of all crop areas
PASTURE_LAND_COVER_CALC(protection, farming_type, location) pasture area (land cover) equals pasture area (land use)
+ AGRIVOLTAICS_LAND_COVER_CALC(location)
LAND_COVER_CHANGE_CALC(land_cover, protection, location) "calc land cover change. landCoverChange(A, A, location) defined as unchanged land cover"
LAND_COVER_CHANGE_CONSTRAINT(land_cover, protection, location) conservation of land area
CONVERSION_COST_CALC(location) cost of land cover conversion
@@ -230,6 +244,9 @@ $gdxin
WOOD_DEMAND_CONSTRAINT(wood_type)
ROUNDWOOD_DEMAND_CONSTRAINT satisfy roundwood demand (can only by obtained from older trees)
TIMBER_FOREST_COST_CALC(location)
+
+ ENERGY_COST_CALC
+ ENERGY_DEMAND_CONSTRAINT
CARBON_CREDIT_CALC(location)
CARBON_CREDIT_CONSTRAINT
@@ -254,17 +271,15 @@ $gdxin
(yieldIrrigOnly(crop, location) - yieldNone(crop, location)) * (1 - exp(-irrigParam(crop, location) * irrigI(crop, farming_type, location))) +
(yieldBoth(crop, location) + yieldNone(crop, location) - yieldFertOnly(crop, location) - yieldIrrigOnly(crop, location)) *
(1 - exp(-fertParam(crop, location) * fertI(crop, farming_type, location))) * (1 - exp(-irrigParam(crop, location) * irrigI(crop, farming_type, location)))
- ) * (1.01 - exp(-otherIntensity(crop, farming_type, location) * otherIParam));
+ ) * (1.01 - exp(-otherIntensity(crop, farming_type, location) * otherIParam)) * agrivoltaicsYieldFactor(crop, farming_type, location);
IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) * suitableLandArea(location) =G= sum((crop, farming_type), irrigMaxRate(crop, location) * irrigI(crop, farming_type, location) * cropArea(crop, farming_type, location));
-
-* FERT_RATE_INCREASE_CONSTRAINT .. sum((crop, location), fertI(crop, location) * cropArea(crop, location)) =L=
-* sum((crop, location), previousFertIntensity(crop, location) * previousCropArea(crop, location)) * (1 + maxFertChange);
-
-* IRRIG_RATE_INCREASE_CONSTRAINT .. sum((crop, location), irrigI(crop, location) * irrigMaxRate(crop, location) * cropArea(crop, location)) =L=
-* sum((crop, location), previousIrrigIntensity(crop, location) * irrigMaxRate(crop, location) * previousCropArea(crop, location)) * (1 + maxIrrigChange);
* Intensities must be <= 1
+ fertI.UP(crop, farming_type, location) = 1;
+ irrigI.UP(crop, farming_type, location) = 1;
+ otherIntensity.UP(crop, farming_type, location) = 1;
+
fertI.UP(crop, 'conventional', location) = fertLimConventional;
irrigI.UP(crop, 'conventional', location) = irrigLimConventional;
otherIntensity.UP(crop, 'conventional', location) = otherLimConventional;
@@ -273,13 +288,17 @@ $gdxin
irrigI.UP(crop, 'restricted', location) = irrigLimRestricted;
otherIntensity.UP(crop, 'restricted', location) = otherLimRestricted;
+ otherIntensity.LO(crop, 'agrivoltaics', location) = 0.2;
+
* Fixing irrigI at 0 where no water to avoid optimisation problems
irrigI.FX(crop, farming_type, location) $ (irrigConstraint(location) = 0) = 0;
*************** Crop supply *************************
NET_SUPPLY_EQ(crop) .. netSupply(crop) =E= sum((location, farming_type), cropArea(crop, farming_type, location) * yield(crop, farming_type, location)) *
- (1 - seedAndWasteRate(crop)) - ruminantFeed(crop) - monogastricFeed(crop) + importAmount(crop) - exportAmount(crop);
+ (1 - seedAndWasteRate(crop)) - ruminantFeed(crop) - monogastricFeed(crop) + importAmount(crop) - exportAmount(crop) - bioenergyFeedstock(crop);
+
+ bioenergyFeedstock.FX(crop)$(not sameAs(crop, 'energycrops')) = 0;
CROP_DEMAND_CONSTRAINT(crop_commodity, crop)$demand_map(crop_commodity, crop) .. netSupply(crop) =G= demand(crop_commodity) * minDemandFraction(crop_commodity, crop);
@@ -301,6 +320,8 @@ $gdxin
PASTURE_LAND_COVER_CALC(protection, farming_type, location)$farming_map(protection, farming_type) .. landCoverArea('pasture', protection, location) =E= cropArea('pasture', farming_type, location);
+ AGRIVOLTAICS_LAND_COVER_CALC(location) .. landCoverArea('agrivoltaics', 'unprotected', location) =E= sum(crop, cropArea(crop, 'agrivoltaics', location));
+
LAND_COVER_CHANGE_CALC(land_cover, protection, location) .. landCoverArea(land_cover, protection, location) =E= previousLandCoverArea(land_cover, protection, location) +
sum(land_cover_before, landCoverChange(land_cover_before, land_cover, protection, location)) -
sum(land_cover_after, landCoverChange(land_cover, land_cover_after, protection, location));
@@ -312,7 +333,7 @@ $gdxin
landCoverChange(land_cover_before, land_cover_after, protection, location) * conversionCost(land_cover_before, land_cover_after, location));
* Cropland area restriction due to slope
- MAX_CROPLAND_CONSTRAINT(location) .. sum(protection, landCoverArea('cropland', protection, location)) =L= maxCroplandArea(location);
+ MAX_CROPLAND_CONSTRAINT(location) .. sum((crop_less_pasture, farming_type), cropArea(crop_less_pasture, farming_type, location)) =L= maxCroplandArea(location);
landCoverChange.UP(land_cover_before, land_cover_after, protection, location) = maxLandCoverChange(land_cover_before, land_cover_after, protection, location);
@@ -345,6 +366,13 @@ $gdxin
CARBON_FOREST_COST_CALC(location) .. forestryCost('carbonForest', location) =E= sum(land_cover$[not sameAs(land_cover, 'carbonForest')], landCoverChange(land_cover, 'carbonForest', 'unprotected', location) * forestManagementCost) +
landCoverArea('carbonForest', 'unprotected', location) * forestBaseCost('carbonForest');
+
+*********** Energy *********************
+
+ ENERGY_COST_CALC .. energyProdCost =E= sum((solar_land, location), solarCost(solar_land, location) * landCoverArea(solar_land, 'unprotected', location)) +
+ bioenergyGenerationCost * energycropsEnergyDensity * bioenergyFeedstock('energycrops');
+
+ ENERGY_DEMAND_CONSTRAINT .. sum((solar_land, location), solarEnergyDensity(solar_land, location) * landCoverArea(solar_land, 'unprotected', location)) + energycropsEnergyDensity * bioenergyFeedstock('energycrops') =G= demand('energy');
*************** Trade *******************************
@@ -366,6 +394,7 @@ $gdxin
sum(location, totalConversionCost(location)) +
sum((managed_forest, location), forestryCost(managed_forest, location)) +
sum(wood_type, woodSupply(wood_type) * loggingCost(wood_type)) +
+ energyProdCost +
sum(import_types, importPrices(import_types) * importAmount(import_types)) -
sum(import_types, exportPrices(import_types) * exportAmount(import_types)) +
sum(import_types, tradeAdjustmentCost(import_types));
@@ -389,7 +418,6 @@ $gdxin
woodYieldRota.L(location) = woodYieldMax(location) * (1 - exp(-woodYieldParam(location) / rotationIntensity.L(location))) * rotationIntensity.L(location);
woodSupply.L(wood_type) = demand(wood_type);
-
* LAND_USE.OptFile = 1;
* display landCoverChange.L
diff --git a/data/agrivoltaic_yield_params.csv b/data/agrivoltaic_yield_params.csv
new file mode 100644
index 00000000..ad890fe0
--- /dev/null
+++ b/data/agrivoltaic_yield_params.csv
@@ -0,0 +1,12 @@
+Crop,ParamA,ParamB
+fruitveg,0.6773239338183652,-1.7756495663799552
+pasture,0.6069465888798169,-1.7756495663799552
+maize,-1.407993616487431,-1.7756495663799552
+oilcropsNFix,-0.9831919372331905,-1.7756495663799552
+wheat,-0.23457877882018402,-1.7756495663799552
+starchyRoots,-0.44668611705014233,-1.7756495663799552
+sugar,-0.5936845126152303,-1.7756495663799552
+pulses,-1.663306897928871,-1.7756495663799552
+rice,-0.23457877882018402,-1.7756495663799552
+energycrops,-1.407993616487431,-1.7756495663799552
+oilcropsOther,-0.9831919372331905,-1.7756495663799552
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index b7f3e01f..d140a83c 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -439,6 +439,7 @@ public class ModelConfig {
public static final double AGRI_LAND_EXPANSION_COST_FACTOR = getDoubleProperty("AGRI_LAND_EXPANSION_COST_FACTOR", 1.0); // for monte carlo
public static final double CROPLAND_CONVERSION_COST = getDoubleProperty("CROPLAND_CONVERSION_COST", 0.3);
public static final double PASTURE_CONVERSION_COST = getDoubleProperty("PASTURE_CONVERSION_COST", 0.2);
+ public static final double SOLAR_CONVERSION_COST = getDoubleProperty("SOLAR_CONVERSION_COST", 0.8);
public static final double FOREST_CONVERSION_COST = getDoubleProperty("FOREST_CONVERSION_COST", 0.5);
public static final double NATURAL_CONVERSION_COST = getDoubleProperty("NATURAL_CONVERSION_COST", 0.02);
@@ -584,4 +585,37 @@ public class ModelConfig {
// How demand is allocated. Options: global (countries can export carbon credits to market), country (global demand allocated by GDP), price (no demand, exogenous price)
public static final String CARBON_DEMAND_METHOD = getProperty("CARBON_DEMAND_METHOD", "global");
+ // Solar
+ public static final boolean IS_PHOTOVOLTAICS_ON = getBooleanProperty("IS_PHOTOVOLTAICS_ON", true);
+ public static final boolean IS_AGRIVOLTAICS_ON = getBooleanProperty("IS_AGRIVOLTAICS_ON", true);
+ public static final String AGRIVOLTAICS_YIELD_PARAM_FILE = getProperty("AGRIVOLTAICS_YIELD_PARAM_FILE", DATA_DIR + File.separator + "agrivoltaic_yield_params.csv");
+ public static final String SOLAR_POTENTIAL_DIR_BASE = getProperty("SOLAR_POTENTIAL_DIR_BASE");
+ public static final String SOLAR_POTENTIAL_DIR_TOP = getProperty("SOLAR_POTENTIAL_DIR_TOP");
+ public static final String SOLAR_POTENTIAL_FILENAME = getProperty("SOLAR_POTENTIAL_FILENAME", "pv_yield.csv");
+ public static final int SOLAR_POTENTIAL_DATA_STEP_SIZE = getIntProperty("SOLAR_POTENTIAL_DATA_STEP_SIZE", 10);
+
+ // PV costs from IRENA, Renewable Power Generation Costs in 2022
+ public static final double PV_INSTALLED_COST = getDoubleProperty("PV_INSTALLED_COST", 0.876); // $1000/kW
+ public static final double PV_MAINTENANCE_COST = getDoubleProperty("PV_MAINTENANCE_COST", 0.0132); // $1000/kW/yr
+ public static final double PV_LIFESPAN = getDoubleProperty("PV_LIFESPAN", 25.0); // years
+ public static final double PV_EFFICIENCY = getDoubleProperty("PV_EFFICIENCY", 0.2);
+
+ // GSR and performance ratio from https://doi.org/10.1016/j.energy.2015.10.108
+ public static final double PV_GSR = getDoubleProperty("PV_GSR", 0.75); // Generator to System Ratio (area of (PV + spacing) / total area)
+ public static final double PV_PERFORMANCE_RATIO = getDoubleProperty("PV_PERFORMANCE_RATIO", 0.65); // Losses e.g. power conversion, PV cell degradation
+
+ // Typical GCR, cost factor and yield factor from https://doi.org/10.1007/s10457-023-00906-3
+ public static final double AV_GCR = getDoubleProperty("AV_GCR", 0.5); // Ground Coverage Ratio of agrivoltaics
+ public static final double AV_COST_FACTOR = getDoubleProperty("AV_COST_FACTOR", 1.3); // Cost multiplier for AV compared to PV
+ public static final double AV_YIELD_FACTOR = getDoubleProperty("AV_YIELD_FACTOR", 0.9); // Yield reduction due to area taken by infrastructure
+
+ // Bioenergy
+
+ // Bioenergy assumptions: installed cost = 2162 $/kW, fixed O&M = 4% of installed cost per year,
+ // variable O&M = 0.005 $/kWh, project lifespan = 20 years.
+ // Source: IRENA, Renewable Power Generation Costs in 2022
+ public static final double BIOENERGY_GENERATION_COST = getDoubleProperty("BIOENERGY_COST", 0.0272); // Cost of electricity generation from biomass $1000/MWh.
+ public static final double ENERGYCROPS_ENERGY_DENSITY = getDoubleProperty("ENERGYCROPS_ENERGY_DENSITY", 5.0); // MWh/t DM
+ public static final double BIOENERGY_GENERATION_EFFICIENCY = getDoubleProperty("BIOENERGY_GENERATION_EFFICIENCY", 0.8); // IRENA, CHP system efficiency
+
}
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 8f3e706e..b439da89 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -13,6 +13,8 @@ import ac.ed.lurg.landuse.*;
import ac.ed.lurg.output.LandCoverChangeOutputer;
import ac.ed.lurg.output.LandUseOutputer;
import ac.ed.lurg.output.LpjgOutputer;
+import ac.ed.lurg.solar.SolarPotentialItem;
+import ac.ed.lurg.solar.SolarPotentialReader;
import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
@@ -46,6 +48,7 @@ public class ModelMain {
private WoodYieldRasterSet woodYieldData;
private CarbonFluxReader carbonFluxReader;
private CarbonFluxRasterSet carbonFluxData;
+ private RasterSet solarPotentialData;
public static void main(String[] args) {
@@ -125,6 +128,8 @@ public class ModelMain {
getWoodYieldData(timestep);
getCarbonFluxData(timestep);
+
+ getSolarPotentialData(timestep);
LogWriter.println("Memory usage 1: " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()) / (1024.0*1024.0*1024.0));
@@ -133,7 +138,7 @@ public class ModelMain {
handleMissingData();
- countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, carbonFluxData, woodYieldData, conversionCosts);
+ countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, carbonFluxData, woodYieldData, solarPotentialData, conversionCosts);
if (ModelConfig.RESET_STOCK_YEAR == timestep.getYear())
internationalMarket.resetStocks();
@@ -882,6 +887,20 @@ public class ModelMain {
}
}
+ /** Get solar potential data */
+ private void getSolarPotentialData(Timestep timestep) {
+ solarPotentialData = new RasterSet<SolarPotentialItem>(desiredProjection) {
+ private static final long serialVersionUID = 2517662660687253540L;
+ protected SolarPotentialItem createRasterData() {
+ return new SolarPotentialItem();
+ }
+ };
+
+ if (ModelConfig.IS_PHOTOVOLTAICS_ON || ModelConfig.IS_AGRIVOLTAICS_ON) {
+ new SolarPotentialReader(solarPotentialData).readData(timestep);
+ }
+ }
+
/** Ugly in situ update of currentIrrigationData, better if IrrigationRasterSets were handled more immutably */
private void getUpdateIrrigationData(Timestep timestep) {
String rootDir = timestep.getYearSubDir(ModelConfig.YIELD_DIR);
@@ -995,7 +1014,8 @@ public class ModelMain {
boolean hasIrrig = currentIrrigationData.containsKey(key);
boolean hasWoodYields = !ModelConfig.IS_FORESTRY_ON || woodYieldData.containsKey(key);
boolean hasCarbonFluxes = !ModelConfig.IS_CARBON_ON || carbonFluxData.containsKey(key);
- boolean isComplete = (hasCropYields && hasIrrig && hasWoodYields && hasCarbonFluxes);
+ boolean hasSolarPotential = !(ModelConfig.IS_AGRIVOLTAICS_ON || ModelConfig.IS_PHOTOVOLTAICS_ON) || solarPotentialData.containsKey(key);
+ boolean isComplete = (hasCropYields && hasIrrig && hasWoodYields && hasCarbonFluxes && hasSolarPotential);
if (isComplete) {
coverage++;
}
@@ -1010,6 +1030,9 @@ public class ModelMain {
if (!hasCarbonFluxes)
carbonFluxData.put(key, CarbonFluxItem.getDefault());
+
+ if(!hasSolarPotential)
+ solarPotentialData.put(key, SolarPotentialItem.getDefault());
}
double coveragePc = (double) coverage / totalCount * 100;
LogWriter.println(String.format("Data coverage: %.2f%%", coveragePc), 1);
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 1d2a809c..b9b1cf05 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -16,6 +16,7 @@ import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.*;
import ac.ed.lurg.shock.MinMaxNetImportManager;
import ac.ed.lurg.shock.MinMaxNetImportManager.LimitType;
+import ac.ed.lurg.solar.SolarPotentialItem;
import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
@@ -103,7 +104,8 @@ public class CountryAgent extends AbstractCountryAgent {
public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData,
Map<CropType, GlobalPrice> worldPrices, RasterSet<CarbonFluxItem> carbonFluxData,
- RasterSet<WoodYieldItem> woodYieldData, Map<LandCoverChangeKey, Double> conversionCosts,
+ RasterSet<WoodYieldItem> woodYieldData, RasterSet<SolarPotentialItem> solarPotentialData,
+ Map<LandCoverChangeKey, Double> conversionCosts,
GlobalPrice carbonPrice, Map<WoodType, GlobalPrice> globalWoodPrices) {
// project demand
@@ -129,7 +131,7 @@ public class CountryAgent extends AbstractCountryAgent {
}
// optimize areas and intensity
- GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, woodYieldData, carbonFluxData, conversionCosts);
+ GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, woodYieldData, carbonFluxData, solarPotentialData, conversionCosts);
GamsRasterOptimiser opti = new GamsRasterOptimiser(input, yieldClusters);
LogWriter.println("Running " + country.getName() + ", currentTimestep " + currentTimestep, 3);
@@ -151,7 +153,8 @@ public class CountryAgent extends AbstractCountryAgent {
}
private GamsRasterInput getGamsRasterInput(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces,
- RasterSet<WoodYieldItem> woodYieldData, RasterSet<CarbonFluxItem> carbonFluxData, Map<LandCoverChangeKey, Double> conversionCosts) {
+ RasterSet<WoodYieldItem> woodYieldData, RasterSet<CarbonFluxItem> carbonFluxData,
+ RasterSet<SolarPotentialItem> solarPotentialData, Map<LandCoverChangeKey, Double> conversionCosts) {
Map<CropType, TradeConstraint> importConstraints = new HashMap<CropType, TradeConstraint>();
Map<CropType, Double> currentExportRestrictions = (exportRestrictions == null) ? null : exportRestrictions.get(currentTimestep.getYear());
@@ -271,7 +274,7 @@ public class CountryAgent extends AbstractCountryAgent {
previousGamsRasterOutput.getCropUsageData(), getTotalDemandFraction(), subsidyRates, currentGen2EcDemand, currentCarbonDemand, currentCarbonPrice, carbonTradeConstraint,
currentWoodDemand, currentWoodPrices, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData());
GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData,
- woodYieldData, carbonFluxData, conversionCosts, countryLevelInputs);
+ woodYieldData, carbonFluxData, solarPotentialData, conversionCosts, countryLevelInputs);
return input;
}
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index 6fd55c28..015f88ea 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -26,6 +26,8 @@ import ac.ed.lurg.landuse.*;
import ac.ed.lurg.shock.ExportRestrictionManager;
import ac.ed.lurg.shock.MinMaxNetImportManager;
import ac.ed.lurg.shock.MinMaxNetImportManager.LimitType;
+import ac.ed.lurg.solar.SolarPotentialItem;
+import ac.ed.lurg.solar.SolarPotentialReader;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.WoodType;
import ac.ed.lurg.utils.LogWriter;
@@ -121,7 +123,8 @@ public class CountryAgentManager {
}
public void determineProductionForAll(Timestep timestep, YieldRaster yieldSurfaces, IrrigationRasterSet currentIrrigationData,
- CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData, Map<LandCoverChangeKey, Double> conversionCosts) {
+ CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData,
+ RasterSet<SolarPotentialItem> solarPotentialData, Map<LandCoverChangeKey, Double> conversionCosts) {
for (AbstractCountryAgent aca : countryAgents) {
aca.setCurrentTimestep(timestep);
@@ -145,7 +148,7 @@ public class CountryAgentManager {
try {
ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), carbonFluxData, woodYieldData,
- conversionCosts, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrices());
+ solarPotentialData, conversionCosts, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrices());
// update global rasters
globalLandUseRaster.putAll(ca.getLandUses());
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index fc619fed..0b9f6627 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -9,6 +9,7 @@ import ac.ed.lurg.forestry.WoodYieldData;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandCoverChangeKey;
import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.solar.SolarPotentialItem;
import ac.ed.lurg.yield.YieldResponsesItem;
public class GamsLocationInput {
@@ -20,6 +21,7 @@ public class GamsLocationInput {
private Map<Integer, ? extends IrrigationItem> irrigationCosts;
private Map<Integer, ? extends CarbonFluxItem> carbonFluxes;
private Map<Integer, ? extends WoodYieldData> woodYields;
+ private Map<Integer, ? extends SolarPotentialItem> solarPotentials;
private Map<LandCoverChangeKey, Double> conversionCosts;
private GamsCountryInput countryInput;
@@ -27,6 +29,7 @@ public class GamsLocationInput {
Map<Integer, ? extends LandUseItem> previousLandUse,
Map<Integer, ? extends IrrigationItem> irrigationCosts, Map<Integer, ? extends CarbonFluxItem> carbonFluxes,
Map<Integer, ? extends WoodYieldData> woodYields,
+ Map<Integer, ? extends SolarPotentialItem> solarPotentials,
Map<LandCoverChangeKey, Double> conversionCosts, GamsCountryInput countryInput) {
super();
this.timestep = timestep;
@@ -36,6 +39,7 @@ public class GamsLocationInput {
this.irrigationCosts = irrigationCosts;
this.carbonFluxes = carbonFluxes;
this.woodYields = woodYields;
+ this.solarPotentials = solarPotentials;
this.conversionCosts = conversionCosts;
this.countryInput = countryInput;
}
@@ -59,7 +63,11 @@ public class GamsLocationInput {
public Map<Integer, ? extends WoodYieldData> getWoodYields() {
return woodYields;
}
-
+
+ public Map<Integer, ? extends SolarPotentialItem> getSolarPotentials() {
+ return solarPotentials;
+ }
+
public Map<LandCoverChangeKey, Double> getConversionCosts() {
return conversionCosts;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 6ee50225..6ad2b85b 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -11,6 +11,8 @@ import java.util.*;
import java.util.Map.Entry;
import ac.ed.lurg.landuse.*;
+import ac.ed.lurg.solar.AgrivoltaicsYieldManager;
+import ac.ed.lurg.solar.SolarPotentialItem;
import ac.ed.lurg.types.*;
import com.gams.api.GAMSDatabase;
import com.gams.api.GAMSException;
@@ -147,8 +149,11 @@ public class GamsLocationOptimiser {
Vector<String> v = new Vector<>();
v.add(crop.getGamsName());
LandProtectionType lpType = landKey.getLpType();
+ LandCoverType lcType = landKey.getLcType();
- if (lpType.equals(LandProtectionType.UNPROTECTED)) {
+ if (lcType.equals(LandCoverType.AGRIVOLTAICS)) {
+ v.add(FarmingType.AGRIVOLTAICS.getName());
+ } else if (lpType.equals(LandProtectionType.UNPROTECTED)) {
v.add(FarmingType.CONVENTIONAL.getName());
} else if (lpType.equals(LandProtectionType.PROTECTED)) {
v.add(FarmingType.RESTRICTED.getName());
@@ -488,6 +493,52 @@ public class GamsLocationOptimiser {
}
}
}
+
+ // Solar
+ GAMSParameter solarEnergyP = inDB.addParameter("solarEnergyDensity", 2);
+ GAMSParameter solarCostP = inDB.addParameter("solarCost", 2);
+ GAMSParameter avYieldFactorP = inDB.addParameter("agrivoltaicsYieldFactor", 3);
+ AgrivoltaicsYieldManager agrivoltParamManager = AgrivoltaicsYieldManager.getInstance();
+
+ for (Entry<Integer, ? extends SolarPotentialItem> entry : inputData.getSolarPotentials().entrySet()) {
+ Integer locationId = entry.getKey();
+ String locString = Integer.toString(locationId);
+ SolarPotentialItem solarPot = entry.getValue();
+ Vector<String> v = new Vector<>();
+
+ v.add(LandCoverType.PHOTOVOLTAICS.getName());
+ v.add(locString);
+ setGamsParamValue(solarEnergyP.addRecord(v), solarPot.getEnergyDensityPV(), -1);
+ setGamsParamValue(solarCostP.addRecord(v), solarPot.getCostPV(), -1);
+
+ v.clear();
+ v.add(LandCoverType.AGRIVOLTAICS.getName());
+ v.add(locString);
+ setGamsParamValue(solarEnergyP.addRecord(v), solarPot.getEnergyDensityAV(), -1);
+ setGamsParamValue(solarCostP.addRecord(v), solarPot.getCostAV(), -1);
+
+ for (CropType crop : CropType.getNonMeatTypes()) {
+ if (crop.equals(CropType.SETASIDE)) {
+ continue;
+ }
+ for (FarmingType farmingType : FarmingType.values()) {
+ v.clear();
+ v.add(crop.getGamsName());
+ v.add(farmingType.getName());
+ v.add(locString);
+ if (farmingType.equals(FarmingType.AGRIVOLTAICS)){
+ setGamsParamValue(avYieldFactorP.addRecord(v), agrivoltParamManager.getYieldFactor(crop, solarPot.getAgrivoltaicsGCR()), -1);
+ } else {
+ // Not agrivoltaics so yields not affected
+ setGamsParamValue(avYieldFactorP.addRecord(v), 1.0, -1);
+ }
+ v.clear();
+ }
+ }
+ }
+
+ addScalar(inDB, "energycropsEnergyDensity", ModelConfig.ENERGYCROPS_ENERGY_DENSITY * ModelConfig.BIOENERGY_GENERATION_EFFICIENCY, -1);
+ addScalar(inDB, "bioenergyGenerationCost", ModelConfig.BIOENERGY_GENERATION_COST, -1);
// Land cover conversion cost
GAMSParameter conversionCostP = inDB.addParameter("conversionCost", 3);
@@ -601,22 +652,24 @@ public class GamsLocationOptimiser {
for (GAMSVariableRecord rec : varAreas) {
String itemName = rec.getKey(0);
- String farmingType = rec.getKey(1);
+ String farmingTypeName = rec.getKey(1);
String locationName = rec.getKey(2);
area = rec.getLevel();
if (area == 0) {
continue;
}
- fertIntensity = varFertIntensities.findRecord(itemName, farmingType, locationName).getLevel();
- irrigIntensity = varIrrigIntensities.findRecord(itemName, farmingType, locationName).getLevel();
- otherIntensity = varOtherIntensities.findRecord(itemName, farmingType, locationName).getLevel();
- yield = varYields.findRecord(itemName, farmingType, locationName).getLevel();
- unitCost = varUnitEnergies.findRecord(itemName, farmingType, locationName).getLevel();
+ fertIntensity = varFertIntensities.findRecord(itemName, farmingTypeName, locationName).getLevel();
+ irrigIntensity = varIrrigIntensities.findRecord(itemName, farmingTypeName, locationName).getLevel();
+ otherIntensity = varOtherIntensities.findRecord(itemName, farmingTypeName, locationName).getLevel();
+ yield = varYields.findRecord(itemName, farmingTypeName, locationName).getLevel();
+ unitCost = varUnitEnergies.findRecord(itemName, farmingTypeName, locationName).getLevel();
int locId = Integer.parseInt(locationName);
CropType cropType = CropType.getForGamsName(itemName);
- LandCoverType lcType = cropType.equals(CropType.PASTURE) ? LandCoverType.PASTURE : LandCoverType.CROPLAND;
- LandProtectionType lpType = farmingType.equals("restricted") ? LandProtectionType.PROTECTED : LandProtectionType.UNPROTECTED;
+ FarmingType farmingType = FarmingType.getForName(farmingTypeName);
+
+ LandCoverType lcType = farmingType.equals(FarmingType.AGRIVOLTAICS) ? LandCoverType.AGRIVOLTAICS : cropType.equals(CropType.PASTURE) ? LandCoverType.PASTURE : LandCoverType.CROPLAND;
+ LandProtectionType lpType = farmingType.equals(FarmingType.RESTRICTED) ? LandProtectionType.PROTECTED : LandProtectionType.UNPROTECTED;
LandKey landKey = new LandKey(lcType, lpType);
LandUseItem landUseItem = landUses.lazyGet(locId); // returns landUseItem for location. If does not exist, creates new one
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index f272fc8c..6919f4d5 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -8,6 +8,7 @@ import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.LandCoverChangeKey;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.IrrigationItem;
+import ac.ed.lurg.solar.SolarPotentialItem;
import ac.ed.lurg.yield.YieldRaster;
import ac.sac.raster.RasterSet;
@@ -19,11 +20,12 @@ public class GamsRasterInput {
private RasterSet<IrrigationItem> irrigationCost;
private RasterSet<WoodYieldItem> woodYields;
private RasterSet<CarbonFluxItem> carbonFluxes;
+ private RasterSet<SolarPotentialItem> solarPotentialData;
private Map<LandCoverChangeKey, Double> conversionCosts;
private GamsCountryInput countryInput;
public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses, RasterSet<IrrigationItem> irrigationCost,
- RasterSet<WoodYieldItem> woodYields, RasterSet<CarbonFluxItem> carbonFluxes,
+ RasterSet<WoodYieldItem> woodYields, RasterSet<CarbonFluxItem> carbonFluxes, RasterSet<SolarPotentialItem> solarPotentialData,
Map<LandCoverChangeKey, Double> conversionCosts, GamsCountryInput countryInput) {
super();
this.timestep = timestep;
@@ -32,6 +34,7 @@ public class GamsRasterInput {
this.irrigationCost = irrigationCost;
this.woodYields = woodYields;
this.carbonFluxes = carbonFluxes;
+ this.solarPotentialData = solarPotentialData;
this.conversionCosts = conversionCosts;
this.countryInput = countryInput;
}
@@ -55,7 +58,11 @@ public class GamsRasterInput {
public RasterSet<CarbonFluxItem> getCarbonFluxes() {
return carbonFluxes;
}
-
+
+ public RasterSet<SolarPotentialItem> getSolarPotentialData() {
+ return solarPotentialData;
+ }
+
public Map<LandCoverChangeKey, Double> getConversionCosts() {
return conversionCosts;
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 58ea51f2..0292d606 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -8,6 +8,7 @@ import ac.ed.lurg.carbon.CarbonFluxItem;
import ac.ed.lurg.forestry.WoodYieldData;
import ac.ed.lurg.forestry.WoodYieldItem;
import ac.ed.lurg.landuse.*;
+import ac.ed.lurg.solar.SolarPotentialItem;
import ac.ed.lurg.types.*;
import ac.ed.lurg.utils.LazyTreeMap;
import ac.ed.lurg.utils.LogWriter;
@@ -231,6 +232,8 @@ public class GamsRasterOptimiser {
RasterSet<WoodYieldItem> woodYieldRaster = rasterInputData.getWoodYields();
+ RasterSet<SolarPotentialItem> solarPotentialRaster = rasterInputData.getSolarPotentialData();
+
// look for inconsistencies
for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
YieldResponsesItem yresp = entry.getValue();
@@ -262,6 +265,11 @@ public class GamsRasterOptimiser {
LazyTreeMap<Integer, WoodYieldData> aggregatedWoodYields = new LazyTreeMap<Integer, WoodYieldData>() {
protected WoodYieldData createValue() { return new WoodYieldData(); }
};
+ LazyTreeMap<Integer, SolarPotentialItem> aggregatedSolarPotentials = new LazyTreeMap<Integer, SolarPotentialItem>() {
+ protected SolarPotentialItem createValue() {
+ return new SolarPotentialItem();
+ }
+ };
for (RasterKey key : landUseRaster.keySet()) {
int clusterId = mapping.get(key).getInt();
@@ -271,12 +279,14 @@ public class GamsRasterOptimiser {
IrrigationItem irrigItem = irrigRaster.get(key);
WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
CarbonFluxItem carbonFluxItem = carbonFluxRaster.get(key);
+ SolarPotentialItem solarPotentialItem = solarPotentialRaster.get(key);
YieldResponsesItem aggYResp = aggregatedYields.lazyGet(clusterId);
LandUseItem aggLandUse = aggregatedAreas.lazyGet(clusterId);
IrrigationItem aggIrrig = aggregatedIrrigCosts.lazyGet(clusterId);
WoodYieldData aggWYield = aggregatedWoodYields.lazyGet(clusterId);
CarbonFluxItem aggCFlux = aggregatedCarbonFluxes.lazyGet(clusterId);
+ SolarPotentialItem aggSolarPot = aggregatedSolarPotentials.lazyGet(clusterId);
double suitableAreaThisTime = landUseItem.getSuitableArea();
double suitableAreaSoFar = aggLandUse.getSuitableArea();
@@ -321,6 +331,23 @@ public class GamsRasterOptimiser {
}
}
+ if (ModelConfig.IS_PHOTOVOLTAICS_ON || ModelConfig.IS_AGRIVOLTAICS_ON) {
+ double areaThisTime = landUseItem.getTotalLandArea();
+ double areaSoFar = landUseItem.getTotalLandArea();
+
+ double pvEnergyAgg = aggregateMean(aggSolarPot.getEnergyDensityPV(), areaSoFar, solarPotentialItem.getEnergyDensityPV(), areaThisTime);
+ double avEnergyAgg = aggregateMean(aggSolarPot.getEnergyDensityAV(), areaSoFar, solarPotentialItem.getEnergyDensityAV(), areaThisTime);
+ double pvCostAgg = aggregateMean(aggSolarPot.getCostPV(), areaSoFar, solarPotentialItem.getCostPV(), areaThisTime);
+ double avCostAgg = aggregateMean(aggSolarPot.getCostAV(), areaSoFar, solarPotentialItem.getCostAV(), areaThisTime);
+ double avGCRAgg = aggregateMean(aggSolarPot.getAgrivoltaicsGCR(), areaSoFar, solarPotentialItem.getAgrivoltaicsGCR(), areaThisTime);
+
+ aggSolarPot.setEnergyDensityPV(pvEnergyAgg);
+ aggSolarPot.setEnergyDensityAV(avEnergyAgg);
+ aggSolarPot.setCostPV(pvCostAgg);
+ aggSolarPot.setCostAV(avCostAgg);
+ aggSolarPot.setAgrivoltaicsGCR(avGCRAgg);
+ }
+
// Crops yields and area fractions
for (CropType crop : CropType.getNonMeatTypes()) {
if (crop.equals(CropType.SETASIDE)) {
@@ -390,7 +417,7 @@ public class GamsRasterOptimiser {
}
return new GamsLocationInput(rasterInputData.getTimestep(), locations, aggregatedYields, aggregatedAreas, aggregatedIrrigCosts,
- aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput());
+ aggregatedCarbonFluxes, aggregatedWoodYields, aggregatedSolarPotentials, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput());
}
private void logWarningWithCoord(String message, RasterKey key, YieldRaster yieldRaster, CropType crop) {
diff --git a/src/ac/ed/lurg/landuse/ConversionCostReader.java b/src/ac/ed/lurg/landuse/ConversionCostReader.java
index 4d1de464..9681936b 100644
--- a/src/ac/ed/lurg/landuse/ConversionCostReader.java
+++ b/src/ac/ed/lurg/landuse/ConversionCostReader.java
@@ -70,31 +70,37 @@ public class ConversionCostReader {
}
LandCoverChangeKey key = new LandCoverChangeKey(fromLc, toLc);
- switch(fromLc) {
- case CROPLAND:
- conversionCosts.put(key, ModelConfig.getAdjParam("CROPLAND_CONVERSION_COST") * costAdjFactor);
- break;
- case PASTURE:
- conversionCosts.put(key, ModelConfig.getAdjParam("PASTURE_CONVERSION_COST") * costAdjFactor);
- break;
- case TIMBER_FOREST:
- conversionCosts.put(key, ModelConfig.getAdjParam("FOREST_CONVERSION_COST") * costAdjFactor);
- break;
- case CARBON_FOREST:
- conversionCosts.put(key, ModelConfig.getAdjParam("FOREST_CONVERSION_COST") * costAdjFactor
- + ModelConfig.CARBON_FOREST_CONVERSION_PENALTY);
- break;
- case NATURAL:
- conversionCosts.put(key, ModelConfig.getAdjParam("NATURAL_CONVERSION_COST") * costAdjFactor);
- break;
- default:
- conversionCosts.put(key, 0.2 * costAdjFactor);
- break;
+ switch (fromLc) {
+ case CROPLAND:
+ conversionCosts.put(key, ModelConfig.getAdjParam("CROPLAND_CONVERSION_COST") * costAdjFactor);
+ break;
+ case PASTURE:
+ conversionCosts.put(key, ModelConfig.getAdjParam("PASTURE_CONVERSION_COST") * costAdjFactor);
+ break;
+ case PHOTOVOLTAICS:
+ conversionCosts.put(key, ModelConfig.getAdjParam("SOLAR_CONVERSION_COST") * costAdjFactor);
+ break;
+ case AGRIVOLTAICS:
+ conversionCosts.put(key, ModelConfig.getAdjParam("SOLAR_CONVERSION_COST") * costAdjFactor);
+ break;
+ case TIMBER_FOREST:
+ conversionCosts.put(key, ModelConfig.getAdjParam("FOREST_CONVERSION_COST") * costAdjFactor);
+ break;
+ case CARBON_FOREST:
+ conversionCosts.put(key, ModelConfig.getAdjParam("FOREST_CONVERSION_COST") * costAdjFactor
+ + ModelConfig.CARBON_FOREST_CONVERSION_PENALTY);
+ break;
+ case NATURAL:
+ conversionCosts.put(key, ModelConfig.getAdjParam("NATURAL_CONVERSION_COST") * costAdjFactor);
+ break;
+ default:
+ conversionCosts.put(key, 0.2 * costAdjFactor);
+ break;
}
// Agricultural expansion cost factor
- if (fromLc.equals(LandCoverType.NATURAL) && (toLc.equals(LandCoverType.CROPLAND) || toLc.equals(LandCoverType.PASTURE))) {
+ if (fromLc.equals(LandCoverType.NATURAL) && (toLc.equals(LandCoverType.CROPLAND) || toLc.equals(LandCoverType.PASTURE) || toLc.equals(LandCoverType.AGRIVOLTAICS))) {
double baseCost = conversionCosts.get(key);
conversionCosts.put(key, baseCost * ModelConfig.getAdjParam("AGRI_LAND_EXPANSION_COST_FACTOR"));
}
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index a2a7f88f..7948ec3d 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -25,10 +25,14 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
setProtectedFraction(landCover.getProtectedFraction(landCoverAreas));
AgricultureDetail cropDetail = (AgricultureDetail) getLandCoverDetail(new LandKey(LandCoverType.CROPLAND, LandProtectionType.UNPROTECTED));
+ AgricultureDetail avDetail = (AgricultureDetail) getLandCoverDetail(new LandKey(LandCoverType.AGRIVOLTAICS, LandProtectionType.UNPROTECTED));
for (CropType crop : CropType.getCropsLessPasture()) {
cropDetail.setCropFraction(crop, landCover.getCropFraction(crop));
cropDetail.setIntensity(crop, new Intensity(0.5, 0.5, 0.5));
+
+ avDetail.setCropFraction(crop, 0);
+ avDetail.setIntensity(crop, new Intensity(0.5, 0.5, 0.5));
}
AgricultureDetail pastureDetail = (AgricultureDetail) getLandCoverDetail(new LandKey(LandCoverType.PASTURE, LandProtectionType.UNPROTECTED));
@@ -71,6 +75,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
return new AgricultureDetail(0);
case PASTURE:
return new AgricultureDetail(0);
+ case AGRIVOLTAICS:
+ return new AgricultureDetail(0);
+ case PHOTOVOLTAICS:
+ return new LandCoverDetail(0);
case TIMBER_FOREST:
return new ForestDetail(0);
case CARBON_FOREST:
@@ -100,6 +108,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
Set<LandKey> keys = new HashSet<>();
keys.addAll(getLandKeysForLandCover(LandCoverType.CROPLAND));
keys.addAll(getLandKeysForLandCover(LandCoverType.PASTURE));
+ keys.addAll(getLandKeysForLandCover(LandCoverType.AGRIVOLTAICS));
return keys;
}
diff --git a/src/ac/ed/lurg/output/LandUseOutputer.java b/src/ac/ed/lurg/output/LandUseOutputer.java
index ff6ae8c1..2613aac6 100644
--- a/src/ac/ed/lurg/output/LandUseOutputer.java
+++ b/src/ac/ed/lurg/output/LandUseOutputer.java
@@ -9,6 +9,7 @@ import java.util.Map.Entry;
import ac.ed.lurg.landuse.*;
import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.types.FarmingType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.LandProtectionType;
import ac.ed.lurg.utils.LogWriter;
@@ -58,8 +59,17 @@ public class LandUseOutputer extends AbstractLandUseOutputer {
cropData.add(intensity==null ? 0.0 : intensity.getOtherIntensity());
cropData.add(intensity==null ? 0.0 : intensity.getYield());
+ String farmingType = "ERROR";
+ if (landKey.getLcType().equals(LandCoverType.AGRIVOLTAICS)) {
+ farmingType = FarmingType.AGRIVOLTAICS.getName();
+ } else if (landKey.getLpType().equals(LandProtectionType.PROTECTED)) {
+ farmingType = FarmingType.RESTRICTED.getName();
+ } else if (landKey.getLpType().equals(LandProtectionType.UNPROTECTED)) {
+ farmingType = FarmingType.CONVENTIONAL.getName();
+ }
+
StringBuilder sbData = new StringBuilder(String.format("%.2f %.2f %s %s", lon, lat,
- landKey.getLpType().getGamsName(), crop.getGamsName()));
+ farmingType, crop.getGamsName()));
for (Double d : cropData) {
String formatStr = d > 0 ? " %.8f" : " %.0f";
@@ -95,7 +105,7 @@ public class LandUseOutputer extends AbstractLandUseOutputer {
try {
String fileName = outputDir.getPath() + File.separator + "LandCover.txt";
writer = new BufferedWriter(new FileWriter(fileName));
- writer.write("Lon Lat Protection TotalArea Cropland Pasture TimberForest CarbonForest UnmanagedForest OtherNatural Barren Urban");
+ writer.write("Lon Lat Protection TotalArea Cropland Pasture TimberForest CarbonForest UnmanagedForest OtherNatural Photovoltaics Agrivoltaics Barren Urban");
writer.newLine();
for (Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) {
@@ -116,6 +126,8 @@ public class LandUseOutputer extends AbstractLandUseOutputer {
double unmanagedForestFraction = detail.getUnmanagedForestFraction();
appendData(sbData, totalNatural * unmanagedForestFraction);
appendData(sbData, totalNatural * (1 - unmanagedForestFraction));
+ appendData(sbData, item.getLandCoverArea(new LandKey(LandCoverType.PHOTOVOLTAICS, lpType)));
+ appendData(sbData, item.getLandCoverArea(new LandKey(LandCoverType.AGRIVOLTAICS, lpType)));
appendData(sbData, item.getLandCoverArea(new LandKey(LandCoverType.BARREN, lpType)));
appendData(sbData, item.getLandCoverArea(new LandKey(LandCoverType.URBAN, lpType)));
diff --git a/src/ac/ed/lurg/solar/AgrivoltaicsYieldManager.java b/src/ac/ed/lurg/solar/AgrivoltaicsYieldManager.java
new file mode 100644
index 00000000..7b42b039
--- /dev/null
+++ b/src/ac/ed/lurg/solar/AgrivoltaicsYieldManager.java
@@ -0,0 +1,42 @@
+package ac.ed.lurg.solar;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.StringTabularReader;
+
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+public class AgrivoltaicsYieldManager {
+
+ private final StringTabularReader reader;
+
+ private static AgrivoltaicsYieldManager apm;
+ public static AgrivoltaicsYieldManager getInstance() {
+ if (apm == null)
+ apm = new AgrivoltaicsYieldManager();
+
+ return apm;
+ }
+
+ private AgrivoltaicsYieldManager() {
+ reader = new StringTabularReader(",", new String[]{"Crop", "ParamA", "ParamB"});
+ reader.read(ModelConfig.AGRIVOLTAICS_YIELD_PARAM_FILE);
+
+ }
+
+ private double getParam(CropType crop, String param) {
+ Map<String, String> query = new HashMap<>();
+ query.put("Crop", crop.getGamsName());
+ List<Map<String, String>> queryResult = reader.query(query);
+ String result = queryResult.get(0).get(param);
+ return Double.parseDouble(result);
+ }
+
+ public double getYieldFactor(CropType crop, double avGCR) {
+ double shading = avGCR + 0.05;
+ return Math.exp(getParam(crop, "ParamA") * shading + getParam(crop, "ParamB") * Math.pow(shading, 2)) * ModelConfig.AV_YIELD_FACTOR;
+ }
+
+}
diff --git a/src/ac/ed/lurg/solar/SolarPotentialItem.java b/src/ac/ed/lurg/solar/SolarPotentialItem.java
new file mode 100644
index 00000000..60a3ba26
--- /dev/null
+++ b/src/ac/ed/lurg/solar/SolarPotentialItem.java
@@ -0,0 +1,97 @@
+package ac.ed.lurg.solar;
+
+import ac.ed.lurg.ModelConfig;
+import ac.sac.raster.RasterItem;
+import ac.sac.raster.RasterKey;
+
+public class SolarPotentialItem implements RasterItem {
+ private double pvPotential; // in MWh/ha
+ private double energyDensityPV = Double.NaN; // MWh/ha
+ private double energyDensityAV = Double.NaN; // MWh/ha
+ private double costPV = Double.NaN; // $1000/ha
+ private double costAV = Double.NaN; // $1000/ha
+ private double optimalGCR; // optimal Ground Coverage Ratio (PV area) / (land area)
+ private double agrivoltaicsGCR;
+
+ public void setPvPotential(double pvPotential) {
+ this.pvPotential = pvPotential;
+ }
+
+ public void calcOptimalGCR(double latitude) {
+ // https://doi.org/10.1016/j.solener.2023.04.038
+ // Optimal panel spacing and therefore Ground Coverage ratio depends on latitude
+ this.optimalGCR = -0.55 / (1 + Math.exp(-0.138 * (latitude - 43.4))) + 0.71;
+ this.agrivoltaicsGCR = optimalGCR * ModelConfig.AV_GCR;
+ }
+
+ public void calcEnergyDensity() {
+ // MWh per ha including spacing between panels (GCR) and additional infrastructure area (GSR)
+ this.energyDensityPV = pvPotential * ModelConfig.PV_EFFICIENCY * optimalGCR * ModelConfig.PV_GSR;
+ // Ground Coverage Ratio (GCR) for AV is fixed.
+ this.energyDensityAV = pvPotential * ModelConfig.PV_EFFICIENCY * agrivoltaicsGCR * ModelConfig.PV_GSR;
+ }
+
+ public void calcCostPerHectarePV() { // photovoltaics cost $1000/ha
+ // Capacity in kW/ha of total land used. Factor 10000 to convert from kW per m2 to hectare.
+ double capacity = ModelConfig.PV_EFFICIENCY * optimalGCR * ModelConfig.PV_GSR * 10000;
+ // Cost
+ this.costPV = capacity * (ModelConfig.PV_INSTALLED_COST / ModelConfig.PV_LIFESPAN + ModelConfig.PV_MAINTENANCE_COST);
+ }
+
+ public void getCostPerHectareAV() { // agrivoltaics cost $1000/ha
+ // Capacity in kW/ha of total land used. Factor 10000 to convert from kW per m2 to hectare.
+ double capacity = ModelConfig.PV_EFFICIENCY * agrivoltaicsGCR * ModelConfig.PV_GSR * 10000;
+ // Cost
+ this.costAV = capacity * (ModelConfig.PV_INSTALLED_COST / ModelConfig.PV_LIFESPAN + ModelConfig.PV_MAINTENANCE_COST) * ModelConfig.AV_COST_FACTOR;
+ }
+
+ public double getEnergyDensityPV() {
+ return ModelConfig.IS_PHOTOVOLTAICS_ON ? energyDensityPV : 0;
+ }
+
+ public void setEnergyDensityPV(double energyDensityPV) {
+ this.energyDensityPV = energyDensityPV;
+ }
+
+ public double getEnergyDensityAV() {
+ return ModelConfig.IS_AGRIVOLTAICS_ON ? energyDensityAV : 0;
+ }
+
+ public void setEnergyDensityAV(double energyDensityAV) {
+ this.energyDensityAV = energyDensityAV;
+ }
+
+ public double getCostPV() {
+ return costPV;
+ }
+
+ public void setCostPV(double costPV) {
+ this.costPV = costPV;
+ }
+
+ public double getCostAV() {
+ return costAV;
+ }
+
+ public void setCostAV(double costAV) {
+ this.costAV = costAV;
+ }
+
+ public double getAgrivoltaicsGCR() {
+ return agrivoltaicsGCR;
+ }
+
+ public void setAgrivoltaicsGCR(double agrivoltaicsGCR) {
+ this.agrivoltaicsGCR = agrivoltaicsGCR;
+ }
+
+ public static SolarPotentialItem getDefault() {
+ SolarPotentialItem item = new SolarPotentialItem();
+ item.setPvPotential(0);
+ item.setEnergyDensityPV(0);
+ item.setEnergyDensityAV(0);
+ item.setCostPV(0);
+ item.setCostAV(0);
+ return item;
+ }
+}
diff --git a/src/ac/ed/lurg/solar/SolarPotentialReader.java b/src/ac/ed/lurg/solar/SolarPotentialReader.java
new file mode 100644
index 00000000..b8967ae6
--- /dev/null
+++ b/src/ac/ed/lurg/solar/SolarPotentialReader.java
@@ -0,0 +1,37 @@
+package ac.ed.lurg.solar;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.Timestep;
+import ac.sac.raster.AbstractTabularRasterReader;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+
+import java.io.File;
+import java.util.Map;
+
+public class SolarPotentialReader extends AbstractTabularRasterReader<SolarPotentialItem> {
+
+ private static final int MIN_COL = 3;
+ private static final String DELIM = ",";
+ private static final double CONVERSION_FACTOR = 10; // convert from kWh/m2 to MWh/ha
+
+ public SolarPotentialReader(RasterSet<SolarPotentialItem> dataset) {
+ super(DELIM, MIN_COL, dataset);
+ }
+
+ public void readData(Timestep timestep) {
+ int dataYear = Math.floorDiv(timestep.getYear(), ModelConfig.SOLAR_POTENTIAL_DATA_STEP_SIZE) * ModelConfig.SOLAR_POTENTIAL_DATA_STEP_SIZE;
+ String filePath = ModelConfig.SOLAR_POTENTIAL_DIR_BASE + File.separator + ModelConfig.SOLAR_POTENTIAL_DIR_TOP +
+ File.separator + dataYear + File.separator + ModelConfig.SOLAR_POTENTIAL_FILENAME;
+ getRasterDataFromFile(filePath);
+ }
+
+ @Override
+ protected void setData(RasterKey key, SolarPotentialItem item, Map<String, Double> rowValues) {
+ item.setPvPotential(getValueForCol(rowValues, "PVYield") * CONVERSION_FACTOR);
+ item.calcOptimalGCR(dataset.getYCoordin(key));
+ item.calcEnergyDensity();
+ item.calcCostPerHectarePV();
+ item.getCostPerHectareAV();
+ }
+}
diff --git a/src/ac/ed/lurg/types/FarmingType.java b/src/ac/ed/lurg/types/FarmingType.java
index e6cc4ce6..8fd854c8 100644
--- a/src/ac/ed/lurg/types/FarmingType.java
+++ b/src/ac/ed/lurg/types/FarmingType.java
@@ -7,7 +7,7 @@ public enum FarmingType {
CONVENTIONAL("conventional"),
RESTRICTED("restricted"),
- AGRIVOLTAIC("agrivoltaic");
+ AGRIVOLTAICS("agrivoltaics");
private String name;
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index bd47294d..1120f9c4 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -14,6 +14,8 @@ public enum LandCoverType {
NATURAL ("natural", true),
CROPLAND ("cropland", true),
PASTURE ("pasture", true),
+ PHOTOVOLTAICS("photovoltaics", true),
+ AGRIVOLTAICS("agrivoltaics", true),
BARREN ("barren", false),
URBAN("urban", false);
--
GitLab