Find a feasible solution with LinearProgramming function?

Questions
,
1

I have got a specific LP optimization problem as attached picture.

LP.png
LP.png563×454 14.3 KB

Given a integer n and \varepsilon,\ p \in (0,1), this is a LP optimization over 2n+3 variables. I am trying to solve this with Mathematica LinearProgramming function. When n is less than 18, everything is fine and it gives an optimal solution. But when n is greater or equal than 19, it gives an warning that “no solution can be found that satisfies the constraint”.

warning.png
warning.png753×131 3.92 KB

Apparently, there is feasible solution that y_i = 0 and r_i = p_i. I believe that Mathematica can handle hundreds of variables without doubt.

But why it does not give a solution when n is over 19 in this case? Btw, I have also posted this problem on Stackexchange Mathematica which can be found here (http://mathematica.stackexchange.com/questions/139976/linearprogramming-function-cannot-find-a-feasible-solution).

My code is as followed:

ClearAll["Global`*"];

n = 3; e = 0.01; p = 2/3;

u[k_,i_]:=If[k>=i,Binomial[k,i]*2^i,0];
c = Table[u[i,j],{i,0,n},{j,0,n}];
zblock = SparseArray[{},{n+1,n+1}];
zcolum = SparseArray[{},{n+1,1}];
ones = ConstantArray[1,{n+1,1}];

(*Constraints*)
c1 = Join[c,zblock,zcolum,2];
c2 = Join[zblock,c,zcolum,2];
c3 = Join[zblock,-c,ones,2];
c4 = {Flatten[Join[SparseArray[{},{1,n+1}],Reverse[c[[n+1]]],{0}]]};
c5 = {Flatten[Join[Reverse[c[[n+1]]],SparseArray[{},{1,n+1}],{0}]]};
c6 = Join[c,-c,zcolum,2];
M = Join[c1,c2,c3,c4,c5,c6];

(*Objective*)
ob = Flatten[Join[SparseArray[{},{1,2*(n+1)}],{{1}},2]];

(*Bounds*)
b1 = ConstantArray[{0,1},3*(n+1)];
b2 = {{1,0},{e,-1}};
(*b31 = {Table[-Sum[Binomial[k,i]*2^i*(p/2)^(n-i)*(1-p)^i,{i,0,k}],{k,0,n}]};*)
(*b31 = {Table[-(p/2)^(n-k)*(2-3*p/2)^k,{k,0,n}]};*)
b31 = {Table[-(1/3)^(n-k),{k,0,n}]};
b3 = Join[Transpose[b31],ones,2];
bounds = Join[b1,b2,b3];
bounds[[All,1]] = Total[M[[All,1;;2*n+2]],{2}]+bounds[[All,1]];
ob = SparseArray[ob];

L = LinearProgramming[ob,M,bounds];
ob.L

Note: this problem has been solved by replacing e = 0.01 to e = 1/100 as suggested by Mr. Lichtblau.

2

I cannot say for certain what the underlying algorithm is, but it is possible that the numbers in your second argument are too big for it when n=19.

For example, even for the small example you gave with n=3,

Det[Transpose[M].M] == 4971003904000
3

Thanks for your reply, Todd. Yes, big numbers (coming from Binomial term in the constraints) and small numbers (coming from p_i) will show up when n goes to large. If this is the case, is there any better way to handle such problem with extreme numbers?

4

According to the documentation, the arguments are supposed to be real numbers. That means Reals in Wolfram Language. If you give it integers or rationals (like your code is doing) the documentation says it tries to give an exact rational result. For example, instead of

Binomial[19, 9] 2^19 (*48432676864*)

you would use

N[Binomial[19, 9] 2^19]

(better to just use N on the arguments ob, M) but you’d want to make sure the Precision is high enough. See the documentation for N and try using the Tolerance option for LinearProgramming.

5

If you use e=1/100 this will deliver a result for the case n=19. But there is something unclear about the description. Is this intended as an integer linear programming (ILP) problem? If so, it was not set up quite right since LinearProgramming needs to be told that. (Point of confusion: positing p is an element of (0,1) seems to mean it lies in the set of two values 0 and 1; if the continuum was intended it should be [0,1].)

6

Thanks for your replying, Mr. Lichtblau! Your suggestion completely solves my problem and I can finally relieve a little bit now. Since I am beginner to Mathematica, I did not realize the difference between 0.01 and 1/100. Really learnt something from you.

Btw, this LP is not intended as an integer linear programming. The parameter n is given before the optimization start. I just want to observe how the optimal value changes when n goes to large. So I can plot a 2-D picture to show this change. As for p, it is just a fixed parameter which lies between zero and one in my problem.

Thanks again!