6.8 Integer Optimization

An optimization problem where one or more of the variables are constrained to integer values is called a (mixed) integer optimization problem. MOSEK supports integer variables in combination with linear, quadratic and quadratically constrtained and conic problems (except semidefinite). See the previous tutorials for an introduction to how to model these types of problems.

6.8.1 Basic linear example

We use the example

(6.27)\[\begin{split}\begin{array}{lccl} \mbox{maximize} & x_0 + 0.64 x_1 & & \\ \mbox{subject to} & 50 x_0 + 31 x_1 & \leq & 250, \\ & 3 x_0 - 2 x_1 & \geq & -4, \\ & x_0, x_1 \geq 0 & & \mbox{and integer} \end{array}\end{split}\]

to demonstrate how to set up and solve a problem with integer variables. It has the structure of a linear optimization problem except for integrality constraints on the variables. Therefore, only the specification of the integer constraints requires something new compared to the linear optimization problem discussed previously.

First, the integrality constraints are imposed using the function MSK_putvartype or one of its bulk analogues:

    for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
      r = MSK_putvartype(task, j, MSK_VAR_TYPE_INT);

Next, the example demonstrates how to set various useful parameters of the mixed-integer optimizer. See Sec. 13.5 (The Mixed-Integer optimizer) for details.

      /* Set max solution time */
      r = MSK_putdouparam(task,
                          MSK_DPAR_MIO_MAX_TIME,
                          60.0);

The complete source for the example is listed Listing 6.13. Please note that when we fetch the solution then the integer solution is requested by using MSK_SOL_ITG. No dual solution is defined for integer optimization problems.

Listing 6.13 Source code implementing problem (6.27). Click here to download.
#include <stdio.h>
#include "mosek.h" /* Include the MOSEK definition file. */

static void MSKAPI printstr(void *handle,
                            const char str[])
{
  printf("%s", str);
} /* printstr */

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

  double       c[]   = {  1.0, 0.64 };
  MSKboundkeye bkc[] = { MSK_BK_UP,    MSK_BK_LO };
  double       blc[] = { -MSK_INFINITY, -4.0 };
  double       buc[] = { 250.0,        MSK_INFINITY };

  MSKboundkeye bkx[] = { MSK_BK_LO,    MSK_BK_LO };
  double       blx[] = { 0.0,          0.0 };
  double       bux[] = { MSK_INFINITY, MSK_INFINITY };


  MSKint32t    aptrb[] = { 0, 2 },
               aptre[] = { 2, 4 },
               asub[] = { 0,    1,   0,    1 };
  double       aval[] = { 50.0, 3.0, 31.0, -2.0 };
  MSKint32t    i, j;

  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;
  MSKrescodee  r;

  /* Create the mosek environment. */
  r = MSK_makeenv(&env, NULL);

  /* Check if return code is ok. */
  if (r == MSK_RES_OK)
  {
    /* Create the optimization task. */
    r = MSK_maketask(env, 0, 0, &task);

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

    /* Optionally add a constant term to the objective. */
    if (r == MSK_RES_OK)
      r = MSK_putcfix(task, 0.0);
    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.*/

    /* Specify integer variables. */
    for (j = 0; j < numvar && r == MSK_RES_OK; ++j)
      r = MSK_putvartype(task, j, MSK_VAR_TYPE_INT);

    if (r == MSK_RES_OK)
      r =  MSK_putobjsense(task,
                           MSK_OBJECTIVE_SENSE_MAXIMIZE);

    if (r == MSK_RES_OK)
      /* Set max solution time */
      r = MSK_putdouparam(task,
                          MSK_DPAR_MIO_MAX_TIME,
                          60.0);

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

      if (r == MSK_RES_OK)
      {
        MSKint32t  j;
        MSKsolstae solsta;
        double     *xx = NULL;

        MSK_getsolsta(task, MSK_SOL_ITG, &solsta);

        xx = calloc(numvar, sizeof(double));
        if (xx)
        {
          switch (solsta)
          {
            case MSK_SOL_STA_INTEGER_OPTIMAL:
              MSK_getxx(task,
                        MSK_SOL_ITG,    /* Request the integer solution. */
                        xx);

              printf("Optimal solution.\n");
              for (j = 0; j < numvar; ++j)
                printf("x[%d]: %e\n", j, xx[j]);
              break;
            case MSK_SOL_STA_PRIM_FEAS:
              /* A feasible but not necessarily optimal solution was located. */
              MSK_getxx(task, MSK_SOL_ITG, xx);

              printf("Feasible solution.\n");
              for (j = 0; j < numvar; ++j)
                printf("x[%d]: %e\n", j, xx[j]);
              break;
            case MSK_SOL_STA_UNKNOWN:
              {
                MSKprostae prosta;
                MSK_getprosta(task, MSK_SOL_ITG, &prosta);
                switch (prosta)
                {
                  case MSK_PRO_STA_PRIM_INFEAS_OR_UNBOUNDED:
                    printf("Problem status Infeasible or unbounded\n");
                    break;
                  case MSK_PRO_STA_PRIM_INFEAS:
                    printf("Problem status Infeasible.\n");
                    break;
                  case MSK_PRO_STA_UNKNOWN:
                    printf("Problem status unknown. Termination code %d.\n", trmcode);
                    break;
                  default:
                    printf("Other problem status.");
                    break;
                }
              }
              break;
            default:
              printf("Other solution status.");
              break;
          }
        }
        else
        {
          r = MSK_RES_ERR_SPACE;
        }
        free(xx);
      }
    }

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

    MSK_deletetask(&task);
  }
  MSK_deleteenv(&env);

  printf("Return code: %d.\n", r);
  return (r);
} /* main */

