OctaveからCPLEXを使う

概要

• David Musicant氏とClaude Tadonki氏が作成した MatlabからCPLEXを使うためのmexファイルを参考にして作成しました。 Octaveもmexファイルをつかえるので， それを使えばよいのですが，liboctaveを使うように書きました。
• 凸2次問題と線形計画問題のために，qp_cplex.ccとlp_cplex.ccを作成しました。
• バージョンが12.2と12.3でソースファイルが1行異なり(ソースファイルのどちらかをコメントアウトしてください)， コンパイルのコマンドも異なりますので，注意してください。
• 内容は無保証です。

qp_cplexの使い方

• 形式

   [x, lambda, info] = qp_cplex(H, c, A, b, lb, ub, x0, nEq)
  最適化  1/2 * (x' * H * x) + c' * x
  条件    A x (= or <=) b    (初めのnEq個が等式，残りが不等式)
          vb <= x <= ub 
  出力
  　　　x      : 最適解
  　　　lambda : ラグランジュ未定乗数
  　　　info   : 解に関する情報
  

• 使用例(Octave testQP.m)

  H = [8 2 3 ; 2 6 1 ; 3 1 7]
  c = [1 -1 1]
  A = [1 1 1]
  b = [3]
  x0 = zeros(1, 3);
  lb = - 1e+50 * ones(1, 3);
  ub = 1e+50 * ones(1, 3);
  [x, lambda, info] = qp_cplex(H, c ,A, b, lb, ub, x0, 1);
  

lp_cplexの使い方

• 形式

   [x, lambda, info] = lp_cplex(c, A, b, lb, ub, x0, ctype)
  最適化  c' * x
  条件    A x (= or <= or <=) b  (ctype[i] = 69, 71, 76のときそれぞれ，>=, =, <=)
          vb <= x <= ub 
  出力
  　　　x      : 最適解
  　　　lambda : ラグランジュ未定乗数
  　　　info   : 解に関する情報
  
  

• 使用例(Octave testLP.m)

  c  = [-1 -1 -1 -1]
  A  = [2 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1 ; 0 1 0 1 ; 1 0 1 0]
  b  = [1 ; 1 ; 1 ; -1 ; 1 ; 1];
  lb = [-2 -2 -2 -2];
  ub = [ 2  2  2  2];
  ctype = [76 76 76 69 76 76];
  [x lambda info] = lp_cplex(c, A, b, lb, ub, ctype)
  

qp_cplexのソースとコンパイル

• ソース

  /*********************************************************************
   * Proguram : qp_cplex_oct.cc                                        *
   * Author: Yukihiko Yamashita                                        *
   * Thanks to                                                         *
   *   Program: qp_cplex_mex.c                                         *
   *   Author  David Musicant & Claude Tadonki                         *
   *********************************************************************/
  
  #include 
  #include 
  #include 
  
  #define Inf 1.0e+30
  
  #define MINIMIZE 1
  #define MAX_ITER_DEFAULT 10000
  
  /* Input arguments */
  #define H_IN   0
  #define c_IN   1
  #define A_IN   2
  #define b_IN   3
  #define lb_IN  4
  #define ub_IN  5
  #define x0_IN  6
  #define nEq_IN 7
  
  /* Output arguments */
  #define x_OUT      0
  #define lambda_OUT 1
  #define STAT_OUT   2
  
  extern "C" {
    /* remove the folloing two backslash for 12.2 */
    // #include "/opt/ILOG/CPLEX_Studio_AcademicResearch122/cplex/include/ilcplex/cplex.h" 
    /* remove the folloing two backslash for 12.3 */
    #include "/opt/ibm/ILOG/CPLEX_Studio_Academic123/cplex/include/ilcplex/cplex.h" 
  }
  
  DEFUN_DLD(qp_cplex, args, ,
             "Return Usage: [x, lambda, info] = qp(H, c ,A, b, lb, ub, x0, nEq)")
  {
  
  
    int nArgs = args.length();
  
    if (nArgs != 8) {
      octave_stdout << "Incorrect number of arguments \n";
      return octave_value(DiagMatrix(1,1,1.0));
    }
  
    CPXENVptr env;
    CPXLPptr  lp = NULL;
    int       maxIter = MAX_ITER_DEFAULT;
    int       status, info;
    char      errorMsg[256];
  
    int dim   = (args(H_IN).matrix_value()).rows();
    int nCond = (args(A_IN).matrix_value()).rows();
  
    NDArray HMat   = args(H_IN).array_value();
    NDArray cVect  = args(c_IN).array_value();
    NDArray AMat   = args(A_IN).array_value();
    NDArray bVect  = args(b_IN).array_value();
    NDArray lbVect = args(lb_IN).array_value();
    NDArray ubVect = args(ub_IN).array_value();
    NDArray x0Vect = args(x0_IN).array_value();
    int     nEq    = args(nEq_IN).int_value();
  
    
    double *H  = (double *) malloc( dim   * dim * sizeof(double));
    double *c  = (double *) malloc( dim         * sizeof(double)); 
    double *A  = (double *) malloc( nCond * dim * sizeof(double)); 
    double *b  = (double *) malloc( nCond       * sizeof(double)); 
    double *lb = (double *) malloc( dim         * sizeof(double));
    double *ub = (double *) malloc( dim         * sizeof(double));
    double *x0 = (double *) malloc( dim         * sizeof(double)); 
  
    
    double *objVal   = (double *) malloc(         sizeof(double));
    double *x        = (double *) malloc( dim   * sizeof(double));
    double *lambda   = (double *) malloc( nCond * sizeof(double));
    double  dStat;
    double *colStat  = (double *) malloc(         sizeof(double));
    int    *iColStat = (int *)    malloc( dim  * sizeof(int));
    double *iterCnt  = (double *) malloc(        sizeof(double));
  
    /* Build the matrix of the quadratic part */
    /* Only dense is assumed */
    int   *HcolLoc = (int *) malloc( dim       * sizeof(int));
    int   *Hcount  = (int *) malloc( dim       * sizeof(int));
    int   *HrowInd = (int *) malloc( dim * dim * sizeof(int));
  
    for (int j = 0 ; j < dim ; ++j) {
      HcolLoc[j] = j * dim;   /* start of each column*/
      Hcount[j]  = dim;       /* number of nonzero value in the column */
      for (int i = 0 ; i < dim ; ++i) {
        HrowInd[j * dim + i] = i; /* row location*/
        H[j * dim + i] = HMat(i,j);
      }
      c[j]  = cVect(j);
      lb[j] = (lbVect(j) < -Inf ? -CPX_INFBOUND : lbVect(j));
      ub[j] = (ubVect(j) >  Inf ?  CPX_INFBOUND : ubVect(j));
      x0[j] =  x0Vect(j);
    }
  
    /* Build the matrix of the lienar part */
    /* Only dense is assumed */
    int   *AcolLoc = (int *) malloc( dim       * sizeof(int));
    int   *Acount  = (int *) malloc( dim       * sizeof(int));
    int   *ArowInd = (int *) malloc( nCond * dim * sizeof(int));
  
    for (int j = 0 ; j < dim ; ++j) {
      AcolLoc[j] = j * nCond;  /* column location */
      Acount[j]  = nCond;       /* number of nonzero value in the column */
      for (int i = 0 ; i < nCond ; ++i) {
        ArowInd[j * nCond + i] = i;  /* row location*/
        A[j * nCond + i] = AMat(i,j);
      }
    }
  
    char *sense = (char *) malloc(nCond * sizeof(char)); 
    int loopEq = 0;
    while (loopEq < nEq)   sense[loopEq++] = 'E';
    while (loopEq < nCond) sense[loopEq++] = 'L';
    
    for(int i = 0 ; i < nCond ; i++) {
      b[i] = bVect(i);
    }
  
    /* Open CPLEX environment */
    /*  env = CPXopenCPLEXdevelop(&status); */
    env = CPXopenCPLEX(&status);
    if (! env) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      octave_stdout << "\nCould not open CPLEX environment.\n";
    }
    
    /* Create CPLEX problem space */
    lp = CPXcreateprob(env, &status, "octave");
    if (! lp) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXcloseCPLEX(& env);
      octave_stdout << "\nCould not create CPLEX problem.\n";
    }
          
    /* Copy LP into CPLEX environment */
    status = CPXcopylp(env, lp, dim, nCond, MINIMIZE, c, b, sense,
                       AcolLoc, Acount, ArowInd, A, lb, ub, NULL);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nFailed to copy LP.\n";
    }
               
   /* Copy QP into CPLEX environment */
    status = CPXcopyquad(env, lp, HcolLoc, Hcount, HrowInd, H);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env, &lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nFailed to copy quadratic matrix.\n";
    }
    
    /* Set iteration limit. */
    status = CPXsetintparam(env, CPX_PARAM_ITLIM, maxIter);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nCould not set number of iterations.\n";
    }
  
    /* Perform optimization  status = CPXbaropt(env,lp); */
    status = CPXbaropt(env,lp);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nOptimization error.\n";
    }
  
    /* Obtain solution */
    double lpstat;
    status = CPXsolution(env, lp, &info, objVal, x, lambda, NULL, NULL);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env, &lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nFailure when retrieving solution.\n";
    }
  
    /* Get status of columns 
    status = CPXgetbase(env, lp, istat, NULL);
    if (status) {
      printf(CPXgeterrorstring(env,status,errorMsg));
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      mexErrMsgTxt("\nUnable to get basis status.");
    }*/
   
    /* Copy int column values to double column values */
    //for (i=0; i < cols; i++)
    //  cstat[i] = istat[i];
  
    // Clean up CPLEX
    CPXfreeprob(env,&lp); 
    CPXcloseCPLEX(&env);
  
  
    /* Return values */
  
    ColumnVector xVect(dim);
    ColumnVector lambdaVect(nCond);
    ColumnVector infoVect(1);
  
    for(int j = 0 ; j < dim ; j++) xVect(j)      = x[j];
    for(int i = 0 ; i < nCond ; i++) lambdaVect(i) = lambda[i];
    infoVect(0) = info;
  
    octave_value_list ret = octave_value_list(octave_value(xVect));
    ret.append(octave_value(lambdaVect));
    ret.append(octave_value(infoVect));
  
    // Free up memory
    free(H);
    free(c);
    free(A);
    free(b);
    free(lb);
    free(ub);
    free(x0);
  
    free(objVal);
    free(x);
    free(lambda);
    free(colStat);
    free(iColStat);
    free(iterCnt);
  
    free(HcolLoc);
    free(Hcount);
    free(HrowInd);
    free(AcolLoc);
    free(Acount);
    free(ArowInd);
  
    return ret;
  }
  

