Mixed Integer Programming

This chapter extends the model developed in Chapter 10 to a Mixed Integer Programming (MIP) problem. It describes how to

• define different types of discrete variables,
• get the MIP solution status and understand the MIP optimization log produced by Xpress Optimizer.

Chapter 6 shows how to formulate and solve the same example with Mosel and in Chapter 16 the problem is input and solved directly with Xpress Optimizer.

Extended problem description

The investor is unwilling to have small share holdings. He looks at the following two possibilities to formulate this constraint:

1. Limiting the number of different shares taken into the portfolio.
2. If a share is bought, at least a minimum amount 10% of the budget is spent on the share.

We are going to deal with these two constraints in two separate models.

MIP model 1: limiting the number of different shares

To be able to count the number of different values we are investing in, we introduce a second set of variables buy_s in the LP model developed in Chapter 2. These variables are indicator variables or binary variables. A variable buy_s takes the value 1 if the share s is taken into the portfolio and 0 otherwise.

We introduce the following constraint to limit the total number of assets to a maximum of MAXNUM. It expresses the constraint that at most MAXNUM of the variables buy_s may take the value 1 at the same time.

We now still need to link the new binary variables buy_s with the variables frac_s, the quantity of every share selected into the portfolio. The relation that we wish to express is `if a share is selected into the portfolio, then it is counted in the total number of values' or `if frac_s > 0 then buy_s = 1'. The following inequality formulates this implication:

If, for some s, frac_s is non-zero, then buy_s must be greater than 0 and hence 1. Conversely, if buy_s is at 0, then frac_s is also 0, meaning that no fraction of share s is taken into the portfolio. Notice that these constraints do not prevent the possibility that buy_s is at 1 and frac_s at 0. However, this does not matter in our case, since any solution in which this is the case is also valid with both variables, buy_s and frac_s, at 0.

Implementation with BCL

We extend the LP model developed in Chapter 10 with the new variables and constraints. The fact that the new variables are binary variables (i.e. they only take the values 0 and 1) is expressed through the type XPRB_BV at their creation.

Another common type of discrete variable is an integer variable, that is, a variable that can only take on integer values between given lower and upper bounds. These variables are defined in BCL with the type XPRB_UI. In the following section (MIP model 2) we shall see yet another example of discrete variables, namely semi-continuous variables.

#include <iostream>
#include "xprb_cpp.h"

using namespace std;
using namespace ::dashoptimization;

#define MAXNUM 4                   // Max. number of shares to be selected

#define NSHARES 10                 // Number of shares
#define NRISK 5                    // Number of high-risk shares
#define NNA 4                      // Number of North-American shares

double RET[] = {5,17,26,12,8,9,7,6,31,21};  // Estimated return in investment
int RISK[] = {1,2,3,8,9};          // High-risk values among shares
int NA[] = {0,1,2,3};              // Shares issued in N.-America

int main(int argc, char **argv)
{
 int s;
 XPRBprob p("FolioMIP1");          // Initialize a new problem in BCL
 XPRBexpr Risk,Na,Return,Cap,Num;
 XPRBvar frac[NSHARES];            // Fraction of capital used per share
 XPRBvar buy[NSHARES];             // 1 if asset is in portfolio, 0 otherwise

// Create the decision variables (including upper bounds for `frac')
 for(s=0;s<NSHARES;s++)
 {
  frac[s] = p.newVar("frac", XPRB_PL, 0, 0.3);
  buy[s] = p.newVar("buy", XPRB_BV);
 }

// Objective: total return
 for(s=0;s<NSHARES;s++) Return += RET[s]*frac[s];
 p.setObj(Return);                 // Set the objective function

// Limit the percentage of high-risk values
 for(s=0;s<NRISK;s++) Risk += frac[RISK[s]];
 p.newCtr(Risk <= 1.0/3);

// Minimum amount of North-American values
 for(s=0;s<NNA;s++) Na += frac[NA[s]];
 p.newCtr(Na >= 0.5);

// Spend all the capital
 for(s=0;s<NSHARES;s++) Cap += frac[s];
 p.newCtr(Cap == 1);

// Limit the total number of assets
 for(s=0;s<NSHARES;s++) Num += buy[s];
 p.newCtr(Num <= MAXNUM);

// Linking the variables
 for(s=0;s<NSHARES;s++) p.newCtr(frac[s] <= buy[s]);

// Solve the problem
 p.setSense(XPRB_MAXIM);
 p.mipOptimize("");

// Solution printing
 cout << "Total return: " << p.getObjVal() << endl;
 for(s=0;s<NSHARES;s++)
  cout << s << ": " << frac[s].getSol()*100 << "% (" << buy[s].getSol()
       << ")" << endl;

 return 0;
}

