tutorial2r.cpp

Go to the documentation of this file.

 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <iostream>
#include "conopt.hpp"
 
class Tut2r_ModelData : public ConoptModelData
{
public:
   /* declaring the parameters */
   double Al;
   double Ak;
   double Ainp;
   double Rho;
   double K;
 
   /* declaring the variable and constraint indices. */
   int varl;
   int varinp;
   int varout;
   int varp;
   int varint;
   int consobj;
   int consprod;
   int consnew;
 
   Tut2r_ModelData()
       : Al(0.16), Ak(2.0), Ainp(0.16), Rho(1.0), K(4.0), varl(-1), varinp(-1), varout(-1),
         varp(-1), consobj(-1), consprod(-1)
   {
   }

   void buildModel()
   {
 
      double initL = 0.5;
      double initInp = 0.5;
      /*                                                               */
      /*   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:    */
      /*                                                               */
      varl = addVariable(0.1, Conopt::Infinity, initL);
      /*                                                               */
      /*   Lower bound on INP = X[1] = 0.1 and initial value = 0.5:    */
      /*                                                               */
      varinp = addVariable(0.1, Conopt::Infinity, initInp);
      /*                                                               */
      /*   Lower bound on OUT = X[2] and P = X[3] are both 0 and the   */
      /*   default initial value of 0 is used:                         */
      /*                                                               */
      varout = addVariable(0., Conopt::Infinity);
      varp = addVariable(0., Conopt::Infinity);
      /*                                                               */
      /*   Lower bound on the new Int is set to 0.01 and the initial   */
      /*   value is defined from row 1                                 */
      /*                                                               */
      varint = addVariable(
            0., Conopt::Infinity, Al * pow(initL, -Rho) + Ak * pow(K, -Rho) + Ainp * pow(initInp, -Rho));
      /*                                                               */
      /*   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                             */
      /*                                                               */
      consobj = addConstraint(ConoptConstraintType::Free, -0.1, {varl, varinp, varout, varp},
            {-1, -1, 0, 0}, {0, 0, 1, 1});
      /*                                                               */
      /*   Constraint 1 (Production Function)                          */
      /*   Rhs = 0 and type Equality                                   */
      /*                                                               */
      consprod = addConstraint(
            ConoptConstraintType::Eq, 0.0, {varl, varinp, varint}, {0, 0, -1}, {1, 1, 0});
      /*                                                               */
      /*   Constraint 2 (Price equation)                               */
      /*   Rhs = 4.0 and type Equality                                 */
      /*                                                               */
      addConstraint(ConoptConstraintType::Eq, 4.0, {varout, varp}, {1, 2}, {0, 0});
      /*                                                               */
      /*   Constraint 3 (New equation)                                 */
      /*   Rhs = 0.0 and type Equality                                 */
      /*                                                               */
      consnew = addConstraint(ConoptConstraintType::Eq, 0.0, {varout, varint}, {-1, 0}, {0, 1});
 
      /* setting the objective constraint */
      setObjectiveElement(ConoptObjectiveElement::Constraint, consobj);
 
      /* setting the optimisation direction */
      setOptimizationSense(ConoptSense::Maximize);
 
      /* setting the structure of the hessian */
      setSDLagrangianStructure({0, 1, 3, 4}, {0, 1, 2, 4});
   }

   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
   {
      /*                                                                   */
      /*   Declare local copies of the optimization variables. This is     */
      /*   just for convenience to make the expressions easier to read.    */
      /*                                                                   */
      double L, Inp, Out, P, Int;
 
      /*                                                                   */
      /*   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[varl];
      Inp = x[varinp];
      Out = x[varout];
      P = x[varp];
      Int = x[varint];
      /*                                                                   */
      /*   Row 0: the objective function is nonlinear                      */
      /*                                                                   */
 
      if (rowno == consobj)
      {
         /*                                                                   */
         /*   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[varout] = P; /* derivative w.r.t. Out = X[2]               */
            jac[varp] = Out; /* derivative w.r.t. P   = X[3]               */
         }
      }
      /*                                                                   */
      /*   Row 1: The first part of the production function is nonlinear   */
      /*          Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho))          */
      /*                                                                   */
      else if (rowno == consprod)
      {
         /*                                                                   */
         /*   Mode = 1 or 3: Function value                                   */
         /*                                                                   */
         if (mode == 1 || mode == 3)
            *g = (Al * pow(L, (-Rho)) + Ak * pow(K, (-Rho)) + Ainp * pow(Inp, (-Rho)));
         /*                                                                   */
         /*   Mode = 2 or 3: Derivatives                                      */
         /*                                                                   */
         if (mode == 2 || mode == 3)
         {
            jac[varl] = Al * (-Rho) * pow(L, (-Rho - 1.));       /* derivative w.r.t. L   = X[0]  */
            jac[varinp] = Ainp * (-Rho) * pow(Inp, (-Rho - 1.)); /* derivative w.r.t. Inp = X[1]  */
         }
      }
      /*                                                                   */
      /*   Row = 2: The row is linear and will not be called.              */
      /*                                                                   */
      /*   Row = 3: The second part of the production function             */
      /*            Int**(-1/Rho)                                          */
      /*                                                                   */
      else if (rowno == consnew)
      {
         /*                                                                   */
         /*   Mode = 1 or 3: Function value                                   */
         /*                                                                   */
         if (mode == 1 || mode == 3)
            *g = pow(Int, (-1.0 / Rho));
         /*                                                                   */
         /*   Mode = 2 or 3: Derivatives                                      */
         /*                                                                   */
         if (mode == 2 || mode == 3)
         {
            jac[varint] = (-1.0 / Rho) * pow(Int, (-1.0 / Rho - 1.0));
         }
      }
 
      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
   {
      /*
        Declare local copies of the optimization variables. This is
        just for convenience to make the expressions easier to read.        */
 
      double L, Inp;
 
      /*  Normal Evaluation mode                                            */
 
      hsvl[2] = u[consobj]; /* Second derivative of constraint 0 is 1       */
      if (u[consprod] != 0.0)
      {
         L = x[varl];
         Inp = x[varinp];
         hsvl[0] = (-Rho) * (-Rho - 1.0) * Al * pow(L, -Rho - 2.0) * u[consprod];
         hsvl[1] = (-Rho) * (-Rho - 1.0) * Ainp * pow(Inp, -Rho - 2.0) * u[consprod];
      };
      if (u[consnew] != 0.0)
         hsvl[3] = (-1.0 / Rho) * (-1.0 / Rho - 1.0) * pow(x[varint], -1.0 / Rho - 2.0);
 
      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);
   Tut2r_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() - 0.572943) > 0.000001)
      cpp_log(conopt, "Incorrect objective returned", -1);
 
   conopt.printStatus();
 
   cpp_log(conopt, "Successful Solve", COI_Error);
}