Using lpsolve from AMPL

AMPL?

AMPL (A Mathematical Programming Language) is a high-level language for describing mathematical programs. AMPL allows a mathematical programming model to be specified independently of the data used for a specific instance of the model. AMPL's language for describing mathematical programs closely follows that used by humans to describe mathematical programs to each other. For this reason, modellers may spend more time improving the model and less time on the tedious details of data manipulation and problem solution.
To start, AMPL needs a mathematical programming model, which describes variables, objectives and relationships without referring to specific data. AMPL also needs an instance of the data, or a particular data set. The model and one (or more) data files are fed into the AMPL program. AMPL works like a compiler: the model and input are put into an intermediate file that can be read by a solver. The solver actually finds an optimal solution to the problem by reading in the intermediate file produced by AMPL and applying an appropriate algorithm. The solver outputs the solution as a text file, which can be viewed directly and cross-referenced with the variables and constraints specified in the model file.

We will not discuss the specifics of AMPL here but instead refer the reader to AMPL (A Mathematical Programming Language) at the University of Michigan Documentation and the AMPL website http://www.ampl.com.

AMPL and lpsolve

One of the possible solvers that can be used by AMPL is lpsolve. However note that AMPL allows also defining non-linear models. These are not solvable by lpsolve because lpsolve can only handle MILP models. A message will be given if the AMPL model results in a model that cannot be handled by lpsolve.

To make this possible, a driver program is needed: lpsolve(.exe). This program must be put in the AMPL directory and AMPL can call the lpsolve solver. That is all.

Solve a model from AMPL via lpsolve

In the following text, ampl: before the AMPL commands is the AMPL prompt. Only the text after the : must be entered.

To select the lpsolve solver, the following command must be executed in AMPL:

ampl: option solver lpsolve;

To solve the model, the following command must be executed in AMPL:

ampl: solve;

Options can be passed to lpsolve. For example to enable scaling, the following command must be executed in AMPL:

ampl: option lpsolve_options 'scale';

Multiple options can be specified by separating them with at least one space:

ampl: option lpsolve_options 'scale scalemode=7 verbose';

A list of all options is given at the end of this document.

An example

ampl: model diet.mod;
ampl: data diet.dat;
ampl: option solver lpsolve;
ampl: option lpsolve_options 'scale scalemode=7';
ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: scale
scalemode=7
LP_SOLVE 5.1.0.6: optimal, objective 88.2
1 simplex iterations

solve_result_num

AMPL has a parameter that is used to indicate the outcome of the optimisation process. It is used as follows

ampl: display solve_result_num

solve_result_num can take the values shown in following table which also presents a short explanation for each value.

Interpretation of solve_result_num.
Value Message
0 the solution is optimal
100 suboptimal solution
200 infeasible problem
300 unbounded problem
500 solution status is unknown
501 failure
502 integer programming failure

Restart

This is also new with version 5. Before, when lpsolve was used from AMPL, the model was each time completely resolved at each AMPL solve command. The last base was not used to do further iterations. This is specific to AMPL because lpsolve is called as a seperate program each time a solve is done and to make restart possible the last base must be provided again to lpsolve to start from there again. From version 5 on this is now implemented by default. This means if a solve is done immediately after a previous solve that no iterations are needed and the solution can be given immediately. This can be shown via an example:
ampl: model diet.mod;
ampl: data diet.dat;
ampl: option solver lpsolve;
ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: optimal, objective 88.2
3 simplex iterations

Note that 3 iterations were needed to solve the model. Now solve the model again:

ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: optimal, objective 88.2
0 simplex iterations

Note now that no iterations are needed to solve the model. This because the result of the previous solve is used as a starting point for the new solve. It is even allowed to add or deleted variables and constraints in the model, lpsolve will still try to start from the last result to continue.

If, for some reason, you don't want this and let lpsolve resolve the model from scratch, there are two possibilities to let lpsolve ignore the starting base. Either via an lpsolve option or an AMPL option.

The lpsolve option that can be set is:

ampl: option lpsolve_options 'coldstart';

The AMPL option that can be set is:

ampl: option send_statuses 0;

One of both is enough.

Now solve the model again:
ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: coldstart
LP_SOLVE 5.1.0.6: optimal, objective 88.2
3 simplex iterations

