tutorial2str.cpp

Go to the documentation of this file.

 
#include <cstdio>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <iostream>
#include <set>
#include <utility>
#include <vector>
#include "conopt.hpp"
 
#include <adolc/adolc.h>

struct PairComp {
   bool operator() (const std::pair<int, int> a, const std::pair<int, int> b) const {
      return a.first < b.first || (a.first == b.first && a.second < b.second);
   }
};
 
class Tut_ModelData : public ConoptModelData
{
public:
   Tut_ModelData()
      : ConoptModelData()
   {
   }

   void evaluateNonlinearExpression(adouble* x, adouble* g, int rowno, int numvar)
   {
      /*                                                                   */
      /*   Declare local copies of the optimization variables. This is     */
      /*   just for convenience to make the expressions easier to read.    */
      /*                                                                   */
      adouble 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;
 
      /* helper variables */
      adouble hold1 = 0, hold2 = 0;
 
      /*                                                                   */
      /*   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 ) {
         *g    = P * Out;
      }
      /*                                                                   */
      /*   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 ));
         *g = hold2;
      }
      /*                                                                   */
      /*   Row = 2: The row is linear and will not be called.              */
      /*                                                                   */
 
   }

   void initialiseAutoDiff()
   {
      int numcons = numCons();
      int numvar = numVar();
 
      for (int c = 0; c < numcons; c++)
      {
         double res;
         /* starting the trace for constraint c */
         trace_on(c);
 
         /* these are variable types for ADOL-C. Both the dependent and independent variables must be declared as
          * adouble
          */
         adouble* ax;
         adouble ag = 0;
 
         /* setting the values of the independent variables */
         ax = new adouble[numvar];
         for (int i = 0; i < numvar; i++)
            ax[i] <<= getVariable(i).curr;
 
         evaluateNonlinearExpression(ax, &ag, c, numvar);
 
         /* setting the results as the value of the independent variable. */
         ag >>= res;
 
         delete[] ax;
 
         /* stopping the trace for constraint c */
         trace_off();
      }
 
   }

   void computeHessianStructure()
   {
      int numcons = numCons();
      int numvar = numVar();
 
      double* x;
      unsigned int** hesspat;
      std::set<std::pair<int, int>, PairComp> hesspairs;
      std::vector<int> rowidx;
      std::vector<int> colidx;
 
      x = new double[numvar];
      hesspat = new unsigned int*[numvar];
      for (int i = 0; i < numvar; i++)
         hesspat[i] = new unsigned int[numvar];
 
      /* assigning initial values to the x vector */
      for (int i = 0; i < numvar; i++)
         x[i] = getVariable(i).curr;
      
      /* computing the structure of the hessian */
      for (int c = 0; c < numcons; c++)
      {
         /* using the hess_pat() method from ADOL-C to compute the structure of the Hesssion */
         hess_pat(c, numvar, x, hesspat, 0);
 
         for (unsigned int i = 0; i < numvar; i++)
         {
            for (unsigned int j = 0; j < hesspat[i][0]; j++)
            {
               /* we are only interested in the lower triangle */
               if (hesspat[i][j + 1] <= i)
                  hesspairs.insert(std::make_pair(i, hesspat[i][j + 1]));
            }
         }
      }
 
      for (const auto &p : hesspairs)
      {
         rowidx.push_back(p.first);
         colidx.push_back(p.second);
      }
 
      for (int i = numvar - 1; i >= 0; i--)
         delete[] hesspat[i];
      delete[] hesspat;
      delete[] x;
 
      /* setting the structure of the hessian */
      setSDLagrangianStructure(rowidx, colidx);
   }

   void buildModel()
   {
      /*                                                               */
      /*   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:    */
      /*                                                               */
      addVariable(0.1, CONOPT_INF, 0.5);
      /*                                                               */
      /*   Lower bound on INP = X[1] = 0.1 and initial value = 0.5:    */
      /*                                                               */
      addVariable(0.1, CONOPT_INF, 0.5);
      /*                                                               */
      /*   Lower bound on OUT = X[2] and P = X[3] are both 0 and the   */
      /*   default initial value of 0 is used:                         */
      /*                                                               */
      addVariable(0., CONOPT_INF);
      addVariable(0., CONOPT_INF);
      /*                                                               */
      /*   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                             */
      /*                                                               */
      addConstraint(ConoptConstraintType::Free, -0.1, {0, 1, 2, 3}, {-1, -1, 0, 0}, {0, 0, 1, 1});
      /*                                                               */
      /*   Constraint 1 (Production Function)                          */
      /*   Rhs = 0 and type Equality                                   */
      /*                                                               */
      addConstraint(ConoptConstraintType::Eq, 0.0, {0, 1, 2}, {0, 0, -1}, {1, 1, 0});
      /*                                                               */
      /*   Constraint 2 (Price equation)                               */
      /*   Rhs = 4.0 and type Equality                                 */
      /*                                                               */
      addConstraint(ConoptConstraintType::Eq, 4.0, {2, 3}, {1, 2}, {0, 0});
 
      /* setting the objective constraint */
      setObjectiveElement(ConoptObjectiveElement::Constraint, 0);
 
      /* setting the optimisation direction */
      setOptimizationSense(ConoptSense::Maximize);
 
      /* initialising the automatic differentiation */
      initialiseAutoDiff();
 
      /* computes and sets the Hessian structure */
      computeHessianStructure();
   }

   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)
   {
      /* using the function() method from ADOL-C to compute the function value from the trace */
      if (mode == 1 || mode == 3)
         function(rowno, 1, numvar, const_cast<double*>(x), g);
 
      /* using the gradient() method from ADOL-C to compute the gradient from the trace */
      if (mode == 2 || mode == 3)
         gradient(rowno, numvar, x, jac);
 
      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)
   {
      double** hessres;
 
      // allocating memory for the hessian result
      hessres = new double*[numvar];
      for(int i = 0; i < numvar; i++)
         hessres[i] = new double[numvar];
 
      for(int c = 0; c < numcon; c++)
      {
         /* using the hessian() method from ADOL-C to compute the hessian value from the trace */
         hessian(c, numvar, const_cast<double*>(x), hessres);
 
         for(int i = 0; i < nhess; i++)
         {
            hsvl[i] += u[c]*hessres[hsrw[i]][hscl[i]];
         }
      }
 
      for (int i = numvar - 1; i >= 0; i--)
         delete[] hessres[i];
      delete[] hessres;
 
      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
   ConoptCpp conopt(pname);
   Tut_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 )
      cpp_log(conopt, "Skipping COI_Solve due to license error. COI_Error = " + COI_Error, COI_Error);
 
   COI_Error =  conopt.solve();           /* Optimize                      */
 
   // checking the statuses and objective value
   if ( conopt.modelStatus() != 2 || conopt.solutionStatus() != 1 )
   {
      cpp_log(conopt, "Incorrect Model or Solver Status", -1);
   }
   else if ( fabs( conopt.objectiveValue() - 0.572943 ) > 0.000001 )
   {
      cpp_log(conopt, "Incorrect objective returned", -1);
   }
 
   // printing the final status of the optimisation
   conopt.printStatus();
 
   cpp_log(conopt, "Successful Solve", COI_Error);
}