6.1 Linear Optimization

The simplest optimization problem is a purely linear problem. A linear optimization problem (see also Sec. 12.1 (Linear Optimization)) is a problem of the following form:

Minimize or maximize the objective function

\[\sum_{j=0}^{n-1} c_j x_j + c^f\]

subject to the linear constraints

\[l_k^c \leq \sum_{j=0}^{n-1} a_{kj} x_j \leq u_k^c,\quad k=0,\ldots ,m-1,\]

and the bounds

\[l_j^x \leq x_j \leq u_j^x, \quad j=0,\ldots ,n-1.\]

The problem description consists of the following elements:

  • \(m\) and \(n\) — the number of constraints and variables, respectively,

  • \(x\) — the variable vector of length \(n\),

  • \(c\) — the coefficient vector of length \(n\)

    \[\begin{split}c = \left[ \begin{array}{c} c_0 \\ \vdots \\ c_{n-1} \end{array} \right],\end{split}\]
  • \(c^f\) — fixed term in the objective,

  • \(A\) — an \(m\times n\) matrix of coefficients

    \[\begin{split}A = \left[ \begin{array}{ccc} a_{0,0} & \cdots & a_{0,(n-1)} \\ \vdots & \cdots & \vdots \\ a_{(m-1),0} & \cdots & a_{(m-1),(n-1)} \end{array} \right],\end{split}\]
  • \(l^c\) and \(u^c\) — the lower and upper bounds on constraints,

  • \(l^x\) and \(u^x\) — the lower and upper bounds on variables.

Please note that we are using \(0\) as the first index: \(x_0\) is the first element in variable vector \(x\).

6.1.1 Example LO1

The following is an example of a small linear optimization problem:

(6.1)\[\begin{split}\begin{array} {lccccccccl} \mbox{maximize} & 3 x_0 & + & 1 x_1 & + & 5 x_2 & + & 1 x_3 & & \\ \mbox{subject to} & 3 x_0 & + & 1 x_1 & + & 2 x_2 & & & = & 30, \\ & 2 x_0 & + & 1 x_1 & + & 3 x_2 & + & 1 x_3 & \geq & 15, \\ & & & 2 x_1 & & & + & 3 x_3 & \leq & 25, \end{array}\end{split}\]

under the bounds

\[\begin{split}\begin{array}{ccccc} 0 & \leq & x_0 & \leq & \infty , \\ 0 & \leq & x_1 & \leq & 10, \\ 0 & \leq & x_2 & \leq & \infty ,\\ 0 & \leq & x_3 & \leq & \infty . \end{array}\end{split}\]

Solving the problem

To solve the problem above we go through the following steps:

  1. (Optionally) Creating an environment.

  2. Creating an optimization task.

  3. Loading a problem into the task object.

  4. Optimization.

  5. Extracting the solution.

Below we explain each of these steps.

Creating an environment.

The user can start by creating a MOSEK environment, but it is not necessary if the user does not need access to other functionalities, license management, additional routines, etc. Therefore in this tutorial we don’t create an explicit environment.

Creating an optimization task.

We create an empty task object. A task object represents all the data (inputs, outputs, parameters, information items etc.) associated with one optimization problem.

    let mut task = match Task::new() {
        Some(e) => e,
        None => return Err("Failed to create task".to_string()),
    }.with_callbacks();

    /* Directs the log task stream to the 'printstr' function. */
    task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;

We also connect a call-back function to the task log stream. Messages related to the task are passed to the call-back function. In this case the stream call-back function writes its messages to the standard output stream. See Sec. 7.4 (Input/Output).

Loading a problem into the task object.

Before any problem data can be set, variables and constraints must be added to the problem via calls to the functions Task.append_cons and Task.append_vars.

    /* Append 'numcon' empty constraints.
     * The constraints will initially have no bounds. */
    task.append_cons(numcon as i32)?;

    /* Append 'numvar' variables.
     * The variables will initially be fixed at zero (x=0). */
    task.append_vars(numvar as i32)?;

New variables can now be referenced from other functions with indexes in \(\idxbeg, \ldots, \idxend{\mathtt{numvar}}\) and new constraints can be referenced with indexes in \(\idxbeg, \ldots , \idxend{\mathtt{numcon}}\). More variables and/or constraints can be appended later as needed, these will be assigned indexes from \(\mathtt{numvar}\)/\(\mathtt{numcon}\) and up. Optionally one can add names.

