Improved MAIDADS optimisation.

GAMS/elasticDemand.gms

@@ -20,7 +20,10 @@ $load i, alpha, beta, kappa, delta, tau, omega
 $gdxin
 
 VARIABLE
-  u                  Utility;
+  u                Utility
+  error            Error term for NLS
+  ressq            Residual squared
+;
 
 POSITIVE VARIABLE
   DiscretionaryC(i)  Discretionary consumption
@@ -49,35 +52,45 @@ EQUATIONS
   EQ_u         "Implicit definition of utility"
   EQ_DiscretionaryC(i)
   EQ_SubsistenceC(j)
+  EQ_u_NLS     "Relaxed constraint for NLS"
+  EQ_ressq_NLS
 ;
 
 EQ_w(j)..
-  w(j) =e= (SubsistenceC(j)+DiscretionaryC(j))*p(j)/m;
-
-EQ_u..
-  sum(i, (alpha(i)+beta(i)*exp(u))/(1+exp(u))*log(DiscretionaryC(i)))-u =e= kappa ;
-
-EQ_DiscretionaryC(i)$p(i).. DiscretionaryC(i) =e=
-    [(alpha(i)+beta(i)*exp(u))/(1+exp(u))]/p(i)*
-    [m-sum(j,p(j)*SubsistenceC(j))];
-
-EQ_SubsistenceC(j)..
-  SubsistenceC(j) =e= (delta(j)+tau(j)*exp(omega*u))/(1+exp(omega*u));
-
-model MAIDADS_Sim "MAIDADS model for simulation" /
-      EQ_w.w
-      EQ_u.u
-      EQ_DiscretionaryC.DiscretionaryC
-      EQ_SubsistenceC.SubsistenceC
- / ;
-
-MAIDADS_Sim.optfile = 1;
+  w(j) =e= (SubsistenceC(j) + DiscretionaryC(j)) * p(j) / m;
+
+EQ_u ..     sum(i, (alpha(i)+ beta(i)* exp(u)) / (1 + exp(u)) * log(DiscretionaryC(i))) - u =e= kappa;
+  
+EQ_u_NLS .. sum(i, (alpha(i)+ beta(i)* exp(u)) / (1 + exp(u)) * log(DiscretionaryC(i))) - u =e= kappa + error;
+
+EQ_DiscretionaryC(i)$p(i).. DiscretionaryC(i) =e= [(alpha(i) + beta(i) * exp(u))/(1 + exp(u))] / p(i) * [m - sum(j, p(j) * SubsistenceC(j))];
+
+EQ_SubsistenceC(j).. SubsistenceC(j) =e= (delta(j) + tau(j) * exp(omega * u)) / (1 + exp(omega * u));
+  
+EQ_ressq_NLS .. ressq =e= sqr(error);
+
+model MAIDADS_Sim_CNS "MAIDADS model for simulation" /
+      EQ_w
+      EQ_u
+      EQ_DiscretionaryC
+      EQ_SubsistenceC
+/ ;
+
+model MAIDADS_Sim_NLP /
+      EQ_w
+      EQ_u_NLS
+      EQ_DiscretionaryC
+      EQ_SubsistenceC
+      EQ_ressq_NLS
+/ ;
+
+MAIDADS_Sim_CNS.optfile = 1;
 option cns=path;
 
 desiredConsumpFactor =  m * maxIncomePropFoodSpend / sum(i, p(i)*delta(i));
 display desiredConsumpFactor;
 
-Scalar ms 'model status', ss 'solve status'; 
+Scalar ms 'model status', ss 'solve status';
 u = previousU;
 
 If (desiredConsumpFactor > 1.0,
@@ -87,14 +100,23 @@ If (desiredConsumpFactor > 1.0,
 	SubsistenceC.LO(j) = 0;
 	w(j) = (SubsistenceC(j)+DiscretionaryC(j))*p(j)/m;
 	desiredConsumpFactor = 1;
-	
-	solve MAIDADS_Sim using cns;
-	
-	ms=MAIDADS_Sim.modelstat;
-	ss=MAIDADS_Sim.solvestat;
+
+*   Initial solution by maximinsing u to get a better start for CNS. If no solution then use least squares
+    solve MAIDADS_Sim_CNS using nlp maximising u;
+    solve MAIDADS_Sim_CNS using cns;
+    
+    ms=MAIDADS_Sim_CNS.modelstat;
+    ss=MAIDADS_Sim_CNS.solvestat;
+ 
+    if ((MAIDADS_Sim_CNS.modelstat ne 16),
+        solve MAIDADS_Sim_NLP using nlp minimising ressq;
+        ms=MAIDADS_Sim_NLP.modelstat;
+        ss=MAIDADS_Sim_NLP.solvestat;
+    );
 
 else
 	SubsistenceC(i)  = delta(i) * desiredConsumpFactor;
 );
 
+
 display SubsistenceC.l, DiscretionaryC.l;

src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java

@@ -17,6 +17,7 @@ import com.gams.api.GAMSVariable;
 import com.gams.api.GAMSVariableRecord;
 import com.gams.api.GAMSWorkspace;
 import com.gams.api.GAMSWorkspaceInfo;
+import com.gams.api.GAMSGlobals.ModelStat;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.types.CommodityType;
@@ -134,15 +135,31 @@ public class GamsDemandOptimiser {
 	}
 		
 	private GamsDemandOutput handleResults(GAMSDatabase outDB) {
-		int modelStatus = (int) outDB.getParameter("ms").findRecord().getValue();
-		String status = GAMSGlobals.ModelStat.lookup(modelStatus).toString();
+		int modelStatusInt = (int) outDB.getParameter("ms").findRecord().getValue();
+		ModelStat modelStatus = GAMSGlobals.ModelStat.lookup(modelStatusInt);
+
 		LogWriter.println(String.format("\nDemand %s: Modelstatus (%d) %s, Solvestatus %s", inputData.getCountry(),
-				modelStatus, status, GAMSGlobals.SolveStat.lookup((int) outDB.getParameter("ss").findRecord().getValue()) ));
+				modelStatusInt, modelStatus.toString(), GAMSGlobals.SolveStat.lookup((int) outDB.getParameter("ss").findRecord().getValue()) ));
 		
 		GAMSVariable varU = outDB.getVariable("u");
 		double utility = varU.getFirstRecord().getLevel();
 		LogWriter.println(String.format("u = %.4f", utility)); 
 		
