Problem n.1 : normal LP

This computes 

min  c x
     A x = b
     lb <= x <= ub
       
through column generation, by finding the convex combination
of the 2^n vertices of the box   lb <= x <= ub, satisfying
A x = b and minimizing c x.

Problem data

Instance n. 1

A = {{ 1,-3,2, 1,-1,1,2, 1, 0,1,0,0, 0},
      { 2, 1,5, 2, 1,0,0, 0, 1,0,1,0, 0},
      { 1, 0,2,-1, 0,1,0, 2, 3,0,0,1, 0},
      {-2,-1,3, 0, 0,0,0,-1,-2,0,0,0,-1}};

m = Length[A]

n = Length[A[[1]]]

inf = 10^10;

lb ={0,0,0,0,0,0,0,0,0,0,0,0,0};
ub={5,5,4,3,3,5,6,2,4,10,10,10,10};

c = {1, -3 , -1, 4,-3, 10,-6,-4,-1,0,0,0,0};

b = {8,5,4,-6};

If[Length[b]!=m, Print["Check b !"]];
If[Length[c]!=n, Print["Check c !"]];
If[Length[lb]!=n, Print["Check lb !"]];
If[Length[ub]!=n, Print["Check ub !"]];

Instance n . 2

A = {{ 1,0,1,1,0,0,0,0,0,0,0},
      { -1,1,0,0,-1,1,1,0,0,0,0},
      { 0,-1,0,0,0,0,0,-1,1,0,0},
      {0,0,-1,0,1,0,0,0,0,1,0},
      {0,0,0,-1,0,-1,0,1,0,-1,1}};

m = Length[A]

n = Length[A[[1]]]

tolerance=0;

inf = 10^10;

lb ={0,0,0,0,0,0,0,0,0,0,0};
ub ={15,8,4,12,7,6,13,18,12,10,6};

c = {5,3,2,1,4,10,7,3,2,1,4};

b = {20,0,0,0,0};

If[Length[b]!=m, Print["Check b !"]];
If[Length[c]!=n, Print["Check c !"]];
If[Length[lb]!=n, Print["Check lb !"]];
If[Length[ub]!=n, Print["Check ub !"]];

Phase One

Initialization

bb=PrependTo[b,1]

sign=Table[If[bb[[i]]<0,-1,1],{i,m+1}]

BT=DiagonalMatrix[sign];

B=BT;

xlist=Table[{},{m+1}];

cb=Table[1,{m+1}];

Iteration

Basis Change ( to be skipped on the first iteration !  )

newcol=A.xt

PrependTo[newcol,1]

y = W . newcol

Block[{i,ratio,min},
      min = 1000000;
      For[i=1, i<=m+1, i++,
          If[ y[[i]]>0,
              ratio = x[[i]]/y[[i]];
              If[ratio < min,
                 min = ratio; q = i]
            ]  
      ]
]

q

BT[[q]]=newcol

cb[[q]]=0

xlist[[q]]=xt;

B=Transpose[BT];

MatrixForm[B];

start here on first iteration !

W=Inverse[B];

u=cb.W

cb

x=W.bb

u0=u[[1]];uu=Drop[u,1];

cc=-uu.A

xt=Table[If[cc[[i]]>=0,lb[[i]],ub[[i]]],{i,Length[cc]}]

-u0+cc.xt

If equal to 0 stop, otherwise run second rightmost bracket

Feasibility check

cb

x

xlist

sol=Transpose[xlist].x

A.sol

Phase Two

Initialization

cb=xlist.c

Iteration

Basis Change ( to be skipped on the first iteration !  )

newcol=A.xt

PrependTo[newcol,1]

y = W . newcol

Block[{i,ratio,min},
      min = 1000000;
      For[i=1, i<=m+1, i++,
          If[ y[[i]]>0,
              ratio = x[[i]]/y[[i]];
              If[ratio < min,
                 min = ratio; q = i]
            ]  
      ]
]

q

BT[[q]]=newcol

cb[[q]]=c.xt

xlist[[q]]=xt;

B=Transpose[BT];

MatrixForm[B];

W=Inverse[B];

start here on first iteration !

u=cb.W

x=W.bb

u0=u[[1]];uu=Drop[u,1];

cc=c-uu.A

xt=Table[If[cc[[i]]>=0,lb[[i]],ub[[i]]],{i,Length[cc]}]

-u0+cc.xt

