6.2 From Linear to Conic Optimization

In Sec. 6.1 (Linear Optimization) we demonstrated setting up the linear part of an optimization problem, that is the objective, linear bounds, linear equalities and inequalities. In this tutorial we show how to define conic constraints. We recommend going through this general conic tutorial before proceeding to examples with specific cone types.

MOSEK accepts conic constraints in the form

\[Fx+g\in \D\]

where

• \(x\in\real^n\) is the optimization variable,

• \(D\subseteq \real^k\) is a conic domain of some dimension \(k\), representing one of the cone types supported by MOSEK,

• \(F\in\real^{k\times n}\) and \(g\in \real^k\) are data which constitute the sequence of \(k\) affine expressions appearing in the rows of \(Fx+g\).

Constraints of this form will be called affine conic constraints, or ACC for short. Therefore in this section we show how to set up a problem 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, \\ & & & Fx+g & \in & \D_1\times\cdots\times \D_p, \end{array}\end{split}\]

with some number \(p\) of affine conic constraints.

Note that conic constraints are a natural generalization of linear constraints to the general nonlinear case. For example, a typical linear constraint of the form

\[Ax+b\geq 0\]

can be also written as membership in the cone of nonnegative real numbers:

\[Ax+b \in \real_{\geq 0}^d,\]

and that naturally generalizes to

\[Fx+g\in \D\]

for more complicated domains \(\D\) from Sec. 15.9 (Supported domains) of which \(\D=\real_{\geq 0}^d\) is a special case.

6.2.1 Running example

In this tutorial we will consider a sample problem of the form

(6.2)\[\begin{split}\begin{array}{ll} \mbox{maximize} & c^T x \\ \mbox{subject to} & \sum_i x_i = 1, \\ & \gamma \geq \| Gx+h \|_2, \end{array}\end{split}\]

where \(x\in \real^n\) is the optimization variable and \(G\in\real^{k\times n}\), \(h\in\real^k\), \(c\in\real^n\) and \(\gamma\in\real\). We will use the following sample data:

\[\begin{split}n=3,\quad k=2,\quad x\in \real^3, \quad c = [2, 3, -1]^T,\quad \gamma=0.03,\quad G = \left[\begin{array}{ccc}1.5 & 0.1 & 0\\0.3 & 0 & 2.1\end{array}\right],\quad h = \left[\begin{array}{c}0 \\ 0.1\end{array}\right].\end{split}\]

To be explicit, the problem we are going to solve is therefore:

(6.3)\[\begin{split}\begin{array}{ll} \mbox{maximize} & 2x_0+3x_1-x_2 \\ \mbox{subject to} & x_0+x_1+x_2 = 1, \\ & 0.03 \geq \sqrt{(1.5x_0+0.1x_1)^2+(0.3x_0+2.1x_2+0.1)^2}. \end{array}\end{split}\]

Consulting the definition of a quadratic cone \(\Q\) we see that the conic form of this problem is:

(6.4)\[\begin{split}\begin{array}{ll} \mbox{maximize} & 2x_0+3x_1-x_2 \\ \mbox{subject to} & x_0+x_1+x_2 = 1, \\ & (0.03,\ 1.5x_0+0.1x_1,\ 0.3x_0+2.1x_2+0.1) \in \Q^3. \end{array}\end{split}\]

The conic constraint has an affine conic representation \(Fx+g\in\D\) as follows:

(6.5)\[\begin{split}\left[\begin{array}{ccc}0 & 0 & 0\\ 1.5 & 0.1 & 0\\0.3 & 0 & 2.1\end{array}\right]x + \left[\begin{array}{c}0.03\\ 0\\ 0.1\end{array}\right] \in \Q^3.\end{split}\]

Of course by the same logic in the general case the conic form of the problem (6.2) would be

(6.6)\[\begin{split}\begin{array}{ll} \mbox{maximize} & c^T x \\ \mbox{subject to} & \sum_i x_i = 1, \\ & (\gamma, Gx+h)\in\Q^{k+1} \end{array}\end{split}\]

and the ACC representation of the constraint \((\gamma, Gx+h)\in\Q^{k+1}\) would be

\[\begin{split}\left[\begin{array}{c}0\\ G\end{array}\right]x + \left[\begin{array}{c}\gamma\\ h\end{array}\right] \in \Q^{k+1}.\end{split}\]

Now we show how to add the ACC (6.5). This involves three steps:

• storing the affine expressions which appear in the constraint,

• creating a domain, and

• combining the two into an ACC.

6.2.2 Step 1: add affine expressions

