tutorial2.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 COI_CALLCONV Tut_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 )
{
 /*                                                               */
 /*   Information about Variables:                                */
 /*   Default: Lower = -Inf, Curr = 0, and Upper = +inf.          */
 /*   Default: the status information in Vsta is not used.        */
 /*                                                               */
 /*   Lower bound on L   = X[0] = 0.1 and initial value = 0.5:    */
 /*                                                               */
   LOWER[0] = 0.1 ;
   CURR[0]  = 0.5 ;
 /*                                                               */
 /*   Lower bound on INP = X[1] = 0.1 and initial value = 0.5:    */
 /*                                                               */
   LOWER[1] = 0.1 ;
   CURR[1]  = 0.5 ;
 /*                                                               */
 /*   Lower bound on OUT = X[2] and P = X[3] are both 0 and the   */
 /*   default initial value of 0 is used:                         */
 /*                                                               */
   LOWER[2] = 0.  ;
   LOWER[3] = 0.  ;
 /*                                                               */
 /*   Information about Constraints:                              */
 /*   Default: Rhs = 0                                            */
 /*   Default: the status information in Esta and the function    */
 /*            value in FV are not used.                          */
 /*   Default: Type: There is no default.                         */
 /*            0 = Equality,                                      */
 /*            1 = Greater than or equal,                         */
 /*            2 = Less than or equal,                            */
 /*            3 = Non binding.                                   */
 /*                                                               */
 /*   Constraint 0 (Objective)                                    */
 /*   Rhs = -0.1 and type Non binding                             */
 /*                                                               */
   RHS[0]  = -0.1 ;
   TYPE[0] = 3    ;
 /*                                                               */
 /*   Constraint 1 (Production Function)                          */
 /*   Rhs = 0 and type Equality                                   */
 /*                                                               */
   TYPE[1] = 0    ;
 /*                                                               */
 /*   Constraint 2 (Price equation)                               */
 /*   Rhs = 4.0 and type Equality                                 */
 /*                                                               */
   RHS[2]  = 4.   ;
   TYPE[2] = 0    ;
 /*                                                               */
 /*   Information about the Jacobian. We use the standard method with */
 /*   Rowno, Value, Nlflag and Colsta and we do not use Colno.    */
 /*                                                               */
 /*   Colsta = Start of column indices (No Defaults):             */
 /*   Rowno  = Row indices                                        */
 /*   Value  = Value of derivative (by default only linear        */
 /*                                 derivatives are used)         */
 /*   Nlflag = 0 for linear and 1 for nonlinear derivative        */
 /*            (not needed for completely linear models)          */
 /*                                                               */
 /*   Indices                                                     */
 /*        X[0]  X[1]  X[2]  X[3]                                 */
 /*   0:    0     2     4     7                                   */
 /*   1:    1     3     5                                         */
 /*   2:                6     8                                   */
 /*                                                               */
   COLSTA[0] =  0 ;
   COLSTA[1] =  2 ;
   COLSTA[2] =  4 ;
   COLSTA[3] =  7 ;
   COLSTA[4] =  9 ;
   ROWNO[0]  =  0 ;
   ROWNO[1]  =  1 ;
   ROWNO[2]  =  0 ;
   ROWNO[3]  =  1 ;
   ROWNO[4]  =  0 ;
   ROWNO[5]  =  1 ;
   ROWNO[6]  =  2 ;
   ROWNO[7]  =  0 ;
   ROWNO[8]  =  2 ;
 /*                                                               */
 /*   Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear */
 /*        X[0]  X[1]  X[2]  X[3]                                 */
 /*   0:    L     L     NL    NL                                  */
 /*   1:    NL    NL    L                                         */
 /*   2:                L     L                                   */
 /*                                                               */
   NLFLAG[0] =  0 ;
   NLFLAG[1] =  1 ;
   NLFLAG[2] =  0 ;
   NLFLAG[3] =  1 ;
   NLFLAG[4] =  1 ;
   NLFLAG[5] =  0 ;
   NLFLAG[6] =  0 ;
   NLFLAG[7] =  1 ;
   NLFLAG[8] =  0 ;
 /*                                                               */
 /*   Value (Linear only)                                         */
 /*        X[0]  X[1]  X[2]  X[3]                                 */
 /*   0:   -1    -1     NL    NL                                  */
 /*   1:    NL    NL    -1                                        */
 /*   2:                 1     2                                  */
 /*                                                               */
   VALUE[0]  = -1. ;
   VALUE[2]  = -1. ;
   VALUE[5]  = -1. ;
   VALUE[6]  =  1. ;
   VALUE[8]  =  2. ;
 
   return 0;
}
 