If equal to 0 stop, otherwise run second rightmost bracket

Solution check

x

xlist

sol=Transpose[xlist].x

A.sol

Problem n. 2 : Bin Packing (relaxation)

This solves a relaxation of the Bin Packing problem.
Object types  j=1,..,m  are defined. For each type
a size  s  and a quantity  b are defined. Each bin
has capacity K. We want to fit all objects in the bins
while minimizing the number of used bins.

Routines

Off[General::spell1]

ZeroOneKnapsack[v_,w_,K_]:=Block[
{val,index,solution,n,pnt},
n=Length[v];
val=Table[0,{K+1}];
pnt=Table[0,{K+1},{n}];
Do[
Do[If[val[[j+1]] < val[[j+1 - w[[i]]]]+v[[i]],
       val[[j+1]]=val[[j+1 - w[[i]]]]+v[[i]];
        pnt[[j+1,i]]=1]
       ,{j,K,w[[i]],-1 }]
,{i,n}];
solution=Table[0,{Length[v]}];
index=K;
Do[If[pnt[[index+1,i]]==1,index=index-w[[i]];solution[[i]]=1]   ,{i,n,1,-1} ];
 solution
]

intknapsack[v_,w_,K_]:=Block[
{mininfinity,val,type,maxval,index,solution,obj},
mininfinity=-1000000;
val=Table[mininfinity,{K+1}];
val[[1]]=0;
type=Table[0,{K+1}];
maxval=0;
Do[
Do[If[i-w[[j]]>=0,
      If[val[[i-w[[j]]+1]]+v[[j]] > val[[i+1]],
         val[[i+1]]=val[[i-w[[j]]+1]]+v[[j]] ;
         type[[i+1]]=j;
         If[maxval<val[[i+1]],
             maxval=val[[i+1]];
             index=i
           ]
         ]
       ]
,{j,Length[v]}]
,{i,K}];
solution=Table[0,{Length[v]}];
While[index>0,
      obj=type[[index+1]];
      solution[[obj]]++;
      index=index-w[[obj]];
      ];
 solution
]

Problem data

data #1

capacity of each bin

cap=15;

type sizes

s = {1,2,3,4,5,7,9,10,11,12};

required type quantities

b = {10,10,5,16,15,8,9,10,22,5};

If[Length[s]!=Length[b],Print["check data !"]]

m=Length[s];

data #2

capacity of each bin

cap=15;

type sizes

s = {1,2,2,2,3,3,4,5,5,5,5,7,7,9,9,9,9,10,11,12};

required type quantities

m = Length[s]

b = Table[1,{m}];

Initialization

BT=IdentityMatrix[m];

B=BT;

xlist=Table[{},{m}];

cb=Table[1,{m}];

Iteration

Basis Change ( to be skipped on the first iteration !  )

newcol=pattern

y = W . newcol

Block[{i,ratio,min},
      min = 1000000;
      For[i=1, i<=m, i++,
          If[ y[[i]]>0,
              ratio = x[[i]]/y[[i]];
              If[ratio < min,
                 min = ratio; q = i]
            ]  
      ]
]

q

BT[[q]]=newcol

B=Transpose[BT];

start here on first iteration !

W=Inverse[B];

u=cb.W

x=W.b

current active patterns

MatrixForm[BT]

current number of bins

Print[IntegerPart[u . b], " + ", FractionalPart[u . b]]

new pattern

pattern=ZeroOneKnapsack[u,s,cap]

cc=1-u.pattern

If equal to 0 stop, otherwise run second rightmost bracket

Solution check

x

opt=Sum[x[[i]],{i,m}]

Ceiling[opt]

randomized rounding

xf=Table[fp=FractionalPart[x[[i]]];If[fp==0,x[[i]],IntegerPart[x[[i]]]+If[Random[]<fp,1,0]],{i,Length[x]}]

Sum[xf[[i]],{i,m}]

the rounding is feasible is next vector is nonnegative

B.xf -b

BT[[Flatten[Position[xf, 1], 1]]]

Problem n. 3 : Dual of TSP

This solves the Lagrangean dual of the TSP (as in the Held-Karp formulation) by solving the dual of the dual (an LP problem) by column generation

Routines

Off[General::spell1]

onetreedegree[onetree_]:=Block[{values},
values=Table[0,{nn}];
ft=Flatten[onetree];
Do[values[[ft[[i]]]]++
 ,{i,Length[ft]}];
 values];

