Handle export and import prices separately - i.e. losses and trade barriers

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;

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);
 		}