This computes a feasible solution of a system
Ax>=b in two dimensions through the Ellipsoid method
and displays all geometric objects.
Then it computes the optimum for the objective function  cx via a sliding objective function

How to use this notebook

First you need run the Routine section.

Then you need input the problem data in the section Problem data. 
You may either 
-- provide a set of points whose convex hull defines the feasible set; or
-- provide a set of linear inequalities Ax>=b via a matrix A and a rhs vector  b.

In the first case there is a routine which computes the defining inequalities and the vertices. In the second case ther e is a routine which computes the vertices.  Alternatively you may use the data already present in the notebook.

Before starting the computation you need initial values for the Ellipsoid. This is done in the section Initialization.

Then you are ready to find a feasible solution by following carefully the instructions within the section Feasibility Iteration.
Once you have found a feasible solution pass to Optimality Iteration to find an optimal solution

Routines

Off[General::spell1];

Off[ParametricPlot::ppcom];

convexhull[p_]:=Block[{ps,slopes,maxslope,minslope},
   ps=Sort[p];
   firstslope=If[(ps[[2,1]]-ps[[1,1]])>0
                ,
                 (ps[[2,2]]-ps[[1,2]])/(ps[[2,1]]-ps[[1,1]])
                ,If[(ps[[2,2]]-ps[[1,2]])>0,
                     Infinity,
                     -Infinity
                    ]
                 ];
   upperqueue={1,2};
   lowerqueue={1,2};
   upperslopes={Infinity,firstslope};
   lowerslopes={-Infinity,firstslope};
   Do[
      If[(ps[[i,1]]-ps[[upperqueue[[-1]],1]])>0,
            newupperslope=(ps[[i,2]]-ps[[upperqueue[[-1]],2]])/
                    (ps[[i,1]]-ps[[upperqueue[[-1]],1]]),
            newupperslope=If[(ps[[i,2]]-ps[[upperqueue[[-1]],2]])>0,
                Infinity,-Infinity];
        ];
        If[(ps[[i,1]]-ps[[lowerqueue[[-1]],1]])>0,
            newlowerslope=(ps[[i,2]]-ps[[lowerqueue[[-1]],2]])/
                    (ps[[i,1]]-ps[[lowerqueue[[-1]],1]]),
            newlowerslope=If[(ps[[i,2]]-ps[[lowerqueue[[-1]],2]])>0,
                Infinity,-Infinity];
        ]
            ;
      While[newupperslope >= upperslopes[[-1]],
            upperqueue=Drop[upperqueue,-1];
            upperslopes=Drop[upperslopes,-1];
            newupperslope=If[(ps[[i,1]]-ps[[upperqueue[[-1]],1]])>0
                  ,(ps[[i,2]]-ps[[upperqueue[[-1]],2]])/
                    (ps[[i,1]]-ps[[upperqueue[[-1]],1]]) 
                  ,If[(ps[[i,2]]-ps[[upperqueue[[-1]],2]])>0,
                     Infinity,
                     -Infinity]
                    ];
            ];
      AppendTo[upperslopes,newupperslope];
      AppendTo[upperqueue,i];
      While[newlowerslope <= lowerslopes[[-1]],
            lowerqueue=Drop[lowerqueue,-1];
            lowerslopes=Drop[lowerslopes,-1];
            newlowerslope=If[(ps[[i,1]]-ps[[lowerqueue[[-1]],1]])>0
                  ,(ps[[i,2]]-ps[[lowerqueue[[-1]],2]])/
                    (ps[[i,1]]-ps[[lowerqueue[[-1]],1]]) 
                  ,If[(ps[[i,2]]-ps[[lowerqueue[[-1]],2]])>0,
                     Infinity,
                     -Infinity]
                    ];
           ];
      AppendTo[lowerslopes,newlowerslope];
      AppendTo[lowerqueue,i];
      ,{i,3,Length[ps]}
     ];
      upperslopes=Drop[upperslopes,1];
     lowerslopes=Drop[lowerslopes,1];
     A={};
     b={};
     Do[If[Abs[upperslopes[[i]]]==Infinity
              ,If[upperslopes[[i]] > 0
                  ,ns=1;ds=0;
                  ,ns=-1;ds=0;
                  ]
              ,ns=Numerator[upperslopes[[i]]];
               ds=Denominator[upperslopes[[i]]];
           ];
        AppendTo[A,{ns,-ds}];
        AppendTo[b,
          -ds ps[[upperqueue[[i]],2]]+ns ps[[upperqueue[[i]],1]]];
         ,{i,Length[upperqueue]-1}];
     Do[If[Abs[lowerslopes[[i]]]==Infinity
              ,If[lowerslopes[[i]] > 0
                  ,ns=1;ds=0;
                  ,ns=-1;ds=0;
                  ]
              ,ns=Numerator[lowerslopes[[i]]];
               ds=Denominator[lowerslopes[[i]]];
           ];
        AppendTo[A,{-ns,ds}];
        AppendTo[b,
          ds ps[[lowerqueue[[i]],2]]-ns ps[[lowerqueue[[i]],1]]];
        ,{i,Length[lowerqueue]-1}];
     {ps[[Join[upperqueue,Reverse[Drop[lowerqueue,-1]]]]],A,b}
    
]

