Implementation with Python

For the implementation of the variable fixing solution heuristic we work with the MIP 1 model from Chapter Mixed Integer Programming. 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 Mixed Integer Programming:

import xpress as xp

# Problem data
MAXNUM = 4
NSHARES = 10
RET = [5, 17, 26, 12, 8, 9, 7, 6, 31, 21]
RISK = [1, 2, 3, 8, 9]
NA = [0, 1, 2, 3]

def printSolution(prob, name):
    # Solution printing.
    print(f"Total return {name}:", prob.attributes.objval)
    sol = prob.getSolution()
    for i in range(NSHARES):
        print(f"{frac[i].name} : {sol[i] * 100:.2f} %")

def solveHeuristic(prob):
    # Disable automatic cuts.
    prob.controls.cutstrategy = 0
    # Switch presolve off.
    prob.controls.presolve = 0
    prob.controls.mippresolve = 0
    # Get feasibility tolerance.
    tol = prob.controls.feastol

    prob.lpoptimize()

    # Save the current basis.
    rowstat = []
    colstat = []
    prob.getbasis(rowstat,colstat)

    # Fix all variables 'buy' for which `frac' is at 0 or at a relatively large value
    fsol = prob.getSolution(frac)      # get the solution values of `frac'
    for i in range(NSHARES):
        if fsol[i] < tol:
            buy[i].lb = 0
            buy[i].ub = 0
        elif fsol[i] > 0.2 - tol:
            buy[i].lb = 1
            buy[i].ub = 1

    prob.mipoptimize()

    print(f"Problem status: \n\t Solve status: {p.attributes.solvestatus.name} \n\t Sol status: \
        {p.attributes.solstatus.name}")

    printSolution(prob, "Heuristic solution")

    # Reset variables to their original bounds.
    for i in range(NSHARES):
        if fsol[i] < tol or fsol[i] > 0.2 - tol:
            idx = prob.getIndex(buy[i])
            buy[i].lb = 0
            buy[i].ub = 1

    # Load basis.
    prob.loadbasis(rowstat, colstat)

    # Set cutoff to the best known solution.
    prob.controls.mipabscutoff = prob.attributes.objval - tol

p = xp.problem(name="Folio")

# VARIABLES.
frac = p.addVariables(NSHARES, ub=0.3, name="frac")
buy = p.addVariables(NSHARES, vartype=xp.binary, name="buy")

# CONSTRAINTS.
# Limit the percentage of high-risk values.
p.addConstraint(xp.Sum(frac[i] for i in RISK) <= 1/3)

# Minimum amount of North-American values.
p.addConstraint(xp.Sum(frac[i] for i in NA) >= 0.5)

# Spend all the capital.
p.addConstraint(xp.Sum(frac) == 1)

# Limit the total number of assets.
p.addConstraint(xp.Sum(buy) <= MAXNUM)

# Linking the variables.
p.addConstraint(frac[i] <= buy[i] for i in range(NSHARES))

# Objective: maximize total return.
p.setObjective(xp.Sum(frac[i] * RET[i] for i in range(NSHARES)), sense=xp.maximize)

# Solve with heuristic.
solveHeuristic(p)

# Solve original problem.
p.optimize()

# Print problem status.
print(f"Problem status: \n\t Solve status: {p.attributes.solvestatus.name} \n\t Sol status: \
    {p.attributes.solstatus.name}")

printSolution(p, "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) 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), because we interact with the problem in the course of its solution, and this can be done correctly if the matrix has not been modified by Xpress.

In addition to the parameter settings we also retrieve the feasibility tolerance used by Xpress: 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 Xpress (e.g. after executing an optimization statement or following an explicit call to method loadbasis) bound changes via lb and ub are passed on directly to Xpress.

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 solver without the heuristic could not find a better solution. The heuristic solution is therefore optimal.