6.3 Conic Quadratic Optimization

The structure of a typical conic optimization problem is

\[\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, \\ & & & Fx+g & \in & \D, \end{array}\end{split}\]

(see Sec. 12 (Problem Formulation and Solutions) for detailed formulations). We recommend Sec. 6.2 (From Linear to Conic Optimization) for a tutorial on how problems of that form are represented in MOSEK and what data structures are relevant. Here we discuss how to set-up problems with the (rotated) quadratic cones.

MOSEK supports two types of quadratic cones, namely:

  • 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, consider the following constraint:

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

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

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

For other types of cones supported by MOSEK, see Sec. 15.10 (Supported domains) and the other tutorials in this chapter. Different cone types can appear together in one optimization problem.

6.3.1 Example CQO1

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

(6.10)\[\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}\]

The two conic constraints can be expressed in the ACC form as shown in (6.11)

(6.11)\[\begin{split}\left[\begin{array}{cccccc}0&0&0&1&0&0\\1&0&0&0&0&0\\0&1&0&0&0&0\\0&0&0&0&1&0\\0&0&0&0&0&1\\0&0&1&0&0&0\end{array}\right] \left[\begin{array}{c}x_1\\x_2\\x_3\\x_4\\x_5\\x_6\end{array}\right] + \left[\begin{array}{c}0\\0\\0\\0\\0\\0\end{array}\right] \in \Q^3 \times \Q_r^3.\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 Sec. 6.1 (Linear Optimization) 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

In order to append the conic constraints we first input the matrix \(\afef\) and vector \(\afeg\) appearing in (6.11). The matrix \(\afef\) is sparse and we input only its nonzeros using MSK_putafefentrylist. Since \(\afeg\) is zero, nothing needs to be done about this vector.

Each of the conic constraints is appended using the function MSK_appendacc. In the first case we append the quadratic cone determined by the first three rows of \(\afef\) and then the rotated quadratic cone depending on the remaining three rows of \(\afef\).

        /* Set the non-zero entries of the F matrix */
        r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val); 

        /* Append quadratic cone domain */
        if (r == MSK_RES_OK)
          r = MSK_appendquadraticconedomain(task, 3, domidx);
        /* Append rotated quadratic cone domain */        
        if (r == MSK_RES_OK)
          r = MSK_appendrquadraticconedomain(task, 3, domidx+1);

        /* Append two ACCs made up of the AFEs and the domains defined above. */
        if (r == MSK_RES_OK)
          r = MSK_appendaccsseq(task, numacc, domidx, numafe, afeidx[0], NULL);

The first argument selects the domain, which must be appended before being used, and must have the dimension matching the number of affine expressions appearing in the constraint. Variants of this method are available to append multiple ACCs at a time. It is also possible to define the matrix \(\afef\) using a variety of methods (row after row, column by column, individual entries, etc.) similarly as for the linear constraint matrix \(A\).

For a more thorough exposition of the affine expression storage (AFE) matrix \(\afef\) and vector \(\afeg\) see Sec. 6.2 (From Linear to Conic Optimization).

Source code

Listing 6.4 Source code solving problem (6.10). 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;
  const MSKint64t numafe = 6,
                  numacc = 2,
                  f_nnz  = 6;

  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};

  MSKint64t   afeidx[] = {0, 1, 2, 3, 4, 5};
  MSKint32t   varidx[] = {3, 0, 1, 4, 5, 2};
  MSKrealt    f_val[]  = {1, 1, 1, 1, 1, 1};

  MSKint64t   domidx[] = {0, 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);

      /* Append 'numafe' affine expressions.
      The affine expressions will initially be empty. */
      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, numafe);


      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)
      {
        /* Set the non-zero entries of the F matrix */
        r = MSK_putafefentrylist(task, f_nnz, afeidx, varidx, f_val); 

        /* Append quadratic cone domain */
        if (r == MSK_RES_OK)
          r = MSK_appendquadraticconedomain(task, 3, domidx);
        /* Append rotated quadratic cone domain */        
        if (r == MSK_RES_OK)
          r = MSK_appendrquadraticconedomain(task, 3, domidx+1);

        /* Append two ACCs made up of the AFEs and the domains defined above. */
        if (r == MSK_RES_OK)
          r = MSK_appendaccsseq(task, numacc, domidx, numafe, afeidx[0], NULL);
      }

      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:
              {
                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:
              printf("Primal or dual infeasibility certificate found.\n");
              break;
            case MSK_SOL_STA_UNKNOWN:
              printf("The status of the solution could not be determined. Termination code: %d.\n", trmcode);
              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 */