Implementation with BCL

For the implementation of the variable fixing solution heuristic we work with the MIP 1 model from Chapter 11. Through the definition of the heuristic in a separate function we only make minimal changes to the model itself: before solving our problem with the standard call to the method mipOptimize we execute our own solution heuristic and the solution printing also has been adapted.

#include <iostream>
#include "xprb_cpp.h"
#include "xprs.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

void solveHeur();

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

XPRBprob p("FolioMIPHeur");        // Initialize a new problem in BCL
XPRBvar frac[NSHARES];             // Fraction of capital used per share
XPRBvar buy[NSHARES];              // 1 if asset is in portfolio, 0 otherwise

int main(int argc, char **argv)
{
 int s;
 XPRBexpr Risk,Na,Return,Cap,Num;

// 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 problem heuristically
 p.setSense(XPRB_MAXIM);
 solveHeur();

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

// Solution printing
 if(p.getMIPStat()==4 || p.getMIPStat()==6)
 {
  cout << "Exact solution: Total return: " << p.getObjVal() << endl;
  for(s=0;s<NSHARES;s++)
   cout << s << ": " << frac[s].getSol()*100 << "%" << endl;
 }
 else
  cout << "Heuristic solution is optimal." << endl;

 return 0;
}

void solveHeur()
{
 XPRBbasis basis;
 int s, ifgsol;
 double solval, bsol[NSHARES],TOL;

 XPRSsetintcontrol(p.getXPRSprob(), XPRS_CUTSTRATEGY, 0);
                                   // Disable automatic cuts
 XPRSsetintcontrol(p.getXPRSprob(), XPRS_HEURSTRATEGY, 0);
                                   // Disable MIP heuristics
 XPRSsetintcontrol(p.getXPRSprob(), XPRS_PRESOLVE, 0);
                                   // Switch presolve off
 XPRSgetdblcontrol(p.getXPRSprob(), XPRS_FEASTOL, &TOL);
                                   // Get feasibility tolerance

 p.mipOptimize("l");               // Solve the LP-relaxation
 basis=p.saveBasis();              // Save the current basis

// Fix all variables `buy' for which `frac' is at 0 or at a relatively
// large value
 for(s=0;s<NSHARES;s++)
 {
  bsol[s]=buy[s].getSol();         // Get the solution values of `frac'
  if(bsol[s] < TOL)  buy[s].setUB(0);
  else if(bsol[s] > 0.2-TOL)  buy[s].setLB(1);
 }

 p.mipOptimize("c");               // Solve the MIP-problem
 ifgsol=0;
 if(p.getMIPStat()==4 || p.getMIPStat()==6)
 {                                 // If an integer feas. solution was found
  ifgsol=1;
  solval=p.getObjVal();            // Get the value of the best solution
  cout << "Heuristic solution: Total return: " << p.getObjVal() << endl;
  for(s=0;s<NSHARES;s++)
   cout << s << ": " << frac[s].getSol()*100 << "%" << endl;
 }

// Reset variables to their original bounds
 for(s=0;s<NSHARES;s++)
  if((bsol[s] < TOL) || (bsol[s] > 0.2-TOL))
  {
   buy[s].setLB(0);
   buy[s].setUB(1);
  }

 p.loadBasis(basis);        /* Load the saved basis: bound changes are
                               immediately passed on from BCL to the
                               Optimizer if the problem has not been modified
                               in any other way, so that there is no need to
                               reload the matrix */
 basis.reset();             // No need to store the saved basis any longer
 if(ifgsol==1)
  XPRSsetdblcontrol(p.getXPRSprob(), XPRS_MIPABSCUTOFF, solval+TOL);
                            // Set the cutoff to the best known solution
}

The implementation of the heuristic certainly requires some explanations.

In this example for the first time we use the direct access to Xpress Optimizer. To do so, we need to include the Optimizer header file xprs.h. The Optimizer functions are applied to the problem representation (of type XPRSprob) held by the Optimizer which can be retrieved with the method getXPRSprob of the BCL problem. For more detail on how to use the BCL and Optimizer libaries in combination the reader is refered to the `BCL Reference Manual'. The complete documentation of all Optimizer functions and parameters is provided in the `Optimizer Reference Manual'.