Note that again 3 iterations are needed to solve the model.

Special considerations for integer models. From the moment that a model contains integer variables, the B&B algorithm must be used to solve this. This algorithm must go trough a tree of possible solutions each time a solve is done. The last base of the best integer solution can not be used as a starting base for a resolve. This is not an lpsolve limitation, but a B&B algorithm property...
For this reason, by default, AMPL doesn't provide a starting to the solver and solve is done from scratch. By default... However, a compromise can be implemented by the solver. And this is also done in lpsolve. Solving integer models is always started with solving the non-integer model. When this is done, the non-integer variables are made integer via the B&B algorithm. lpsolve returns the base of the non-integer model and this base is used as a starting base for the model at the next solve. As said, this is a compromise because it results in a fast solve of the non-integer model, but the B&B algorithm must still be executed. So even if no modifications are done to the model and solve is done again, there will be iterations needed by the B&B algorithm to solve the integers. Because this must be explicitly implemented in the solver, AMPL doesn't provide the starting base by default for integer models. It must explicitly be activated:

ampl: option send_statuses 2;

Only then, lpsolve will get a starting base when there are integers.

Example:

ampl: model multmip3.mod;
ampl: data multmip3.dat;
ampl: option solver lpsolve;
ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: optimal, objective 235625
1986 simplex iterations
592 branch & bound nodes: depth 21

Note that 1986 iterations were needed to solve the model. Now solve the model again:

ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: optimal, objective 235625
1986 simplex iterations
592 branch & bound nodes: depth 21

The same as above. So no restart was done. Now activate the starting base for integer models and solve again:

ampl: option send_statuses 2;
ampl: solve;

This gives as result:

LP_SOLVE 5.1.0.6: optimal, objective 235625
1922 simplex iterations
592 branch & bound nodes: depth 21

Note that less iterations are needed, but not 0 as in non-integer models. These are the iterations needed for the B&B algorithm.

lpsolve options

lpsolve accepts a lot of options. Some options are just on/off switches and for there you just specify the option keyword (for example scale). Other options need a value and this is specified by an equal (=) sign after the option keyword and then the value (for example scalemode=7). Multiple options can be specified by separating them by a space (for example scale scalemode=7).

Here is a list of the possible lpsolve options:

bb=...         branch-and-bound rule: one of
  		    0 (default) for lowest indexed non-integer variable
  		    1 for selection based on distance from the current bounds
  		    2 for selection based on the largest current bound
  		    3 for selection based on largest fractional value
  		    4 for simple, unweighted pseudo-cost of a variable
  		    5 this is an extended pseudo-costing strategy based on
                      minimizing the number of integer infeasibilities
  		    6 this is an extended pseudo-costing strategy based on
                      maximizing the normal pseudo-cost divided by the
                      number of infeasibilities.
  	       plus
  		   32 for greedy mode
  		   64 for pseudo cost mode
  		  128 select the node that has already been selected before
                      the most number of times
  		 4096 select the node that has been selected before the
                      fewest number of times or not at all
bfp=...        set basis factorization package
cauto          in IPs, algorithm decides which branch being taken first
cfirst         in IPs, take ceiling branch first
coldstart      ignore starting base
crash=...      determines a starting base: one of
  		    0 (default) none
  		    2 most feasible basis
debug          debug mode
degen          perturb degeneracies
degenx=...     anti-degen handling: one of
  		    0 (default) no anti-degeneracy handling
  		    1 check if there are equality slacks in the basis and
                      try to drive them out in order to reduce chance of
                      degeneracy in Phase 1
                    2 ColumnCheck
                    4 Stalling
                    8 NumFailure
                   16 LostFeas
                   32 Infeasible
                   64 Dynamic
                  128 During BB
depth=...      set branch-and-bound depth limit
dual           prefer the dual simplex for both phases
eps=...        tolerance for rounding to integer
epsb=...       minimum tolerance for the RHS
epsd=...       minimum tolerance for reduced costs
epsel=...      minimum tolerance for rounding values to zero
epsp=...       value that is used as perturbation scalar for
               degenerative problems
f              specifies that branch-and-bound algorithm stops at first
               found solution