Besides the additional variables and constraints, the choice of optimization algorithm needs to be adapted to the problem type: we now wish to solve a MIP problem via Branch-and-Bound, and we therefore use the method mipOptimize.

Just as with the LP problem in the previous chapter, it is usually helpful to check the solution status before accessing the MIP solution—only if the MIP status is `unfinished (solution found)' or `optimal' will BCL print out a meaningful solution:

 char *MIPSTATUS[] = {"not loaded", "not optimized", "LP optimized",
                      "unfinished (no solution)",
                      "unfinished (solution found)", "infeasible",
                      "optimal", "unbounded"};

 cout << "Problem status: " << MIPSTATUS[p.getMIPStat()] << endl;

Analyzing the solution

As the result of the execution of our program we obtain the following output:

Reading Problem FolioMIP1
Problem Statistics
          14 (    514 spare) rows
          20 (      0 spare) structural columns
          49 (   5056 spare) non-zero elements
Global Statistics
          10 entities        0 sets        0 set members
Maximizing MILP FolioMIP1
Original problem has:
        14 rows           20 cols           49 elements        10 globals
Presolved problem has:
        13 rows           19 cols           46 elements         9 globals
LP relaxation tightened
Will try to keep branch and bound tree memory usage below 14.8Gb
Starting concurrent solve with dual

 Concurrent-Solve,   0s
            Dual
    objective   dual inf
 D  14.066667   .0000000
------- optimal --------
Concurrent statistics:
      Dual: 4 simplex iterations, 0.00s
Optimal solution found

   Its         Obj Value      S   Ninf  Nneg   Sum Dual Inf  Time
     4         14.066667      D      0     0        .000000     0
Dual solved problem
  4 simplex iterations in 0s

Final objective                         : 1.406666666666667e+01
  Max primal violation      (abs / rel) : 5.551e-17 / 5.551e-17
  Max dual violation        (abs / rel) :       0.0 /       0.0
  Max complementarity viol. (abs / rel) :       0.0 /       0.0
All values within tolerances

Starting root cutting & heuristics

 Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time
c           13.100000    14.066667      1                  6.87
   1  K     13.100000    13.908571      1      1      0    5.81      0
   2  K     13.100000    13.580000      1     12      0    3.53      0
 *** Search completed ***     Time:     0 Nodes:          1
Number of integer feasible solutions found is 1
Best integer solution found is    13.100000
Best bound is    13.100014
Uncrunching matrix
Problem status: optimal
Total return: 13.1
0: 20% (1)
1: 0% (0)
2: 30% (1)
3: 0% (0)
4: 20% (1)
5: 30% (1)
6: 0% (0)
7: 0% (0)
8: 0% (0)
9: 0% (0)