SetOptions[ParametricPlot,ColorOutput->RGBColor,
          DefaultColor->RGBColor[0,0,1],
          Axes->False];

signvertex[x_,a_,b_]:=Block[{c},
      c=a.x-b;
      Which[c>0,1,
            c<0,-1,
            c==0,0]
  ];

showpol[A_,b_,box_]:=Block[{},
boxsequence={{box[[1,1]],box[[2,1]]},
             {box[[1,2]],box[[2,1]]},
             {box[[1,2]],box[[2,2]]},
             {box[[1,1]],box[[2,2]]},
             {box[[1,1]],box[[2,1]]}};
drawing={RGBColor[0.8,0.8,1]};
Do[signsequence=Table[signvertex[boxsequence[[i]],A[[j]],b[[j]]]
     ,{i,5}];
   newsequence={};
   Do[If[signsequence[[i]]!=signsequence[[i+1]]
        ,lambda=(b[[j]]-A[[j]].boxsequence[[i+1]])/
             (A[[j]].boxsequence[[i]]-A[[j]].boxsequence[[i+1]]);
         AppendTo[newsequence,
          {0,lambda boxsequence[[i]]+(1-lambda) boxsequence[[i+1]]}]
         ,
         ];
        AppendTo[newsequence,{signsequence[[i+1]],boxsequence[[i+1]]}]
      ,{i,4}];
     infeasiblearea={};
     Do[If[newsequence[[i,1]]<=0,
          AppendTo[infeasiblearea,newsequence[[i,2]]]
          ]
        ,{i,Length[newsequence]}];
     AppendTo[drawing,Polygon[infeasiblearea]];

   ,{j,Length[A]}
   ];
  bbox=Graphics[{RGBColor[1,0.2,0.2],Polygon[boxsequence]}];
  aaa=Graphics[drawing];
Show[aaa,PlotRange->box,AspectRatio->Automatic,
Axes->Automatic];
  ];

getvertices[AA_,b_]:=Block[{},
  vertices={};
  firstvertex[AA,b];
  AppendTo[vertices,v];
  firstbase=base;
  end=False;
  While[!end,
         nextvertex[AA,b,base[[1]],base[[2]]];
         AppendTo[vertices,v];
         end=base==firstbase
        ]
];

firstvertex[A_,b_]:=Block[{found,m},
  m=Length[A];
  totset=Table[i,{i,m}];
  found=False;
  ik=1;jk=2;
  While[!found && (ik < m-1),
        base={ik,jk};
        nonbase=Complement[totset,base];
        If[Det[A[[base]]]!=0
          ,v=LinearSolve[A[[base]],b[[base]]];
           found=True;
           Do[found=found &&
                    A[[nonbase[[j]]]].v >= b[[nonbase[[j]]]];
             ,{j,Length[nonbase]}
             ];
          ];
         jk++;
         If[jk>m, ik++;jk=ik+1];
       ];
   If[!found,Print["no feasible vertices !"]];
];

