pinopt.c

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "coiheader.h"
 
#include "comdecl.h"
 
 
 
int T    = 16;   /* Number of time periods */
int Vpp  =  7;   /* Variables per period */
int Epp  =  6;   /* Equations per period */
 

int COI_CALLCONV Pin_ReadMatrix( double LOWER[], double CURR[], double UPPER[], int VSTA[], int TYPE[],
                             double RHS[], int ESTA[], int COLSTA[], int ROWNO[], double VALUE[],
                             int NLFLAG[], int NUMVAR, int NUMCON, int NUMNZ, void* USRMEM )
{
   int it, is, i, icol, iz;
/*
  Define the information for the columns.
 
  We should not supply status information, VSTA.
 
  We order the variables as follows:
  0: td, 1: cs, 2: s, 3: d, 4: r, 5: p, and 6: rev. All variables for
  period 1 appears first followed by all variables for period 2 etc.
 
  td, cs, s, and d have lower bounds of 0, r and p have lower
  bounds of 1, and rev has no lower bound.
  All have infinite upper bounds (default).
  The initial value of td is 18, s is 7, cs is 7*t, d is td-s,
  p is 14, and r is r(t-1)-d. No initial value for rev.           */
 
   for ( it=1; it<=T; it++ ) {
      is = Vpp*(it-1);
      LOWER[is  ] = 0;
      LOWER[is+1] = 0;
      LOWER[is+2] = 0;
      LOWER[is+3] = 0;
      LOWER[is+4] = 1;
      LOWER[is+5] = 1;
      CURR[is  ]  = 18;
      CURR[is+1]  = 7*it;
      CURR[is+2]  = 7;
      CURR[is+3]  = CURR[is] - CURR[is+2];
      if ( it > 1 )
         CURR[is+4] = CURR[is+4-Vpp] - CURR[is+3];
      else
         CURR[is+4] = 500 - CURR[is+3];
      CURR[is+5]  = 14;
      }
/*
  Define the information for the rows
 
  We order the constraints as follows: The objective is first,
  followed by tddef, sdef, csdef, ddef, rdef, and revdef for
  the first period, the same for the second period, etc.
 
  The objective is a nonbinding constraint:                      */
 
   TYPE[0] = 3;
 
/*  All others are equalities:                                   */
 
   for ( i = 1; i < NUMCON; i++ )
      TYPE[i] = 0;
/*
  Right hand sides: In all periods except the first, only tddef
  has a nonzero right hand side of 1+2.3*1.015**(t-1).
  In the initial period there are contributions from lagged
  variables in the constraints that have lagged variables.       */
 
   for ( it = 1; it<=T; it++ ) {
      is      = 1 + 6*(it-1);
      RHS[is] = 1.0+2.3*pow(1.015,(it-1.));
      }
 
/* tddef: + 0.87*td(0)                                           */
 
   RHS[1] = RHS[1]+ 0.87*18.;
 
/* sdef: +0.75*s(0)                                              */
 
   RHS[2] =  0.75*6.5;
 
/* csdef: +1*cs(0)                                               */
 
   RHS[3] =  0;
 
/*  rdef: +1*r(0)                                                */
 
   RHS[5] = 500;
 
/*
  Define the structure and content of the Jacobian:
  To help define the Jacobian pattern and values it can be useful to
  make a picture of the Jacobian. We describe the variables for one
  period and the constraints they are part of:
 
                 td   cs   s    d    r    p    rev
                 0    1    2    3    4    5    6
 0: Obj                                      (1+r)**(1-t)
 Period t:
 0: tddef       1.0                      0.13
 1: sdef             NL   1.0            NL
 2: csdef            1.0 -1.0
 3: ddef       -1.0       1.0  1.0
 4: rdef                       1.0  1.0
 5: revdef                     NL   NL   NL   1.0
 Period t+1:
 0: tddef      -0.87
 1: sdef                 -0.75
 2: csdef           -1.0
 3: ddef
 4: rdef                           -1.0
 5: revdef
 
  The Jacobian has to be sorted column-wise so we will just define
  the elements column by column according to the table above:         */
 
   iz   = 0;  /* Running index for Jacobian in C convention 0 to nz-1 */
   icol = 0;  /* Running column index in C convention 0 to N-1        */
   for ( it = 1; it <= T; it++ )
      {
 
/* is points to the position before the first equation for the period */
 
      is = 1 + 6*(it-1);  /* 1 correspond to the objective in row 0   */
 
/* Column td:                                                         */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = is;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+3;
      VALUE[iz]  = -1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      if ( it < T ) {
         ROWNO[iz]  = is+Epp;
         VALUE[iz]  = -0.87;
         NLFLAG[iz] = 0;
         iz         = iz + 1;
         }
 
/*  Column cs                                                         */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = is+1;
      NLFLAG[iz] = 1;
      iz         = iz + 1;
      ROWNO[iz]  = is+2;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      if ( it < T ) {
         ROWNO[iz]  = is+2+Epp;
         VALUE[iz]  = -1.0;
         NLFLAG[iz] = 0;
         iz         = iz + 1;
      }
 
/*  Column s                                                          */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = is+1;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+2;
      VALUE[iz]  = -1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+3;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      if ( it < T ) {
         ROWNO[iz]  = is+1+Epp;
         VALUE[iz]  = -0.75;
         NLFLAG[iz] = 0;
         iz         = iz + 1;
      }
 
/*  Column d:                                                         */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = is+3;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+4;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+5;
      NLFLAG[iz] = 1;
      iz         = iz + 1;
 
/*  Column r:                                                         */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = is+4;
      VALUE[iz]  = +1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+5;
      NLFLAG[iz] = 1;
      iz         = iz + 1;
      if ( it < T ) {
         ROWNO[iz]  = is+4+Epp;
         VALUE[iz]  = -1.0;
         NLFLAG[iz] = 0;
         iz         = iz + 1;
      }
 
/*  Column p:                                                         */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = is;
      VALUE[iz]  = +0.13;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+1;
      NLFLAG[iz] = 1;
      iz         = iz + 1;
      ROWNO[iz]  = is+5;
      NLFLAG[iz] = 1;
      iz         = iz + 1;
 
/*  Column rev:                                                       */
      COLSTA[icol] = iz;
      icol       = icol + 1;
      ROWNO[iz]  = +0;                   /* Objective -- no is term   */
      VALUE[iz]  = pow(1.05,(1.0-it));
      NLFLAG[iz] = 0;
      iz         = iz + 1;
      ROWNO[iz]  = is+5;
      VALUE[iz]  = 1.0;
      NLFLAG[iz] = 0;
      iz         = iz + 1;
   }
   COLSTA[icol] = iz;  /* First elements past last column */
 
   return 0;
}
/*
 =====================================================================
 Compute nonlinear terms and non-constant Jacobian elements           */
 