To store affine expressions (AFE for short) MOSEK provides a matrix \(\afef\) and a vector \(\afeg\) with the understanding that every row of

\[\afef x + \afeg\]

defines one affine expression. The API functions with infix afe are used to operate on \(\afef\) and \(\afeg\), add rows, add columns, set individual elements, set blocks etc. similarly to the methods for operating on the \(A\) matrix of linear constraints. The storage matrix \(\afef\) is a sparse matrix, therefore only nonzero elements have to be explicitly added.

Remark: the storage \(\afef,\afeg\) may, but does not have to be, equal to the pair \(F,g\) appearing in the expression \(Fx+g\). It is possible to store the AFEs in different order than the order they will be used in \(F,g\), as well as store some expressions only once if they appear multiple times in \(Fx+g\). In this first turorial, however, we will for simplicity store all expressions in the same order we will later use them, so that \((\afef,\afeg)=(F,g)\).

In our example we create only one conic constraint (6.5) with three (in general \(k+1\)) affine expressions

\[\begin{split}\begin{array}{l} 0.03, \\ 1.5x_0+0.1x_1,\\ 0.3x_0 +2.1x_2 +0.1. \end{array}\end{split}\]

Given the previous remark, we initialize the AFE storage as:

(6.7)\[\begin{split}\afef = \left[\begin{array}{ccc}0 & 0 & 0\\ 1.5 & 0.1 & 0\\0.3 & 0 & 2.1\end{array}\right],\quad \afeg=\left[\begin{array}{c}0.03\\ 0\\ 0.1\end{array}\right].\end{split}\]

Initially \(\afef\) and \(\afeg\) are empty (have 0 rows). We construct them as follows. First, we append a number of empty rows:

    // Append empty AFE rows for affine expression storage
    task.append_afes(k + 1)?;

We now have \(\afef\) and \(\afeg\) with 3 rows of zeros and we fill them up to obtain (6.7).

    // G matrix in sparse form
    let Gsubi : &[i64] = &[0, 0, 1, 1];
    let Gsubj : &[i32] = &[0, 1, 0, 2];
    let Gval           = &[1.5, 0.1, 0.3, 2.1];
    // Other data
    let h     = &[0.0, 0.1];
    let gamma = 0.03;

    // Construct F matrix in sparse form
    let Fsubi : Vec<i64> = Gsubi.iter().map(|i| *i+1).collect(); // G will be placed from row number 1 in F
    let Fsubj = Gsubj;
    let Fval  = Gval;

    // Fill in F storage
    task.put_afe_f_entry_list(Fsubi.as_slice(), Fsubj, Fval)?;
    // Fill in g storage
    task.put_afe_g(0, gamma)?;
    task.put_afe_g_slice(1, k+1, h)?;

We have now created the matrices from (6.7). Note that at this point we have not defined any ACC yet. All we did was define some affine expressions and place them in a generic AFE storage facility to be used later.

6.2.3 Step 2: create a domain

Next, we create the domain to which the ACC belongs. Domains are created with functions with infix domain. In the case of (6.5) we need a quadratic cone domain of dimension 3 (in general \(k+1\)), which we create with:

    // Define a conic quadratic domain
    let quadDom = task.append_quadratic_cone_domain(k + 1)?;

The function returns a domain index, which is just the position in the list of all domains (potentially) created for the problem. At this point the domain is just stored in the list of domains, but not yet used for anything.

6.2.4 Step 3: create the actual constraint

We are now in position to create the affine conic constraint. ACCs are created with functions with infix acc. The most basic variant, Task.append_acc will append an affine conic constraint based on the following data:

• the list afeidx of indices of AFEs to be used in the constraint. These are the row numbers in \(\afef,\afeg\) which contain the required affine expressions.

• the index domidx of the domain to which the constraint belongs.

Note that number of AFEs used in afeidx must match the dimension of the domain.

In case of (6.5) we have already arranged \(\afef,\afeg\) in such a way that their (only) three rows contain the three affine expressions we need (in the correct order), and we already defined the quadratic cone domain of matching dimension 3. The ACC is now constructed with the following call:

    // Create the ACC
    task.append_acc(quadDom,    // Domain index
                    (0..k+1).collect::<Vec<i64>>().as_slice(), // Indices of AFE rows [0,...,k]
                    vec![0.0; (k+1) as usize].as_slice())?;       // Ignored

This completes the setup of the affine conic constraint.

6.2.5 Example ACC1