Parameters: The solution heuristic starts with parameter settings for the Xpress Optimizer. Switching off the automated cut generation (parameter XPRS_CUTSTRATEGY) and the MIP heuristics (parameter XPRS_HEURSTRATEGY) is optional, whereas it is required in our case to disable the presolve mechanism (a treatment of the matrix that tries to reduce its size and improve its numerical properties, set with parameter XPRS_PRESOLVE), because we interact with the problem in the Optimizer in the course of its solution and this is only possible correctly if the matrix has not been modified by the Optimizer.

In addition to the parameter settings we also retrieve the feasibility tolerance used by Xpress Optimizer: the Optimizer works with tolerance values for integer feasibility and solution feasibility that are typically of the order of 10^-6 by default. When evaluating a solution, for instance by performing comparisons, it is important to take into account these tolerances.

The fine tuning of output printing mentioned in Chapter 10 can be obtained by setting the parameters XPRS_LPLOG and XPRS_MIPLOG (both to be set with function XPRSsetintcontrol).

Optimization calls: We use the optimization method mipOptimize with the argument "l", indicating that we only want to solve the top node LP relaxation (and not yet the entire MIP problem). To continue with MIP solving from the point where we have stopped the algorithm we use the argument "c".

Saving and loading bases: To speed up the solution process, we save (in memory) the current basis of the Simplex algorithm after solving the initial LP relaxation, before making any changes to the problem. This basis is loaded again at the end, once we have restored the original problem. The MIP solution algorithm then does not have to re-solve the LP problem from scratch, it resumes the state where it was `interrupted' by our heuristic.

Bound changes: When a problem has already been loaded into the Optimizer (e.g. after executing an optimization statement or following an explicit call to method loadMat) bound changes via setLB and setUB are passed on directly to the Optimizer. Any other changes (addition or deletion of constraints or variables) always lead to a complete reloading of the problem.

The program produces the following output. As can be seen, when solving the original problem for the second time the Simplex algorithm performs 0 iterations because it has been started with the basis of the optimal solution saved previously.

Reading Problem FolioMIPHeur
Problem Statistics
          14 (      0 spare) rows
          20 (      0 spare) structural columns
          49 (      0 spare) non-zero elements
Global Statistics
          10 entities        0 sets        0 set members
Maximizing MILP FolioMIPHeur
Original problem has:
        14 rows           20 cols           49 elements        10 globals
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: 5 simplex iterations, 0.00s
Optimal solution found

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

Final objective                         : 1.406666666666667e+01
  Max primal violation      (abs / rel) :       0.0 /       0.0
  Max dual violation        (abs / rel) :       0.0 /       0.0
  Max complementarity viol. (abs / rel) :       0.0 /       0.0
All values within tolerances
 *** Search unfinished ***    Time:     0 Nodes:          0
Number of integer feasible solutions found is 0
Best bound is    14.066667

Starting root cutting & heuristics

 Its Type    BestSoln    BestBound   Sols    Add    Del     Gap     GInf   Time
a           13.100000    14.066667      1                  6.87
Heuristic search started
Heuristic search stopped
 *** 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
Heuristic solution: Total return: 13.1
0: 20%
1: 0%
2: 30%
3: 0%
4: 20%
5: 30%
6: 0%
7: 0%
8: 0%
9: 0%
Maximizing MILP FolioMIPHeur
Original problem has:
        14 rows           20 cols           49 elements        10 globals
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: 0 simplex iterations, 0.00s
Optimal solution found

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

Final objective                         : 1.406666666666667e+01
  Max primal violation      (abs / rel) :       0.0 /       0.0
  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
Heuristic search started
Heuristic search stopped

Starting tree search.
Deterministic mode with up to 8 running threads and up to 16 tasks.

    Node     BestSoln    BestBound   Sols Active  Depth     Gap     GInf   Time
 *** Search completed ***     Time:     0 Nodes:          5
Problem is integer infeasible
Number of integer feasible solutions found is 0
Best bound is    13.100001
Heuristic solution is optimal.

Observe that the heuristic found a solution of 13.1, and that the MIP optimizer without the heuristic could not find a better solution (hence the infeasible message). The heuristic solution is therefore optimal.