Message Boards Message Boards

Find a feasible solution with LinearProgramming function?

GROUPS:

I have got a specific LP optimization problem as attached picture.enter image description here

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".

enter image description here

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.

POSTED BY: Kun Fang
Answer
9 months ago

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
POSTED BY: Todd Rowland
Answer
9 months ago

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?

POSTED BY: Kun Fang
Answer
9 months ago

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.

POSTED BY: Todd Rowland
Answer
9 months ago

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].)

POSTED BY: Daniel Lichtblau
Answer
9 months ago

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!

POSTED BY: Kun Fang
Answer
9 months ago

Group Abstract Group Abstract