int COI_CALLCONV Tut_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 )
{
 /*                                                                   */
 /*   Declare local copies of the optimization variables. This is     */
 /*   just for convenience to make the expressions easier to read.    */
 /*                                                                   */
   double L, Inp, Out, P;
 /*                                                                   */
 /*   Declare parameters and their data values.                       */
 /*                                                                   */
   double Al   = 0.16;
   double Ak   = 2.0;
   double Ainp = 0.16;
   double Rho  = 1.0;
   double K    = 4.0;
   double hold1, hold2, hold3;
 
 /*                                                                   */
 /*   Move the optimization variables from the X vector to a set      */
 /*   of local variables with the same names as the variables in      */
 /*   the model description. This is not necessary, but it should make*/
 /*   the equations easier to recognize.                              */
 /*   This time we work with the C numbering convention               */
 /*                                                                   */
   L   = X[0];
   Inp = X[1];
   Out = X[2];
   P   = X[3];
 /*                                                                   */
 /*   Row 0: the objective function is nonlinear                      */
 /*                                                                   */
 
   if ( ROWNO == 0 ) {
 /*                                                                   */
 /*   Mode = 1 or 3. Function value: G = P * Out                      */
 /*                                                                   */
      if ( MODE == 1 || MODE == 3 )
         *G    = P * Out;
 /*                                                                   */
 /*   Mode = 2 or 3: Derivative values:                               */
 /*                                                                   */
      if ( MODE == 2 || MODE == 3 ) {
         JAC[2] = P;    /* derivative w.r.t. Out = X[2]               */
         JAC[3] = Out;  /* derivative w.r.t. P   = X[3]               */
         }
      }
 /*                                                                   */
 /*   Row 1: The production function is nonlinear                     */
 /*                                                                   */
   else if ( ROWNO == 1 ) {
 /*                                                                   */
 /*   Compute some common terms                                       */
 /*                                                                   */
      hold1 = (Al*pow(L,(-Rho)) + Ak*pow(K,(-Rho)) + Ainp*pow(Inp,(-Rho)));
      hold2 = pow(hold1,( -1./Rho ));
 /*                                                                   */
 /*   Mode = 1 or 3: Function value                                   */
 /*                                                                   */
      if ( MODE == 1 || MODE == 3 )
         *G = hold2;
 /*                                                                   */
 /*   Mode = 2 or 3: Derivatives                                      */
 /*                                                                   */
      if ( MODE == 2 || MODE == 3 ) {
         hold3  = hold2 / hold1;
         JAC[0] = hold3 * Al   * pow(L  ,(-Rho-1.)); /* derivative w.r.t. L   = X[0]  */
         JAC[1] = hold3 * Ainp * pow(Inp,(-Rho-1.)); /* derivative w.r.t. Inp = X[1]  */
         }
      }
 /*                                                                   */
 /*   Row = 2: The row is linear and will not be called.              */
 /*                                                                   */
 
   return 0;
}

int COI_CALLCONV Tut_2DLagrStr( int HSRW[], int HSCL[], int* NODRV,
                                int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
/*  There are two nonlinear constraints, P * Out in row 0 and
  Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho ) in
  row 1 and they have the following Hessians:
 
  Row 0: P * Out
  Hessian (diagonal and lower triangle): A total of 1 nonzero element
      Out    P
       2     3
  2
  3    1
 
  Row 1: Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho )
         = hold1**(-1.d0/Rh0) = hold2
  with L and Inp being the variables.
  Hessian (diagonal and lower triangle): A total of 3 nonzero element
       L    Inp
       0     1
   0   d1
   1   d2   d3
 
  where d1, d2, and d3 are derived below.
 
  The two Hessian matrices do not overlap so the total number of nonzeros
  is 4. The structure using numerical indices for rows and columns and
  running indices for the Hessian element is:
 
       0   1   2   3
  0    0
  1    1   2
  2
  3            3
 
  and the same with names of variables and numbers of contributing
  equations:
 
       L  Inp Out  P
  L    1
  Inp  1   1
  Out
  P            0
                                                                      */
 
/*  Define structure of Hessian                                       */
 
   HSRW[0] = 0; HSCL[0] = 0;
   HSRW[1] = 1; HSCL[1] = 0;
   HSRW[2] = 1; HSCL[2] = 1;
   HSRW[3] = 3; HSCL[3] = 2;
 
   return 0;
}