6.8.2 Specifying an initial solution (hot-start)

It is a common strategy to provide a starting feasible point (if one is known in advance) to the mixed-integer solver. This is known as hot-start or warm-start. The feasible point may come from the user’s prior knowledge of the model, a heuristic, or a solution from a preceding solve if the problem was modified so that the solution remains feasible. Using hot-start allows the solver to skip worse solutions and potentially get closer to the optimum faster.

There are two modes for MOSEK to utilize an initial solution.

  • A complete solution. MOSEK will first try to check if the current value of the primal variable solution is a feasible point. The solution can either come from a previous solver call or can be entered by the user, however the full solution with values for all variables (both integer and continuous) must be provided. This check is always performed and does not require any extra action from the user. The outcome of this process can be inspected via information items MSK_IINF_MIO_INITIAL_FEASIBLE_SOLUTION and MSK_DINF_MIO_INITIAL_FEASIBLE_SOLUTION_OBJ, and via the Initial feasible solution objective entry in the log.

  • A partial integer solution. MOSEK can also try to construct a feasible solution by fixing integer variables to the values provided by the user (rounding if necessary) and optimizing over the remaining continuous variables. In this setup the user must provide initial values for all integer variables. This action is only performed if the parameter MSK_IPAR_MIO_CONSTRUCT_SOL is switched on. The outcome of this process can be inspected via information items MSK_IINF_MIO_CONSTRUCT_SOLUTION and MSK_DINF_MIO_CONSTRUCT_SOLUTION_OBJ, and via the Construct solution objective entry in the log.

In the following example we focus on inputting a partial integer solution.

(6.28)\[\begin{split}\begin{array} {ll} \mbox{maximize} & 7 x_0 + 10 x_1 + x_2 + 5 x_3 \\ \mbox{subject to} & x_0 + x_1 + x_2 + x_3 \leq 2.5\\ & x_0,x_1,x_2 \in \integral \\ & x_0,x_1,x_2,x_3 \geq 0 \end{array}\end{split}\]

Solution values can be set using MSK_putxx, MSK_putxxslice or similar . If the solution to be used as hot-start happens to come from a previous solve of the same model/task then it will be used by the solver automatically and does not have to be input explicitly through the API again.

Listing 6.14 Implementation of problem (6.28) specifying an initial solution. Click here to download.
  if (r == MSK_RES_OK)
  {
    /* Assign values to integer variables 
       (we only set a slice of xx) */
    double       xxInit[] = {1.0, 1.0, 0.0};
    r = MSK_putxxslice(task, MSK_SOL_ITG, 0, 3, xxInit);
  }
  if (r == MSK_RES_OK)
  {
    /* Request constructing the solution from integer variable values */
    r = MSK_putintparam(task, MSK_IPAR_MIO_CONSTRUCT_SOL, MSK_ON);
  }