+		double leastSquareError = Math.abs(outDB.getVariable("error").getFirstRecord().getLevel());
+
+		switch (modelStatus) {
+			case SOLVED: break;
+			case UNDEFINED_STAT: break;
+			case OPTIMAL_LOCAL:
+				if (leastSquareError < 1e-5)
+					LogWriter.printlnWarning("No exact solution. Approximate solution found. Residual = " + leastSquareError);
+				else
+					LogWriter.printlnError("Critical! No exact solution. Approximate solution found with high error. Residual = " + leastSquareError);
+				break;
+			default:
+				LogWriter.printlnError("Critical! Solution failed.");			
+		}
+		
 		GAMSVariable varSubs = outDB.getVariable("SubsistenceC");
 		GAMSVariable varDisc = outDB.getVariable("DiscretionaryC");
 
@@ -157,7 +174,7 @@ public class GamsDemandOptimiser {
 				LogWriter.println(String.format("Can't find commodity from Gams demand, so ignoring: " + commodityName)); 
 			else {
 				double subsGams = rec.getLevel();
-				double discGams = modelStatus==0 ? 0 : varDisc.findRecord(commodityName).getLevel();  //modelStatus==0 is optimiser not run, i.e. too low GDP
+				double discGams = modelStatusInt==0 ? 0 : varDisc.findRecord(commodityName).getLevel();  //modelStatus==0 is optimiser not run, i.e. too low GDP
 				GamsCommodityDemand demand = new GamsCommodityDemand(commodity, subsGams, discGams,inputData.getKcalPerT(commodity));
 				LogWriter.println(demand.toString());
 				demandMap.add(demand);
@@ -174,6 +191,6 @@ public class GamsDemandOptimiser {
 		double desiredConsumpFactor = parConsumpFact.getFirstRecord().getValue();
 		LogWriter.println(String.format("desiredConsumpFactor=%.4f", desiredConsumpFactor));
 		
-		return new GamsDemandOutput(status, demandMap, utility, desiredConsumpFactor);
+		return new GamsDemandOutput(modelStatus.toString(), demandMap, utility, desiredConsumpFactor);
 	}
 }

src/ac/ed/lurg/country/gams/GamsDemandOutput.java

@@ -46,10 +46,4 @@ public class GamsDemandOutput {
 	public String getStatus() {
 		return status;
 	}
-
-	public void resetForNewAttempt(int i, double inputUtility) {
-		
-		if ((i & 1) == 0) utility = inputUtility * -1.0;
-		else utility = inputUtility + 0.8*i;
-	}
 }

src/ac/ed/lurg/demand/ElasticDemandManager.java

@@ -72,18 +72,8 @@ public class ElasticDemandManager extends AbstractSSPDemandManager {
 		GamsDemandInput inputData = new GamsDemandInput(c, year, gdpPc, consumerPrices,kcalPerT, usaGdpPc, previousGamsDemands.get(c));
 
 		// Do the projection
-		GamsDemandOutput gamsOutput = null;
-		for (int i = 1; i < 50; ++i) {
-			gamsOutput = new GamsDemandOptimiser(inputData).getDemandPc();
-			if (gamsOutput.getStatus().equals("SOLVED") || gamsOutput.getStatus().equals("UNDEFINED_STAT"))
-				break;
-			else {
-				LogWriter.printlnError(i + ": Problem solving " + c + ", " + year + " got " + gamsOutput.getStatus());
-				double inputUtility = inputData.getPreviousDemands() == null ? 0 : inputData.getPreviousDemands().getUtility();
-				gamsOutput.resetForNewAttempt(i, inputUtility);
-				inputData = new GamsDemandInput(c, year, gdpPc, consumerPrices, kcalPerT,usaGdpPc, gamsOutput);
-			}
-		}
+		GamsDemandOptimiser gamsDemandOptimiser = new GamsDemandOptimiser(inputData);
+		GamsDemandOutput gamsOutput = gamsDemandOptimiser.getDemandPc();
 		
 		previousGamsDemands.put(c, gamsOutput);
 		Map<CommodityType, Double> foodPc = gamsOutput.getPlumDemands();