nextvertex[A_,b_,indexin_,indexout_]:=Block[{},
    nonbase=Complement[totset,{indexin,indexout}];
    found=False;
    ik=1;
    While[!found && (ik <= Length[nonbase]),
          base={nonbase[[ik]],indexin};
          If[Det[A[[base]]]!=0
            ,v=LinearSolve[A[[base]],b[[base]]];
             found=A[[indexout]].v>=b[[indexout]];
             If[found
                ,Do[found=found &&
                    (A[[nonbase[[j]]]].v >= b[[nonbase[[j]]]]);
                   ,{j,Length[nonbase]}
                   ];
                ];
            ];
            ik++;
           ]
];

Data

Input as convex of hull of a given point set

points={{15,19},{18,18},{14,17},{15,17},
     {15,18},{12,18},{20,20}};

{vertices,A,b}=convexhull[points]

G4=Graphics[{RGBColor[1,0,0],Line[vertices]}]

Show[G4 ,Axes->Automatic,AspectRatio->Automatic]

Input as linear inequalities Ax>=b

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

b = {-1., -1., -1,1, 0.,0};

getvertices[A,b];

vertices

G4=Graphics[{RGBColor[1,0,0],Line[vertices]}]

Show[G4,AspectRatio->Automatic,Axes->Automatic]

Initialization

This computes the initial radius R which guarantees the inclusion of the feasible set in the first ellipsoid (a circle). Alternatively you may use your own value.

n=2

K=Max[ Abs[{A,b}]]

R=K^2  Sqrt[8]

The ellipsoids are defined as  (x-c).Inverse[Q].(x-c)<=1 The first Q is computed from  the initial radius. The first center c is the origin.

Q = R^2 IdentityMatrix[n]

X=IdentityMatrix[n]
lam={R^2,R^2}

c= {0,0}

centerhistory={c};

G1=ParametricPlot[{R Sin[t],
   R Cos[t]},{t,0,2Pi},
     AspectRatio->Automatic,DisplayFunction->Identity];

C1=Graphics[{RGBColor[0,0,1],PointSize[0.01],Point[c]}];

Show[G1,G4,C1,DisplayFunction->$DisplayFunction,
      AspectRatio->Automatic,
      Axes->Automatic,AxesStyle->RGBColor[0,1,0],
      Ticks->None]

Feasibility Iteration

The procedure is interactive and you are required to provide at each iteration the most violated inequality in the section Feasibility check and cut depth

Feasibility check and cut depth

Run these espressions. You will see the indices of the violated inequalities together with the cut depth. If no output appears then the center c is feasible and you may go to the section Optimality Iteration

G2=G1;
C2=C1;

res=A.c - b

Do[If[res[[i]]<0,
     p=i;
     a = A[[p]];
     bp= b[[p]];
     f = a.Q.a;
     Delta =N[( bp-a.c )/Sqrt[f]];
     Print[i,"  ",Delta]
     ]
  ,{i,Length[res]}
  ]

Provide for p one violated index

p=1

a = A[[p]];
bp= b[[p]];

Next various quantities are computed. The display of yhe new ellipsoid  is in the section Graphical data computation

segment computation

bb=-res[[p]];
DD=DiagonalMatrix[Table[Sqrt[lam[[i]]],{i,2}]];

aa=a.X.DD;

aaa={aa[[2]],-aa[[1]]};

z1=(1/Sqrt[aaa.aaa])X.DD.aaa;

New ellipsoid computation

f = a.Q.a;

z0=(bb/f) Q.a;

Clear[segment]

segment={c+z0-z1,c+z0+z1};

G3=Graphics[{RGBColor[0,1,0],Line[segment]}];

Delta =N[( bp-a.c )/Sqrt[f]]

alfa = n^2(1-Delta^2)/(n^2-1)

