pindyck.cpp

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <iostream>
#include <map>
#include <utility>
#include "conopt.hpp"
#include <cassert>
 
#define CONS_SEQ 0
#define CONS_DREV 1
 
class PD_ModelData : public ConoptModelData
{
public:
   int T = 16;  /* Number of time periods */
   int Vpp = 7; /* Variables per period */
   int Epp = 6; /* Equations per period */
 
   /* the demand vector */
   std::vector<double> demand;
 
   /* index arrays for the variables */
   std::vector<int> vartd;
   std::vector<int> varcs;
   std::vector<int> vars;
   std::vector<int> vard;
   std::vector<int> varr;
   std::vector<int> varp;
   std::vector<int> varrev;
 
   std::map<int, std::pair<int, int>> consmapping;

   void buildModel()
   {
      int varidx;
 
      /* defining the demand for each time period. */
      for (int t = 0; t < T; t++)
         demand.push_back(1.0 + 2.3 * pow(1.015, t));

 
      for (int t = 0; t < T; t++)
      {
         /* td:  Lower=0, Upper=inf, Curr=18 */
         varidx = addVariable(0, Conopt::Infinity, 18);
         vartd.push_back(varidx);
 
         /* cs:  Lower=0, Upper=inf, Curr=7*t */
         varidx = addVariable(0, Conopt::Infinity, 7 * (t + 1));
         varcs.push_back(varidx);
 
         /* s:   Lower=0, Upper=inf, Curr=7 */
         varidx = addVariable(0, Conopt::Infinity, 7);
         vars.push_back(varidx);
 
         /* d:   Lower=0, Upper=inf, Curr=td-s */
         varidx =
               addVariable(0, Conopt::Infinity, getVariable(vartd[t]).curr - getVariable(vars[t]).curr);
         vard.push_back(varidx);
 
         /* r:   Lower=1, Upper=inf, Curr=r(t-1)-d */
         if (t > 0)
            varidx = addVariable(
                  1, Conopt::Infinity, getVariable(varr[t - 1]).curr - getVariable(vard[t]).curr);
         else
            varidx = addVariable(1, Conopt::Infinity, 500 - getVariable(vard[t]).curr);
         varr.push_back(varidx);
 
         /* p:   Lower=1, Upper=inf, Curr=14 */
         varidx = addVariable(1, Conopt::Infinity, 14);
         varp.push_back(varidx);
 
         /* rev: Lower=-inf, Upper=inf, Curr=0 */
         varidx = addVariable(-Conopt::Infinity, Conopt::Infinity);
         varrev.push_back(varidx);
      }
 
      /* adding the objective */
      std::vector<double> objcoeff;
      std::vector<int> objnlflag(T, 0);
      for (int t = 0; t < T; ++t)
      {
         objcoeff.push_back(pow(1.05, 1 - (t + 1)));
      }
      int objidx = addConstraint(ConoptConstraintType::Free, 0.0, varrev, objcoeff, objnlflag);
 
      /* adding the constraint information */
      int considx;
      for (int t = 0; t < T; t++)
      {
         /* adding the tdeq equations
          * for the initial time period, the fixed value for td_{t - 1} is applied to the RHS
          *
          * tdeq(t): td(t) = 0.87*td(t-1) - 0.13*p(t) + demand(t)
          *
          * In GAMS, td("1974") is fixed (td.fx = 18), and the first equation is for 1975.
          * So here for t=0 (1974), we apply the fixed value 18.0 directly on the RHS.
          *
          * That is why below we check if t == 0 and apply the fixed value of 18.0 directly
          * to the RHS. Otherwise, we use the value of td(t-1) from the previous period.
          *
          * The same applies to the remaining equations.
          */
         if (t == 0)
         {
            addConstraint(ConoptConstraintType::Eq, demand[t] + 0.87 * 18.0,
                  {vartd[t], varp[t]}, // the variables
                  {1.0, 0.13},         // the coefficients
                  {0, 0}               // nlflags
            );
         }
         else
         {
            addConstraint(ConoptConstraintType::Eq, demand[t], {vartd[t], vartd[t - 1], varp[t]},
                  {1.0, -0.87, 0.13}, {0, 0, 0});
         }
 
         /* adding the seq Equations
          * for the initial time period, the fixed value of s(1974)=6.5 is applied to the RHS
          * GAMS: s(t) = 0.75*s(t-1) + (1.1+0.1*p(t))*1.02^(-cs(t)/7)
          */
         if (t == 0)
         {
            considx = addConstraint(ConoptConstraintType::Eq, 0.75 * 6.5,
                  {vars[t], varp[t], varcs[t]}, {1.0, 0.0, 0.0}, {0, 1, 1});
         }
         else
         {
            considx = addConstraint(ConoptConstraintType::Eq, 0.0,
                  {vars[t], vars[t - 1], varp[t], varcs[t]}, {1.0, -0.75, 0.0, 0.0}, {0, 0, 1, 1});
         }
 
         /* storing the constraint index for the seq constraint, so that we can identify then
          * constraint in FDEval
          */
         consmapping.emplace(considx, std::make_pair(CONS_SEQ, t));
 
         /* adding the cseq Equations
          * for the initial time period, cs(1974)=0.0
          * cs(t) = cs(t-1) + s(t)
          * So for t=0, we use cs(0) = s(0) and drop cs(t-1)
          */
         if (t == 0)
         {
            addConstraint(ConoptConstraintType::Eq, 0.0, {varcs[t], vars[t]}, {1.0, -1.0}, {0, 0});
         }
         else
         {
            addConstraint(ConoptConstraintType::Eq, 0.0, {varcs[t], varcs[t - 1], vars[t]},
                  {1.0, -1.0, -1.0}, {0, 0, 0});
         }
 
         /* adding the deq equations */
         addConstraint(ConoptConstraintType::Eq, 0.0, {vard[t], vartd[t], vars[t]},
               {1.0, -1.0, 1.0}, {0, 0, 0});
 
         /* adding the req equations
          * for the initial time period, the value of r(1974)=500 is applied to the RHS
          * r(t) = r(t-1) - d(t)
          * So for t=0, we use r(0) + d(0) = 500
          */
         if (t == 0)
         {
            addConstraint(ConoptConstraintType::Eq, 500.0, {varr[t], vard[t]}, {1.0, 1.0}, {0, 0});
         }
         else
         {
            addConstraint(ConoptConstraintType::Eq, 0.0, {varr[t], varr[t - 1], vard[t]},
                  {1.0, -1.0, 1.0}, {0, 0, 0});
         }
 
         /* adding the drev equations */
         considx = addConstraint(ConoptConstraintType::Eq, 0.0,
               {varrev[t], vard[t], varp[t], varr[t]}, {1.0, 0.0, 0.0, 0.0}, {0, 1, 1, 1});
 
         /* storing the constraint index for the drev constraint, so that we can identify then
          * constraint in FDEval
          */
         consmapping.emplace(considx, std::make_pair(CONS_DREV, t));
      }
 
      /* setting the objective constraint */
      setObjectiveElement(ConoptObjectiveElement::Constraint, objidx);
 
      /* setting the optimisation direction */
      setOptimizationSense(ConoptSense::Maximize);
   }