basecheck:=Do[If[BT[[i]]==tourvect,
   Print["1-tree #",i," is a tour! 
   Change on BT and cb must be done manually"] ],
   {i,nn-2}];

costcompute[tree_]:=Sum[cost[[tree[[i,1]],tree[[i,2]]]],
             {i,Length[tree]}]

show[a_,b_,c_,d_]:=Block[{b1,dd1,dd2,dd3},
           Which[b==1,b1=0.005,b==2,b1=0.008,b==3,b1=0.015];
           Which[d=="R",dd1=1;dd2=0;dd3=0;,
                 d=="B",dd1=0;dd2=0;dd3=1;,
                 d=="G",dd1=0;dd2=1;dd3=0;,
                 d=="Y",dd1=1;dd2=1;dd3=0.5;,
                 d=="Bk",dd1=0;dd2=0;dd3=0;
                 ];
           If[c==0,
          Graphics[{Thickness[b1],RGBColor[dd1,dd2,dd3],
          Table[Line[{nodes[[a[[i,1]]]],
                         nodes[[a[[i,2]]]]}
                       ]
                   ,{i,Length[a]}]}],
          Graphics[{Thickness[b1],Dashing[{0.01}],RGBColor[dd1,dd2,dd3],
          Table[Line[{nodes[[a[[i,1]]]],
                         nodes[[a[[i,2]]]]}
                       ]
                   ,{i,Length[a]}]}]
                   ]]

minspantree[arcs_]:=Block[{p,r,kk,lt,ki,kj,ri,rj},
root[k_]:=Block[{kp},kp=k;While[kp!=p[[kp]],kp=p[[kp]]];kp];
p=Table[i,{i,nn}];
r=Table[1,{nn}];
For[kk=1;lt=0;tree={},
    (lt<nn-2)&&(kk<=Length[arcs]),
    kk++,
    ki=arcs[[kk,1]];
    kj=arcs[[kk,2]];
    ri=root[ki];
    rj=root[kj];
    If[ri!=rj,
       Which[r[[ri]]>r[[rj]]
            ,p[[rj]]=ri
            ,r[[rj]]>r[[ri]]
            ,p[[ri]]=rj
            ,r[[ri]]==r[[rj]]
            ,p[[rj]]=ri;
             r[[ri]]++
             ];
        lt++;
        AppendTo[tree,arcs[[kk]]];
       ];
     ];
   If[lt<nn-2,Print["disconnected graph !"];
             Do[pp=ii;
                While[!p[[pp]]==pp,pp=p[[pp]]];
                Print[ii," -> ",pp]
                ,{ii,n}]
     ];

]

onetreefind[costmatrix_]:=Block[
{specarcs,arcs,first,second,val},
specarcs={};
Do[AppendTo[specarcs,{costmatrix[[i,j]],{i,j}}],
{i,2,nn-1},{j,i+1,nn}];
arcs=Transpose[Sort[specarcs]][[2]];
minspantree[arcs];
first={10000,{0,0}};
second={10000,{0,0}};
Do[val=costmatrix[[1,i]];
   If[val <second[[1]],
     If[val >= first[[1]],
          second={val,{1,i}},
          second=first;
          first={val,{1,i}}]
      ]
   ,{i,2,nn}];
onetree=Append[tree,first[[2]]];
AppendTo[onetree,second[[2]]];
   ];

How to use this notebook

First run the section Routines. 

Then generate the data as described in the section Data. You have alternative ways of generating the data.

Data

The nodes are randomly generated within a square.
Then the arc costs are taken as Euclidean distances.

(* nice seeds=40 (nn=10) -> many iterations
              10 (nn=10) -> conv. comb. of two 1-trees,
                            many iterations
              50 (nn=10) -> solution at 48th iteration,
                            nice graphs.
              40 (nn= 8) -> optimal tour !
              40 (nn= 7) -> optimal tour !
              25 (nn= 5) -> very fast with duality gap
              35 (nn= 5) -> fast with no gap
              25 (nn= 6) -> no few iterations, no gap*)

seed=55;nn=5;

na=nn(nn-1)/2;

SeedRandom[seed]

nodes={};
Do[AppendTo[nodes,
 {Random[Integer,{1,99}],Random[Integer,{1,99}]}],{nn}]

nodes  (* coordinates *)

bb=Flatten[Table[Floor[N[Sqrt[(nodes[[i,1]]-nodes[[j,1]])^2+
                    (nodes[[i,2]]-nodes[[j,2]])^2]]],
           {i,nn-1},{j,i+1,nn}]];

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];

MatrixForm[costs]

cost=Transpose[Drop[Transpose[Drop[costs,1]],1]];

rangeplot={{-5,105},{-5,105}};

aaa=Graphics[{PointSize[0.015],Table[{Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],{-2,0}]},{i,Length[nodes]}]}];