int COI_CALLCONV Pin_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, void* USRMEM )
{
 
   int it, is;
   double  h1, h2;
 
/*  Compute the number of the period, it                              */
 
   it = (ROWNO+Epp-1) / Epp;
   is = Vpp*(it-1);
   if ( 1+(it-1)*Epp+1 == ROWNO ) {
 
/*  sdef equation. Nonlinear term = -(1.1+0.1*p)*1.02**(-cs/7)        */
 
      h1 = (1.1+0.1*X[is+5]);
      h2 = pow(1.02,-X[is+1]/7.0);
      if ( 1 == MODE || 3 == MODE )
         *G = -h1*h2;
      if ( 2 == MODE || 3 == MODE ) {
         JAC[is+1] = h1*h2*log(1.02)/7.0;
         JAC[is+5] = -h2*0.1;
         }
      }
   else if ( 1+(it-1)*Epp+5 == ROWNO ) {
 
/*  revdef equation. Nonlinear term = -d*(p-250/r)                    */
 
      if ( 1 == MODE || 3 == MODE )
         *G = -X[is+3]*(X[is+5]-250./X[is+4]);
      if ( 2 == MODE || 3 == MODE ) {
         JAC[is+3] = -(X[is+5]-250./X[is+4]);
         JAC[is+4] = -X[is+3]*250./pow(X[is+4],2);
         JAC[is+5] = -X[is+3];
         }
      }
   else {
 
/*  Error - this equation is not nonlinear                            */
 
      printf("\nError. Pin_FDEval called with ROWNO = %d. MODE = %d\n\n", ROWNO, MODE);
      return 1;
      };
   return 0;
}
 

