6.3 Conic Quadratic Optimization

Conic optimization is a generalization of linear optimization, allowing constraints of the type

\[x^t \in \K_t,\]

where \(x^t\) is a subset of the problem variables and \(\K_t\) is a convex cone. Since the set \(\real^n\) of real numbers is also a convex cone, we can simply write a compound conic constraint \(x\in \K\) where \(\K=\K_1\times\cdots\times\K_l\) is a product of smaller cones and \(x\) is the full problem variable.

MOSEK can solve conic quadratic optimization problems of the form

\[\begin{split}\begin{array}{lccccl} \mbox{minimize} & & & c^T x + c^f & & \\ \mbox{subject to} & l^c & \leq & A x & \leq & u^c, \\ & l^x & \leq & x & \leq & u^x, \\ & & & x \in \K, & & \end{array}\end{split}\]

where the domain restriction, \(x \in \K\), implies that all variables are partitioned into convex cones

\[x = (x^0, x^1, \ldots , x^{p-1}),\quad \mbox{with } x^t \in \K_t \subseteq \real^{n_t}.\]

For convenience, a user defining a conic quadratic problem only needs to specify subsets of variables \(x^t\) belonging to quadratic cones. These are:

  • Quadratic cone:

    \[\Q^n = \left\lbrace x \in \real^n: x_0 \geq \sqrt{\sum_{j=1}^{n-1} x_j^2} \right\rbrace.\]
  • Rotated quadratic cone:

    \[\Qr^n = \left\lbrace x \in \real^n: 2 x_0 x_1 \geq \sum_{j=2}^{n-1} x_j^2,\quad x_0\geq 0,\quad x_1 \geq 0 \right\rbrace.\]

For example, the following constraint:

\[(x_4, x_0, x_2) \in \Q^3\]

describes a convex cone in \(\real^3\) given by the inequality:

\[x_4 \geq \sqrt{x_0^2 + x_2^2}.\]

Furthermore, each variable may belong to one cone at most. The constraint \(x_i - x_j = 0\) would however allow \(x_i\) and \(x_j\) to belong to different cones with same effect.

6.3.1 Example CQO1

Consider the following conic quadratic problem which involves some linear constraints, a quadratic cone and a rotated quadratic cone.

(1)\[\begin{split}\begin{array} {lccc} \mbox{minimize} & x_4 + x_5 + x_6 & & \\ \mbox{subject to} & x_1+x_2+ 2 x_3 & = & 1, \\ & x_1,x_2,x_3 & \geq & 0, \\ & x_4 \geq \sqrt{x_1^2 + x_2^2}, & & \\ & 2 x_5 x_6 \geq x_3^2 & & \end{array}\end{split}\]

Setting up the linear part

The linear parts (constraints, variables, objective) are set up using exactly the same methods as for linear problems, and we refer to Section 6.1 for all the details. The same applies to technical aspects such as defining an optimization task, retrieving the solution and so on.

Setting up the conic constraints

A cone is defined using the function MSK_appendcone:

        csub[0] = 3;
        csub[1] = 0;
        csub[2] = 1;

        r = MSK_appendcone(task,
                           MSK_CT_QUAD,
                           0.0, /* For future use only, can be set to 0.0 */
                           3,
                           csub);

The first argument selects the type of quadratic cone, in this case either MSK_CT_QUAD for a quadratic cone or MSK_CT_RQUAD for a rotated quadratic cone. The second parameter is currently ignored and passing 0.0 will work.

The next argument denotes the number of variables in the cone, in this case 3, and the last argument is a list of indexes of the variables appearing in the cone.

Variants of this method are available to append multiple cones at a time.

Source code