Show[aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Starting base

First a minimal spanning 1-tree is found. onetree is the list of
arcs of the 1-tree. tree is the minimal spamming tree over the nodes
{2,..,nn}

onetreefind[cost]

onetree

tree

treeshow=show[onetree,2,0,"B"];

Show[aaa,treeshow,AspectRatio->Automatic];

b is the rhs vector

b=Table[2,{nn}];

b[[1]]=1;

The first nn-2 columns are computed from tree by adding in turn
the arcs {1,2} and {1,i}  for  i=3 to nn. The following table gives
the costs of these 1-trees.

cb=costcompute[tree]+
         Table[cost[[1,2]]+cost[[1,i]],{i,3,nn}]

The (nn-1)-th column is given  by any tour. In this case it is the following

tour={};
Do[AppendTo[tour,{i,i+1}],{i,nn-1}];
AppendTo[tour,{nn,1}]

ct=costcompute[tour]

The nn-th column is the dummy column {-1,0,...0} with zero
cost attached. Now cb is the cost of the base elements.

cb=Join[cb,{ct,0}]

Next the basis matrix is computed. Each column has 1 as 
its first element. The remaining entries are the degrees of
the nodes  {2,3,...,nn} for each one of the (nn-2) 1-trees
and the tour. onetrees records the 1-trees in base.

cv=onetreedegree[tree]

cv[[1]]=1

cv

onetrees={};

addmat=Table[0,{i,nn},{j,nn-2}];

Do[addmat[[2,j]]=1;
   addmat[[j+2,j]]=1;
   thisonetree=Join[tree,{{1,2},{1,j+2}}];
   AppendTo[onetrees,thisonetree];
   ,{j,1,nn-2}]

MatrixForm[addmat]

BT is the transpose of the base matrix B. At the moment it
refers to the first nn-2 rows.

BT=Table[cv+Transpose[addmat][[i]],{i,nn-2}]

Adding the row associated to the tour

tourvect=Table[If[i==1,1,2],{i,nn}];

AppendTo[BT,tourvect]

AppendTo[onetrees,tour];

Adding the dummy row

AppendTo[BT,Table[If[i==1,-1,0],{i,nn}]]

basecheck

Iteration

Basis Change ( to be skipped on the first iteration !  )

newcol=onetreedegree[onetree];newcol[[1]]=1;

newcol

y = W . newcol

Block[{i,ratio,min},
      min = 1000000;
      For[i=1, i<=nn, i++,
          If[ y[[i]]>0,
              ratio = x[[i]]/y[[i]];
              If[ratio < min,
                 min = ratio; q = i]
            ]  
      ]
]

q

BT[[q]]=newcol

cb[[q]]=cx

onetrees[[q]]=onetree;

start here on first iteration !

B=Transpose[BT];

MatrixForm[B];

W = Inverse[B];
MatrixForm[W];

x = W. b

u = cb . W

v=N[u.b]

dual=u;dual[[1]]=0;

dual

newcost=Table[cost[[i,j]]-dual[[i]]-dual[[j]],{i,nn},{j,nn}];

onetreefind[newcost]

treeshow=show[onetree,1,0,"B"];

Show[aaa,treeshow,AspectRatio->Automatic];

onetreedegree[onetree]

cx=costcompute[onetree]

uAx=dual.onetreedegree[onetree]

N[cx-uAx+2 Sum[dual[[i]],{i,nn}]]  (*lower bound on TSP*)

cx-u[[1]]-uAx

If equal to 0 stop, otherwise run second rightmost bracket

Solution Display

posel={};
Do[If[x[[i]]>0,AppendTo[posel,i]],{i,nn-1}];posel

colortable={"B","R"};

treeshow=Table[show[onetrees[[posel[[i]]]],1,0,colortable[[i]]],
{i,Length[posel]}];

Show[GraphicsArray[treeshow],AspectRatio->Automatic];

Show[aaa,treeshow,AspectRatio->Automatic,
PlotRange->rangeplot]