diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 9ae48a424bf663cb24c893b93a44d65f31014547..5d8dca6383d9d316fee0560574c1f4082ed9348a 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -19,7 +19,8 @@
PARAMETER fertParam(crop, location) yield response to fertilizer parameter;
PARAMETER irrigParam(crop, location) yield response to irrigation parameter;
PARAMETER demand(all_types) in t;
- PARAMETER worldInputPrices(import_crop) average input energy from world exports used to determine if we should import or export;
+ PARAMETER worldExportPrices(import_crop) prices for exports;
+ PARAMETER worldImportPrices(import_crop) prices for imports;
PARAMETER maxNetImport(import_crop) maximum net import for each crop based on world market;
PARAMETER minNetImport(import_crop) minimum net import for each crop based on world market;
PARAMETER irrigCost(location) irrigation cost;
@@ -31,17 +32,14 @@
SCALAR landChangeEnergy energy required to add ha of agricultural land;
SCALAR minFeedRate minimum rate of feed for producing animal products (why is this needed?);
- SCALAR fitToPreviousAreas controlling parameter 0 means dont fit and 1 is to calibrate the cropAdj in the calibration year
-
*$gdxin "/Users/peteralexander/Documents/R_Workspace/temp/GamsTest/_gams_java_1091455539.gdx"
$gdxin %gdxincname%
$load location, suitableLandArea, previousArea, demand, landChangeEnergy
$load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
-$load fertParam, irrigParam, worldInputPrices, maxNetImport, minNetImport
-$load meatEfficency, minFeedRate, irrigCost, irrigMaxRate, irrigConstraint
-$load cropAdj, fitToPreviousAreas
-$gdxin
-
+$load fertParam, irrigParam, worldExportPrices, worldImportPrices, maxNetImport, minNetImport
+$load meatEfficency, minFeedRate, irrigCost, irrigMaxRate, irrigConstraint, cropAdj
+$gdxin
+
SCALAR delta use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf / 0.00000000001 /
SCALAR minDemandPerCereal / 0.1 /;
@@ -73,7 +71,8 @@ $gdxin
irrigI(crop, location) irrigation intensity for each crop - factor between 0 and 1
otherIntensity(crop, location)
feedAmount(crop) amount of feed for each crop - Mt
- netImportAmount(all_types) net imports of crops and meat - Mt
+ importAmount(all_types) imports of crops and meat - Mt
+ exportAmount(all_types) exports of crops and meat - Mt
yield(crop, location) yield per area for each crop - t per ha
unitEnergy(crop, location) energy per area for each crop - energy
net_supply(crop) supply after exports and feed
@@ -84,7 +83,8 @@ $gdxin
pastureDecrease(location)
energy total input energy - energy;
- POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount, agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease;
+ POSITIVE VARIABLE area, fertI, irrigI, otherIntensity, feedAmount, importAmount, exportAmount,
+ agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease;
EQUATIONS
UNIT_ENERGY_EQ(crop, location) energy per area
@@ -131,7 +131,7 @@ $gdxin
) * cropAdj(crop) * otherIntensity(crop, location);
- NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + netImportAmount(crop);
+ NET_SUPPLY_EQ(crop) .. net_supply(crop) =E= sum(location, area(crop, location) * yield(crop, location)) - feedAmount(crop) + importAmount(crop) - exportAmount(crop);
NON_CEREAL_DEMAND_CONSTRAINT(non_cereal_crop) .. net_supply(non_cereal_crop) =G= demand(non_cereal_crop) ;
@@ -139,7 +139,7 @@ $gdxin
TOTAL_CEREAL_DEMAND_CONSTRAINT .. sum(cereal_crop, net_supply(cereal_crop)) =G= demand('cereals');
- MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= demand('meat') - netImportAmount('meat');
+ MEAT_DEMAND_CONSTRAINT .. meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount(feed_crop)) =G= demand('meat') - importAmount('meat') - exportAmount('meat');
MAX_FERT_INTENSITY_CONSTRAINT(crop, location) .. fertI(crop, location) =L= 1;
MAX_IRRIG_INTENSITY_CONSTRAINT(crop, location) .. irrigI(crop, location) =L= 1;
@@ -149,11 +149,11 @@ $gdxin
NON_FEED_CROP_CONSTRAINT(not_feed_crop) .. feedAmount(not_feed_crop) =E= 0;
- MAX_NET_IMPORT_CONSTRAINT(import_crop) .. netImportAmount(import_crop) =L= maxNetImport(import_crop);
- MIN_NET_IMPORT_CONSTRAINT(import_crop) .. netImportAmount(import_crop) =G= minNetImport(import_crop);
- PASTURE_IMPORT_CONSTRAINT .. netImportAmount('pasture') =E= 0;
+ MAX_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =L= maxNetImport(import_crop);
+ MIN_NET_IMPORT_CONSTRAINT(import_crop) .. importAmount(import_crop) - exportAmount(import_crop) =G= minNetImport(import_crop);
+ PASTURE_IMPORT_CONSTRAINT .. importAmount('pasture') - exportAmount('pasture') =E= 0;
- MIN_FEED_CONSTRAINT .. sum(feed_crop_less_pasture, feedDM(feed_crop_less_pasture) * feedAmount(feed_crop_less_pasture)) =G= minFeedRate * (demand('meat') - netImportAmount('meat'));
+ MIN_FEED_CONSTRAINT .. sum(feed_crop_less_pasture, feedDM(feed_crop_less_pasture) * feedAmount(feed_crop_less_pasture)) =G= minFeedRate * (demand('meat') - importAmount('meat') + exportAmount('meat'));
IRRIGATION_CONSTRAINT(location) .. irrigConstraint(location) =G= sum(crop, irrigMaxRate(crop, location) * irrigI(crop, location) * area(crop, location)) / suitableLandArea(location);
@@ -168,90 +168,47 @@ $gdxin
+ (sum(location, agriLandExpansion(location)) * landChangeEnergy)
+ (sum(location, (0.2 * SQR((cropIncrease(location) + cropDecrease(location)) / suitableLandArea(location)) + (cropIncrease(location) + cropDecrease(location)) * 0.1) * (cropIncrease(location) + cropDecrease(location)) * landChangeEnergy))
+ (sum(location, (0.2 * SQR((pastureIncrease(location) + pastureDecrease(location)) / suitableLandArea(location)) + (pastureIncrease(location) + pastureDecrease(location)) * 0.1) * (pastureIncrease(location) + pastureDecrease(location)) * landChangeEnergy))
- + sum(import_crop, (netImportAmount(import_crop)) * worldInputPrices(import_crop))) / 1000000;
+ + sum(import_crop, (importAmount(import_crop)) * worldImportPrices(import_crop))
+ - sum(import_crop, (exportAmount(import_crop)) * worldExportPrices(import_crop))) / 1000000;
MODEL LAND_USE /ALL/ ;
-
- scalar basePasture;
- scalar baseCropland;
- scalar count / 1 /;
- scalar newPasture / 0 /;
- scalar newCropland / 0 /;
- scalar cropAdjSingle / 0 /;
- scalar pastureAdjSingle / 0 /;
-
- parameter temp(location);
-
-* if fitToPreviousAreas is zero then we just optimise
- if (fitToPreviousAreas eq 0,
- fertI.L(crop, location) = 0.5;
- irrigI.L(crop, location) = 0.5;
- otherIntensity.L(crop, location) = 0.1;
- netImportAmount.L(import_crop) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
- area.L(crop, location) = previousArea(crop, location)
- SOLVE LAND_USE USING NLP MINIMIZING energy;
-
-
+
+ fertI.L(crop, location) = 0.5;
+ irrigI.L(crop, location) = 0.5;
+ otherIntensity.L(crop, location) = 0.1;
+
+ importAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) gt 0) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
+ exportAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) lt 0) = -(maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
+
+ area.L(crop, location) = previousArea(crop, location)
+ SOLVE LAND_USE USING NLP MINIMIZING energy;
-* if fitToPreviousAreas is not zero then we try to find a cropAdj factors that give the previous areas ie we attempt to calibrate
- else
- basePasture = sum(location, previousArea("pasture", location));
- baseCropland = sum((location, crop_less_pasture), previousArea(crop_less_pasture, location));
-
- while((count le 10 and not (abs(basePasture-newPasture) le 0.01 and abs(baseCropland-newCropland) le 0.1)),
- fertI.L(crop, location) = 0.5;
- irrigI.L(crop, location) = 0.5;
- area.L(crop, location) = previousArea(crop, location);
- feedAmount.L(crop) = 0.0;
- netImportAmount.L(import_crop) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
- SOLVE LAND_USE USING NLP MINIMIZING energy;
-
- display pastureIncrease.l, pastureDecrease.l, cropDecrease.l, cropIncrease.l, unitEnergy.l;
-
- newPasture = sum(location, area.l("pasture", location));
- newCropland = sum((location, crop_less_pasture), area.l(crop_less_pasture, location));
-
- cropAdjSingle = newCropland / baseCropland;
- cropAdjSingle = min(max(cropAdjSingle, 0.5), 2);
- cropAdj(crop_less_pasture) = cropAdj(crop_less_pasture) * cropAdjSingle;
- pastureAdjSingle = newPasture / basePasture;
- pastureAdjSingle = min(max(pastureAdjSingle, 0.5), 2);
- cropAdj("pasture") = cropAdj("pasture") * pastureAdjSingle;
-
- cropAdj(crop) = min(max(cropAdj(crop), 0.2), 5);
-
- display area.l, count, basePasture, baseCropland, newPasture, newCropland, cropAdj;
- count = count+1 ;
- );
- );
-
- display agriLandExpansion.l, previousArea, area.l, net_supply.l, demand, feedAmount.l, yield.l, netImportAmount.l;
+* display agriLandExpansion.l, previousArea, area.l, net_supply.l, demand, feedAmount.l, yield.l, netImportAmount.l;
area.l(crop_less_pasture, location) = area.l(crop_less_pasture, location) / (1.0 - unhandledCropArea);
parameter totalProd(all_types);
- totalProd(crop) = sum(location, area.l(crop, location) * yield.l(crop, location));
- totalProd('meat') = meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount.l(feed_crop));
-
parameter totalProdCost(all_types);
- parameter totalImportEnergy(all_types);
+ parameter totalCropland(location);
+ parameter netImportAmount(all_types);
+ parameter netImportEnergy(all_types);
parameter feedEnergy(all_types);
+ totalProd(crop) = sum(location, area.l(crop, location) * yield.l(crop, location));
+ totalProd('meat') = meatEfficency*sum(feed_crop, feedDM(feed_crop) * feedAmount.l(feed_crop));
totalProdCost(crop) = sum(location, unitEnergy.l(crop, location) * area.l(crop, location));
- totalImportEnergy(import_crop) = netImportAmount.l(import_crop) * worldInputPrices(import_crop);
-
- feedEnergy(feed_crop)$(netImportAmount.l(feed_crop) gt 0) = (totalProdCost(feed_crop) + totalImportEnergy(feed_crop)) * feedAmount.l(feed_crop) / (totalProd(feed_crop) + netImportAmount.l(feed_crop));
- feedEnergy(feed_crop)$(netImportAmount.l(feed_crop) le 0 and totalProd(feed_crop) gt 0) = totalProdCost(feed_crop) * feedAmount.l(feed_crop) / totalProd(feed_crop);
-
- totalProdCost('meat') = sum(feed_crop, feedEnergy(feed_crop))
-
- parameter totalCropland(location);
totalCropland(location) = sum(crop_less_pasture, area.l(crop_less_pasture, location));
- display totalProdCost, totalProd, totalImportEnergy, feedEnergy, totalProdCost
+ netImportEnergy(import_crop) = importAmount.l(import_crop) * worldImportPrices(import_crop) - exportAmount.l(import_crop) * worldExportPrices(import_crop);
+ netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
+
+ feedEnergy(feed_crop)$(netImportAmount(feed_crop) gt 0) = (totalProdCost(feed_crop) + netImportEnergy(feed_crop)) * feedAmount.l(feed_crop) / (totalProd(feed_crop) + netImportAmount(feed_crop));
+ feedEnergy(feed_crop)$(netImportAmount(feed_crop) le 0 and totalProd(feed_crop) gt 0) = totalProdCost(feed_crop) * feedAmount.l(feed_crop) / totalProd(feed_crop);
+ totalProdCost('meat') = sum(feed_crop, feedEnergy(feed_crop));
+* display totalProdCost, totalProd, totalImportEnergy, feedEnergy, totalProdCost
Scalar ms 'model status', ss 'solve status';
ms=LAND_USE.modelstat;
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 8c95adb7c5e030f4dbe820eb4aeb1209bd7a96e1..0c71372cc1d76f166ef77751dee19af52ef3531b 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -164,14 +164,19 @@ public class GamsLocationOptimiser {
}
}
- if (DEBUG) LogWriter.println("\nWorld input prices");
- addItemMapParm(inDB.addParameter("worldInputPrices", 1), countryInput.getWorldInputPrices());
-
- if (DEBUG) LogWriter.println("\nMax Net Import");
- addItemMapParm(inDB.addParameter("maxNetImport", 1), countryInput.getMaxNetImport());
-
- if (DEBUG) LogWriter.println("\nMin Net Import");
- addItemMapParm(inDB.addParameter("minNetImport", 1), countryInput.getMinNetImport());
+ if (DEBUG) {
+ LogWriter.println("\nImport-export, min net imports, max net imports, import price, export price");
+ for (CropType crop : CropType.getImportedTypes()) {
+ LogWriter.println(String.format(" %15s, \t %5.1f, \t %5.1f, \t %5.3f, \t %5.3f",
+ crop.getGamsName(), countryInput.getMinNetImport().get(crop), countryInput.getMaxNetImport().get(crop),
+ countryInput.getWorldInputPrices().get(crop), countryInput.getWorldInputPrices().get(crop)));
+ }
+ }
+
+ addItemMapParm(inDB.addParameter("minNetImport", 1), countryInput.getMinNetImport(), false);
+ addItemMapParm(inDB.addParameter("maxNetImport", 1), countryInput.getMaxNetImport(), false);
+ addItemMapParm(inDB.addParameter("worldImportPrices", 1), countryInput.getWorldInputPrices(), false);
+ addItemMapParm(inDB.addParameter("worldExportPrices", 1), countryInput.getWorldInputPrices(), false);
addScalar(inDB.addParameter("meatEfficency", 0), countryInput.getMeatEfficiency());
if (DEBUG) LogWriter.println("\nMeatEfficiency: " + countryInput.getMeatEfficiency());
@@ -182,12 +187,8 @@ public class GamsLocationOptimiser {
addScalar(inDB.addParameter("minFeedRate", 0), countryInput.getMinFeedRate());
if (DEBUG) LogWriter.println("\nMinFeedRate: " + countryInput.getMinFeedRate());
- double fitToPreviousAreas = countryInput.isCalibrateToObserved() ? 1 : 0;
- if (DEBUG) LogWriter.println("\nfitToPreviousAreas: " + fitToPreviousAreas);
- addScalar(inDB.addParameter("fitToPreviousAreas", 0), fitToPreviousAreas);
-
if (DEBUG) LogWriter.println("\ncropAdj");
- addItemMapParm(inDB.addParameter("cropAdj", 1), countryInput.getCropAdjustments());
+ addItemMapParm(inDB.addParameter("cropAdj", 1), countryInput.getCropAdjustments(), DEBUG);
}
private void setPreviousArea(GAMSParameter prevCropP, String locString, String cropTypeString, double d) {
@@ -211,7 +212,7 @@ public class GamsLocationOptimiser {
GAMSVariable varIrrigIntensities = outDB.getVariable("irrigI");
GAMSVariable varOtherIntensities = outDB.getVariable("otherIntensity");
GAMSVariable varFeedAmount = outDB.getVariable("feedAmount");
- GAMSVariable varNetImports = outDB.getVariable("netImportAmount");
+ GAMSParameter parmNetImports = outDB.getParameter("netImportAmount");
GAMSVariable varYields = outDB.getVariable("yield");
GAMSVariable varUnitEnergies = outDB.getVariable("unitEnergy");
GAMSParameter parmCropAdj = outDB.getParameter("cropAdj");
@@ -247,7 +248,7 @@ public class GamsLocationOptimiser {
if (!cropUsageData.containsKey(cropType)) { // then we must not have seen this crop type before, so need to do all non location specific stuff
feedAmount = varFeedAmount.findRecord(itemName).getLevel();
- netImport = cropType.isImportedCrop() ? varNetImports.findRecord(itemName).getLevel() : 0;
+ netImport = cropType.isImportedCrop() ? getParmValue(parmNetImports, itemName) : 0;
cropAdj = getParmValue(parmCropAdj, itemName);
prod = getParmValue(parmProd, itemName);
prodCost = getParmValue(parmProdCost, itemName);
@@ -278,7 +279,7 @@ public class GamsLocationOptimiser {
landUseItem.setCropFraction(cropType, croplandArea > 0 ? area/croplandArea : 0);
}
- netImport = varNetImports.findRecord(CropType.MEAT.getGamsName()).getLevel();
+ netImport = getParmValue(parmNetImports, CropType.MEAT.getGamsName());
prod = getParmValue(parmProd, CropType.MEAT.getGamsName());
prodCost = getParmValue(parmProdCost, CropType.MEAT.getGamsName());
@@ -307,10 +308,10 @@ public class GamsLocationOptimiser {
param.addRecord().setValue(val);
}
- private void addItemMapParm(GAMSParameter parm, Map<CropType, Double> itemMap) {
+ private void addItemMapParm(GAMSParameter parm, Map<CropType, Double> itemMap, boolean printOutput) {
for (Map.Entry<CropType, Double> entry : itemMap.entrySet()) {
double d = entry.getValue();
- if (DEBUG) LogWriter.println(String.format(" %15s,\t %.3f", entry.getKey().getGamsName(), d));
+ if (printOutput) LogWriter.println(String.format(" %15s,\t %.3f", entry.getKey().getGamsName(), d));
if (!Double.isNaN(d))
parm.addRecord(entry.getKey().getGamsName()).setValue(d);
}