ga=...         specifies the absolute MIP gap for branch-and-bound
gr=...         specifies the relative MIP gap for branch-and-bound
improve=...    the iterative improvement level: one of
  		    0 (default) improve none
  		    1 improve FTRAN
  		    2 improve BTRAN
  		    3 improve FTRAN + BTRAN
                    4 triggers automatic inverse accuracy control in the dual
                      simplex, and when an error gap is exceeded the basis is
                      reinverted
n=...          specify which solution number to return
o=...          specifies that branch-and-bound algorithm stops when
               objective value is better than value
objno=...      objective number: 0 = none, 1 (default) = first
obound=...     a lower bound for the objective function.
               may speed up the calculations
parse_only     parse input file but do not solve
piv=...        simplex pivot rule: one of
  		    0 select first
  		    1 select according to Dantzig
  		    2 (default) select Devex pricing from Paula Harris
  		    3 select steepest edge
piva           temporarily use first index if cycling is detected
pivf           in case of Steepest Edge, fall back to DEVEX in primal
pivla          scan entering/leaving columns alternatingly left/right
pivll          scan entering/leaving columns left rather than right
pivm           multiple pricing
pivr           adds a small randomization effect to the selected pricer
presolve       presolve problem before start optimizing
presolvel      also eliminate linearly dependent rows
presolver      if the phase 1 solution process finds that a constraint is
               redundant then this constraint is deleted
prim           prefer the primal simplex for both phases
printsol=...   print solution: one of
  		    0 (default) print nothing
  		    1 only objective value
  		    2 obj value+variables
  		    3 obj value+variables+constraints
  		    4 obj value+variables+constraints+duals
  		    5 obj value+variables+constraints+duals+lp model
  		    6 obj value+variables+constraints+duals+lp model+scales
  		    7 obj value+variables+constraints+duals+lp model+scales+
                      lp tableau
prlp           print the LP
psols          print (intermediate) feasible solutions
psolsa         print (intermediate) feasible solutions (non-zeros)
r=...          max nbr of pivots between a re-inversion of the matrix
scale          scale the problem
scalemode=...  scale mode:  one of
  		    1 for scale to convergence using largest absolute value,
  		    2 for scale based on the simple numerical range, or
  		    3 (default) for numerical range-based scaling, or
  		    4 for geometric scaling, or
  		    7 for Curtis & Reid scaling
  	       plus
  		   16 for scale to convergence using logarithmic mean
                      of all values
  		   32 for also do Power scaling
  		   64 to make sure that no scaled number is above 1
  		  128 to scale integer variables
simplexdd      set Phase1 Dual, Phase2 Dual
simplexdp      set Phase1 Dual, Phase2 Primal
simplexpd      set Phase1 Primal, Phase2 Dual
simplexpp      set Phase1 Primal, Phase2 Primal
timeout=...    timeout after sec seconds when not solution found
trace          trace pivot selections
trej=...       minimum pivot value
verbose        verbose mode
version        report version details
wafter         write model after solve (usefull if presolve used)
wantsol=...    solution report without -AMPL: sum of
		1 ==> write .sol file
		2 ==> print primal variable values
		4 ==> print dual variable values
		8 ==> do not print solution message
wfmps=...      write to MPS file in free format
wlp=...        write to LP filename
wmps=...       write to MPS file in fixed format
wxli=...       write file with xli library
wxliname=...   xli library
wxliopt=...    options for xli library

lpsolve command line options

Normally, the lpsolve program is called from AMPL. However the program can also be called stand-alone. When the command is invoked without an option or with the option -?, a list of the possible options is shown:

usage: lpsolve [options] stub [-AMPL] [<assignment> ...]

Options:
        --  {end of options}
        -=  {show name= possibilities}
        -?  {show usage}
        -e  {suppress echoing of assignments}
        -s  {write .sol file (without -AMPL)}
        -v  {just show version}

stub and the -AMPL option is used when the program is called from AMPL.
The -v option shows the version of lpsolve.
The -= option shows all the options that can be passed to lpsolve from within AMPL.
These options can also be specified in an environment variable lpsolve_options.

To call the lpsolve command to solve a model, you first must have a stub file that can be read by lpsolve. This stub file can be created by AMPL as follows:

ampl: model models\diet.mod;
ampl: data models\diet.dat;
ampl: write bdiet;
ampl: quit

This created a binary file diet.nl. As an alternative the command write gdiet; can be used. This creates an ascii file that also can be read. However it is somewhat slower to read, especially when the models are larger.

This creates a file diet.nl in the current directory. This file is the stub file needed by lpsolve:

lpsolve lpsolve diet.nl

This gives:

LP_SOLVE 5.1.0.6: optimal, objective 88.19999999999999
3 simplex iterations

Options can be passed to lpsolve via the environment variable lpsolve_options:

set lpsolve_options=scale scalemode=7 verbose

Then the result of the previous lpsolve command is:

scale
scalemode=7
verbose
Model name:  '' - run #1
Objective:   Minimize(R0)

Submitted:
Model size:        4 constraints,       8 variables,           39 non-zeros.
Constraints:       0 equality,          0 GUB,                  0 SOS.
Variables:         0 integer,           0 semi-cont.,           0 SOS.

Using DUAL simplex for phase 1 and PRIMAL simplex for phase 2.

Optimal solution with dual simplex at iteration          1

lp_solve solution         88.2 final at iteration        1,       0 nodes explor
ed

Excellent numeric accuracy ||*|| = 1.13687e-013

 Memo: Largest [etaPFI v1.0] inv(B) had 0 NZ entries, 0.0x largest basis.
      In the total iteration count 1, 0 (0.0%) were minor/bound swaps.
      There were 0 refactorizations, 0 triggered by time and 0 by density.
             ... on average 1.0 major pivots per refactorization.
      Total solver time was 0.000 seconds.
LP_SOLVE 5.1.0.6: optimal, objective 88.19999999999999
1 simplex iterations

In this mode, the lpsolve option wantsol shows the solution:

set lpsolve_options=wantsol=2
lpsolve diet.nl

This gives:

wantsol=2
LP_SOLVE 5.1.0.6: optimal, objective 88.19999999999999
3 simplex iterations

variable  value
_svar[1]  0
_svar[2]  0
_svar[3]  0
_svar[4]  0
_svar[5]  46.666666666666664
_svar[6]  0
_svar[7]  -3.552713678800501e-15
_svar[8]  0

The -e option results in not echoing the options passed to lpsolve:

set lpsolve_options=wantsol=2
lpsolve -e diet.nl

This gives:

LP_SOLVE 5.1.0.6: optimal, objective 88.19999999999999
3 simplex iterations

variable  value
_svar[1]  0
_svar[2]  0
_svar[3]  0
_svar[4]  0
_svar[5]  46.666666666666664
_svar[6]  0
_svar[7]  -3.552713678800501e-15
_svar[8]  0

lpsolve options can also be passed at the command line, they don't overrule the lpsolve_options environment variable, they are added:

set lpsolve_options=wantsol=2
lpsolve diet.nl "verbose scale"

This gives:

wantsol=2
verbose
scale
Model name:  '' - run #1
Objective:   Minimize(R0)

Submitted:
Model size:        4 constraints,       8 variables,           39 non-zeros.
Constraints:       0 equality,          0 GUB,                  0 SOS.
Variables:         0 integer,           0 semi-cont.,           0 SOS.

Using DUAL simplex for phase 1 and PRIMAL simplex for phase 2.

Optimal solution with dual simplex at iteration          1

lp_solve solution         88.2 final at iteration        1,       0 nodes explor
ed

Excellent numeric accuracy ||*|| = 1.13687e-013

 Memo: Largest [etaPFI v1.0] inv(B) had 0 NZ entries, 0.0x largest basis.
      In the total iteration count 1, 0 (0.0%) were minor/bound swaps.
      There were 0 refactorizations, 0 triggered by time and 0 by density.
             ... on average 1.0 major pivots per refactorization.
      Total solver time was 0.000 seconds.
LP_SOLVE 5.1.0.6: optimal, objective 88.19999999999997
1 simplex iterations

variable  value
_svar[1]  0
_svar[2]  0
_svar[3]  0
_svar[4]  0
_svar[5]  46.66666666666666
_svar[6]  0
_svar[7]  0
_svar[8]  0