int COI_CALLCONV Pin_Option( int NCALL, double* RVAL, int* IVAL, int* LVAL, char* NAME, void* USRMEM )
{
   fprintf(fd,"Inside Pin_Option. NCALL = %d\n",NCALL);
   if (NCALL == 0) {
      strcpy(NAME, "Doprep");
      *LVAL = 0;
        }
   else if (NCALL == 1) {
      strcpy(NAME, "Lotria");
      *IVAL = 1;
        }
   else {
      strcpy(NAME,"");
        };
   return 0;
}
#include "std.c"

int main(int argc, char** argv)
{
   int major, minor, patch;
 
   c_log( "Starting to execute", START );
 
/* Write which version of CONOPT we are using.                        */
 
   COIGET_Version( &major, &minor, &patch );
   printf("\nSolving Pindyck Model using CONOPT version %d.%d.%d:\n", major, minor, patch );
 
/*  Tell CONOPT about the sizes in the model                          */
/*  Number of variables (excl. slacks): 7 per period                  */
 
   COI_Error +=  COIDEF_NumVar    ( CntVect, 7*T  );
 
/*  Number of equations: 1 objective + 6 per period                   */
 
   COI_Error +=  COIDEF_NumCon    ( CntVect, 1+6*T  );
 
/*  Number of nonzeros in the Jacobian. See the counting in           */
/*  ReadMatrix above. For each period there is 1 in the objective, 16 */
/*  for unlagged variables and 4 for lagged variables.                */
 
   COI_Error +=  COIDEF_NumNz     ( CntVect, 17*T+4*(T-1) );
 
/*  Number of nonlinear nonzeros. 5 unlagged for each period.         */
 
   COI_Error +=  COIDEF_NumNlNz   ( CntVect, 5*T );
 
/*  Objective: Maximize Constraint no 0                               */
 
   COI_Error +=  COIDEF_OptDir    ( CntVect, 1  );
   COI_Error +=  COIDEF_ObjCon    ( CntVect, 0  );
 
/* Turn debugging of FDEval on (1) or off (0)                         */
 
   COI_Error +=  COIDEF_DebugFV   ( CntVect, 0  );
 
/* Register the options file                                          */
 
   COI_Error +=  COIDEF_Optfile   ( CntVect, "pindyck.opt");
 /*
   Register the necessary callback routines with CONOPT
 */
   COI_Error +=  COIDEF_Message   ( CntVect, &Std_Message  );  /* Register the callback Message     */
   COI_Error +=  COIDEF_ErrMsg    ( CntVect, &Std_ErrMsg   );  /* Register the callback ErrMsg      */
   COI_Error +=  COIDEF_Status    ( CntVect, &Std_Status   );  /* Register the callback Status      */
   COI_Error +=  COIDEF_Solution  ( CntVect, &Std_Solution );  /* Register the callback Solution    */
   COI_Error +=  COIDEF_ReadMatrix( CntVect, &Pin_ReadMatrix); /* Register the callback ReadMatrix  */
   COI_Error +=  COIDEF_FDEval    ( CntVect, &Pin_FDEval);     /* Register the callback FDEval      */
   COI_Error +=  COIDEF_Option    ( CntVect, &Pin_Option);     /* Register the callback Option      */
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
   COI_Error +=  COIDEF_License   ( CntVect, LICENSE_INT_1, LICENSE_INT_2, LICENSE_INT_3, LICENSE_TEXT);
#endif
 
   if ( COI_Error ) {
      printf("Skipping COI_Solve due to setup errors. COI_Error = %d\n",COI_Error);
      c_log( "Skipping Solve due to setup errors", COI_Error);
      }
   COI_Error =  COI_Solve        ( CntVect );           /* Optimize                      */
   printf("After solving. COI_Error =%d\n",COI_Error);
   if ( COI_Error ) {
      c_log( "Errors encountered during solution", COI_Error); }
   else if ( stacalls == 0 || solcalls == 0 ) {
      c_log( "Status or Solution routine was not called", -1); }
   else if ( mstat != 2 || sstat != 1 ) {
      c_log( "Incorrect Model or Solver Status", -1); }
   else if ( fabs( OBJ-1170.486 ) > 0.001 ) {
      c_log( "Incorrect objective returned", -1); }
 
   c_log( "Successful Solve", OK );
}