Implementation with Java

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. In the code snippet below we highlight the main changes in relation to the MIP model in Chapter 11:

public class FolioHeuristic {

    ... // declare and initialize data objects

    public XpressProblem prob;
    /* Fraction of capital used per share */
    public Variable[] frac;
    /* 1 if asset is in portfolio, 0 otherwise */
    public Variable[] buy;

    public FolioHeuristic(XpressProblem p) {
        prob = p;
        /**** VARIABLES ****/
        frac = prob.addVariables(NSHARES)
                /* Fraction of capital used per share */
                .withName(i -> String.format("frac_%d", i))
                /* Upper bounds on the investment per share */
                .withUB(0.3).toArray();

        buy = prob.addVariables(NSHARES)
                /* Fraction of capital used per share */
                .withName(i -> String.format("buy_%d", i))
                .withType(ColumnType.Binary).toArray();

    ... // print functions and model definition as in Chapter 11

    private static void solveHeuristic(FolioHeuristic folio) {
        XpressProblem p = folio.prob;
        /* Disable automatic cuts */
        p.controls().setCutStrategy(XPRSconstants.CUTSTRATEGY_NONE);
        // Switch presolve off
        p.controls().setPresolve(XPRSconstants.PRESOLVE_NONE);
        p.controls().setMIPPresolve(0);
        /* Get feasibility tolerance */
        double tol = p.controls().getFeasTol();

        /* Solve the LP-problem */
        p.lpOptimize();

        /* Get Solution */
        double[] sol = p.getSolution();

        /* Basis information */
        int[] rowstat = new int[p.attributes().getRows()];
        int[] colstat = new int[p.attributes().getCols()];
        /* Save the current basis */
        p.getBasis(rowstat, colstat);

        /*
         * Fix all variables `buy' for which `frac' is at 0 or at a relatively large
         * value
         */
        double[] fsol = new double[NSHARES];
        range(0, NSHARES).forEach(i -> {
            /* Get the solution values of `frac' */
            fsol[i] = folio.frac[i].getValue(sol);
            if (fsol[i] < tol)
                folio.buy[i].fix(0);
            else if (fsol[i] > 0.2 - tol)
                folio.buy[i].fix(1);
        });

        /* Solve with the new bounds on 'buy' */
        p.mipOptimize();

        printProblemStatus(p);
        printProblemSolution(folio, "Heuristic solution");

        /* Reset variables to their original bounds */
        range(0, NSHARES).forEach(i -> {
            if ((fsol[i] < tol) || (fsol[i] > 0.2 - tol)) {
                folio.buy[i].setLB(0);
                folio.buy[i].setUB(1);
            }
        });

        /* Load basis */
        p.loadBasis(rowstat, colstat);
    }

    public static void main(String[] args) {
        try (XpressProblem prob = new XpressProblem()) {
            /* Solve with heuristic */
            FolioHeuristic folio = new FolioHeuristic(prob);
            model(folio);
            solveHeuristic(folio);
        }

        try (XpressProblem prob = new XpressProblem()) {
            FolioHeuristic folio = new FolioHeuristic(prob);
            model(folio);
            // Solve
            folio.prob.optimize();
            // Solution printing
            printProblemStatus(prob);
            printProblemSolution(folio, "Exact Solve");
        }
    }
}

The implementation of the heuristic certainly requires some explanation.

Parameters: The solution heuristic starts with parameter settings for the Xpress Optimizer. Switching off the automated cut generation (parameter CUTSTRATEGY_NONE) is optional. However, 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 PRESOLVE_NONE), because we interact with the problem in the Optimizer in the course of its solution, and this can be done 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 example, by performing comparisons), it is important to take into account these tolerances.

Optimization calls: We use the optimization method lpOptimize, indicating that we only want to solve the top node LP relaxation (and not yet the entire MIP problem).

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, after 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 loadBasis) 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:

Maximizing LP  using up to 20 threads and up to 31GB memory, with these control settings:

...

Optimal solution found
Dual solved problem
  5 simplex iterations in 0.00 seconds at time 0

Final objective                       : 1.406666666666666e+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

...

Maximizing MILP  using up to 20 threads and up to 31GB memory, with these control settings:

...

 *** Search completed ***
Final MIP objective                   : 1.310000000000000e+01
Final MIP bound                       : 1.310001310000000e+01
  Solution time / primaldual integral :      0.01s/ 32.322812%
  Number of solutions found / nodes   :         1 /         3
  Max primal violation      (abs/rel) :       0.0 /       0.0
  Max integer violation     (abs    ) :       0.0
Problem status:
	Solve status: Completed
	LP status: Optimal
	MIP status: Optimal
	Sol status: Optimal
Total return (Heuristic solution): 13.099999999999998

...

Maximizing MILP  using up to 20 threads and up to 31GB memory, with default controls

...

 *** Search completed ***
Uncrunching matrix
Final MIP objective                   : 1.310000000000000e+01
Final MIP bound                       : 1.310001310000000e+01
  Solution time / primaldual integral :      0.00s/ 41.952984%
  Number of solutions found / nodes   :         1 /         1
  Max primal violation      (abs/rel) : 5.551e-17 / 5.551e-17
  Max integer violation     (abs    ) :       0.0
Problem status:
	Solve status: Completed
	LP status: CutOffInDual
	MIP status: Optimal
	Sol status: Optimal
Total return (Exact Solve): 13.1

This output shows that the heuristic found a solution of 13.1, and that the MIP optimizer without the heuristic could not find a better solution. The heuristic solution is therefore optimal.