Listing 6.4 Source code solving problem (1). Click here to download.
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  MSKrescodee  r;

  const MSKint32t numvar = 6,
                  numcon = 1;

  MSKboundkeye bkc[] = { MSK_BK_FX };
  double       blc[] = { 1.0 };
  double       buc[] = { 1.0 };

  MSKboundkeye bkx[] = {MSK_BK_LO,
                        MSK_BK_LO,
                        MSK_BK_LO,
                        MSK_BK_FR,
                        MSK_BK_FR,
                        MSK_BK_FR
                       };
  double       blx[] = {0.0,
                        0.0,
                        0.0,
                        -MSK_INFINITY,
                        -MSK_INFINITY,
                        -MSK_INFINITY
                       };
  double       bux[] = { +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY,
                         +MSK_INFINITY
                       };

  double       c[]   = {0.0,
                        0.0,
                        0.0,
                        1.0,
                        1.0,
                        1.0
                       };

  MSKint32t   aptrb[] = {0, 1, 2, 3, 3, 3},
              aptre[] = {1, 2, 3, 3, 3, 3},
              asub[]  = {0, 0, 0, 0};
  double      aval[]  = {1.0, 1.0, 2.0};


  MSKint32t   i, j, csub[3];

  MSKenv_t    env  = NULL;
  MSKtask_t   task = NULL;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  if ( r == MSK_RES_OK )
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    if ( r == MSK_RES_OK )
    {
      MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

      /* Append 'numcon' empty constraints.
      The constraints will initially have no bounds. */
      if ( r == MSK_RES_OK )
        r = MSK_appendcons(task, numcon);

      /* Append 'numvar' variables.
      The variables will initially be fixed at zero (x=0). */
      if ( r == MSK_RES_OK )
        r = MSK_appendvars(task, numvar);

      for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
      {
        /* Set the linear term c_j in the objective.*/
        if (r == MSK_RES_OK)
          r = MSK_putcj(task, j, c[j]);

        /* Set the bounds on variable j.
        blx[j] <= x_j <= bux[j] */
        if (r == MSK_RES_OK)
          r = MSK_putvarbound(task,
                              j,           /* Index of variable.*/
                              bkx[j],      /* Bound key.*/
                              blx[j],      /* Numerical value of lower bound.*/
                              bux[j]);     /* Numerical value of upper bound.*/

        /* Input column j of A */
        if (r == MSK_RES_OK)
          r = MSK_putacol(task,
                          j,                 /* Variable (column) index.*/
                          aptre[j] - aptrb[j], /* Number of non-zeros in column j.*/
                          asub + aptrb[j],   /* Pointer to row indexes of column j.*/
                          aval + aptrb[j]);  /* Pointer to Values of column j.*/

      }

      /* Set the bounds on constraints.
       for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
      for (i = 0; i < numcon && r == MSK_RES_OK; ++i)
        r = MSK_putconbound(task,
                            i,           /* Index of constraint.*/
                            bkc[i],      /* Bound key.*/
                            blc[i],      /* Numerical value of lower bound.*/
                            buc[i]);     /* Numerical value of upper bound.*/

      if ( r == MSK_RES_OK )
      {
        /* Append the first cone. */
        csub[0] = 3;
        csub[1] = 0;
        csub[2] = 1;

        r = MSK_appendcone(task,
                           MSK_CT_QUAD,
                           0.0, /* For future use only, can be set to 0.0 */
                           3,
                           csub);
      }

      if ( r == MSK_RES_OK )
      {
        /* Append the second cone. */
        csub[0] = 4;
        csub[1] = 5;
        csub[2] = 2;

        r = MSK_appendcone(task,
                           MSK_CT_RQUAD,
                           0.0,
                           3,
                           csub);
      }

      if ( r == MSK_RES_OK )
      {
        MSKrescodee trmcode;

        /* Run optimizer */
        r = MSK_optimizetrm(task, &trmcode);


        /* Print a summary containing information
           about the solution for debugging purposes*/
        MSK_solutionsummary (task, MSK_STREAM_MSG);

        if ( r == MSK_RES_OK )
        {
          MSKsolstae solsta;

          MSK_getsolsta (task, MSK_SOL_ITR, &solsta);

          switch (solsta)
          {
            case MSK_SOL_STA_OPTIMAL:
            case MSK_SOL_STA_NEAR_OPTIMAL:
              {
                double *xx = NULL;

                xx = calloc(numvar, sizeof(double));
                if ( xx )
                {
                  MSK_getxx (task,
                             MSK_SOL_ITR,    /* Request the interior solution. */
                             xx);

                  printf("Optimal primal solution\n");
                  for (j = 0; j < numvar; ++j)
                    printf("x[%d]: %e\n", j, xx[j]);
                }
                else
                {
                  r = MSK_RES_ERR_SPACE;
                }
                free(xx);
              }
              break;
            case MSK_SOL_STA_DUAL_INFEAS_CER:
            case MSK_SOL_STA_PRIM_INFEAS_CER:
            case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
            case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined.\n");
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          printf("Error while optimizing.\n");
        }
      }

      if (r != MSK_RES_OK)
      {
        /* In case of an error print error code and description. */
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];

        printf("An error occurred while optimizing.\n");
        MSK_getcodedesc (r,
                         symname,
                         desc);
        printf("Error %s - '%s'\n", symname, desc);
      }
    }
    /* Delete the task and the associated data. */
    MSK_deletetask(&task);
  }

  /* Delete the environment and the associated data. */
  MSK_deleteenv(&env);

  return ( r );
} /* main */