Inputting and solving a Linear Programming problem
In this chapter we take the example formulated in Chapter 2 and show how to implement the model with BCL. With some extensions to the initial formulation we also introduce input and output functionalities of BCL:
- writing an LP model with BCL,
- data input from file using index sets,
- output facilities of BCL,
- exporting a problem to a matrix file.
Chapter 3 shows how to formulate and solve the same example with Mosel and in Chapter 15 the problem is input and solved directly with Xpress Optimizer.
Implementation with BCL
All BCL examples in this book are written with C++. Due to the possibility of overloading arithmetic operators in the C++ programming language, this interface provides the most convenient way of stating models in a close to algebraic form. The same models can also be implemented using the C, Java, or C# interfaces of BCL.
The following BCL program implements the LP example introduced in Chapter 2:
#include <iostream>
#include "xprb_cpp.h"
using namespace std;
using namespace ::dashoptimization;
#define NSHARES 10 // Number of shares
#define NRISK 5 // Number of high-risk shares
#define NNA 4 // Number of North-American shares
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
int main(int argc, char **argv)
{
int s;
XPRBprob p("FolioLP"); // Initialize a new problem in BCL
XPRBexpr Risk,Na,Return,Cap;
XPRBvar frac[NSHARES]; // Fraction of capital used per share
// Create the decision variables
for(s=0;s<NSHARES;s++) frac[s] = p.newVar("frac");
// 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", Risk <= 1.0/3);
// Minimum amount of North-American values
for(s=0;s<NNA;s++) Na += frac[NA[s]];
p.newCtr("NA", Na >= 0.5);
// Spend all the capital
for(s=0;s<NSHARES;s++) Cap += frac[s];
p.newCtr("Cap", Cap == 1);
// Upper bounds on the investment per share
for(s=0;s<NSHARES;s++) frac[s].setUB(0.3);
// Solve the problem
p.setSense(XPRB_MAXIM);
p.lpOptimize("");
// Solution printing
cout << "Total return: " << p.getObjVal() << endl;
for(s=0;s<NSHARES;s++)
cout << s << ": " << frac[s].getSol()*100 << "%" << endl;
return 0;
}Let us now have a closer look at what we have just written.
Initialization
To use the BCL C++ interface you need to include the header file xprb_cpp.h. We also define the namespace to which the BCL classes belong.
If the software has not been initialized previously, BCL is initialized automatically when the first problem is created, that is by the line
XPRBprob p("FolioLP");which creates a new problem with the name `FolioLP'.
General structure
The definition of the model itself starts with the creation of the decision variables (method newVar), followed by the definition of the objective function and the constraints. In C++ (and Java, C#) constraints may be created starting with linear expressions as shown in the example. Equivalently, they may be constructed termwise, for instance the constraint limiting the percentage of high-risk shares:
XPRBctr CRisk;
CRisk = p.newCtr("Risk");
for(s=0;s<NRISK;s++) CRisk.addTerm(frac[RISK[s]], 1);
CRisk.setType(XPRB_L);
CRisk.addTerm(1.0/3);This second type of constraint definition is common to all BCL interfaces and is the only method of defining constraints in C where overloading is not available.
Notice that in the definition of equality constraints (here the constraint stating that we wish to spend all the capital) we need to employ a double equality sign ==.
The method setUB is used to set the upper bounds on the decision variables frac. Alternatively to this separate function call, we may also specify the bounds directly at the creation of the variables, but in this case we need to provide the full information including the name, variable type (XPRB_PL for continuous), lower and upper bound values:
for(s=0;s<NSHARES;s++) frac[s] = p.newVar("frac", XPRB_PL, 0, 0.3);Giving string names to modeling objects (decision variables, constraints, etc.) as shown in our example program is optional. If the user does not specify any name, BCL will generate a default name. However, user-defined names may be helpful for debugging and for the interpretation of output produced by the Optimizer.
Solving
Prior to launching the solver, the optimization direction is set to maximization with a call to setSense. With the method lpOptimize we then call Xpress Optimizer to maximize the objective function (Return) set with the method setObj, subject to all constraints that have been defined. The empty string argument of lpOptimize indicates that the default LP algorithm is to be used. Other possible values are "p" for primal, "d" for dual Simplex, and "b" for Newton-Barrier.
Output printing
The last few lines print out the value of the optimal solution and the solution values for all variables.
Compilation and program execution
If you have followed the standard installation procedure of Xpress Optimizer and BCL, you may compile this file with the following command under Windows (note that it is important to use the flag /MD):
cl /MD /I%XPRESSDIR%\include %XPRESSDIR%\lib\xprb.lib foliolp.cpp
For Linux or Solaris use
cc -D_REENTRANT -I${XPRESSDIR}/include -L${XPRESSDIR}/lib foliolp.C -o foliolp -lxprbFor other systems please refer to the example makefile provided with the corresponding distribution.
Running the resulting program will generate the following output:
Reading Problem FolioLP
Problem Statistics
3 ( 0 spare) rows
10 ( 0 spare) structural columns
19 ( 0 spare) non-zero elements
Global Statistics
0 entities 0 sets 0 set members
Maximizing LP FolioLP
Original problem has:
3 rows 10 cols 19 elements
Presolved problem has:
3 rows 10 cols 19 elements
Its Obj Value S Ninf Nneg Sum Inf Time
0 42.600000 D 2 0 .000000 0
5 14.066667 D 0 0 .000000 0
Uncrunching matrix
Optimal solution found
Dual solved problem
5 simplex iterations in 0s
Final objective : 1.406666666666667e+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
All values within tolerances
Problem status: optimal
Total return: 14.0667
0: 30%
1: 0%
2: 20%
3: 0%
4: 6.66667%
5: 30%
6: 0%
7: 0%
8: 13.3333%
9: 0%The upper half of this display is the log of Xpress Optimizer: the size of the matrix, 3 rows (i.e. constraints) and 10 columns (i.e. decision variables), and the log of the LP solution algorithm (here: `D' for dual Simplex). The lower half is the output produced by our program: the maximum return of 14.0667 is obtained with a portfolio consisting of shares 1, 3, 5, 6, and 9. 30% of the total amount are spent in shares 1 and 6 each, 20% in 3, 13.3333% in 9 and 6.6667% in 5. It is easily verified that all constraints are indeed satisfied: we have 50% of North-American shares (1 and 3) and 33.33% of high-risk shares (3 and 9).
It is possible to modify the amount of output printing by BCL and Xpress Optimizer by adding the following line before the start of the optimization:
p.setMsgLevel(1);
This setting will disable all output (including warnings) from BCL and Xpress Optimizer, with the exception of error messages. The possible values for the printing level range from 0 to 4. In Chapter 13 we show how to access the Optimizer control parameters directly which, for instance, allows fine tuning the message display.
Data input from file
Instead of simply numbering the decision variables, we may wish to use more meaningful indices in our model. For instance, the problem data may be given in file(s) using string indices such as the file foliocpplp.dat we now wish to read:
! Return "treasury" 5 "hardware" 17 "theater" 26 "telecom" 12 "brewery" 8 "highways" 9 "cars" 7 "bank" 6 "software" 31 "electronics" 21
We modify our previous model as follows to work with this data file:
#include <iostream>
#include "xprb_cpp.h"
using namespace std;
using namespace ::dashoptimization;
#define DATAFILE "foliocpplp.dat"
#define NSHARES 10 // Number of shares
#define NRISK 5 // Number of high-risk shares
#define NNA 4 // Number of North-American shares
double RET[NSHARES]; // Estimated return in investment
char RISK[][100] = {"hardware", "theater", "telecom", "software",
"electronics"}; // High-risk values among shares
char NA[][100] = {"treasury", "hardware", "theater", "telecom"};
// Shares issued in N.-America
XPRBindexSet SHARES; // Set of shares
XPRBprob p("FolioLP"); // Initialize a new problem in BCL
void readData(void)
{
double value;
int s;
FILE *datafile;
char name[100];
SHARES=p.newIndexSet("Shares",NSHARES); // Create the `SHARES' index set
// Read `RET' data from file
datafile=fopen(DATAFILE,"r");
for(s=0;s<NSHARES;s++)
{
XPRBreadlinecb(XPRB_FGETS, datafile, 200, "T g", name, &value);
RET[SHARES+=name]=value;
}
fclose(datafile);
SHARES.print(); // Print out the set contents
}
int main(int argc, char **argv)
{
int s;
XPRBexpr Risk,Na,Return,Cap;
XPRBvar frac[NSHARES]; // Fraction of capital used per share
// Read data from file
readData();
// Create the decision variables
for(s=0;s<NSHARES;s++) frac[s] = p.newVar("frac");
// 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[SHARES[RISK[s]]];
p.newCtr("Risk", Risk <= 1.0/3);
// Minimum amount of North-American values
for(s=0;s<NNA;s++) Na += frac[SHARES[NA[s]]];
p.newCtr("NA", Na >= 0.5);
// Spend all the capital
for(s=0;s<NSHARES;s++) Cap += frac[s];
p.newCtr("Cap", Cap == 1);
// Upper bounds on the investment per share
for(s=0;s<NSHARES;s++) frac[s].setUB(0.3);
// Solve the problem
p.setSense(XPRB_MAXIM);
p.lpOptimize("");
// Solution printing
cout << "Total return: " << p.getObjVal() << endl;
for(s=0;s<NSHARES;s++)
cout << SHARES[s] << ": " << frac[s].getSol()*100 << "%" << endl;
return 0;
}The arrays RISK and NA now store indices in the form of strings and we have added a new object, the index set SHARES that is defined while the return-on-investment values RET are read from the data file. In this example we have initialized the index set with exactly the right size. This is not really necessary since index sets may grow dynamically if more entries are added to them than the initially allocated space. The actual set size can be obtained with the method getSize.
For reading the data we use the function XPRBreadlinecb. It will skip comments preceded by ! and any empty lines in the data file. The format string "T g" indicates that we wish to read a text string (surrounded by single or double quotes if it contains blanks) followed by a real number, the two separated by spaces (including tabulation). If the data file used another separator sign such as `,' then the format string could be changed accordingly (e.g. "T,g").
In the model itself, the definition of the linear expressions Risk and Na has been adapted to the new indices.
Another modification concerns the solution printing: we print the name of every share, not simply its sequence number and hence the solution now gets displayed as follows:
Total return: 14.0667 treasury: 30% hardware: 0% theater: 20% telecom: 0% brewery: 6.66667% highways: 30% cars: 0% bank: 0% software: 13.3333% electronics: 0%
Output functions and error handling
Most BCL modeling objects (XPRBprob, XPRBvar, XPRBctr, XPRBsos, and XPRBindexSet) have a method print. For variables, depending on where the method is invoked either their bounds or their solution value is printed: adding the line
frac[2].print();
before the call to the optimization will print the name of the variable and its bounds,
frac2: [0,0.3]
whereas after the problem has been solved its solution value gets displayed:
frac2: 0.2
Whenever BCL detects an error it stops the program execution with an error message so that it will usually not be necessary to test the return value of every operation. If a BCL program is embedded into some larger application, it may be helpful to use the explicit initialization to check early on that the software can be accessed correctly, for instance:
if(XPRB::init() != 0)
{
cout << "Initialization failed." << endl;
return 1;
}The method getLPStat may be used to test the LP problem status. Only if the LP problem has been solved successfully BCL will return or print out meaningful solution values:
char *LPSTATUS[] = {"not loaded", "optimal", "infeasible",
"worse than cutoff", "unfinished", "unbounded",
"cutoff in dual", "unsolved", "nonconvex"};
cout << "Problem status: " << LPSTATUS[p.getLPStat()] << endl;Exporting matrices
If the optimization process with Xpress Optimizer is started from within a BCL program (methods lpOptimize or mipOptimize of XPRBprob), then the problem matrix is loaded in memory into the solver without writing it out to a file (which would be expensive in terms of running time). However, in certain cases it may still be required to be able to produce a matrix. With Xpress, the user has the choice between two matrix formats: extended MPS and extended LP format, the latter being in general more easily human-readable since constraints are printed in algebraic form.
To export a matrix in MPS format add the following line to your BCL program, immediately before or instead of the optimization statement; this will create the file Folio.mat in your working directory:
p.exportProb(XPRB_MPS, "Folio");
For an LP format matrix use the following:
p.setSense(XPRB_MAXIM); p.exportProb(XPRB_LP, "Folio");
The LP format contains information about the sense of optimization. Since the default is to minimize, for maximization we first need to reset the sense. The resulting matrix file will have the name Folio.lp.