At the beginning we see the log of the execution of Xpress Optimizer: the problem statistics (we now have 14 constraints and 20 variables, out of which 10 are MIP variables, refered to as `entities'), the log of the execution of the LP algorithm (concurrently solving with primal and dual simplex on a multi-core processor), the log of the built-in MIP heuristics (a solution with the value 13.1 has been found) and the automated cut generation (a total of 13 cuts of type `K' = knapsack have been generated). Since this problem is very small, it is solved by the MIP heuristics and the addition of cuts (additional constraints that cut off parts of the LP solution space, but no MIP solution) tightens the LP formulation in such a way that the solution to the LP relaxation becomes integer feasible. The Branch-and-Bound search therefore stops at the first node and no log of the Branch-and-Bound search gets displayed.

The output printed by our program tells us that the problem has been solved to optimality (i.e. the MIP search has been completed and at least one integer feasible solution has been found). The maximum return is now lower than in the original LP problem due to the additional constraint. As required, only four different shares are selected to form the portfolio.

MIP model 2: imposing a minimum investment in each share

To formulate the second MIP model, we start again with the LP model from Chapters 2 and 10. The new constraint we wish to formulate is `if a share is bought, at least a minimum amount 10% of the budget is spent on the share.' Instead of simply constraining every variable frac_s to take a value between 0 and 0.3, it now must either lie in the interval between 0.1 and 0.3 or take the value 0. This type of variable is known as semi-continuous variable. In the new model, we replace the bounds on the variables frac_s by the following constraint:

Implementation with BCL

The following program implements the MIP model 2. The semi-continuous variables are defined by the type XPRB_SC. By default, BCL assumes a continuous limit of 1, se we need to set this value to 0.1 with the method setLim.

A similar type is available for integer variables that take either the value 0 or an integer value between a given limit and their upper bound (so-called semi-continuous integers): XPRB_SI. A third composite type is a partial integer which takes integer values from its lower bound to a given limit value and is continuous beyond this value (marked by XPRB_PI).

#include <iostream>
#include "xprb_cpp.h"

using namespace std;
using namespace ::dashoptimization;

#define NSHARES 10                   // Number of shares
#define NRISK 5                      // Number of high-risk shares
#define NNA 4                        // Number of North-American shares

double RET[] = {5,17,26,12,8,9,7,6,31,21};  // Estimated return in investment
int RISK[] = {1,2,3,8,9};            // High-risk values among shares
int NA[] = {0,1,2,3};                // Shares issued in N.-America

int main(int argc, char **argv)
{
 int s;
 XPRBprob p("FolioSC");              // Initialize a new problem in BCL
 XPRBexpr Risk,Na,Return,Cap;
 XPRBvar frac[NSHARES];              // Fraction of capital used per share

// Create the decision variables
 for(s=0;s<NSHARES;s++)
 {
  frac[s] = p.newVar("frac", XPRB_SC, 0, 0.3);
  frac[s].setLim(0.1);
 }

// Objective: total return
 for(s=0;s<NSHARES;s++) Return += RET[s]*frac[s];
 p.setObj(Return);                   // Set the objective function

// Limit the percentage of high-risk values
 for(s=0;s<NRISK;s++) Risk += frac[RISK[s]];
 p.newCtr(Risk <= 1.0/3);

// Minimum amount of North-American values
 for(s=0;s<NNA;s++) Na += frac[NA[s]];
 p.newCtr(Na >= 0.5);

// Spend all the capital
 for(s=0;s<NSHARES;s++) Cap += frac[s];
 p.newCtr(Cap == 1);

// Solve the problem
 p.setSense(XPRB_MAXIM);
 p.mipOptimize("");

// Solution printing
 cout << "Total return: " << p.getObjVal() << endl;
 for(s=0;s<NSHARES;s++)
  cout << s << ": " << frac[s].getSol()*100 << "%" << endl;

 return 0;
}

When executing this program we obtain the following output (leaving out the part printed by the Optimizer):

Total return: 14.0333
0: 30%
1: 0%
2: 20%
3: 0%
4: 10%
5: 26.6667%
6: 0%
7: 0%
8: 13.3333%
9: 0%

Now five securities are chosen for the portfolio, each forming at least 10% and at most 30% of the total investment. Due to the additional constraint, the optimal MIP solution value is again lower than the initial LP solution value.