We refer to Sec. 6.1 (Linear Optimization) for instructions how to set up the objective and linear constraint \(x_0+x_1+x_2=1\). All else that remains is to set up the MOSEK environment, task, add variables, call the solver with Task.optimize and retrieve the solution with Task.get_xx. Since our problem contains a nonlinear constraint we fetch the interior-point solution. The full code solving problem (6.3) is shown below.

Listing 6.2 Full code of example ACC1. Click here to download.

extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Streamtype,Solsta,Soltype,Boundkey};

// Define problem data
#[allow(non_upper_case_globals)]
const n : i32 = 3;
#[allow(non_upper_case_globals)]
const k : i64 = 2;

#[allow(non_snake_case)]
fn main() -> Result<(),String> {
    // Create a task
    let mut task = match Task::new() {
        Some(e) => e,
        None => return Err("Failed to create task".to_string()),
    }.with_callbacks();
    // Attach a printer to the task
    task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
    // Create n free variables
    task.append_vars(n)?;

    let x : Vec<i32> = (0..n).collect();
    task.put_var_bound_slice_const(0, n, Boundkey::FR, 0.0, 0.0)?;

    // Set up the objective
    let c = &[2.0, 3.0, -1.0];
    task.put_obj_sense(Objsense::MAXIMIZE)?;
    task.put_c_list(x.as_slice(), c)?;

    // One linear constraint - sum(x) = 1
    task.append_cons(1)?;
    task.put_a_row(0,x.as_slice(), vec![1.0; n as usize].as_slice())?;
    task.put_con_bound(0, Boundkey::FX, 1.0, 1.0)?;

    // Append empty AFE rows for affine expression storage
    task.append_afes(k + 1)?;

    // G matrix in sparse form
    let Gsubi : &[i64] = &[0, 0, 1, 1];
    let Gsubj : &[i32] = &[0, 1, 0, 2];
    let Gval           = &[1.5, 0.1, 0.3, 2.1];
    // Other data
    let h     = &[0.0, 0.1];
    let gamma = 0.03;

    // Construct F matrix in sparse form
    let Fsubi : Vec<i64> = Gsubi.iter().map(|i| *i+1).collect(); // G will be placed from row number 1 in F
    let Fsubj = Gsubj;
    let Fval  = Gval;

    // Fill in F storage
    task.put_afe_f_entry_list(Fsubi.as_slice(), Fsubj, Fval)?;
    // Fill in g storage
    task.put_afe_g(0, gamma)?;
    task.put_afe_g_slice(1, k+1, h)?;

    // Define a conic quadratic domain
    let quadDom = task.append_quadratic_cone_domain(k + 1)?;

    // Create the ACC
    task.append_acc(quadDom,    // Domain index
                    (0..k+1).collect::<Vec<i64>>().as_slice(), // Indices of AFE rows [0,...,k]
                    vec![0.0; (k+1) as usize].as_slice())?;       // Ignored

    // Solve and retrieve solution
    let _ = task.optimize()?;
    task.write_data("acc1.ptf")?;
    let mut xx = vec![0.0; n as usize];
    task.get_xx(Soltype::ITR,xx.as_mut_slice())?;

    assert!(task.get_sol_sta(Soltype::ITR)? == Solsta::OPTIMAL);
    println!("Solution: {:?}",xx);

    // Demonstrate retrieving activity of ACC
    let mut activity = vec![0.0; (k+1) as usize];
    let mut doty     = vec![0.0; (k+1) as usize];

    task.evaluate_acc(Soltype::ITR,0,activity.as_mut_slice())?;
    println!("Activity of ACC:: {:?}",activity);

    // Demonstrate retrieving the dual of ACC
    task.get_acc_dot_y(Soltype::ITR,0,doty.as_mut_slice())?;
    println!("Dual of ACC:: {:?}",doty);


    Ok(())
}

fn maxgap(a : &[f64], b : &[f64]) -> f64 {
    a.iter().zip(b.iter()).map(|(&a,&b)| (a - b).abs()).max_by(|a,b| if a < b { std::cmp::Ordering::Less } else if b < a { std::cmp::Ordering::Greater } else {std::cmp::Ordering::Equal} ).unwrap()
}
fn dot(a : &[f64], b : &[f64]) -> f64 {
    a.iter().zip(b.iter()).map(|(&a,&b)| (a * b)).sum()
}

The answer is

[-0.07838011145615721, 1.1289128998004547, -0.0505327883442975]

The dual values \(\doty\) of an ACC can be obtained with Task.get_acc_dot_y if required.

    // Demonstrate retrieving the dual of ACC
    task.get_acc_dot_y(Soltype::ITR,0,doty.as_mut_slice())?;
    println!("Dual of ACC:: {:?}",doty);