• コンパイル方法
  12.2

  mkoctfile qp_cplex.cc -I/opt/ILOG/CPLEX_Studio_AcademicResearch122/cplex/include/ilcplex -L/opt/ILOG/CPLEX_Studio_AcademicResearch122/cplex/lib/x86-64_sles10_4.1/static_pic -lpthread -lcplex
  

  12.3

  mkoctfile qp_cplex.cc -I/opt/ibm/ILOG/CPLEX_Studio_Academic123/cplex/include/ilcplex -L/opt/ibm/ILOG/CPLEX_Studio_Academic123/cplex/lib/x86-64_sles10_4.1/static_pic -lpthread -lcplex
  

lp_cplexのソースとコンパイル

• ソース

  /*********************************************************************
   * Proguram : lq_cplex_oct.cc                                        *
   * Author: Yukihiko Yamashita                                        *
   * Thanks to                                                         *
   *   Program: lp_cplex_mex.c                                         *
   *   Author  David Musicant & Claude Tadonki                         *
   *********************************************************************/
  
  #include 
  #include 
  #include 
  
  #define Inf 1.0e+30
  
  #define MINIMIZE 1
  #define MAX_ITER_DEFAULT 10000
  
  // Large and small
  #define LARGE 1.0e+8
  #define SMALL 1.0e-8
  
  /* Input arguments */
  #define c_IN     0
  #define A_IN     1
  #define b_IN     2
  #define lb_IN    3
  #define ub_IN    4
  #define ctype_IN 5
  
  /* Output arguments */
  #define x_OUT      0
  #define lambda_OUT 1
  #define STAT_OUT   2
  
  extern "C" {
    /* remove the folloing two backslash for 12.2 */
    // #include "/opt/ILOG/CPLEX_Studio_AcademicResearch122/cplex/include/ilcplex/cplex.h" 
    /* remove the folloing two backslash for 12.3 */
    #include "/opt/ibm/ILOG/CPLEX_Studio_Academic123/cplex/include/ilcplex/cplex.h" 
  }
  
  DEFUN_DLD(lp_cplex, args, ,
             "Return Usage: [x, lambda, info] = lp_cplex(c, A, b, lb, ub, ctype)")
  {
  
    int nArgs = args.length();
  
    if (nArgs != 6) {
      octave_stdout << "Incorrect number of arguments \n" <<
        " Usage: [x, lambda, info] = lp_cplex(c, A, b, lb, ub, ctype)\n";
      return octave_value(DiagMatrix(1,1,1.0));
    }
  
    CPXENVptr env;
    CPXLPptr  lp = NULL;
    int       maxIter = MAX_ITER_DEFAULT;
    int       status, info;
    char      errorMsg[256];
  
    int dim   = (args(A_IN).matrix_value()).columns();
    int nCond = (args(A_IN).matrix_value()).rows();
  
    NDArray cVect  = args(c_IN).array_value();
    NDArray AMat   = args(A_IN).array_value();
    NDArray bVect  = args(b_IN).array_value();
    NDArray lbVect = args(lb_IN).array_value();
    NDArray ubVect = args(ub_IN).array_value();
    NDArray ctypeV = args(ctype_IN).array_value(); // condition type
  
    double *objVal   = (double *) malloc(         sizeof(double));
    double *x        = (double *) malloc( dim   * sizeof(double));
    double *lambda   = (double *) malloc( nCond * sizeof(double));
  
    double *colStat  = (double *) malloc(        sizeof(double));
    int    *iColStat = (int *)    malloc( dim  * sizeof(int));
    double *iterCnt  = (double *) malloc(        sizeof(double));
  
    // Build optimum condition, ub, lb, initial vector 
    // count nonzero value in A
    double *c  = (double *) malloc( dim         * sizeof(double)); 
    double *b  = (double *) malloc( nCond       * sizeof(double)); 
    double *lb = (double *) malloc( dim         * sizeof(double));
    double *ub = (double *) malloc( dim         * sizeof(double));
    char   *sense = (char *) malloc(nCond * sizeof(char)); 
  
    int    AnNZ     = 0;
    int   *AcolLoc  = (int *) malloc( dim       * sizeof(int));
    int   *Acount   = (int *) malloc( dim       * sizeof(int));
  
    for (int j = 0 ; j < dim ; ++j) {
      c[j]       = cVect(j);
      lb[j]      = (lbVect(j) < -LARGE ? -CPX_INFBOUND : lbVect(j));
      ub[j]      = (ubVect(j) >  LARGE ?  CPX_INFBOUND : ubVect(j));
      Acount[j]  = 0;
      AcolLoc[j] = AnNZ;
      x[j]       = 0.0;
      for (int i = 0 ; i < nCond ; ++i) {
        if (fabs(AMat(i,j)) > SMALL) {
  	Acount[j] += 1;
  	//octave_stdout << "i,j = " << i << " " << j <<"\n";
  	++AnNZ;
        }
      }
    }
  
    // Build the matrix of the lienar part
    // Only dense is assumed
    double *Aval    = (double *) malloc( AnNZ * sizeof(double)); 
    int    *ArowInd = (int *) malloc( AnNZ  * sizeof(int));
  
    int pos = 0;
    for (int j = 0 ; j < dim ; ++j) {
      for (int i = 0 ; i < nCond ; ++i) {
        if (fabs(AMat(i,j)) > SMALL) {
  	ArowInd[pos] = i;
  	Aval[pos++]  = AMat(i,j);
        }
      }
    }
  
    /*
    pos = 0;
    for (int j = 0 ; j < dim ; ++j) {
      octave_stdout << j << " " << c[j] << " " << lb[j] << " " << ub[j] << " " << Acount[j] << " " << AcolLoc[j] <<"\n" ;
      for (int i = 0 ; i < nCond ; ++i) {
        if (fabs(AMat(i,j)) > SMALL) {
  	octave_stdout << ArowInd[pos] << " " << Aval[pos] <<"\n";
  	++pos;
        }
      }
    }
    */
  
    // Set b and ctype
    for (int i = 0 ; i < nCond ; ++i){
      sense[i]  = (char) (ctypeV(i) + 0.5);
      b[i]      = bVect(i);
      lambda[i] = 0.0;
      // octave_stdout << sense[i] << " " << b[i] << "  ";
    }
    octave_stdout << "\n";
  
  
    /* Open CPLEX environment */
    /* env = CPXopenCPLEXdevelop(&status); */
    env = CPXopenCPLEX(&status);
    if (! env) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      octave_stdout << "\nCould not open CPLEX environment.\n";
    }
  
    /* Create CPLEX problem space */
    lp = CPXcreateprob(env, &status, "octave");
    if (! lp) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXcloseCPLEX(& env);
      octave_stdout << "\nCould not create CPLEX problem.\n";
    }
          
    /* Copy LP into CPLEX environment */
    status = CPXcopylp(env, lp, dim, nCond, MINIMIZE, c, b, sense,
                       AcolLoc, Acount, ArowInd, Aval, lb, ub, NULL);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nFailed to copy LP.\n";
    }
               
    status = CPXsetintparam(env, CPX_PARAM_ITLIM, maxIter);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nCould not set number of iterations.\n";
    }
  
    // ************* Perform optimization ***************
    status = CPXprimopt(env,lp);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nOptimization error.\n";
    }
  
    // Obtain Solution
    status = CPXsolution(env, lp, &info, objVal, x, lambda, NULL, NULL);
    if (status) {
      octave_stdout << CPXgeterrorstring(env,status,errorMsg);
      CPXfreeprob(env, &lp); 
      CPXcloseCPLEX(&env);
      octave_stdout << "\nFailure when retrieving solution.\n";
    }
  
    /* Get status of columns 
    status = CPXgetbase(env, lp, istat, NULL);
    if (status) {
      printf(CPXgeterrorstring(env,status,errorMsg));
      CPXfreeprob(env,&lp); 
      CPXcloseCPLEX(&env);
      mexErrMsgTxt("\nUnable to get basis status.");
    }
   
     // Copy int column values to double column values 
    for (i=0; i < cols; i++)  cstat[i] = istat[i];
    */
  
    /* Return values */
  
    ColumnVector xVect(dim);
    ColumnVector lambdaVect(nCond);
    ColumnVector infoVect(1);
  
    for(int j = 0 ; j < dim ; j++)   xVect(j) = x[j];
    for(int i = 0 ; i < nCond ; i++) lambdaVect(i) = lambda[i];
    infoVect(0) = info;
  
    // Clean up CPLEX
    CPXfreeprob(env,&lp); 
    CPXcloseCPLEX(&env);
  
  
    octave_value_list ret = octave_value_list(octave_value(xVect));
    ret.append(octave_value(lambdaVect));
    ret.append(octave_value(infoVect));
  
    // Free up memory
    free(c);
    free(Aval);
    free(b);
    free(lb);
    free(ub);
  
    free(objVal);
    free(x);
    free(lambda);
    free(colStat);
    free(iColStat);
    free(iterCnt);
  
    free(sense);
    free(AcolLoc);
    free(Acount);
    free(ArowInd);
  
    return ret;
  }
  

• コンパイル方法
  12.2

  mkoctfile lp_cplex.cc -I/opt/ILOG/CPLEX_Studio_AcademicResearch122/cplex/include/ilcplex -L/opt/ILOG/CPLEX_Studio_AcademicResearch122/cplex/lib/x86-64_sles10_4.1/static_pic -lpthread -lcplex
  

  12.2

  mkoctfile lp_cplex.cc -I/opt/ibm/ILOG/CPLEX_Studio_Academic123/cplex/include/ilcplex -L/opt/ibm/ILOG/CPLEX_Studio_Academic123/cplex/lib/x86-64_sles10_4.1/static_pic -lpthread -lcplex