Question: Trouble with the simplex package?

While teaching a linear programming course I put together a worksheet to illustrate finding the largest disk inside a conve polygon as in section 2.6 of Understanding and Using Linear Programming by Jiří Matoušek and Bernd Gärtner.  I used both the Optimization[LPSolve] and the simplex[maximize] routines with the same objective and the same constraints.  Optimization[LPSolve] gives the correct answer but simplex[maximize] does not.  Is this a bug or did I do something wrong?

Below is the worksheet.

This is a worksheet I put together to illustrate finding the largest disk inside a convex polygon

as in section 2.6 of Understanding and Using Linear Programming by Jiří Matoušek and
Bernd Gärtner.

 

Load some packages

restart;
with(plots):
with(LinearAlgebra):
with(simplex);
with(Optimization);
with(plottools):

[basis, convexhull, cterm, define_zero, display, dual, feasible, maximize, minimize, pivot, pivoteqn, pivotvar, ratio, setup, standardize]

 

[ImportMPS, Interactive, LPSolve, LSSolve, Maximize, Minimize, NLPSolve, QPSolve]

(1)

This is the definition of the polygon with a plot.

constr := [x>=0,y>=0,y-x<=1,x+y<=5,2*x-y<=6];
Vector(%);
P1 := inequal(constr,x=-1..5,y=-1..5,optionsfeasible=[color=yellow]);

These are the outward normals for the polygon constraints and the "b" values.

C := [<-1,0>,<0,-1>,<-1,1>,<1,1>,<2,-1>];
b:= [0,0,1,5,6];

This is a check that a point is inside the polygon with plots of the circle that show the distance

from the pont to each defining line.

Pt := <1,.5>;
Pta := convert(Pt,list);
map(is,subs({x=2,y=1},constr));

DistV := [seq((b[i]-DotProduct(C[i],Pt))/Norm(C[i],2),i=1..5)]

P2 := NULL:
for i from 1 to 5 do
   P2 := P2,plottools:-circle(Pta,DistV[i],color=cat("Bright ",i),thickness=2);
   end do:
display([P2,P1],scaling=constrained);
display(map(display,[seq([P1,[P2][i]],i=1..5)]),insequence=true,scaling=constrained);

Now the distance from a point in the polygon to the boundary of the polygon is defined

and plotted over the polygon.

f := proc(x,y)
     local i;
      uses LinearAlgebra;
     min([seq(((b[i]-DotProduct(C[i],<x,y>))/Norm(C[i],2)),i=1..5)]);
     end proc;

plot3d(f(x,y),x=-1..5,y=-1..5,view=[0..4,0..4,0..2],grid=[151,151]);

 

This is the definition of the LP to solve the largest disk problem.

 

We maximize z, a variable less than or equal to the distance from the point x, y to each of the

constraining lines where x, yis in the polygon.

 

C := [<-1,0>,<0,-1>,<-1,1>,<1,1>,<2,-1>];
b:= [0,0,1,5,6];

Here are the constraints.

obj := z;
constr1 := [x>=0,y>=0,y-x<=1,x+y<=5,2*x-y<=6];
constr2 := [seq(z<=(b[i]-DotProduct(C[i],<x,y>))/Norm(C[i],2),i=1..5)];

This is the solution using LPSolve.

soln1 := LPSolve(z,[op(constr1),op(constr2)],maximize);

This is the same problem solved with simplex[maximize].

soln2 := simplex:-maximize(z,[op(constr1),op(constr2)]);
evalf(%);

The answer from the simplex package is wrong as the following plots demonstrate.

 

What, if anything, did I do wrong?

constr := [x>=0,y>=0,y-x<=1,x+y<=5,2*x-y<=6]:
P1 := inequal(constr,x=-1..5,y=-1..5,optionsfeasible=[color=yellow]):
Pt1 := subs(soln1[2],[x,y]);
C1 := disk(Pt1,soln1[1],color=red):
display(C1,P1);
Pt2 := subs(soln2,[x,y]);
r2 := subs(soln2,z);
C2 := disk(Pt2,r2,color=red):
display(C2,P1);

 

 

Download inscribe_circle.mw

Please Wait...