From d4e90b92410c0ba09a1160e0162a22618ab54c48 Mon Sep 17 00:00:00 2001
From: Bart Arendarczyk <s1924442@ed.ac.uk>
Date: Fri, 8 Jul 2022 13:10:31 +0100
Subject: [PATCH] Merged demand optimisation improvements.
---
GAMS/elasticDemand.gms | 76 ++++++++++++-------
.../country/gams/GamsDemandOptimiser.java | 27 +++++--
.../lurg/country/gams/GamsDemandOutput.java | 6 --
.../ed/lurg/demand/ElasticDemandManager.java | 14 +---
4 files changed, 73 insertions(+), 50 deletions(-)
diff --git a/GAMS/elasticDemand.gms b/GAMS/elasticDemand.gms
index 9d62dd08..b5b40b36 100755
--- a/GAMS/elasticDemand.gms
+++ b/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;
diff --git a/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java b/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
index 287a6054..80d55c7d 100755
--- a/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsDemandOptimiser.java
@@ -21,6 +21,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;
@@ -157,15 +158,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");
@@ -180,7 +197,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);
@@ -197,7 +214,7 @@ 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);
}
diff --git a/src/ac/ed/lurg/country/gams/GamsDemandOutput.java b/src/ac/ed/lurg/country/gams/GamsDemandOutput.java
index 06a2c353..5e793e53 100644
--- a/src/ac/ed/lurg/country/gams/GamsDemandOutput.java
+++ b/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;
- }
}
diff --git a/src/ac/ed/lurg/demand/ElasticDemandManager.java b/src/ac/ed/lurg/demand/ElasticDemandManager.java
index e59eb112..eb88edfe 100755
--- a/src/ac/ed/lurg/demand/ElasticDemandManager.java
+++ b/src/ac/ed/lurg/demand/ElasticDemandManager.java
@@ -59,18 +59,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();
--
GitLab