Setting the objective.

Next step is to set the problem data. We first set the objective coefficients \(c_j = \mathtt{c[j]}\). This can be done with functions such as Task.put_c_j or Task.put_c_list.

        task.put_c_j(j as i32,c[j])?;

Setting bounds on variables

For every variable we need to specify a bound key and two bounds according to Table 6.1.

Table 6.1 Bound keys as defined in the enum Boundkey.

Bound key

Type of bound

Lower bound

Upper bound

Boundkey::FX

\(u_j = l_j\)

Finite

Identical to the lower bound

Boundkey::FR

Free

\(-\infty\)

\(+\infty\)

Boundkey::LO

\(l_j \leq \cdots\)

Finite

\(+\infty\)

Boundkey::RA

\(l_j \leq \cdots \leq u_j\)

Finite

Finite

Boundkey::UP

\(\cdots \leq u_j\)

\(-\infty\)

Finite

For instance bkx[0]= Boundkey::LO means that \(x_0 \geq l_0^x\). Finally, the numerical values of the bounds on variables are given by

\[l_j^x = \mathtt{blx[j]}\]

and

\[u_j^x = \mathtt{bux[j]}.\]

Let us assume we have the bounds on variables stored in the arrays

    let bkx = vec![ Boundkey::LO, Boundkey::RA, Boundkey::LO, Boundkey::LO ];
    let blx = vec![ 0.0,       0.0,       0.0,       0.0       ];
    let bux = vec![ INF,      10.0,       INF,       INF       ];

Then we can set them using various functions such Task.put_var_bound, Task.put_var_bound_slice, Task.put_var_bound_list, depending on what is most convenient in the given context. For instance:

        task.put_var_bound(j as i32,    /* Index of variable.*/
                           bkx[j],      /* Bound key.*/
                           blx[j],      /* Numerical value of lower bound.*/
                           bux[j])?;     /* Numerical value of upper bound.*/

Defining the linear constraint matrix.

Recall that in our example the \(A\) matrix is given by

\[\begin{split}A = \left[ \begin{array}{cccc} 3 & 1 & 2 & 0 \\ 2 & 1 & 3 & 1 \\ 0 & 2 & 0 & 3 \end{array} \right].\end{split}\]

This matrix is stored in sparse format:

    let aptrb = vec![ 0, 2, 5, 7 ];
    let aptre = vec![ 2, 5, 7, 9 ];
    let asub  = vec![ 0, 1,
                      0, 1, 2,
                      0, 1,
                      1, 2 ];
    let aval  = vec![ 3.0, 2.0,
                      1.0, 1.0, 2.0,
                      2.0, 3.0,
                      1.0, 3.0 ];

The array aval[j] contains the non-zero values of column \(j\) and asub[j] contains the row indices of these non-zeros.

We now input the linear constraint matrix into the task. This can be done in many alternative ways, row-wise, column-wise or element by element in various orders. See functions such as Task.put_a_row, Task.put_a_row_list, Task.put_aij_list, Task.put_a_col and similar.

        task.put_a_col(j as i32,          /* Variable (column) index.*/
                       & asub[aptrb[j]..aptre[j]],     /* Pointer to row indexes of column j.*/
                       & aval[aptrb[j]..aptre[j]])?;    /* Pointer to Values of column j.*/

Setting bounds on constraints

Finally, the bounds on each constraint are set similarly to the variable bounds, using the bound keys as in Table 6.1. This can be done with one of the many functions Task.put_con_bound, Task.put_con_bound_slice, Task.put_con_bound_list, depending on the situation.

    for i in 0..numcon {
      task.put_con_bound(i as i32,    /* Index of constraint.*/
                         bkc[i],      /* Bound key.*/
                         blc[i],      /* Numerical value of lower bound.*/
                         buc[i])?;     /* Numerical value of upper bound.*/
    }

Optimization

After the problem is set-up the task can be optimized by calling the function Task.optimize.

    let _trmcode = task.optimize()?;

Extracting the solution.

After optimizing the status of the solution is examined with a call to Task.get_sol_sta.

    let solsta = task.get_sol_sta(Soltype::BAS)?;

If the solution status is reported as Solsta::OPTIMAL the solution is extracted:

            task.get_xx(Soltype::BAS,    /* Request the basic solution. */
                        & mut xx[..])?;

