pinadd2ddir.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 PADD2DIR_ModelData : public ConoptModelData
{
public:
   int T = 0;
 
   /* 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;
 
   /* index arrays for the constraints */
   std::vector<int> consttdeq;
   std::vector<int> constseq;
   std::vector<int> constcseq;
   std::vector<int> constdeq;
   std::vector<int> constreq;
   std::vector<int> constdrev;
 
   std::map<int, std::pair<int, int>> consmapping;
 
   std::vector<int> hessianstart;

   void buildModel(int T, const std::vector<double> &xkeep, const std::vector<int> &xstat,
         const std::vector<int> &estat)
   {
      int varidx;
      this->T = T;
 
      /* 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);
      }

      if (xkeep.size() > 0)
      {
         for (size_t i = 0; i < xkeep.size(); i++)
         {
            getVariable(i).curr = xkeep[i];
         }
 
         /* linear extrapolation for the last period of variable td */
         getVariable(vartd[T - 1]).curr =
               2 * getVariable(vartd[T - 2]).curr - getVariable(vartd[T - 3]).curr;
         /* cs */
         getVariable(varcs[T - 1]).curr =
               2 * getVariable(varcs[T - 2]).curr - getVariable(varcs[T - 3]).curr;
         /* s */
         getVariable(vars[T - 1]).curr =
               2 * getVariable(vars[T - 2]).curr - getVariable(vars[T - 3]).curr;
         /* d */
         getVariable(vard[T - 1]).curr =
               2 * getVariable(vard[T - 2]).curr - getVariable(vard[T - 3]).curr;
         /* r */
         getVariable(varr[T - 1]).curr =
               2 * getVariable(varr[T - 2]).curr - getVariable(varr[T - 3]).curr;
         /* p */
         getVariable(varp[T - 1]).curr =
               2 * getVariable(varp[T - 2]).curr - getVariable(varp[T - 3]).curr;
         /* rev */
         getVariable(varrev[T - 1]).curr =
               2 * getVariable(varrev[T - 2]).curr - getVariable(varrev[T - 3]).curr;

         for (size_t i = 0; i < xstat.size(); i++)
         {
            getVariable(i).varstatus = xstat[i];
         }
 
         /* td */
         getVariable(vartd[T - 1]).varstatus = getVariable(vartd[T - 2]).varstatus;
         /* cs */
         getVariable(varcs[T - 1]).varstatus = getVariable(varcs[T - 2]).varstatus;
         /* s */
         getVariable(vars[T - 1]).varstatus = getVariable(vars[T - 2]).varstatus;
         /* d */
         getVariable(vard[T - 1]).varstatus = getVariable(vard[T - 2]).varstatus;
         /* r */
         getVariable(varr[T - 1]).varstatus = getVariable(varr[T - 2]).varstatus;
         /* p */
         getVariable(varp[T - 1]).varstatus = getVariable(varp[T - 2]).varstatus;
         /* rev */
         getVariable(varrev[T - 1]).varstatus = getVariable(varrev[T - 2]).varstatus;
      }
 
      /* 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)
         {
            considx = 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
         {
            considx = addConstraint(ConoptConstraintType::Eq, demand[t],
                  {vartd[t], vartd[t - 1], varp[t]}, {1.0, -0.87, 0.13}, {0, 0, 0});
         }
         consttdeq.push_back(considx);
 
         /* 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});
         }
         constseq.push_back(considx);
 
         /* 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)
         {
            considx = addConstraint(
                  ConoptConstraintType::Eq, 0.0, {varcs[t], vars[t]}, {1.0, -1.0}, {0, 0});
         }
         else
         {
            considx = addConstraint(ConoptConstraintType::Eq, 0.0,
                  {varcs[t], varcs[t - 1], vars[t]}, {1.0, -1.0, -1.0}, {0, 0, 0});
         }
         constcseq.push_back(considx);
 
         /* adding the deq equations */
         considx = addConstraint(ConoptConstraintType::Eq, 0.0, {vard[t], vartd[t], vars[t]},
               {1.0, -1.0, 1.0}, {0, 0, 0});
         constdeq.push_back(considx);
 
         /* 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)
         {
            considx = addConstraint(
                  ConoptConstraintType::Eq, 500.0, {varr[t], vard[t]}, {1.0, 1.0}, {0, 0});
         }
         else
         {
            considx = addConstraint(ConoptConstraintType::Eq, 0.0, {varr[t], varr[t - 1], vard[t]},
                  {1.0, -1.0, 1.0}, {0, 0, 0});
         }
         constreq.push_back(considx);
 
         /* 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});
         constdrev.push_back(considx);
 
         /* 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));
      }
 
      /* Restore values from previous runs */
      if (estat.size() > 0)
      {
         /* the equation status */
         for (size_t i = 0; i < estat.size(); i++)
         {
            getConstraint(i).slackstatus = estat[i];
         }
 
         /* tdeq */
         getConstraint(consttdeq[T - 1]).slackstatus = getConstraint(consttdeq[T - 2]).slackstatus;
         /* seq */
         getConstraint(constseq[T - 1]).slackstatus = getConstraint(constseq[T - 2]).slackstatus;
         /* cseq */
         getConstraint(constcseq[T - 1]).slackstatus = getConstraint(constcseq[T - 2]).slackstatus;
         /* deq */
         getConstraint(constdeq[T - 1]).slackstatus = getConstraint(constdeq[T - 2]).slackstatus;
         /* req */
         getConstraint(constreq[T - 1]).slackstatus = getConstraint(constreq[T - 2]).slackstatus;
         /* drev */
         getConstraint(constdrev[T - 1]).slackstatus = getConstraint(constdrev[T - 2]).slackstatus;
      }
 
      /* setting the objective constraint */
      setObjectiveElement(ConoptObjectiveElement::Constraint, objidx);
 
      /* setting the optimisation direction */
      setOptimizationSense(ConoptSense::Maximize);
 
      /*
        There are two nonlinear constraints, Sdef and RevDef and they have
        the following Hessians (for each time period):
 
        Revdef equation. Nonlinear term = -d*(p-250/r)
        Hessian (diagonal and lower triangle): A total of 3 nonzero elements per period
                      3           4           5
                      d           r           p
        3 : d
        4 : r     -250/r**2   500*d/r**3
        5 : p        -1
 
        Sdef: Nonlinear term = -(1.1+0.1p)*1.02**(-cs/7)
        Hessian (diagonal and lower triangle): A total of 2 nonzero elements per period
                                 1                                  5
                                cs                                  p
        1 : cs  -(1.1+0.1p)*1.02**(-cs/7)*(log(1.02)/7)**2
        5 : p           0.1*1.02**(-cs/7)*(log(1.02)72)
 
        The two Hessian matrices do not overlap so the total number of nonzeros is
        5 per period. The structure using numerical indices for rows and columns and
        running indices for the Hessian element is:
 
             1   3   4   5
        1    0
        3
        4        2   4
        5    1   3
 
        and the same with names of variables and contributing equations
             cs      d       r     p
        cs  sdef
        d
        r         revdef  revdef
        p   sdef  revdef
       */
      /* setting the structure of the hessian.
       * NOTE: the hessian must be in increasing order of columns, and then
       * rows. Care must be taken to ensure that the variables are referenced
       * in the correct order.
       */
      std::vector<int> hsrow;
      std::vector<int> hscol;
 
      for (int t = 0; t < T; ++t)
      {
         /* storing the start index for the hessian in each time period.
          * This will provide a mapping that can be used in the hessian
          * evaluation
          */
         hessianstart.push_back(hscol.size());
 
         hscol.push_back(varcs[t]);
         hsrow.push_back(varcs[t]);
         hscol.push_back(varcs[t]);
         hsrow.push_back(varp[t]);
         hscol.push_back(vard[t]);
         hsrow.push_back(varr[t]);
         hscol.push_back(vard[t]);
         hsrow.push_back(varp[t]);
         hscol.push_back(varr[t]);
         hsrow.push_back(varr[t]);
      }
      setSDLagrangianStructure(hsrow, hscol);
   }

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

   int SDDir(const double x[], const double dx[], double d2g[], int rowno, const int jacnum[],
         int *nodrv, 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;
 
      if (consset == CONS_SEQ)
      {
         /*  seq equation. Nonlinear term = -(1.1+0.1*p)*1.02**(-cs/7)        */
         /*  Hessian (diagonal and lower triangle): A total of 2 nonzero elements per period
                                 1                                  5
                                 cs                                  p
         1 : cs  -(1.1+0.1p)*1.02**(-cs/7)*(log(1.02)/7)**2
         5 : p           0.1*1.02**(-cs/7)*(log(1.02)72)                     */
 
         double cs = x[varcs[t]];
         double p = x[varp[t]];
         double h1 = 1.1 + 0.1 * p;
         double h2 = pow(1.02, -cs / 7.0);
         double h11 = -h1 * h2 * pow(log(1.02) / 7.0, 2);
         double h15 = 0.1 * h2 * (log(1.02) / 7.);
         d2g[varcs[t]] = h11 * dx[varcs[t]] + h15 * dx[varp[t]];
         d2g[varp[t]] = h15 * dx[varcs[t]];
      }
      else if (consset == CONS_DREV)
      {
         /*  drev equation. Nonlinear term = -d*(p-250/r)                    */
         /*  Hessian (diagonal and lower triangle): A total of 3 nonzero elements per period
                         3           4           5
                         d           r           p
         3 : d
         4 : r     -250/r**2   500*d/r**3
         5 : p        -1                                                     */
         double d = x[vard[t]];
         double r = x[varr[t]];
         double h34 = -250.0 / pow(r, 2);
         double h35 = -1.0;
         double h44 = 500.0 * d / pow(r, 3);
         d2g[vard[t]] = h34 * dx[varr[t]] + h35 * dx[varp[t]];
         d2g[varr[t]] = h34 * dx[vard[t]] + h44 * dx[varr[t]];
         d2g[varp[t]] = h35 * dx[vard[t]];
      }
      else
      {
         /*  Error - this equation is not nonlinear                            */
         printf("\nError. Pin_2DDir called with rowno = %d.\n\n", rowno);
         return 1;
      }
 
      return 0;
   }

   int SDLagrVal(const double x[], const double u[], const int hsrw[], const int hscl[],
         double hsvl[], int *nodrv, int numvar, int numcon, int nhess) override
   {
      double cs, d, r, p;
      double h1, h2;
 
      /*  Normal Evaluation mode                                            */
      for (int t = 0; t < T; t++)
      {
         int hessidx = hessianstart[t];
         cs = x[varcs[t]];
         d = x[vard[t]];
         r = x[varr[t]];
         p = x[varp[t]];
 
         h1 = 1.1 + 0.1 * p;
         h2 = pow(1.02, -cs / 7.0);
         /* the second derivative of the SEQ equations */
         hsvl[hessidx] = -h1 * h2 * pow(log(1.02) / 7.0, 2) * u[constseq[t]];
         hsvl[hessidx + 1] = 0.1 * h2 * (log(1.02) / 7.) * u[constseq[t]];
 
         /* the second derivative of the DREV equations */
         hsvl[hessidx + 2] = -250.0 / pow(r, 2) * u[constdrev[t]];
         hsvl[hessidx + 3] = -1.0 * u[constdrev[t]];
         hsvl[hessidx + 4] = 500.0 * d / pow(r, 3) * u[constdrev[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);
   Tut_MessageHandler msghandler(pname);
 
   // adding the message handler to the conopt interface
   conopt.setMessageHandler(msghandler);
 
   for (int i = 16; i <= 20; i++)
   {
 
      PADD2DIR_ModelData modeldata;
 
      auto xkeep = conopt.getVariableValues();
      auto xstat = conopt.getVariableStatus();
      auto estat = conopt.getConstraintStatus();
 
      // building the model
      modeldata.buildModel(i, xkeep, xstat, estat); // Add T, and other vectors as parameters
 
      // 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);
 
      conopt.printStatus();
   }
   cpp_log(conopt, "Successful Solve", COI_Error);
}