From 2b5e9e64c32a29c83162ceb007d36af14898fcdd Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Wed, 11 Jan 2017 12:23:47 +0000
Subject: [PATCH] Changes for new LPJ data version. Also, some demand work.
---
GAMS/IntExtOpt.gms | 4 +-
GAMS/demand.gms | 64 ++++----
debug_config.properties | 14 +-
src/ac/ed/lurg/ModelConfig.java | 6 +-
src/ac/ed/lurg/Timestep.java | 2 +-
.../country/gams/GamsDemandOptimiser.java | 143 ++++++++++++++++++
6 files changed, 193 insertions(+), 40 deletions(-)
create mode 100644 src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index ee3f1546..46da5614 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -119,7 +119,8 @@ $gdxin
CROP_INCREASE_CALC(location)
CROP_DECREASE_CALC(location)
PASTURE_INCREASE_CONV_CALC(location)
- PASTURE_DECREASE_CONV_CALC(location)
+ PASTURE_DECREASE_CONV_CALC(location)
+ PASTURE_TOTAL_CHANGE_CONSTRAINT(location)
COST_EQ total cost objective function;
UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E= ( baseCost(crop) * 0.5 +
@@ -177,6 +178,7 @@ $gdxin
CROP_DECREASE_CALC(location) .. cropDecrease(location) =G= -sum(crop_less_pasture, area(crop_less_pasture, location) - previousArea(crop_less_pasture, location));
PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= area('pasture', location) - previousArea('pasture', location);
PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(area('pasture', location) - previousArea('pasture', location));
+ PASTURE_TOTAL_CHANGE_CONSTRAINT(location) .. pastureIncrease(location) -pastureDecrease(location) =G= area('pasture', location) - previousArea('pasture', location);
COST_EQ .. total_cost =E=
(
diff --git a/GAMS/demand.gms b/GAMS/demand.gms
index 5f5f83f1..16580822 100644
--- a/GAMS/demand.gms
+++ b/GAMS/demand.gms
@@ -3,49 +3,53 @@
ALIAS (com, com2);
ALIAS (com, com3);
- SCALAR income / 100 /;
+
+ SCALAR income;
+
+$gdxin %demandgdx%
+$load income
+$gdxin
TABLE subElasticities(com, com2)
ruminants wheat maize rice oilcrops pulses starchyRoots
monogastrics 0.7 0.15 0.1 0.1 0.1 0.1 0.1
ruminants 0.1 0.1 0.1 0.1 0.1 0.1
- wheat 0.6 0.6 0.4 0.4 0.4
- maize 0.6 0.4 0.4 0.4
- rice 0.4 0.4 0.4
+ wheat 0.6 0.6 0.3 0.3 0.3
+ maize 0.6 0.3 0.3 0.3
+ rice 0.3 0.3 0.3
oilcrops 0.5 0.5
pulses 0.5 ;
subElasticities(com, com2)$subElasticities(com2, com) = subElasticities(com2, com);
PARAMETER calories(com)
- / monogastrics 3,
- ruminants 5,
- wheat 1,
- maize 1,
+ / monogastrics 1.5,
+ ruminants 1.7,
+ wheat 1.02,
+ maize 0.95,
rice 1,
- oilcrops 2,
- pulses 1,
- starchyRoots 1 / ;
- calories(com) = 1;
+ oilcrops 1.3,
+ pulses 1.05,
+ starchyRoots 0.9 / ;
PARAMETER idealDiet(com)
- / monogastrics 1,
+ / monogastrics 1.1,
ruminants 1,
- wheat 1,
- maize 1,
- rice 1,
- oilcrops 1,
- pulses 1,
- starchyRoots 1 / ;
+ wheat 0.9,
+ maize 0.4,
+ rice 0.5,
+ oilcrops 0.6,
+ pulses 0.5,
+ starchyRoots 0.8 / ;
PARAMETER price(com)
- / monogastrics 4,
- ruminants 5,
- wheat 2,
- maize 2,
- rice 2,
- oilcrops 2,
- pulses 1,
+ / monogastrics 5,
+ ruminants 7.7,
+ wheat 1.4,
+ maize 1.25,
+ rice 1.5,
+ oilcrops 3,
+ pulses 1.5,
starchyRoots 1 / ;
PARAMETER comUtility(com);
@@ -105,7 +109,7 @@
total_cost_calc .. totalCost =e= sum(com, consumption(com) * price(com));
- net_utility_calc .. netUtility =e= ((totalCost/income)) * 25 + foodUtility;
+ net_utility_calc .. netUtility =e= ((totalCost/income)) * 30 + foodUtility;
MODEL DEMAND /ALL/ ;
@@ -121,4 +125,8 @@
PARAMETER check(com);
check(com) = 0;
-*display check;
\ No newline at end of file
+*display check;
+
+ Scalar ms 'model status', ss 'solve status';
+ ms=DEMAND.modelstat;
+ ss=DEMAND.solvestat;
\ No newline at end of file
diff --git a/debug_config.properties b/debug_config.properties
index 394a0d91..1e41da7c 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -1,12 +1,12 @@
BASE_DIR=/Users/peteralexander/Documents/R_Workspace/UNPLUM
#YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/yield_anpp_rcps-2016-01-29/rcp8p5_past_fert/output-2016-01-29
-YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/sam-03Nov2016
+YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/10Jan2017/combined
SPATIAL_DIR_NAME=halfdeg
CELL_SIZE_X=0.5
-ANPP_FILENAME=anpp.out.2011-2015
-YIELD_FILENAME=yield.out.2011-2015
-IRRIG_MAX_WATER_FILENAME=gsirrigation.out.2011-2015
+#ANPP_FILENAME=anpp.out
+#YIELD_FILENAME=yield.out
+#IRRIG_MAX_WATER_FILENAME=gsirrigation.out
MIN_FERT_AMOUNT=0
MID_FERT_AMOUNT=200
MAX_FERT_AMOUNT=1000
@@ -16,7 +16,7 @@ PASTURE_FERT_RESPONSE_FROM_LPJ=false
CLEANUP_GAMS_DIR=false
# Properties for testing
-CHANGE_YIELD_DATA_YEAR=false
+CHANGE_YIELD_DATA_YEAR=true
CHANGE_DEMAND_YEAR=false
DEBUG_LIMIT_COUNTRIES=false
DEBUG_COUNTRY_NAME=Spain
@@ -25,7 +25,7 @@ SSP_SCENARIO=SSP2_v9_130325
MAX_IMPORT_CHANGE=0.0
MARKET_ADJ_PRICE=false
-IS_CALIBRATION_RUN = true
+IS_CALIBRATION_RUN = false
#LAND_CHANGE_COST=.6
@@ -36,7 +36,7 @@ IS_CALIBRATION_RUN = true
#NUM_CEREAL_CATEGORIES = 10
#NUM_PASTURE_CATEGORIES = 3
-END_TIMESTEP=28
+END_TIMESTEP=18
TIMESTEP_SIZE=5
INTERPOLATE_OUTPUT_YEARS = false
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 43766901..e81c70d8 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -102,8 +102,8 @@ public class ModelConfig {
// yield data
public static final String YIELD_DIR = getProperty("YIELD_DIR");
public static final int LPJG_MONITOR_TIMEOUT_SEC = getIntProperty("LPJG_MONITOR_TIMEOUT", 60*60*2);
- public static final String ANPP_FILENAME = getProperty("ANPP_FILENAME", "anpp_plum.out");
- public static final String YIELD_FILENAME = getProperty("YIELD_FILENAME", "yield_plum.out");
+ public static final String ANPP_FILENAME = getProperty("ANPP_FILENAME", "anpp.out");
+ public static final String YIELD_FILENAME = getProperty("YIELD_FILENAME", "yield.out");
public static final boolean PASTURE_FERT_RESPONSE_FROM_LPJ = getBooleanProperty("PASTURE_FERT_RESPONSE_FROM_LPJ", true);;
public static final double CALIB_FACTOR_WHEAT = getDoubleProperty("CALIB_FACTOR_WHEAT", 0.966);
@@ -122,7 +122,7 @@ public class ModelConfig {
public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries.asc";
public static final String IRRIGATION_COST_FILE = SPATIAL_DATA_DIR + File.separator + "irrigation_cost.asc";
public static final String IRRIGATION_CONSTRAINT_FILE = SPATIAL_DATA_DIR + File.separator + "blue_water_available_pseudoCRU_rcp8p5_2004_2013_grid_allhdyro_mm.txt";
- public static final String IRRIG_MAX_WATER_FILENAME = getProperty("IRRIG_MAX_WATER_FILENAME", "gsirrigation_plum.out");
+ public static final String IRRIG_MAX_WATER_FILENAME = getProperty("IRRIG_MAX_WATER_FILENAME", "gsirrigation.out");
public static final String PROTECTED_AREAS_FILE = SPATIAL_DATA_DIR + File.separator + "protected_areas.txt";
// Output
diff --git a/src/ac/ed/lurg/Timestep.java b/src/ac/ed/lurg/Timestep.java
index 496697c7..50274b6c 100644
--- a/src/ac/ed/lurg/Timestep.java
+++ b/src/ac/ed/lurg/Timestep.java
@@ -88,6 +88,6 @@ public class Timestep {
if (ModelConfig.CHANGE_YIELD_DATA_YEAR)
return rootDir + File.separator + getYieldYear();
else
- return rootDir;
+ return rootDir + File.separator + ModelConfig.BASE_YEAR;
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java b/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
new file mode 100644
index 00000000..967b3ecf
--- /dev/null
+++ b/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
@@ -0,0 +1,143 @@
+package ac.ed.lurg.country.gams;
+
+import java.io.BufferedWriter;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.HashMap;
+import java.util.Map;
+
+import com.gams.api.GAMSDatabase;
+import com.gams.api.GAMSGlobals;
+import com.gams.api.GAMSJob;
+import com.gams.api.GAMSOptions;
+import com.gams.api.GAMSParameter;
+import com.gams.api.GAMSParameterRecord;
+import com.gams.api.GAMSVariable;
+import com.gams.api.GAMSVariableRecord;
+import com.gams.api.GAMSWorkspace;
+import com.gams.api.GAMSWorkspaceInfo;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.LogWriter;
+
+public class GamsDemandOptimiser {
+ private static final boolean DEBUG = true;
+
+ private double income;
+
+
+ public static void main(String[] args) {
+ BufferedWriter demandWriter = null;
+
+ try {
+ demandWriter = new BufferedWriter(new FileWriter(ModelConfig.OUTPUT_DIR + "/demand.out", false));
+
+ StringBuffer sbHeader = new StringBuffer("income");
+ for (CropType crop : CropType.getImportedTypes())
+ sbHeader.append("," + crop.getGamsName());
+
+ demandWriter.write(sbHeader.toString());
+ demandWriter.newLine();
+
+ for (int i=0; i<50; i++) {
+ double income = 20.0 + i * 2.0;
+ GamsDemandOptimiser theModel = new GamsDemandOptimiser(income);
+ Map<CropType, Double> results = theModel.run();
+
+ StringBuffer sbData = new StringBuffer(String.format("%.6f", income));
+ for (CropType crop : CropType.getImportedTypes()) {
+ Double d = results.get(crop);
+ sbData.append(String.format(",%.6f", d));
+ }
+
+ demandWriter.write(sbData.toString());
+ demandWriter.newLine();
+ }
+ }
+ catch (IOException e) {
+ LogWriter.print(e);
+ }
+ finally {
+ if (demandWriter != null) {
+ try {
+ demandWriter.close();
+ }
+ catch (IOException e) {
+ LogWriter.print(e);
+ }
+ }
+
+ }
+ }
+
+
+ public GamsDemandOptimiser(double income) {
+ this.income = income;
+ }
+
+ public Map<CropType, Double> run() {
+ File workingDirectory = new File(ModelConfig.TEMP_DIR, "GamsTmp");
+ workingDirectory.mkdir();
+
+ GAMSWorkspaceInfo wsInfo = new GAMSWorkspaceInfo();
+ wsInfo.setWorkingDirectory(workingDirectory.getAbsolutePath());
+ // wsInfo.setDebugLevel(DebugLevel.VERBOSE);
+
+ GAMSWorkspace ws = new GAMSWorkspace(wsInfo);
+ GAMSDatabase inDB = ws.addDatabase();
+
+ setupInDB(inDB);
+
+ GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.BASE_DIR + File.separator + "GAMS/demand.gms");
+ GAMSOptions opt = ws.addOptions();
+ // opt.setAllModelTypes("conopt");
+ opt.defines("demandgdx", inDB.getName());
+
+ long startTime = System.currentTimeMillis();
+ gamsJob.run(opt, inDB);
+ LogWriter.println("Took " + (System.currentTimeMillis() - startTime) + " ms to run");
+
+ return handleResults(gamsJob.OutDB());
+ }
+
+ private void setupInDB(GAMSDatabase inDB) {
+ addScalar(inDB, "income", income, 3);
+ }
+
+ private void addScalar(GAMSDatabase gamsDb, String recordName, double val, int places) {
+ GAMSParameter param = gamsDb.addParameter(recordName, 0);
+ setGamsParamValue(param.addRecord(), val, places);
+ if (DEBUG) LogWriter.println(recordName + ": " + val);
+ }
+
+ private void setGamsParamValue(GAMSParameterRecord param, double val, int places) {
+ double dOut = Math.round(val * Math.pow(10,places)) / Math.pow(10,places);
+ param.setValue(dOut);
+ }
+
+ private Map<CropType, Double> handleResults(GAMSDatabase outDB) {
+ int modelStatus = (int) outDB.getParameter("ms").findRecord().getValue();
+ LogWriter.println(String.format("\nModelstatus %s, Solvestatus %s",
+ GAMSGlobals.ModelStat.lookup( modelStatus ),
+ GAMSGlobals.SolveStat.lookup((int) outDB.getParameter("ss").findRecord().getValue()) ));
+
+ GAMSVariable varConsump = outDB.getVariable("consumption");
+
+ double consumpValue;
+ Map<CropType, Double> consumptionData = new HashMap<CropType, Double>();
+
+ for (GAMSVariableRecord rec : varConsump) {
+ String itemName = rec.getKeys()[0];
+ consumpValue = rec.getLevel();
+
+ CropType cropType = CropType.getForGamsName(itemName);
+
+ consumptionData.put(cropType, consumpValue);
+ if (DEBUG) LogWriter.println(String.format("\n%s:\tconsumption= %.3f", itemName, consumpValue));
+ }
+
+ return consumptionData;
+ }
+}
--
GitLab