The Task.get_xx function obtains the solution. MOSEK may compute several solutions depending on the optimizer employed. In this example the basic solution is requested by setting the first argument to Soltype::BAS. For details about fetching solutions see Sec. 7.2 (Accessing the solution).

Source code

The complete source code lo1.rs of this example appears below. See also lo2.rs for a version where the \(A\) matrix is entered row-wise.

Listing 6.1 Linear optimization example. Click here to download.
extern crate mosek;

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

const INF : f64 = 0.0;

fn main() -> Result<(),String> {
    let numvar = 4;
    let numcon = 3;

    let c = vec![3.0, 1.0, 5.0, 1.0];

    /* Below is the sparse representation of the A
     * matrix stored by column. */
    let aptrb = vec![ 0, 2, 5, 7 ];
    let aptre = vec![ 2, 5, 7, 9 ];
    let asub  = vec![ 0, 1,
                      0, 1, 2,
                      0, 1,
                      1, 2 ];
    let aval  = vec![ 3.0, 2.0,
                      1.0, 1.0, 2.0,
                      2.0, 3.0,
                      1.0, 3.0 ];

    /* Bounds on constraints. */
    let bkc = vec![ Boundkey::FX, Boundkey::LO, Boundkey::UP ];
    let blc = vec![ 30.0,      15.0,      -INF      ];
    let buc = vec![ 30.0,      INF,       25.0      ];
    /* Bounds on variables. */
    let bkx = vec![ Boundkey::LO, Boundkey::RA, Boundkey::LO, Boundkey::LO ];
    let blx = vec![ 0.0,       0.0,       0.0,       0.0       ];
    let bux = vec![ INF,      10.0,       INF,       INF       ];

    /* Create the optimization task. */
    let mut task = match Task::new() {
        Some(e) => e,
        None => return Err("Failed to create task".to_string()),
    }.with_callbacks();

    /* Directs the log task stream to the 'printstr' function. */
    task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;

    /* Append 'numcon' empty constraints.
     * The constraints will initially have no bounds. */
    task.append_cons(numcon as i32)?;

    /* Append 'numvar' variables.
     * The variables will initially be fixed at zero (x=0). */
    task.append_vars(numvar as i32)?;

    for j in 0..numvar
    {
        /* Set the linear term c_j in the objective.*/
        task.put_c_j(j as i32,c[j])?;

        /* Set the bounds on variable j.
         * blx[j] <= x_j <= bux[j] */
        task.put_var_bound(j as i32,    /* 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 */
        task.put_a_col(j as i32,          /* Variable (column) index.*/
                       & asub[aptrb[j]..aptre[j]],     /* Pointer to row indexes of column j.*/
                       & aval[aptrb[j]..aptre[j]])?;    /* Pointer to Values of column j.*/
    }

    /* Set the bounds on constraints.
     * for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
    for i in 0..numcon {
      task.put_con_bound(i as i32,    /* Index of constraint.*/
                         bkc[i],      /* Bound key.*/
                         blc[i],      /* Numerical value of lower bound.*/
                         buc[i])?;     /* Numerical value of upper bound.*/
    }

    /* Maximize objective function. */
    task.put_obj_sense(Objsense::MAXIMIZE)?;

    /* Run optimizer */
    let _trmcode = task.optimize()?;

    /* Print a summary containing information
     * about the solution for debugging purposes. */

    task.solution_summary(Streamtype::LOG)?;

    let solsta = task.get_sol_sta(Soltype::BAS)?;

    match solsta
    {
        Solsta::OPTIMAL =>
        {
            let mut xx = vec![0.0,0.0,0.0,0.0];
            task.get_xx(Soltype::BAS,    /* Request the basic solution. */
                        & mut xx[..])?;
            println!("Optimal primal solution");
            for j in 0..numvar as usize
            {
                println!("x[{}]: {}",j,xx[j]);
            }
          }

        Solsta::DUAL_INFEAS_CER |
        Solsta::PRIM_INFEAS_CER =>
        {
            println!("Primal or dual infeasibility certificate found.");
        }

        Solsta::UNKNOWN =>
        {
            /* If the solutions status is unknown, print the termination code
             * indicating why the optimizer terminated prematurely. */

            println!("The solution status is unknown.");
            println!("The optimizer terminitated with code: {}",solsta);
          }
        _ =>
        {
            println!("Other solution status.");
        }
    }
    Ok(())

}