6.2.6 Example ACC2 - more conic constraints

Now that we know how to enter one affine conic constraint (ACC) we will demonstrate a problem with two ACCs. From there it should be clear how to add multiple ACCs. To keep things familiar we will reuse the previous problem, but this time cast it into a conic optimization problem with two ACCs as follows:

(6.8)\[\begin{split}\begin{array}{ll} \mbox{maximize} & c^T x \\ \mbox{subject to} & (\sum_i x_i - 1,\ \gamma,\ Gx+h) \in \{0\} \times \Q^{k+1} \end{array}\end{split}\]

or, using the data from the example:

\[\begin{split}\begin{array}{lll} \mbox{maximize} & 2x_0+3x_1-x_2 & \\ \mbox{subject to} & x_0+x_1+x_2 - 1 & \in \{0\}, \\ & (0.03, 1.5x_0+0.1x_1, 0.3x_0+2.1x_2+0.1) & \in \Q^3 \end{array}\end{split}\]

In other words, we transformed the linear constraint into an ACC with the one-point zero domain.

As before, we proceed in three steps. First, we add the variables and create the storage \(\afef\), \(\afeg\) containing all affine expressions that appear throughout all off the ACCs. It means we will require 4 rows:

(6.9)\[\begin{split}\afef = \left[\begin{array}{ccc}1& 1 & 1 \\ 0 & 0 & 0\\ 1.5 & 0.1 & 0\\0.3 & 0 & 2.1\end{array}\right],\quad \afeg=\left[\begin{array}{c}-1 \\0.03\\ 0\\ 0.1\end{array}\right].\end{split}\]

    // Set AFE rows representing the linear constraint
    task.append_afes(1)?;
    task.put_afe_f_row(0, x.as_slice(), vec![1.0; n as usize].as_slice())?;
    task.put_afe_g(0, -1.0)?;

    // Set AFE rows representing the quadratic constraint
    task.append_afes(k + 1)?;
    task.put_afe_f_row(2,             // afeidx, row number
                    &[0i32, 1],    // varidx, column numbers
                    &[1.5, 0.1])?; // values
    task.put_afe_f_row(3,          // afeidx, row number
                    &[0i32, 2],     // varidx, column numbers
                    &[0.3, 2.1])?; // values

    let h = &[0.0, 0.1];
    let gamma = 0.03;
    task.put_afe_g(1, gamma)?;
    task.put_afe_g_slice(2, k+2, h)?;

Next, we add the required domains: the zero domain of dimension 1, and the quadratic cone domain of dimension 3.

    // Define domains
    let zero_dom = task.append_rzero_domain(1)?;
    let quad_dom = task.append_quadratic_cone_domain(k + 1)?;

Finally, we create both ACCs. The first ACCs picks the 0-th row of \(\afef,\afeg\) and places it in the zero domain:

    task.append_acc(zero_dom,    // Domain index
                    &[0i64],        // Indices of AFE rows
                    &[0.0])?;       // Ignored

The second ACC picks rows \(1,2,3\) in \(\afef,\afeg\) and places them in the quadratic cone domain:

    task.append_acc(quad_dom,    // Domain index
                   &[1i64,2,3],    // Indices of AFE rows
                   &[0.0,0.0,0.0])?; // Ignored

The completes the construction and we can solve the problem like before:

Listing 6.3 Full code of example ACC2. Click here to download.

extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Streamtype,Solsta,Soltype,Boundkey};

