Implementation with Mosel

Topics covered in this section:

For the implementation (file folioheur.mos) of the variable fixing solution heuristic we work with the MIP 1 model from Chapter 6. Through the definition of the heuristic in the form of a subroutine (more precisely, a procedure) we only make minimal changes to the model itself: at the beginning we declare the procedure using the keyword forward, and before solving our problem with the standard call to the maximization function we execute our own solution heuristic. The solution printing also has been adapted.

model "Portfolio optimization solved heuristically"
 uses "mmxprs"                      ! Use Xpress Optimizer

 parameters
  MAXRISK = 1/3                     ! Max. investment into high-risk values
  MAXVAL = 0.3                      ! Max. investment per share
  MINAM = 0.5                       ! Min. investment into N.-American values
  MAXNUM = 4                        ! Max. number of assets
 end-parameters

 forward procedure solve_heur       ! Heuristic solution procedure

 declarations
  SHARES: set of string             ! Set of shares
  RISK: set of string               ! Set of high-risk values among shares
  NA: set of string                 ! Set of shares issued in N.-America
  RET: array(SHARES) of real        ! Estimated return in investment
 end-declarations

 initializations from "folio.dat"
  RISK RET NA
 end-initializations

 declarations
  frac: array(SHARES) of mpvar      ! Fraction of capital used per share
  buy: array(SHARES) of mpvar       ! 1 if asset is in portfolio, 0 otherwise
 end-declarations

! Objective: total return
 Return:= sum(s in SHARES) RET(s)*frac(s)

! Limit the percentage of high-risk values
 sum(s in RISK) frac(s) <= MAXRISK

! Minimum amount of North-American values
 sum(s in NA) frac(s) >= MINAM

! Spend all the capital
 sum(s in SHARES) frac(s) = 1

! Upper bounds on the investment per share
 forall(s in SHARES) frac(s) <= MAXVAL

! Limit the total number of assets
 sum(s in SHARES) buy(s) <= MAXNUM

 forall(s in SHARES) do
  buy(s) is_binary
  frac(s) <= buy(s)
 end-do

! Solve problem heuristically
 solve_heur

! Solve the problem
 maximize(Return)

! Solution printing
 if getprobstat=XPRS_OPT then
  writeln("Exact solution: Total return: ", getobjval)
  forall(s in SHARES) writeln(s, ": ", getsol(frac(s))*100, "%")
 else
  writeln("Heuristic solution is optimal.")
 end-if

!-----------------------------------------------------------------

 procedure solve_heur
  declarations
   TOL: real                       ! Solution feasibility tolerance
   fsol: array(SHARES) of real     ! Solution values for `frac' variables
   bas: basis                      ! LP basis
  end-declarations

  setparam("XPRS_VERBOSE",true)    ! Enable message printing in mmxprs
  setparam("XPRS_CUTSTRATEGY",0)   ! Disable automatic cuts
  setparam("XPRS_HEUREMPHASIS",0)  ! Disable automatic MIP heuristics
  setparam("XPRS_PRESOLVE",0)      ! Switch off presolve
  TOL:=getparam("XPRS_FEASTOL")    ! Get feasibility tolerance
  setparam("ZEROTOL",TOL)          ! Set comparison tolerance

  maximize(XPRS_LPSTOP,Return)     ! Solve the LP problem
  savebasis(bas)                   ! Save the current basis

 ! Fix all variables `buy' for which `frac' is at 0 or at a relatively
 ! large value
  forall(s in SHARES) do
   fsol(s):= getsol(frac(s))       ! Get the solution values of `frac'
   if (fsol(s) = 0) then
    setub(buy(s), 0)
   elif (fsol(s) >= 0.2) then
    setlb(buy(s), 1)
   end-if
  end-do

  maximize(XPRS_CONT,Return)       ! Solve the MIP problem
  ifgsol:=false
  if getprobstat=XPRS_OPT then     ! If an integer feas. solution was found
   ifgsol:=true
   solval:=getobjval               ! Get the value of the best solution
   writeln("Heuristic solution: Total return: ", solval)
   forall(s in SHARES) writeln(s, ": ", getsol(frac(s))*100, "%")
  end-if

 ! Reset variables to their original bounds
  forall(s in SHARES)
   if ((fsol(s) = 0) or (fsol(s) >= 0.2)) then
    setlb(buy(s), 0)
    setub(buy(s), 1)
   end-if

  loadbasis(bas)                   ! Load the saved basis

  if ifgsol then                   ! Set cutoff to the best known solution
   setparam("XPRS_MIPABSCUTOFF", solval+TOL)
  end-if
 end-procedure

end-model

This model certainly requires some more detailed explanations.

Subroutines

A subroutine in Mosel has a similar structure as the model itself: a procedure starts with the keyword procedure, followed by the name of the procedure, and terminates with end-procedure. Similarly, a function starts with the keyword function, followed by its name, and terminates with end-function. Both types of subroutines may take a list of arguments and for functions in addition the return type must be indicated, for example:

 function myfunc(myint: integer, myarray: array(range) of string): real

for a function that returns a real and takes as input arguments an integer and an array of string.

As shown in our example, a subroutine may contain one (or several) declarations blocks. The objects defined in a subroutine are only valid locally and are deleted at the end of the subroutine.

Subroutine definitions may be overloaded, that is, a single subroutine may take different combinations of arguments. It is possible to overload any subroutines defined by Mosel and its modules, provided that the new definition differs from the existing one(s) in at least one argument.

For more detail and further examples of subroutine definition see the `Mosel User Guide'.

Optimizer parameters and functions

Parameters: The solution heuristic starts with parameter settings for Xpress Optimizer. For a detailed explanation of all Optimizer parameters the reader is refered to the `Optimizer Reference Manual'. All parameters are accessed through the Mosel subroutines setparam and getparam. In the example, we first enable the output printing by the module mmxprs. As a result, more information than what is printed by our model will be displayed in the logging pane:

Figure 9.1: Optimizer output display

Switching off the automated cut generation (parameter XPRS_CUTSTRATEGY) and the MIP heuristics (parameter XPRS_HEUREMPHASIS) 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.

Optimization statement: We use a new version of the maximization procedure with an additional argument, XPRS_LPSTOP, 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 XPRS_CONT. This is an example of an overloaded subroutine definition.

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 loadprob) 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.

For more detail on the Optimizer functionality used in this example see the documentation of the module mmxprs in the `Mosel Language Reference Manual'.

Comparison tolerance

After retrieving the feasibility tolerance of the Optimizer we set the comparison tolerance of Mosel (ZEROTOL) to this value TOL. This means that the test fsol(s) = 0 evaluates to true if fsol(s) lies between -TOL and TOL, and fsol(s) >= 0.2 is satisfied if the value of fsol(s) is at least 0.2-TOL.

Comparisons in Mosel always use a tolerance, with a very small default value. By resetting this parameter to the Optimizer feasibility tolerance Mosel evaluates solution values just like the Optimizer.