6.1 Linear Optimization

The simplest optimization problem is a purely linear problem. A linear optimization problem 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:

(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. Create an environment.
  2. Create an optimization task.
  3. Load a problem into the task object.
  4. Optimization.
  5. Extracting the solution.

Below we explain each of these steps.

Create an environment.

Before setting up the optimization problem, a MOSEK environment must be created. All tasks in the program should share the same environment.

  r = MSK_makeenv(&env, NULL);

Create an optimization task.

Next, an empty task object is created:

    /* Create the optimization task. */
    r = MSK_maketask(env, numcon, numvar, &task);

    /* Directs the log task stream to the 'printstr' function. */
    if ( r == MSK_RES_OK )
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

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.

Load 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 MSK_appendcons and MSK_appendvars.

    /* 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);

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

Next step is to set the problem data. We loop over each variable index \(j=0, \ldots, \mathtt{numvar}-1\) calling functions to set problem data. We first set the objective coefficient \(c_j = \mathtt{c[j]}\) by calling the function MSK_putcj.

      /* Set the linear term c_j in the objective.*/
      if (r == MSK_RES_OK)
        r = MSK_putcj(task, j, c[j]);

Setting bounds on variables

The bounds on variables are stored in the arrays

  const MSKboundkeye bkx[]  = {MSK_BK_LO,     MSK_BK_RA, MSK_BK_LO,     MSK_BK_LO     };
  const double       blx[]  = {0.0,           0.0,       0.0,           0.0           };
  const double       bux[]  = { +MSK_INFINITY, 10.0,      +MSK_INFINITY, +MSK_INFINITY };

and are set with calls to MSK_putvarbound.

      /* 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.*/

The Bound key stored in bkx specifies the type of the bound according to Table 6.1.

Table 6.1 Bound keys as defined in the enum MSKboundkeye.
Bound key Type of bound Lower bound Upper bound
MSK_BK_FX \(u_j = l_j\) Finite Identical to the lower bound
MSK_BK_FR Free \(-\infty\) \(+\infty\)
MSK_BK_LO \(l_j \leq \cdots\) Finite \(+\infty\)
MSK_BK_RA \(l_j \leq \cdots \leq u_j\) Finite Finite
MSK_BK_UP \(\cdots \leq u_j\) \(-\infty\) Finite

For instance bkx[0]= MSK_BK_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]}.\]

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 in the arrays:

  const MSKint32t    aptrb[] = {0, 2, 5, 7},
                     aptre[] = {2, 5, 7, 9},
                     asub[]  = { 0, 1,
                                 0, 1, 2,
                                 0, 1,
                                 1, 2
                               };
  const double       aval[]  = { 3.0, 2.0,
                                 1.0, 1.0, 2.0,
                                 2.0, 3.0,
                                 1.0, 3.0
                               };

The ptrb, ptre, asub, and aval arguments define the constraint matrix \(A\) in the column ordered sparse format (for details, see Section 16.1.4.2).

Using the function MSK_putacol we set column \(j\) of \(A\)

        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.*/

There are many alternative formats for entering the \(A\) matrix. See functions such as MSK_putarow, MSK_putarowlist, MSK_putaijlist and similar.

Finally, the bounds on each constraint are set by looping over each constraint index \(i= 0, \ldots,\mathtt{numcon}-1\)

    /* 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.*/

Optimization

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

      r = MSK_optimizetrm(task, &trmcode);

Extracting the solution.

After optimizing the status of the solution is examined with a call to MSK_getsolsta. If the solution status is reported as MSK_SOL_STA_OPTIMAL or MSK_SOL_STA_NEAR_OPTIMAL the solution is extracted in the lines below:

                MSK_getxx(task,
                          MSK_SOL_BAS,    /* Request the basic solution. */
                          xx);

The MSK_getxx 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 MSK_SOL_BAS.

Source code

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

Listing 6.1 Linear optimization example. Click here to download.
#include <stdio.h>
#include "mosek.h"

/* This function prints log output from MOSEK to the terminal. */
static void MSKAPI printstr(void       *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

int main(int argc, const char *argv[])
{
  const MSKint32t    numvar = 4,
                     numcon = 3;

  const double       c[]     = {3.0, 1.0, 5.0, 1.0};
  /* Below is the sparse representation of the A
     matrix stored by column. */
  const MSKint32t    aptrb[] = {0, 2, 5, 7},
                     aptre[] = {2, 5, 7, 9},
                     asub[]  = { 0, 1,
                                 0, 1, 2,
                                 0, 1,
                                 1, 2
                               };
  const double       aval[]  = { 3.0, 2.0,
                                 1.0, 1.0, 2.0,
                                 2.0, 3.0,
                                 1.0, 3.0
                               };

  /* Bounds on constraints. */
  const MSKboundkeye bkc[]  = {MSK_BK_FX, MSK_BK_LO,     MSK_BK_UP    };
  const double       blc[]  = {30.0,      15.0,          -MSK_INFINITY};
  const double       buc[]  = {30.0,      +MSK_INFINITY, 25.0         };
  /* Bounds on variables. */
  const MSKboundkeye bkx[]  = {MSK_BK_LO,     MSK_BK_RA, MSK_BK_LO,     MSK_BK_LO     };
  const double       blx[]  = {0.0,           0.0,       0.0,           0.0           };
  const double       bux[]  = { +MSK_INFINITY, 10.0,      +MSK_INFINITY, +MSK_INFINITY };
  MSKenv_t           env  = NULL;
  MSKtask_t          task = NULL;
  MSKrescodee        r;
  MSKint32t          i, j;

  /* 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);

    /* Directs the log task stream to the 'printstr' function. */
    if ( r == MSK_RES_OK )
      r = 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);

    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.*/

    /* Maximize objective function. */
    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MAXIMIZE);

    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_LOG);

      if ( r == MSK_RES_OK )
      {
        MSKsolstae solsta;

        if ( r == MSK_RES_OK )
          r = MSK_getsolsta (task,
                             MSK_SOL_BAS,
                             &solsta);
        switch (solsta)
        {
          case MSK_SOL_STA_OPTIMAL:
          case MSK_SOL_STA_NEAR_OPTIMAL:
            {
              double *xx = (double*) calloc(numvar, sizeof(double));
              if ( xx )
              {
                MSK_getxx(task,
                          MSK_SOL_BAS,    /* Request the basic solution. */
                          xx);

                printf("Optimal primal solution\n");
                for (j = 0; j < numvar; ++j)
                  printf("x[%d]: %e\n", j, xx[j]);

                free(xx);
              }
              else
                r = MSK_RES_ERR_SPACE;

              break;
            }
          case MSK_SOL_STA_DUAL_INFEAS_CER:
          case MSK_SOL_STA_PRIM_INFEAS_CER:
          case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
          case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:
            printf("Primal or dual infeasibility certificate found.\n");
            break;
          case MSK_SOL_STA_UNKNOWN:
            {
              char symname[MSK_MAX_STR_LEN];
              char desc[MSK_MAX_STR_LEN];

              /* If the solutions status is unknown, print the termination code
                 indicating why the optimizer terminated prematurely. */

              MSK_getcodedesc(trmcode,
                              symname,
                              desc);

              printf("The solution status is unknown.\n");
              printf("The optimizer terminitated with code: %s\n", symname);
              break;
            }
          default:
            printf("Other solution status.\n");
            break;
        }
      }
    }

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