FICO
FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious exampleNext example

Delivery - Data input from file; infeasibility analysis

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.

xbdelvrcpp.zip[download all files]

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.
xbdelvr.cxx[download]
xbdlvriis.cxx[download]
xbdlvriis2iso.cxx[download]
xbdlvriis2rep.cxx[download]

Data Files





xbdlvriis2rep.cxx

/********************************************************
  Xpress-BCL C++ Example Problems
  ===============================

  file xbdlvriis2rep.cxx
  ``````````````````````
  Transportation problem (infeasible data).
  Repairing infeasibility.
  - Using Optimizer functions -

  (c) 2008-2024 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(XPRBprob &p, 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

/***********************************************************************/

void modDelivery(XPRBprob &p)
{
 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(p, 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(XPRBprob &p, 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)
{
 XPRBprob p("Delivery");          // Initialize a new problem in BCL
 readData();                      // Data input from file
 modDelivery(p);                  // Problem formulation and solution

 return 0;
} 

Back to examples browserPrevious exampleNext example