The log output from the optimizer will in this case indicate that the inputted values were used to construct an initial feasible solution:

Construct solution objective       : 1.950000000000e+01

The same information can be obtained from the API:

Listing 6.15 Retrieving information about usage of initial solution Click here to download.
    MSK_getintinf(task, MSK_IINF_MIO_CONSTRUCT_SOLUTION, &constr);
    MSK_getdouinf(task, MSK_DINF_MIO_CONSTRUCT_SOLUTION_OBJ, &constr_obj);
    printf("Construct solution utilization: %d\nConstruct solution objective: %.3f\n", constr, constr_obj);

6.8.3 Basic conic example

Integer variables can also be used arbitrarily in conic problems (except semidefinite). We refer to the previous tutorials for how to set up a conic optimization problem. Here we present sample code that sets up a simple optimization problem:

(6.29)\[\begin{split}\begin{array}{ll} \mbox{minimize} & x^2+y^2 \\ \mbox{subject to} & x \geq e^y+3.8, \\ & x, y \ \mbox{integer}. \end{array}\end{split}\]

The canonical conic formulation of (6.29) suitable for Optimizer API for C is

(6.30)\[\begin{split}\begin{array}{llr} \mbox{minimize} & t & \\ \mbox{subject to} & (t,x,y)\in\Q^3 & (t\geq\sqrt{x^2+y^2}) \\ & (x-3.8, 1, y) \in\EXP & (x-3.8\geq e^y) \\ & x, y \ \mbox{integer}, & \\ & t\in\real. \end{array}\end{split}\]
Listing 6.16 Implementation of problem (6.30). Click here to download.
int main(int argc, char *argv[])
{
  MSKint32t         numvar = 3;  /* x, y, t */

  MSKvariabletypee  vart[] = { MSK_VAR_TYPE_INT, MSK_VAR_TYPE_INT };
  MSKint32t       intsub[] = { 0, 1 }; 

  MSKint32t    i, j;

  MSKenv_t     env = NULL;
  MSKtask_t    task = NULL;
  MSKrescodee  r, trm;

  r = MSK_makeenv(&env, NULL);

  if (r == MSK_RES_OK)
  {
    r = MSK_maketask(env, 0, 0, &task);

    if (r == MSK_RES_OK)
      r = MSK_linkfunctotaskstream(task, MSK_STREAM_LOG, NULL, printstr);

    if (r == MSK_RES_OK)
      r = MSK_appendvars(task, numvar);
    if (r == MSK_RES_OK)
      r = MSK_putvarboundsliceconst(task, 0, numvar, MSK_BK_FR, -0.0, 0.0);    

    /* Integrality constraints */
    if (r == MSK_RES_OK)
      r = MSK_putvartypelist(task, 2, intsub, vart);        

    /* Objective */
    if (r == MSK_RES_OK)
      r = MSK_putobjsense(task, MSK_OBJECTIVE_SENSE_MINIMIZE);
    if (r == MSK_RES_OK)
      r = MSK_putcj(task, 2, 1.0);    /* Minimize t */

    /* Conic part of the problem */
    if (r == MSK_RES_OK) 
    {
      /* Set up the affine expressions */
      /* x, x-3.8, y, t, 1.0 */
      MSKint64t afeidx[] = {0, 1, 2, 3};
      MSKint32t varidx[] = {0, 0, 1, 2};
      MSKrealt     val[] = {1.0, 1.0, 1.0, 1.0};
      MSKrealt       g[] = {0.0, -3.8, 0.0, 0.0, 1.0};
      MSKint64t domExp, domQuad;
      MSKint64t afeidxExp[]  = {1, 4, 2};
      MSKint64t afeidxQuad[] = {3, 0, 2};

      if (r == MSK_RES_OK)
        r = MSK_appendafes(task, 5);
      if (r == MSK_RES_OK)
        r = MSK_putafefentrylist(task,
                                 4,
                                 afeidx,
                                 varidx,
                                 val);
      if (r == MSK_RES_OK)
        r = MSK_putafegslice(task, 0, 5, g);

      // Add constraint (x-3.8, 1, y) \in \EXP
      if (r == MSK_RES_OK)
        r = MSK_appendprimalexpconedomain(task, &domExp);
      if (r == MSK_RES_OK)
        r = MSK_appendacc(task, domExp, 3, afeidxExp, NULL);

      // Add constraint (t, x, y) \in \QUAD
      if (r == MSK_RES_OK)
        r = MSK_appendquadraticconedomain(task, 3, &domQuad);
      if (r == MSK_RES_OK)
        r = MSK_appendacc(task, domQuad, 3, afeidxQuad, NULL);
    }

    /* Optimize the problem */
    if (r == MSK_RES_OK)
      r = MSK_optimizetrm(task, &trm);

    if (r == MSK_RES_OK)
      r = MSK_solutionsummary(task, MSK_STREAM_MSG);

    if (r == MSK_RES_OK)
    {
      MSKrealt xx[] = {0, 0};

      r = MSK_getxxslice(task, MSK_SOL_ITG, 0, 2, xx);
      
      if (r == MSK_RES_OK)
        printf("x = %.2f, y = %.2f\n", xx[0], xx[1]);
    }

    if (task) MSK_deletetask(&task);
  }

  if (env) MSK_deleteenv(&env);
  return r;
} 

