Heuristics
In this chapter we show a simple binary variable fixing solution heuristic that involves a heuristic solution procedure interacting with Xpress-Optimizer through
- parameter settings,
- saving and recovering bases, and
- modifications of variable bounds.
Chapter 8 shows how to implement the same heuristic with Mosel.
Binary variable fixing heuristic
The heuristic we wish to implement should perform the following steps:
- Solve the LP relaxation and save the basis of the optimal solution
- Rounding heuristic: Fix all variables `buy' to 0 if they are close to 0, and to 1 if they have a relatively large value.
- Solve the resulting MIP problem.
- If an integer feasible solution was found, save the value of the best solution.
- Restore the original problem by resetting all variables to their original bounds, and load the saved basis.
- Solve the original MIP problem, using the heuristic solution as cutoff value.
Step 2: Since the fraction variables frac have an upper bound of 0.3, as a `relatively large value' in this case we may choose 0.2. In other applications, for binary
variables a more suitable choice may be 1-
, where is a very small value such as 10-5.
Step 6: Setting a cutoff value means that we only search for solutions that are better than this value. If the LP relaxation of a node is worse than this value it gets cut off, because this node and its descendants can only lead to integer feasible solutions that are even worse than the LP relaxation.
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 maxim 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;
XPRBlinExp 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
solveHeur();
// Solve the problem
p.maxim("g");
// 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.maxim(""); // Solve the LP-problem
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.maxim("g"); // 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;
}
XPRSinitglobal(p.getXPRSprob()); // Re-initialize the global search
// 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).
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
Its Obj Value S Ninf Nneg Sum Inf Time
0 42.600000 D 12 0 6.166667 0
11 14.066667 D 0 0 .000000 0
Optimal solution found
Its Obj Value S Ninf Nneg Sum Inf Time
0 14.066667 D 3 0 2.200000 0
3 14.066667 D 0 0 .000000 0
Optimal solution found
Branch Parent Solution Entity Value/Bound Active GInf Time
*** Solution found **** Time: 0
1 0 13.100000 0 0 0
*** Search completed *** Time: 0 Nodes: 1
Number of integer feasible solutions found is 1
Best integer solution found is 13.100000
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%
Its Obj Value S Ninf Nneg Sum Inf Time
0 14.066667 D 0 0 .000000 0
Optimal solution found
Branch Parent Solution Entity Value/Bound Active GInf Time
*** Search completed *** Time: 0 Nodes: 7
Problem is integer infeasible
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.
If you have any comments or suggestions about these pages, please send mail to docs@dashoptimization.com.