   int FDEval(const double x[], double *g, double jac[], int rowno, const int jacnum[], int mode,
         int ignerr, int *errcnt, int numvar, int numjac, int thread) override
   {
      /* returning an iterator based on the row number.
       * NOTE: a map is just an array of pairs. So for the iterator, the
       * first element is the key and the second is the value. In our case,
       * then value is also a pair.
       */
      auto itr = consmapping.find(rowno);
      if (itr == consmapping.end())
         return 0;
 
      /* storing the value (a pair) as a separate variable to make the
       * extraction easier
       */
      auto conspair = itr->second;
      /* I think that this will be a .first and .second. It could be ->first
       * and ->second
       */
      int consset = conspair.first;
      int t = conspair.second;
      assert(consset == CONS_SEQ || consset == CONS_DREV);
 
      if (consset == CONS_SEQ)
      {
         double h1, h2;
         /*  seq equation. Nonlinear term = -(1.1+0.1*p(t))*1.02**(-cs(t)/7)        */
 
         h1 = (1.1 + 0.1 * x[varp[t]]);
         h2 = pow(1.02, -x[varcs[t]] / 7.0);
         if (1 == mode || 3 == mode)
            *g = -h1 * h2;
         if (2 == mode || 3 == mode)
         {
            jac[varcs[t]] = h1 * h2 * log(1.02) / 7.0;
            jac[varp[t]] = -h2 * 0.1;
         }
      }
      else if (consset == CONS_DREV)
      {
         /*  drev equation. Nonlinear term = -d(to)*(p(to)-250/r(to))        */
 
         if (1 == mode || 3 == mode)
            *g = -x[vard[t]] * (x[varp[t]] - 250. / x[varr[t]]);
         if (2 == mode || 3 == mode)
         {
            jac[vard[t]] = -(x[varp[t]] - 250. / x[varr[t]]);
            jac[varr[t]] = -x[vard[t]] * 250. / pow(x[varr[t]], 2);
            jac[varp[t]] = -x[vard[t]];
         }
      }
 
      return 0;
   }
};
 
#include "std.cpp"

int main(int argc, char **argv)
{
   int COI_Error = 0;
 
   // getting the program name from the executable path
   std::string pname = getProgramName(argv[0]);
 
   // initialising the Conopt Object
   Conopt conopt(pname);
   PD_ModelData modeldata;
   Tut_MessageHandler msghandler(pname);
 
   // adding the message handler to the conopt interface
   conopt.setMessageHandler(msghandler);
 
   // building the model
   modeldata.buildModel();
 
   // loading the model in the conopt object
   conopt.loadModel(modeldata);
 
#if defined(CONOPT_LICENSE_INT_1) && defined(CONOPT_LICENSE_INT_2) && defined(CONOPT_LICENSE_INT_3) && defined(CONOPT_LICENSE_TEXT)
   std::string license = CONOPT_LICENSE_TEXT;
   COI_Error += conopt.setLicense(CONOPT_LICENSE_INT_1, CONOPT_LICENSE_INT_2, CONOPT_LICENSE_INT_3, license);
#endif
 
   if (COI_Error)
   {
      conopt.sendMessage(
            "Skipping COI_Solve due to setup errors. COI_Error = " + std::to_string(COI_Error));
      cpp_log(conopt, "Skipping Solve due to setup errors", COI_Error);
   }
 
   COI_Error = conopt.solve(); /* Optimize                      */
 
   conopt.sendMessage("After solving. COI_Error = " + std::to_string(COI_Error));
   if (conopt.modelStatus() != 2 || conopt.solutionStatus() != 1)
      cpp_log(conopt, "Incorrect Model or Solver Status", -1);
   else if (fabs(conopt.objectiveValue() - 1170.486) > 0.001)
      cpp_log(conopt, "Incorrect objective returned", -1);
 
   conopt.printStatus();
 
   cpp_log(conopt, "Successful Solve", COI_Error);
}