qp4.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"
 
/*  QP Model: Min (x-target)*Q*(x-target)/2 s.t.
              sum( x ) = 1;
    where x is a vector of non-negative variables of length nvar
          target is a vector of +10, and
          Q is a matrix with +1 on the diagonal and 0.1 on the
            first upper and lower bi-diagonal                         */
 
class QP4_ModelData : public ConoptModelData
{
public:
   std::vector<double> Q;      // Q matrix values (diagonal and off-diagonal)
   std::vector<int> Qrow;      // Row indices of Q matrix
   std::vector<int> Qcol;      // Column indices of Q matrix
   std::vector<double> Target; // Target vector
   int nvar;
   int nqnz;

   QP4_ModelData(int numvar) : nvar(numvar), nqnz(2 * numvar - 1)
   {
      Q.resize(nqnz);
      Qrow.resize(nqnz);
      Qcol.resize(nqnz);
      Target.resize(nvar, 10.0); // all targets set to 10.0
 
      int j = 0;
      for (int i = 0; i < nvar; ++i)
      {
         Q[j] = 1.0;
         Qrow[j] = i;
         Qcol[j] = i;
         ++j;
 
         if (i < nvar - 1)
         {
            Q[j] = 0.1;
            Qrow[j] = i + 1;
            Qcol[j] = i;
            ++j;
         }
      }
   }

   void buildModel()
   {
      /*                                                               */
      /*   Information about Variables:                                */
      /*   Default: Lower = -Inf, Curr = 0, and Upper = +inf.          */
      /*   Default: the status information in Vsta is not used.        */
      /*                                                               */
      /*   Lower bound = 0 on all variables:                           */
      /*                                                               */
      for (int i = 0; i < nvar; ++i)
         addVariable(0.0, Conopt::Infinity); // lower = 0, upper = +∞
 
      /*   Build common index/coeff lists                              */
      std::vector<int> varIdx(nvar);
      std::vector<double> ones(nvar, 1.0);    // coefficients
      std::vector<int> nlflags_free(nvar, 1); // nonlinear for objective
      std::vector<int> nlflags_sum(nvar, 0);  // linear for sum-to-one constraint
 
      for (int i = 0; i < nvar; ++i)
         varIdx[i] = i;
 
      /*                                                               */
      /*   Constraint 0: Objective                                     */
      /*   Rhs = 0 and type Non binding                                */
      /*                                                               */
      addConstraint(ConoptConstraintType::Free, 0.0, varIdx, ones, nlflags_free);
 
      /*                                                               */
      /*   Constraint 1: Objective                                     */
      /*   Rhs = 1 and type Equality                                   */
      /*                                                               */
      addConstraint(ConoptConstraintType::Eq, 1.0, varIdx, ones, nlflags_sum);
 
      /* setting the objective constraint */
      setObjectiveElement(ConoptObjectiveElement::Constraint, 0);
 
      /* setting the optimisation direction */
      setOptimizationSense(ConoptSense::Minimize);
 
      /* setting the structure of the hessian */
      setSDLagrangianStructure(Qrow, Qcol);
   }

   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
   {
      if (rowno == 0)
      {
         // Mode 1 or 3: Evaluate function value
         if (mode == 1 || mode == 3)
         {
            double sum = 0.0;
            for (int k = 0; k < nqnz; ++k)
            {
               int i = Qrow[k], j = Qcol[k];
               double qi = Q[k];
               if (i == j)
                  sum += (x[i] - Target[i]) * qi * (x[i] - Target[i]);
               else
                  sum += 2.0 * (x[i] - Target[i]) * qi * (x[j] - Target[j]);
            }
            *g = sum / 2.0;
         }
 
         // Mode 2 or 3: Evaluate gradient
         if (mode == 2 || mode == 3)
         {
            for (int i = 0; i < numvar; ++i)
               jac[i] = 0.0;
 
            for (int k = 0; k < nqnz; ++k)
            {
               int i = Qrow[k], j = Qcol[k];
               double qi = Q[k];
               if (i == j)
                  jac[i] += qi * (x[i] - Target[i]);
               else
               {
                  jac[i] += qi * (x[j] - Target[j]);
                  jac[j] += qi * (x[i] - Target[i]);
               }
            }
         }
      }
 
      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
   {
      int i, j, k;
      /*                                                                    */
      /*   Is only called for rowno = 0, the nonlinear objective            */
      /*   Return H*Dx = Q*Dx in d2g                                        */
      /*                                                                    */
      if (rowno == 0)
      {
         for (i = 0; i < numvar; i++)
            d2g[i] = 0.0;
         for (k = 0; k < nqnz; k++)
         {
            i = Qrow[k];
            j = Qcol[k];
            if (i == j)
               d2g[i] += Q[k] * dx[i];
            else
            {
               d2g[i] += Q[k] * dx[j];
               d2g[j] += Q[k] * dx[i];
            }
         }
      }
      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
   {
      int k;
 
      for (k = 0; k < nqnz; k++)
         hsvl[k] = Q[k] * u[0];
      return 0;
   }
};
 
#include "std.cpp"

int main(int argc, char **argv)
{
   int COI_Error = 0;
 
   // the number of variables in the problem
   int numvar = 1000;
 
   // getting the program name from the executable path
   std::string pname = getProgramName(argv[0]);
 
   // initialising the Conopt Object
   Conopt conopt(pname);
   QP4_ModelData modeldata(numvar);
   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() - 59978.0) > 0.001)
      cpp_log(conopt, "Incorrect objective returned", -1);
 
   conopt.printStatus();
 
   cpp_log(conopt, "Successful Solve", COI_Error);
}