#[allow(non_upper_case_globals)]
fn main() -> Result<(),String> {
    // Define problem data
    const n : i32 = 3;
    const k : i64 = 2;

    // Create a task
    let mut task = match Task::new() {
        Some(e) => e,
        None => return Err("Failed to create task".to_string()),
    }.with_callbacks();
    // Attach a printer to the task
    task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
    // Create n free variables
    task.append_vars(n)?;
    let x : Vec<i32> = (0..n).collect();

    task.put_var_bound_slice_const(0, n, Boundkey::FR, 0.0,0.0)?;

    // Set up the objective
    let c = &[2.0, 3.0, -1.0];
    task.put_obj_sense(Objsense::MAXIMIZE)?;
    task.put_c_list(x.as_slice(), c)?;

    // Set AFE rows representing the linear constraint
    task.append_afes(1)?;
    task.put_afe_f_row(0, x.as_slice(), vec![1.0; n as usize].as_slice())?;
    task.put_afe_g(0, -1.0)?;

    // Set AFE rows representing the quadratic constraint
    task.append_afes(k + 1)?;
    task.put_afe_f_row(2,             // afeidx, row number
                    &[0i32, 1],    // varidx, column numbers
                    &[1.5, 0.1])?; // values
    task.put_afe_f_row(3,          // afeidx, row number
                    &[0i32, 2],     // varidx, column numbers
                    &[0.3, 2.1])?; // values

    let h = &[0.0, 0.1];
    let gamma = 0.03;
    task.put_afe_g(1, gamma)?;
    task.put_afe_g_slice(2, k+2, h)?;

    // Define domains
    let zero_dom = task.append_rzero_domain(1)?;
    let quad_dom = task.append_quadratic_cone_domain(k + 1)?;

    // Append affine conic constraints
    task.append_acc(zero_dom,    // Domain index
                    &[0i64],        // Indices of AFE rows
                    &[0.0])?;       // Ignored
    task.append_acc(quad_dom,    // Domain index
                   &[1i64,2,3],    // Indices of AFE rows
                   &[0.0,0.0,0.0])?; // Ignored

    // Solve and retrieve solution
    let _ = task.optimize()?;
    task.write_data("acc2.ptf")?;
    let mut xx = vec![0.0; n as usize];
    task.get_xx(Soltype::ITR,xx.as_mut_slice())?;
    assert! (task.get_sol_sta(Soltype::ITR)? == Solsta::OPTIMAL);
    println!("Solution: {:?}",xx);

    // Demonstrate retrieving activity of ACC
    let mut activity = vec![0.0; 3];
    let mut doty     = vec![0.0; 3];
    task.evaluate_acc(Soltype::ITR,1,activity.as_mut_slice())?;
    println!("Activity of quadratic ACC:: {:?}",activity);

    // Demonstrate retrieving the dual of ACC
    task.get_acc_dot_y(Soltype::ITR,1,doty.as_mut_slice())?;
    println!("Dual of quadratic ACC:: {:?}",doty);


    Ok(())
}

fn maxgap(a : &[f64], b : &[f64]) -> f64 {
    a.iter().zip(b.iter()).map(|(&a,&b)| (a - b).abs()).max_by(|a,b| if a < b { std::cmp::Ordering::Less } else if b < a { std::cmp::Ordering::Greater } else {std::cmp::Ordering::Equal} ).unwrap()
}
fn dot(a : &[f64], b : &[f64]) -> f64 {
    a.iter().zip(b.iter()).map(|(&a,&b)| (a * b)).sum()
}

We obtain the same result:

[-0.07838011145615721, 1.1289128998004547, -0.0505327883442975]

6.2.7 Summary and extensions

In this section we presented the most basic usage of the affine expression storage \(\afef,\afeg\) to input affine expressions used together with domains to create affine conic constraints. Now we briefly point out additional features of his interface which can be useful in some situations for more demanding users. They will be demonstrated in various examples in other tutorials and case studies in this manual.

• It is important to remember that \(\afef,\afeg\) has only a storage function and during the ACC construction we can pick an arbitrary list of row indices and place them in a conic domain. It means for example that:

  • It is not necessary to store the AFEs in the same order they will appear in ACCs.

  • The same AFE index can appear more than once in one and/or more conic constraints (this can be used to reduce storage if the same affine expression is used in multiple ACCs).

  • The \(\afef,\afeg\) storage can even include rows that are not presently used in any ACC.

• Domains can be reused: multiple ACCs can use the same domain. On the other hand the same type of domain can appear under many domidx positions. In this sense the list of created domains also plays only a storage role: the domains are only used when they enter an ACC.

• Affine expressions can also contain semidefinite terms, ie. the most general form of an ACC is in fact

  \[Fx + \langle \bar{F},\barX\rangle + g \in \D\]

  These terms are input into the rows of AFE storage using the functions with infix afebarf, creating an additional storage structure \(\bar{\afef}\).

• The same affine expression storage \(\afef,\afeg\) is shared between affine conic and disjunctive constraints (see Sec. 6.9 (Disjunctive constraints)).

• If, on the other hand, the user chooses to always store the AFEs one by one sequentially in the same order as they appear in ACCs then sequential functions such as Task.append_acc_seq and Task.append_accs_seq make it easy to input one or more ACCs by just specifying the starting AFE index and dimension.

• It is possible to add a number of ACCs in one go using Task.append_accs.

• When defining an ACC an additional constant vector \(b\) can be provided to modify the constant terms coming from \(\afeg\) but only for this particular ACC. This could be useful to reduce \(\afef\) storage space if, for example, many expressions \(f^Tx-b_i\) with the same linear part \(f^Tx\), but varying constant terms \(b_i\), are to be used throughout ACCs.