7.2 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). 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. 14.8 (Supported domains) and the other tutorials in this chapter. Different cone types can appear together in one optimization problem.

7.2.1 Example CQO1

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

(7.2)\[\begin{split} \begin{array}{ll} \minimize & y_1 + y_2 + y_3 \\ \st & x_1 + x_2 + 2.0 x_3 = 1.0,\\ & x_1,x_2,x_3 \geq 0.0,\\ & (y_1,x_1,x_2) \in \Q^3,\\ & (y_2,y_3,x_3) \in \Qr^3. \end{array}\end{split}\]

We start by creating the optimization model:

  Model::t M = new Model("cqo1"); auto _M = finally([&]() { M->dispose(); });

We then define variables x and y. Two logical variables (aliases) z1 and z2 are introduced to model the quadratic cones. These are not new variables, but map onto parts of x and y for the sake of convenience.

  Variable::t x  = M->variable("x", 3, Domain::greaterThan(0.0));
  Variable::t y  = M->variable("y", 3, Domain::unbounded());

  // Create the aliases
  //      z1 = [ y[0],x[0],x[1] ]
  //  and z2 = [ y[1],y[2],x[2] ]
  Variable::t z1 = Var::vstack(y->index(0),  x->slice(0, 2));
  Variable::t z2 = Var::vstack(y->slice(1, 3), x->index(2));

The linear constraint is defined using the dot product:

  // Create the constraint
  //      x[0] + x[1] + 2.0 x[2] = 1.0
  auto aval = new_array_ptr<double, 1>({1.0, 1.0, 2.0});
  M->constraint("lc", Expr::dot(aval, x), Domain::equalsTo(1.0));

The conic constraints are defined using the logical views z1 and z2 created previously. Note that this is a basic way of defining conic constraints, and that in practice they would have more complicated structure.

  // Create the constraints
  //      z1 belongs to C_3
  //      z2 belongs to K_3
  // where C_3 and K_3 are respectively the quadratic and
  // rotated quadratic cone of size 3, i.e.
  //                 z1[0] >= sqrt(z1[1]^2 + z1[2]^2)
  //  and  2.0 z2[0] z2[1] >= z2[2]^2
  Constraint::t qc1 = M->constraint("qc1", z1, Domain::inQCone());
  Constraint::t qc2 = M->constraint("qc2", z2, Domain::inRotatedQCone());

We only need the objective function:

  // Set the objective function to (y[0] + y[1] + y[2])
  M->objective("obj", ObjectiveSense::Minimize, Expr::sum(y));

Calling the Model.solve method invokes the solver:

  M->solve();

The primal and dual solution values can be retrieved using Variable.level, Constraint.level and Variable.dual, Constraint.dual, respectively:

  // Get the linear solution values
  ndarray<double, 1> xlvl   = *(x->level());
  ndarray<double, 1> ylvl   = *(y->level());
  // Get conic solution of qc1
  ndarray<double, 1> qc1lvl = *(qc1->level());
  ndarray<double, 1> qc1dl  = *(qc1->dual());
Listing 7.2 Fusion implementation of model (7.2). Click here to download.
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model("cqo1"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x  = M->variable("x", 3, Domain::greaterThan(0.0));
  Variable::t y  = M->variable("y", 3, Domain::unbounded());

  // Create the aliases
  //      z1 = [ y[0],x[0],x[1] ]
  //  and z2 = [ y[1],y[2],x[2] ]
  Variable::t z1 = Var::vstack(y->index(0),  x->slice(0, 2));
  Variable::t z2 = Var::vstack(y->slice(1, 3), x->index(2));

  // Create the constraint
  //      x[0] + x[1] + 2.0 x[2] = 1.0
  auto aval = new_array_ptr<double, 1>({1.0, 1.0, 2.0});
  M->constraint("lc", Expr::dot(aval, x), Domain::equalsTo(1.0));

  // Create the constraints
  //      z1 belongs to C_3
  //      z2 belongs to K_3
  // where C_3 and K_3 are respectively the quadratic and
  // rotated quadratic cone of size 3, i.e.
  //                 z1[0] >= sqrt(z1[1]^2 + z1[2]^2)
  //  and  2.0 z2[0] z2[1] >= z2[2]^2
  Constraint::t qc1 = M->constraint("qc1", z1, Domain::inQCone());
  Constraint::t qc2 = M->constraint("qc2", z2, Domain::inRotatedQCone());

  // Set the objective function to (y[0] + y[1] + y[2])
  M->objective("obj", ObjectiveSense::Minimize, Expr::sum(y));

  // Solve the problem
  M->solve();

  // Get the linear solution values
  ndarray<double, 1> xlvl   = *(x->level());
  ndarray<double, 1> ylvl   = *(y->level());
  // Get conic solution of qc1
  ndarray<double, 1> qc1lvl = *(qc1->level());
  ndarray<double, 1> qc1dl  = *(qc1->dual());

  std::cout << "x1,x2,x2 = " << xlvl << std::endl;
  std::cout << "y1,y2,y3 = " << ylvl << std::endl;
  std::cout << "qc1 levels = " << qc1lvl << std::endl;
  std::cout << "qc1 dual conic var levels = " << qc1dl << std::endl;
}