7.3 Power Cone Optimization

The structure of a typical conic optimization problem is

$$\begin{array}{r}\begin{array}{lccccl}\text{minimize}& & & c^{T}x+c^f& & \\ \text{subject to}& l^c& \le & Ax& \le & u^c,\\ & l^x& \le & x& \le & u^x,\\ & & & Fx+g& \in & \mathcal{D},\end{array}\end{array}$$

(see Sec. 12 (Problem Formulation and Solutions) for detailed formulations). Here we discuss how to set-up problems with the primal/dual power cones.

MOSEK supports the primal and dual power cones, defined as below:

• Primal power cone:

  $${\mathcal{P}}_n^{{\alpha }_k}=\{x\in {\mathbb{R}}^n:\prod _{i=0}^{n_{\ell }-1}{x}_i^{{\beta }_i}\ge \sqrt{\sum _{j=n_{\ell }}^{n-1}{x}_j^2},x_0\dots ,x_{n_{\ell }-1}\ge 0\}$$

  where $s=\sum _i{\alpha }_i$ and ${\beta }_i={\alpha }_i/s$, so that $\sum _i{\beta }_i=1$.

• Dual power cone:

  $$({\mathcal{P}}_n^{{\alpha }_k})=\{x\in {\mathbb{R}}^n:\prod _{i=0}^{n_{\ell }-1}{\left(\frac{x_i}{{\beta }_i}\right)}^{{\beta }_i}\ge \sqrt{\sum _{j=n_{\ell }}^{n-1}{x}_j^2},x_0\dots ,x_{n_{\ell }-1}\ge 0\}$$

  where $s=\sum _i{\alpha }_i$ and ${\beta }_i={\alpha }_i/s$, so that $\sum _i{\beta }_i=1$.

Perhaps the most important special case is the three-dimensional power cone family:

$${\mathcal{P}}_3^{\alpha ,1-\alpha }=\{x\in {\mathbb{R}}^3:x_0^{\alpha }{x}_1^{1-\alpha }\ge |x_2|,x_0,x_1\ge 0\}.$$

which has the corresponding dual cone:

For example, the conic constraint $(x,y,z)\in {\mathcal{P}}_3^{0.25,0.75}$ is equivalent to $x^{0.25}{y}^{0.75}\ge |z|$, or simply $x{y}^3\ge z^4$ with $x,y\ge 0$.

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.3.1 Example POW1

Consider the following optimization problem which involves powers of variables:

(7.3)$$\begin{array}{r}\begin{array}{lrcl}\text{maximize}& x_0^{0.2}{x}_1^{0.8}+x_2^{0.4}-x_0& & \\ \text{subject to}& x_0+x_1+\frac{1}{2}{x}_2& =& 2,\\ & x_0,x_1,x_2& \ge & 0.\end{array}\end{array}$$

We convert (7.3) into affine conic form using auxiliary variables as bounds for the power expressions:

(7.4)$$\begin{array}{r}\begin{array}{lrcl}\text{maximize}& x_3+x_4-x_0& & \\ \text{subject to}& x_0+x_1+\frac{1}{2}{x}_2& =& 2,\\ & (x_0,x_1,x_3)& \in & {\mathcal{P}}_3^{0.2,0.8},\\ & (x_2,1.0,x_4)& \in & {\mathcal{P}}_3^{0.4,0.6}.\end{array}\end{array}$$

The two conic constraints shown in (7.4) can be expressed in the ACC form as shown in (7.5):

(7.5)$$\begin{array}{r}\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ x_4\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 1\\ 0\end{array}\right]\in {\mathcal{P}}_3^{0.2,0.8}\times {\mathcal{P}}_3^{0.4,0.6}.\end{array}$$

We start by creating the optimization model:

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

We then define the variable x corresponding to the original problem (7.3), and auxiliary variables appearing in the conic reformulation (7.4).

  Variable::t x  = M->variable("x", 3, Domain::unbounded());
  Variable::t x3 = M->variable();
  Variable::t x4 = M->variable();

The linear constraint is defined using the dot product operator Expr.dot:

  // Create the linear constraint
  auto aval = new_array_ptr<double, 1>({1.0, 1.0, 0.5});
  M->constraint(Expr::dot(x, aval), Domain::equalsTo(2.0));

The primal power cone is referred to via Domain.inPPowerCone with an appropriate list of variables or expressions in each case.

  // Create the conic constraints
  M->constraint(Var::vstack(x->slice(0,2), x3), Domain::inPPowerCone(0.2));
  M->constraint(Expr::vstack(x->index(2), 1.0, x4), Domain::inPPowerCone(0.4));     

We only need the objective function:

  auto cval = new_array_ptr<double, 1>({1.0, 1.0, -1.0});
  M->objective(ObjectiveSense::Maximize, Expr::dot(cval, Var::vstack(x3, x4, x->index(0))));

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. Here we just display the primal solution

  ndarray<double, 1> xlvl   = *(x->level());
  std::cout << "x,y,z = " << xlvl << std::endl;

which is

[ 0.06389298  0.78308564  2.30604283 ]

Listing 7.3 Fusion implementation of model (7.3). 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("pow1"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x  = M->variable("x", 3, Domain::unbounded());
  Variable::t x3 = M->variable();
  Variable::t x4 = M->variable();

  // Create the linear constraint
  auto aval = new_array_ptr<double, 1>({1.0, 1.0, 0.5});
  M->constraint(Expr::dot(x, aval), Domain::equalsTo(2.0));

  // Create the conic constraints
  M->constraint(Var::vstack(x->slice(0,2), x3), Domain::inPPowerCone(0.2));
  M->constraint(Expr::vstack(x->index(2), 1.0, x4), Domain::inPPowerCone(0.4));     

  auto cval = new_array_ptr<double, 1>({1.0, 1.0, -1.0});
  M->objective(ObjectiveSense::Maximize, Expr::dot(cval, Var::vstack(x3, x4, x->index(0))));

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

  // Get the linear solution values
  ndarray<double, 1> xlvl   = *(x->level());
  std::cout << "x,y,z = " << xlvl << std::endl;
}