Delivery - Data input from file; infeasibility analysis
|
|
|
| Type: | Network flow |
| Rating: | 2 (easy-medium) |
| Description: | A simple supply and demand network example showing data input from file and the use of "views": incremental definition of arrays of variables. Also uses constraint templates with the arrays of variables. A second version of this model (file xbdlvriis) has modified data making the problem infeasible. This example shows how to analyze infeasibility with the help of IIS (irreducible infeasible sets), it retrieves the IIS and prints out their contents. It is possible to retrieve more detailed information on the IIS, such as isolation rows or bounds, using Xpress Optimizer functions (file xbdlvriis2iso) or to use the infeasibility repair functionality of the Optimizer (file xbdlvriis2rep) with models defined in BCL. |
| File(s): | xbdelvr.cxx, xbdlvriis.cxx, xbdlvriis2iso.cxx, xbdlvriis2rep.cxx |
| Data file(s): | ifvan.dat, cost.dat |
|
|
|
| xbdelvr.cxx |
/********************************************************
Xpress-BCL C++ Example Problems
===============================
file xbdelvr.cxx
`````````````````
Transportation problem.
(c) 2008 Fair Isaac Corporation
author: S.Heipcke, Jan. 2000, rev. Mar. 2011
********************************************************/
#include <iostream>
#include <cstdio>
#include "xprb_cpp.h"
using namespace std;
using namespace ::dashoptimization;
#define NSupp 10 /* Number of suppliers */
#define NCust 7 /* Number of customers */
#define MaxArcs 100 /* Max. num. of non-zero cost values */
#define VANFILE XPRBDATAPATH "/delivery/ifvan.dat" /* Van data file */
#define COSTFILE XPRBDATAPATH "/delivery/cost.dat" /* Cost data file */
/****DATA****/
/* Supplier: London Luton B'ham Bristl Derby Stckpt York */
double SUPPLY[] = {140.0, 600.0, 50.0, 10.0, 400.0, 200.0, 20.0,
/* Supplier: Derby Soton Scnthp */
90.0, 30.0, 12.0};
/* Customer: London Livpol Doncst York Hull Manchr Shffld */
double DEMAND[] = {123.3, 56.4, 17.1, 192.8, 310.0, 47.0, 86.0};
double COST[NSupp][NCust]; /* Cost per supplier-customer pair */
double IFVAN[NSupp][NCust]; /* Non-zero if route uses vans instead
of lorries */
double VANCAP=40.0; /* Capacity on routes that use vans */
XPRBprob p("Delivery"); /* Initialize a new problem in BCL */
/***********************************************************************/
void modDelivery()
{
XPRBexpr lobj, lc;
int s,c;
XPRBvar x[NSupp][NCust];
/****VARIABLES****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
x[s][c] = p.newVar(XPRBnewname("x_s%d",s));
/****OBJECTIVE****/
for(s=0;s<NSupp;s++) /* Objective: Minimize total cost */
for(c=0; c<NCust; c++)
lobj += COST[s][c]*x[s][c];
p.setObj(p.newCtr("OBJ", lobj)); /* Set objective function */
/****CONSTRAINTS****/
for(c=0; c<NCust; c++) /* Satisfy demand of each customer */
{
lc=0;
for(s=0;s<NSupp;s++) lc += x[s][c];
p.newCtr("Demand", lc >= DEMAND[c]);
}
for(s=0;s<NSupp;s++) /* Keep within supply at each supplier*/
{
lc=0;
for(c=0; c<NCust; c++) lc+= x[s][c];
p.newCtr("Supply", lc <= SUPPLY[s]);
}
/****BOUNDS****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0) x[s][c].setUB(VANCAP);
/****SOLVING + OUTPUT****/
p.exportProb(XPRB_MPS,"delivery"); /* Write out an MPS file */
p.setSense(XPRB_MINIM); /* Set objective sense to minimization */
p.lpOptimize(""); /* Solve the LP-problem */
cout << "Objective: " << p.getObjVal() << endl; /* Get objective value */
for(s=0;s<NSupp;s++) /* Print out the solution values */
for(c=0; c<NCust; c++)
cout << x[s][c].getName() << ":" << x[s][c].getSol() << " ";
cout << endl;
}
/***********************************************************************/
/**** Read data from files ****/
void readData()
{
FILE *datafile;
int s,c;
/* Initialize data tables to 0: in the van data file some entries that are
* zero are simply left out, but the function XPRBreadarrline only initializes
* those elements to 0 that have been read, e.g. a data line ",,," results
* in the first 4 elements of the array to be set to zero, even if the
* maximum number of elements to be read (=last parameter of XPRBreadarrline)
* has a much larger value */
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
{
COST[s][c] = 0;
IFVAN[s][c] = 0;
}
/* Read the demand data file */
datafile=fopen(COSTFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", COST[s], NCust);
fclose(datafile);
/* Read the van data file */
datafile=fopen(VANFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99,"g,", IFVAN[s], NCust);
fclose(datafile);
}
/***********************************************************************/
int main(int argc, char **argv)
{
readData(); /* Data input from file */
modDelivery(); /* Problem formulation and solution */
return 0;
}
|
| xbdlvriis.cxx |
/********************************************************
Xpress-BCL C++ Example Problems
===============================
file xbdlvriis.cxx
``````````````````
Transportation problem (infeasible data).
Retrieving and printing IIS.
(c) 2008 Fair Isaac Corporation
author: S.Heipcke, 2005, rev. Mar. 2011
********************************************************/
#include <iostream>
#include <cstdio>
#include <cstring>
#include "xprb_cpp.h"
#include "xprs.h"
using namespace std;
using namespace ::dashoptimization;
#define NSupp 10 // Number of suppliers
#define NCust 7 // Number of customers
#define MaxArcs 100 // Max. num. of non-zero cost values
#define VANFILE XPRBDATAPATH "/delivery/ifvan.dat" // Van data file
#define COSTFILE XPRBDATAPATH "/delivery/cost.dat" // Cost data file
/****DATA****/
// Supplier: London Luton B'ham Bristl Derby Stckpt York
double SUPPLY[] = {140.0, 200.0, 50.0, 10.0, 400.0, 200.0, 20.0,
// Supplier: Derby Soton Scnthp
90.0, 30.0, 12.0};
// Customer: London Livpol Doncst York Hull Manchr Shffld
double DEMAND[] = {1230.3, 560.4, 117.1, 592.8, 310.0, 1247.0, 86.0};
double COST[NSupp][NCust]; // Cost per supplier-customer pair
double IFVAN[NSupp][NCust]; // Non-zero if route uses vans instead
// of lorries
double VANCAP=40.0; // Capacity on routes that use vans
XPRBprob p("Delivery"); // Initialize a new problem in BCL
/***********************************************************************/
void modDelivery()
{
XPRBexpr lobj, lc;
int s,c,i;
XPRBvar x[NSupp][NCust];
int numv, numc, numiis, len, ncol, nrow, scode, ct;
char *vnames, *cnames;
int *viis,*ciis,*vsizes,*csizes;
char **vindex,**cindex;
double *suminfeas = NULL;
XPRSprob op;
/****VARIABLES****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
x[s][c] = p.newVar(XPRBnewname("x_s%d",s));
/****OBJECTIVE****/
for(s=0;s<NSupp;s++) // Objective: Minimize total cost
for(c=0; c<NCust; c++)
lobj += COST[s][c]*x[s][c];
p.setObj(p.newCtr("OBJ", lobj)); // Set objective function
/****CONSTRAINTS****/
for(c=0; c<4; c++) // Satisfy demand of each customer
{
lc=0;
for(s=0;s<5;s++) lc += x[s][c];
p.newCtr("Demand", lc >= DEMAND[c]);
}
for(c=4; c<NCust; c++) // Satisfy demand of each customer
{
lc=0;
for(s=5;s<NSupp;s++) lc += x[s][c];
p.newCtr("Demand", lc >= DEMAND[c]);
}
for(s=0;s<5;s++) // Keep within supply at each supplier
{
lc=0;
for(c=0; c<4; c++) lc+= x[s][c];
p.newCtr("Supply", lc <= SUPPLY[s]);
}
for(s=5;s<NSupp;s++) // Keep within supply at each supplier
{
lc=0;
for(c=4; c<NCust; c++) lc+= x[s][c];
p.newCtr("Supply", lc <= SUPPLY[s]);
}
/****BOUNDS****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0) x[s][c].setUB(VANCAP);
/****SOLVING + OUTPUT****/
p.setSense(XPRB_MINIM); // Set objective sense to minimization
p.lpOptimize(""); // Solve the LP-problem
cout << "LP status: " << p.getLPStat() << endl;
if (p.getLPStat()==XPRB_LP_OPTIMAL)
{
cout << "Objective: " << p.getObjVal() << endl; // Get objective value
}
else if (p.getLPStat()==XPRB_LP_INFEAS) // Get the IIS
{
op = p.getXPRSprob(); // Retrieve the Optimizer problem
/**** Get all IIS ****/
numiis = p.getNumIIS(); // Generate IIS, get number of independent IIS
// Alternatively: use Optimizer functions to generate IIS
/*
XPRSiisall(op);
XPRSgetintattrib(op, XPRS_NUMIIS, &numiis);
*/
cout << "Number of IIS: " << numiis << endl;
/**** IIS summary: Retrieve set sizes, including for the approximation ****/
vsizes = new int[numiis+1];
csizes = new int[numiis+1];
suminfeas = new double[numiis+1];
XPRSiisstatus(op, &numiis, csizes, vsizes, suminfeas, NULL);
cout << "Size of subproblem Sum of infeasibilities"<< endl;
cout << " (bounds/constr.)" << endl;
cout << " Approx: (" << vsizes[0] << "/" << csizes[0] << ") ";
cout << suminfeas[0] << endl;
for(i=1; i<=numiis;i++) {
cout << " IIS " << i << ": (" << vsizes[i] << "/" << csizes[i];
cout << ") " << suminfeas[i] << endl;
}
delete [] vsizes; delete [] csizes; delete [] suminfeas;
/**** Alternatively: Enumeration with 'next' ****/
/* Call 'iisfirst' it with argument 0 to obtain the information on the
IIS approximation before the first 'iisnext' call */
XPRSiisfirst(op, 0, &scode);
ct=-1;
while (scode==0)
{
ct++;
XPRSgetiisdata(op, ct, &numc, &numv, NULL, NULL, NULL, NULL,
NULL, NULL, NULL, NULL);
if (ct==0) {
/* IIS number zero denotes the initial IIS approximation. 'getiisdata'
returns the [numbers of] columns and rows with nonzero reduced
costs/duals. The subproblem defined by these nonzero reduced costs
and duals is infeasible by itself.
The sizes returned by 'iisstatus' for the IIS appoximation may be
larger since it returns all columns and rows of the approximation,
independent of their reduced cost/dual values.
*/
cout << "IIS approximation: " << numv << " columns and " << numc;
cout << " constraints with nonzero duals" << endl;
}
else {
/* It's a valid IIS, so the sizes are returned */
cout << "IIS " << ct << ": " << numv << " bounds and " << numc;
cout << " constraints" << endl;
}
XPRSiisnext(op, &scode);
}
/**** Obtain variable and constraint names for later use in printout ****/
// Retrieve variable names
XPRSgetintattrib(op, XPRS_ORIGINALCOLS, &ncol);
XPRSgetnamelist(op, 2, NULL, 0, &len, 0, ncol-1);
// Get number of bytes required for retrieving names
vnames = new char[len];
vindex = new char*[ncol];
XPRSgetnamelist(op, 2, vnames, len, NULL, 0, ncol-1);
vindex[0]=vnames;
for(i=1; i<ncol; i++) vindex[i] =vindex[i-1]+strlen(vindex[i-1])+1;
// Retrieve constraint names
XPRSgetintattrib(op, XPRS_ORIGINALROWS, &nrow);
XPRSgetnamelist(op, 1, NULL, 0, &len, 0, nrow-1);
cnames = new char[len];
cindex = new char*[nrow];
XPRSgetnamelist(op, 1, cnames, len, NULL, 0, nrow-1);
cindex[0]=cnames;
for(i=1; i<nrow; i++) cindex[i] =cindex[i-1]+strlen(cindex[i-1])+1;
/**** Retrieve variables and constraints contained in IIS ****/
for(s=1;s<=numiis;s++)
{
XPRSgetiisdata(op, s, &numc, &numv, NULL, NULL, NULL, NULL,
NULL, NULL, NULL, NULL);
ciis = new int[numc];
viis = new int[numv];
XPRSgetiisdata(op, s, &numc, &numv, ciis, viis, NULL, NULL,
NULL, NULL, NULL, NULL);
cout << "IIS " << s << ": " << numv << " variables, ";
cout << numc << " constraints" << endl;
if (numv>0)
{
cout << " Variables: "; // Print all variables in the IIS
for(i=0;i<numv;i++) cout << vindex[viis[i]] << " ";
cout << endl;
delete [] viis; // Free the array of variables
}
if (numc>0)
{
cout << " Constraints: "; // Print all constraints in the IIS
for(i=0;i<numc;i++) cout << cindex[ciis[i]] << " ";
cout << endl;
delete [] ciis; // Free the array of constraints
}
}
delete [] vnames;
delete [] cnames;
delete [] vindex;
delete [] cindex;
}
}
/***********************************************************************/
/**** Read data from files ****/
void readData()
{
FILE *datafile;
int s,c;
// Initialize data tables to 0
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
{
COST[s][c] = 0;
IFVAN[s][c] = 0;
}
// Read the demand data file
datafile=fopen(COSTFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", COST[s], NCust);
fclose(datafile);
// Read the van data file
datafile=fopen(VANFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", IFVAN[s], NCust);
fclose(datafile);
}
/***********************************************************************/
int main(int argc, char **argv)
{
readData(); // Data input from file
modDelivery(); // Problem formulation and solution
return 0;
}
|
| xbdlvriis2iso.cxx |
/********************************************************
Xpress-BCL C++ Example Problems
===============================
file xbdlvriis2iso.cxx
``````````````````````
Transportation problem (infeasible data).
Retrieving and printing IIS.
- Using Optimizer functions to retrieve detailed
IIS information including isolation rows/bounds -
(c) 2008 Fair Isaac Corporation
author: S.Heipcke, Jan. 2008, rev. Mar. 2011
********************************************************/
#include <iostream>
#include <cstdio>
#include <cstring>
#include "xprb_cpp.h"
#include "xprs.h"
using namespace std;
using namespace ::dashoptimization;
#define NSupp 10 // Number of suppliers
#define NCust 7 // Number of customers
#define MaxArcs 100 // Max. num. of non-zero cost values
#define VANFILE XPRBDATAPATH "/delivery/ifvan.dat" // Van data file
#define COSTFILE XPRBDATAPATH "/delivery/cost.dat" // Cost data file
/****DATA****/
// Supplier: London Luton B'ham Bristl Derby Stckpt York
double SUPPLY[] = {140.0, 200.0, 50.0, 10.0, 400.0, 200.0, 20.0,
// Supplier: Derby Soton Scnthp
90.0, 30.0, 12.0};
// Customer: London Livpol Doncst York Hull Manchr Shffld
double DEMAND[] = {1230.3, 560.4, 117.1, 592.8, 310.0, 1247.0, 86.0};
double COST[NSupp][NCust]; // Cost per supplier-customer pair
double IFVAN[NSupp][NCust]; // Non-zero if route uses vans instead
// of lorries
double VANCAP=40.0; // Capacity on routes that use vans
XPRBprob p("Delivery"); // Initialize a new problem in BCL
/***********************************************************************/
void modDelivery()
{
XPRBexpr lobj, lc;
int s,c,i;
XPRBvar x[NSupp][NCust];
int numv, numc, numiis, len, ncol, nrow;
char *vnames, *cnames;
int *viis,*ciis;
char **vindex,**cindex;
double bnd, rhs;
char *ctrtype = NULL;
char *bndtype = NULL;
double *duals = NULL;
double *rdcs = NULL;
char *isolationrows = NULL;
char *isolationbnds = NULL;
const char *isotype[] = {"N/A", "No ", "Yes"};
XPRSprob op;
/****VARIABLES****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
x[s][c] = p.newVar(XPRBnewname("x_s%d",s));
/****OBJECTIVE****/
for(s=0;s<NSupp;s++) // Objective: Minimize total cost
for(c=0; c<NCust; c++)
lobj += COST[s][c]*x[s][c];
p.setObj(p.newCtr("OBJ", lobj)); // Set objective function
/****CONSTRAINTS****/
for(c=0; c<4; c++) // Satisfy demand of each customer
{
lc=0;
for(s=0;s<5;s++) lc += x[s][c];
p.newCtr("Demand", lc >= DEMAND[c]);
}
for(c=4; c<NCust; c++) // Satisfy demand of each customer
{
lc=0;
for(s=5;s<NSupp;s++) lc += x[s][c];
p.newCtr("Demand", lc >= DEMAND[c]);
}
for(s=0;s<5;s++) // Keep within supply at each supplier
{
lc=0;
for(c=0; c<4; c++) lc+= x[s][c];
p.newCtr("Supply", lc <= SUPPLY[s]);
}
for(s=5;s<NSupp;s++) // Keep within supply at each supplier
{
lc=0;
for(c=4; c<NCust; c++) lc+= x[s][c];
p.newCtr("Supply", lc <= SUPPLY[s]);
}
/****BOUNDS****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0) x[s][c].setUB(VANCAP);
/****SOLVING + OUTPUT****/
p.setSense(XPRB_MINIM); // Set objective sense to minimization
p.lpOptimize(""); // Solve the LP-problem
cout << "LP status: " << p.getLPStat() << endl;
if (p.getLPStat()==XPRB_LP_OPTIMAL)
{
cout << "Objective: " << p.getObjVal() << endl; // Get objective value
}
else if (p.getLPStat()==XPRB_LP_INFEAS) // Problem is infeasible
{
op = p.getXPRSprob(); // Retrieve the Optimizer problem
/**** Get all IIS ****/
numiis = p.getNumIIS(); // Generate IIS, get number of independent IIS
// Alternatively: use Optimizer functions to generate IIS
/*
XPRSiisall(op);
XPRSgetintattrib(op, XPRS_NUMIIS, &numiis);
*/
cout << "Number of IIS: " << numiis << endl;
/**** Obtain variable and constraint names for later use in printout ****/
// Retrieve variable names
XPRSgetintattrib(op, XPRS_ORIGINALCOLS, &ncol);
XPRSgetnamelist(op, 2, NULL, 0, &len, 0, ncol-1);
// Get number of bytes required for retrieving names
vnames = new char[len];
vindex = new char*[ncol];
XPRSgetnamelist(op, 2, vnames, len, NULL, 0, ncol-1);
vindex[0]=vnames;
for(i=1; i<ncol; i++) vindex[i] =vindex[i-1]+strlen(vindex[i-1])+1;
// Retrieve constraint names
XPRSgetintattrib(op, XPRS_ORIGINALROWS, &nrow);
XPRSgetnamelist(op, 1, NULL, 0, &len, 0, nrow-1);
cnames = new char[len];
cindex = new char*[nrow];
XPRSgetnamelist(op, 1, cnames, len, NULL, 0, nrow-1);
cindex[0]=cnames;
for(i=1; i<nrow; i++) cindex[i] = cindex[i-1]+strlen(cindex[i-1])+1;
/**** Retrieve detailed IIS info (incl. isolations) ****/
for(s=1;s<=numiis;s++)
{
XPRSgetiisdata(op, s, &numc, &numv, NULL, NULL, NULL, NULL,
NULL, NULL, NULL, NULL);
XPRSiisisolations(op, s); // Find isolations
ciis = new int[numc];
viis = new int[numv];
ctrtype = new char[numc];
bndtype = new char[numv];
duals = new double[numc];
rdcs = new double[numv];
isolationrows = new char[numc];
isolationbnds = new char[numv];
XPRSgetiisdata(op, s, &numc, &numv, ciis, viis, ctrtype, bndtype,
duals, rdcs, isolationrows, isolationbnds);
cout << "IIS " << s << ": " << numv << " variables, " << numc;
cout << " constraints" << endl;
cout << " Name Type Sense Bound Dual values In iso. " << endl;
if (numv>0)
{ // Print all variables in the IIS
for(i=0;i<numv;i++)
{
if (bndtype[i] == 'L') XPRSgetlb(op, &bnd, viis[i], viis[i]);
else XPRSgetub(op, &bnd, viis[i], viis[i]);
cout << " " << vindex[viis[i]] << " column " << bndtype[i];
cout << " " << bnd << " " << rdcs[i] << " ";
cout << isotype[1+(int)isolationbnds[i]] << endl;
}
delete [] viis; // Free the array of variables
delete [] bndtype;
delete [] rdcs;
delete [] isolationbnds;
}
if (numc>0)
{ // Print all constraints in the IIS
for(i=0;i<numc;i++)
{
XPRSgetrhs(op, &rhs, ciis[i], ciis[i]);
cout << " " << cindex[ciis[i]]<< " row " << ctrtype[i];
cout << " " << rhs << " " << duals[i] << " ";
cout << isotype[1+(int)isolationrows[i]] << endl;
}
delete [] ciis; // Free the array of constraints
}
}
delete [] vnames;
delete [] cnames;
delete [] vindex;
delete [] cindex;
}
}
/***********************************************************************/
/**** Read data from files ****/
void readData()
{
FILE *datafile;
int s,c;
// Initialize data tables to 0
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
{
COST[s][c] = 0;
IFVAN[s][c] = 0;
}
// Read the demand data file
datafile=fopen(COSTFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", COST[s], NCust);
fclose(datafile);
// Read the van data file
datafile=fopen(VANFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", IFVAN[s], NCust);
fclose(datafile);
}
/***********************************************************************/
int main(int argc, char **argv)
{
readData(); // Data input from file
modDelivery(); // Problem formulation and solution
return 0;
}
|
| xbdlvriis2rep.cxx |
/********************************************************
Xpress-BCL C++ Example Problems
===============================
file xbdlvriis2rep.cxx
``````````````````````
Transportation problem (infeasible data).
Repairing infeasibility.
- Using Optimizer functions -
(c) 2008 Fair Isaac Corporation
author: S.Heipcke, Jan. 2008, rev. Mar. 2011
********************************************************/
#include <iostream>
#include <cstdio>
#include <cstring>
#include "xprb_cpp.h"
#include "xprs.h"
using namespace std;
using namespace ::dashoptimization;
#define NSupp 10 // Number of suppliers
#define NCust 7 // Number of customers
#define MaxArcs 100 // Max. num. of non-zero cost values
#define VANFILE XPRBDATAPATH "/delivery/ifvan.dat" // Van data file
#define COSTFILE XPRBDATAPATH "/delivery/cost.dat" // Cost data file
void printsolution(int scode, XPRBvar x[][NCust]);
void printsolution2(XPRSprob op, int scode, XPRBvar x[][NCust]);
/****DATA****/
// Supplier: London Luton B'ham Bristl Derby Stckpt York
double SUPPLY[] = {140.0, 200.0, 50.0, 10.0, 400.0, 200.0, 20.0,
// Supplier: Derby Soton Scnthp
90.0, 30.0, 12.0};
// Customer: London Livpol Doncst York Hull Manchr Shffld
double DEMAND[] = {1230.3, 560.4, 117.1, 592.8, 310.0, 1247.0, 86.0};
double COST[NSupp][NCust]; // Cost per supplier-customer pair
double IFVAN[NSupp][NCust]; // Non-zero if route uses vans instead
// of lorries
double VANCAP=40.0; // Capacity on routes that use vans
XPRBprob p("Delivery"); // Initialize a new problem in BCL
/***********************************************************************/
void modDelivery()
{
XPRBexpr lobj, lc;
XPRBctr ctr, CSupply[NSupp], CDemand[NCust];
int s,c;
XPRBvar x[NSupp][NCust];
int ncol, nrow, scode;
double *lrp, *grp, *lbp, *ubp;
XPRSprob op;
/****VARIABLES****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
x[s][c] = p.newVar(XPRBnewname("x_s%d",s));
/****OBJECTIVE****/
for(s=0;s<NSupp;s++) // Objective: Minimize total cost
for(c=0; c<NCust; c++)
lobj += COST[s][c]*x[s][c];
p.setObj(p.newCtr("OBJ", lobj)); // Set objective function
/****CONSTRAINTS****/
for(c=0; c<4; c++) // Satisfy demand of each customer
{
lc=0;
for(s=0;s<5;s++) lc += x[s][c];
CDemand[c] = p.newCtr("Demand", lc >= DEMAND[c]);
}
for(c=4; c<NCust; c++) // Satisfy demand of each customer
{
lc=0;
for(s=5;s<NSupp;s++) lc += x[s][c];
CDemand[c] = p.newCtr("Demand", lc >= DEMAND[c]);
}
for(s=0;s<5;s++) // Keep within supply at each supplier
{
lc=0;
for(c=0; c<4; c++) lc+= x[s][c];
CSupply[s] = p.newCtr("Supply", lc <= SUPPLY[s]);
}
for(s=5;s<NSupp;s++) // Keep within supply at each supplier
{
lc=0;
for(c=4; c<NCust; c++) lc+= x[s][c];
CSupply[s] = p.newCtr("Supply", lc <= SUPPLY[s]);
}
/****BOUNDS****/
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0) x[s][c].setUB(VANCAP);
/****SOLVING + OUTPUT****/
p.setSense(XPRB_MINIM); // Set objective sense to minimization
p.lpOptimize(""); // Solve the LP-problem
cout << "LP status: " << p.getLPStat() << endl;
if (p.getLPStat()==XPRB_LP_OPTIMAL)
{
cout << "Objective: " << p.getObjVal() << endl; // Get objective value
}
else if (p.getLPStat()==XPRB_LP_INFEAS) // Get the IIS
{
op = p.getXPRSprob(); // Retrieve the Optimizer problem
/**** Trying to fix infeasibilities ****/
/*
lrp: (affects = and <= rows)
ax - aux_var = b
ax - aux_var <= b
grp: (affects = and >= rows)
ax + aux_var = b
ax + aux_var >= b
lbp:
x_i + aux_var >= l
ubp:
x_i - aux_var <= u
*/
/**** Simplified infeasibility repair:
specifying preferences per constraint/bound type ****/
cout << "\n**** Repair infeasibility:" << endl;
XPRSrepairinfeas(op, &scode, 'c', 'o', ' ', 10, 9, 0, 20, 0.001);
printsolution2(op, scode, x); /* Print out the solution values */
/**** Weighted infeasibility repair:
specifying preferences for every constraint/bound separately ****/
XPRSgetintattrib(op, XPRS_ORIGINALCOLS, &ncol);
XPRSgetintattrib(op, XPRS_ORIGINALROWS, &nrow);
lrp = new double[nrow];
grp = new double[nrow];
lbp = new double[ncol];
ubp = new double[ncol];
memset(lrp, 0, nrow*sizeof(double));
memset(grp, 0, nrow*sizeof(double));
memset(lbp, 0, ncol*sizeof(double));
memset(ubp, 0, ncol*sizeof(double));
/* Relax bounds due to van capacity */
/* Repairweightedinfeas for upper bounds of concerned flow variables */
for(s=0; s<NSupp; s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0) ubp[x[s][c].getColNum()] = 20;
cout << "\n**** Relax van capacity:" << endl;
XPRSrepairweightedinfeas(op, &scode, lrp, grp, lbp, ubp, 'd', 0.001, "");
printsolution2(op, scode, x);
/* Relax supply limits (may buy in additional quantities) */
/* Repairinfeas for 'less or equal' side of Supply constraints */
for(s=0; s<NSupp; s++) lrp[CSupply[s].getRowNum()] = 10;
cout << "\n**** Relax supply limits:" << endl;
XPRSrepairweightedinfeas(op, &scode, lrp, grp, lbp, ubp, 'd', 0.001, "");
printsolution(scode, x);
/* Relax demand constraints (may not satisfy all customers) */
/* Repairinfeas for 'greater or equal' side of Demand constraints */
for(c=0; c<NCust; c++) grp[CDemand[c].getRowNum()] = 9;
cout << "\n**** Relax demand constraints:" << endl;
XPRSrepairweightedinfeas(op, &scode, lrp, grp, lbp, ubp, 'd', 0.001, "");
printsolution2(op, scode, x);
delete [] lrp;
delete [] grp;
delete [] lbp;
delete [] ubp;
}
}
/***********************************************************************/
/**** Read data from files ****/
void readData()
{
FILE *datafile;
int s,c;
// Initialize data tables to 0
for(s=0;s<NSupp;s++)
for(c=0; c<NCust; c++)
{
COST[s][c] = 0;
IFVAN[s][c] = 0;
}
// Read the demand data file
datafile=fopen(COSTFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", COST[s], NCust);
fclose(datafile);
// Read the van data file
datafile=fopen(VANFILE,"r");
for(s=0;s<NSupp;s++)
XPRBreadarrlinecb(XPRB_FGETS, datafile, 99, "g,", IFVAN[s], NCust);
fclose(datafile);
}
/**** Print out the solution values ****/
void printsolution(int scode, XPRBvar x[][NCust])
{
int s,c;
double sup,dem;
const char *rstat[] = {"relaxed optimum found", "relaxed problem infeasible",
"relaxed problem unbounded", "solution nonoptimal for original objective",
"error", "numerical instability"};
cout << "Status: " << rstat[scode] << endl;
if(scode==0)
{
p.sync(XPRB_XPRS_SOL);
for(s=0; s<NSupp; s++)
for(c=0; c<NCust; c++)
if(x[s][c].getSol()>0.01)
cout << x[s][c].getName() << ":" << x[s][c].getSol() << " ";
cout << endl;
cout << "Violations:" << endl;
for(c=0; c<NCust; c++)
{
sup=0;
for(s=0; s<NSupp; s++) sup+=x[s][c].getSol();
if(sup<DEMAND[c])
cout << " Customer " << c << ": " << DEMAND[c]-sup << endl;
}
for(s=0; s<NSupp; s++)
{
dem=0;
for(c=0; c<NCust; c++) dem+=x[s][c].getSol();
if(dem>SUPPLY[s])
cout << " Supplier " << s << ": " << dem-SUPPLY[s] << endl;
}
for(s=0; s<NSupp; s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0 && VANCAP<x[s][c].getSol())
cout << " Van " << s << "-" << c << ": " << x[s][c].getSol()-VANCAP << endl;
}
}
void printsolution2(XPRSprob op, int scode, XPRBvar x[][NCust])
{
int s,c,ncol;
double sup,dem;
const char *rstat[] = {"relaxed optimum found", "relaxed problem infeasible",
"relaxed problem unbounded", "solution nonoptimal for original objective",
"error", "numerical instability"};
double *sol;
cout << "Status: " << rstat[scode] << endl;
if(scode==0)
{
XPRSgetintattrib(op, XPRS_ORIGINALCOLS, &ncol);
sol = new double[ncol];
XPRSgetlpsol(op, sol, NULL, NULL, NULL); // Get the solution values
for(s=0; s<NSupp; s++)
for(c=0; c<NCust; c++)
if(x[s][c].getColNum()>=0 && sol[x[s][c].getColNum()]>0.01)
cout << x[s][c].getName() << ":" << sol[x[s][c].getColNum()] << " ";
cout << endl;
cout << "Violations:" << endl;
for(c=0; c<NCust; c++)
{
sup=0;
for(s=0; s<NSupp; s++)
if(x[s][c].getColNum()>=0) sup+=sol[x[s][c].getColNum()];
if(sup<DEMAND[c])
cout << " Customer " << c << ": " << DEMAND[c]-sup << endl;
}
for(s=0; s<NSupp; s++)
{
dem=0;
for(c=0; c<NCust; c++)
if(x[s][c].getColNum()>=0) dem+=sol[x[s][c].getColNum()];
if(dem>SUPPLY[s])
cout << " Supplier " << s << ": " << dem-SUPPLY[s] << endl;
}
for(s=0; s<NSupp; s++)
for(c=0; c<NCust; c++)
if(IFVAN[s][c]!=0 && x[s][c].getColNum()>=0 &&
VANCAP<sol[x[s][c].getColNum()]) {
cout << " Van " << s << "-" << c << ": ";
cout << sol[x[s][c].getColNum()]-VANCAP << endl;
}
}
}
/***********************************************************************/
int main(int argc, char **argv)
{
readData(); // Data input from file
modDelivery(); // Problem formulation and solution
return 0;
}
|