beta = (2 +2 n Delta)/((1+Delta)(n+1))

gamma = (1+Delta n)/(n+1)

New matrix

QQ = alfa(Q - (beta/f) Q.Outer[Times,a,a].Q);
MatrixForm[QQ]

New center

cc = c + N[(gamma/Sqrt[f]) Q.a]

AppendTo[centerhistory,cc]

Graphical data computation

Q=QQ;
c=cc;

{lam,XX}=Eigensystem[Q]

r=Table[1/Sqrt[XX[[i]].XX[[i]]],{i,n}]

X=Transpose[XX].DiagonalMatrix[r]

G1=ParametricPlot[c+X.{Sqrt[lam[[1]]]Sin[t],
   Sqrt[lam[[2]]]Cos[t]},{t,0,2Pi},
     AspectRatio->Automatic,DisplayFunction->Identity,
     DefaultColor->RGBColor[0,0,1]];

C1=Graphics[{RGBColor[0,0,1],PointSize[0.01],Point[c]}];

Show[G3,G1,G2,G4,C1,C2,DisplayFunction->$DisplayFunction,
     AspectRatio->Automatic,
     DefaultColor->RGBColor[0,0,1]]

Optimality Iteration

cost coefficients (maximum problem)

cc={1,1}

Feasibility check and cut depth

G2=G1;
C2=C1;

res=A.c - b

Do[If[res[[i]]<0,
     p=i;
     a = A[[p]];
     bp= b[[p]];
     f = a.Q.a;
     Delta =N[( bp-a.c )/Sqrt[f]];
     Print[i,"  ",Delta]
     ]
  ,{i,Length[res]}
  ]

Give violated index or p=0 if no violation (sliding objective function)

p=0

If[p==0
   ,a=cc;bp=cc.c;
   ,a = A[[p]];
    bp= b[[p]];
    ]

segment computation

If[p==0,bb=0,bb=-res[[p]]];
DD=DiagonalMatrix[Table[Sqrt[lam[[i]]],{i,2}]];

aa=a.X.DD;

aaa={aa[[2]],-aa[[1]]};

z1=(1/Sqrt[aaa.aaa])X.DD.aaa;

New ellipsoid computation

f = a.Q.a;

z0=(bb/f) Q.a;

Clear[segment]

segment={c+z0-z1,c+z0+z1};

G3=Graphics[{RGBColor[0,1,0],Line[segment]}];

If[p==0,G5=Graphics[{RGBColor[1,0,0],Line[segment]}]];

Delta =N[( bp-a.c )/Sqrt[f]]

alfa = n^2(1-Delta^2)/(n^2-1)

beta = (2 +2 n Delta)/((1+Delta)(n+1))

gamma = (1+Delta n)/(n+1)

New matrix

QQ = alfa(Q - (beta/f) Q.Outer[Times,a,a].Q);
MatrixForm[QQ]

New center

cc = c + N[(gamma/Sqrt[f]) Q.a]

AppendTo[centerhistory,cc]

Graphical data computation

Q=QQ;
c=cc;

{lam,XX}=Eigensystem[Q]

r=Table[1/Sqrt[XX[[i]].XX[[i]]],{i,n}]

X=Transpose[XX].DiagonalMatrix[r]

G1=ParametricPlot[c+X.{Sqrt[lam[[1]]]Sin[t],
   Sqrt[lam[[2]]]Cos[t]},{t,0,2Pi},
     AspectRatio->Automatic,DisplayFunction->Identity,
     DefaultColor->RGBColor[0,0,1]];

C1=Graphics[{RGBColor[0,0,1],PointSize[0.01],Point[c]}];

Show[G3,G1,G2,G4,G5,C1,C2,DisplayFunction->$DisplayFunction,
     AspectRatio->Automatic,
     DefaultColor->RGBColor[0,0,1]]

Animation

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]

Show[G1,G2,G3,G4,C1,C2,DisplayFunction->$DisplayFunction,
PlotRange->{{-R 1.5,R 1.5},{-R 1.5,R 1.5}}]