int COI_CALLCONV Tut_2DLagrVal( const double X[], const double U[], const int HSRW[], const int HSCL[], double HSVL[],
                                int* NODRV, int NUMVAR, int NUMCON, int NHESS, void* USRMEM )
{
/*
  Declare local copies of the optimization variables. This is
  just for convenience to make the expressions easier to read.        */
 
   double L, Inp, Out, P;
 /*                                                                   */
 /*   Declare parameters and their data values.                       */
 /*                                                                   */
   double Al   = 0.16;
   double Ak   = 2.0;
   double Ainp = 0.16;
   double Rho  = 1.0;
   double K    = 4.0;
   double hold1, hold2, hold3, hold4;
   int    i;
 
/*  Normal Evaluation mode                                            */
 
   HSVL[3] = U[0]; /* Second derivative of constraint 0 is 1       */
   if ( U[1] != 0.0 ) {
      L   = X[0];
      Inp = X[1];
      Out = X[2];
      P   = X[3];
      hold1 = (Al*pow(L,-Rho) + Ak*pow(K,-Rho) + Ainp*pow(Inp,-Rho));
      hold2 = pow( hold1,-1.0/Rho);
      hold3 = hold2 / hold1;
/*
  The code for the first derivative was:
        jac(1) = hold3 * Al   * L   ** (-Rho-1.0)
        jac(2) = hold3 * Ainp * Inp ** (-Rho-1.0)
  Define Hold4 as the derivative of Hold3 with respect to Hold1       */
 
     hold4 = hold3 / hold1 * (-1.0/Rho-1.0);
/*
  The second derivatives are computed as:
  (1) The derivative of Hold3 multiplied by the other terms PLUS
  (2) Hold3 multiplied by the derivative of the other terms
  where the derivative of Hold3 is Hold4 * the derivative of Hold1
  with respect to the variable in question.                           */
 
      HSVL[0] = hold4 * (-Rho) * pow(Al*pow(L,-Rho-1.0),2) +
                hold3 * Al   * (-Rho-1.0)*pow(L,-Rho-2.0);
      HSVL[1] = hold4 * (-Rho) * (Al*pow(L,-Rho-1.0)) *
                                 (Ainp*pow(Inp,-Rho-1.0));
      HSVL[2] = hold4 * (-Rho) * pow(Ainp*pow(Inp,-Rho-1.0),2) +
                hold3 * Ainp * (-Rho-1.0)*pow(Inp,-Rho-2.0);
 
/*  Scale the derivatives with the Lagrange multiplier, U[1]:         */
 
      for ( i = 0; i < 3; i++ )
         HSVL[i] = HSVL[i] * U[1];
   }
 
   (void)P;
   (void)Out;
 
   return 0;
}
 
#include "std.c"

int main(int argc, char** argv)
{
   c_log( "Starting to execute", START );
 
 /*
   Tell CONOPT about the sizes in the model
 */
   COI_Error +=  COIDEF_NumVar    ( CntVect, 4  );       /* 4 variables                       */
   COI_Error +=  COIDEF_NumCon    ( CntVect, 3  );       /* 3 constraints                     */
   COI_Error +=  COIDEF_NumNz     ( CntVect, 9  );       /* 9 nonzeros in the Jacobian        */
   COI_Error +=  COIDEF_NumNlNz   ( CntVect, 4  );       /* 4 of which are nonlinear          */
   COI_Error +=  COIDEF_NumHess   ( CntVect, 4  );       /* 4 Hessian nonzeros                */
   COI_Error +=  COIDEF_OptDir    ( CntVect, 1  );       /* Maximize                          */
   COI_Error +=  COIDEF_ObjCon    ( CntVect, 0  );       /* Objective is constraint 0         */
   COI_Error +=  COIDEF_DebugFV   ( CntVect, 0  );       /* 0 means no debugging, 1 each iter */
   COI_Error +=  COIDEF_Debug2D   ( CntVect, 0  );       /* 0 means no debugging, 1 each iter */
   COI_Error +=  COIDEF_StdOut    ( CntVect, 0  );       /* 1 means Allow output to StdOut    */
   COI_Error +=  COIDEF_Optfile   ( CntVect, "tutorial.opt"); /* Register the Options file    */
 /*
   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, &Tut_ReadMatrix); /* Register the callback ReadMatrix     */
   COI_Error +=  COIDEF_FDEval    ( CntVect, &Tut_FDEval);     /* Register the callback FDEval         */
   COI_Error +=  COIDEF_2DLagrStr ( CntVect, &Tut_2DLagrStr);  /* Register the callback Tut_2DLagrStr  */
   COI_Error +=  COIDEF_2DLagrVal ( CntVect, &Tut_2DLagrVal);  /* Register the callback Tut_2DLagrVal  */
 
#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-0.572943 ) > 0.000001 ) {
      c_log( "Incorrect objective returned", -1); }
 
   c_log( "Successful Solve", OK );
}