Error and solution status handling were omitted for readability.

6.8.4 Fixed problem and dual values

The dual solution is not defined for mixed-integer problems, but in some cases the user may want to obtain some dual information (shadow prices). One typical strategy is to compute shadow price information under the assumption that the combinatorial decisions (integer variable values) do not change, that is:

  • solve the mixed-integer model (to some feasible solution, not necessarily optimal),

  • fix all integer variables to their values in the solution,

  • solve the fixed model as a continuous problem and extract the dual values.

Optimizer API for C facilitates the construction of the fixed model with with MSK_getfixedproblem.

As an example we consider a toy production planning model with two plants with maximum capacities 70, 80 units and a demand of 100 units:

(6.31)\[\begin{split}\begin{array}{ll} \mbox{minimize} & 10x_1 + 14x_2 + 500s_1 + 300s_2 \\ \mbox{subject to} & x_1+x_2\geq 100, \\ & 0 \leq x_1 \leq 70s_1, \\ & 0 \leq x_2 \leq 80s_2, \\ & s_1,s_2 \in \{0,1\}. \end{array}\end{split}\]

We begin by solving the mixed-integer problem and verifying that it has a feasible solution with \(s_1=s_2=1\) i.e. both plants active with production levels \((x_1,x_2)=(70,30)\). Then, assuming no modifications are made to the problem structure or numerical data, we can immediately construct and solve the fixed integer model and retrieve its solution as shown below.

Listing 6.17 Solving the fixed integer model for (6.31). Click here to download.
  /* Create and solve the problem with fixed variables */
  MOSEKCALL(r, MSK_getfixedproblem(task, &fixTask));
  MOSEKCALL(r, MSK_linkfunctotaskstream(fixTask, MSK_STREAM_LOG, NULL, printstr));

  MOSEKCALL(r, MSK_optimize(fixTask));

  MOSEKCALL(r, MSK_getprosta(fixTask, MSK_SOL_BAS, &fixProsta));
  if (r == MSK_RES_OK && fixProsta != MSK_PRO_STA_PRIM_AND_DUAL_FEAS) {
    printf("Unsuitable problem status, exiting\n");
    return -1;
  }

  if (r == MSK_RES_OK) {
    double xfix[2], y[3];
    MOSEKCALL(r, MSK_getxxslice(fixTask, MSK_SOL_BAS, 0, 2, xfix));
    MOSEKCALL(r, MSK_getyslice(fixTask, MSK_SOL_BAS, 0, 3, y));
    printf("xfix = %.2f, %.2f\n", xfix[0], xfix[1]);
    printf("demand dual = %.2f\n", y[0]);
    printf("production dual = %.2f, %.2f\n", y[1], y[2]);
  }

For instance, the shadow price for the first plant’s production constraint is \(-4\), corresponding to the fact that increasing the capacity of the first plant by 1 unit would allow shifting one production unit form plant 2 to plant 1, reducing the objective cost by \(14-10=4\).

Note that the solution of the mixed-integer problem will typically have small violations, which implies that the fixed model may be declared (borderline) infeasible. To avoid this effect, the fixing algorithm will by default introduce small perturbations of bounds in order to make the fixed problem strictly feasible. If needed this can be controlled by setting the parameter MSK_IPAR_FIXING_METHOD in the original task prior to fixing.