Inequalities implied by a system of linear inequalities/equalities in Matlab: numerical arguments or counterexample?

Question

I have a system of linear inequalities/equalities to solve in Matlab and I use linprog. As some of the inequalities are strict I use a very small costant eps to get get the strict inclusion as explained here

The function solve below solves the system after having provided a value for eps.

function pj=solve(eps)

%Inequalities
%x(1)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)<=0;
%x(2)-x(6)-x(9)-x(12)-x(13)-x(15)-x(17)-x(18)-x(19)<=0;
%x(3)-x(7)-x(10)-x(12)-x(14)-x(16)-x(17)-x(18)-x(19)<=0;
%x(4)-x(8)-x(11)-x(13)-x(14)-x(15)-x(16)-x(18)-x(19)<=0;

%x(1)+x(2)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)-...
%              x(6)-x(12)-x(13)-x(18)<=0;
%x(1)+x(3)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)-...
%              x(7)-x(12)-x(14)-x(18)<=0;
%x(1)+x(4)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)-...
%             x(8)-x(13)-x(14)-x(18)<=0;
%x(2)+x(3)-x(6)-x(9)-x(12)-x(13)-x(15)-x(17)-x(18)-x(19)-...
%              x(7)-x(10)-x(14)-x(16)<=0;
%x(2)+x(4)-x(6)-x(9)-x(12)-x(13)-x(15)-x(17)-x(18)-x(19)-...
%              x(8)-x(11)-x(14)-x(16)<=0;
%x(3)+x(4)-x(7)-x(10)-x(12)-x(14)-x(16)-x(17)-x(18)-x(19)-...
%              x(8)-x(11)-x(13)-x(15)<=0;


%x(1)+x(2)+x(3)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)-...
%                  x(6)-x(12)-x(13)-x(18)-...
%                  x(7)-x(14)<=0;  
%x(1)+x(2)+x(4)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)-...
%                  x(6)-x(12)-x(13)-x(18)-...
%                  x(8)-x(14)<=0;   
%x(1)+x(3)+x(4)-x(5)-x(9)-x(10)-x(11)-x(15)-x(16)-x(17)-x(19)-...
%                  x(7)-x(12)-x(14)-x(18)-...
%                  x(8)-x(13)<=0; 
%x(2)+x(3)+x(4)-x(6)-x(9)-x(12)-x(13)-x(15)-x(17)-x(18)-x(19)-...
%                  x(7)-x(10)-x(14)-x(16)-...
%                  x(8)-x(11)<=0; 


%Equalities
%x(1)+x(2)+x(3)+x(4)=1;
%x(5)+x(6)+x(7)+x(8)+x(9)+x(10)+x(11)+x(12)+x(13)+x(14)+x(15)+x(16)+x(17)+x(18)+x(19)=1;


%I also want each component of x to be different from 1 and 0 (strictly included between 1 and 0 given the equalities constraints)
%x(1)>0 ---> x(1)>=eps ---> -x(1)<=-eps
%...
%x(19)>0
%x(1)<1 ---> x(1)<=1-eps
%...
%x(19)<1

%52 inequalities (14+19+19)
%2 equalities
%19 unknowns

A=[1  0  0  0  -1  0  0  0  -1  -1  -1  0  0  0  -1  -1  -1  0  -1;...
   0  1  0  0   0 -1  0  0  -1   0   0 -1 -1  0  -1   0  -1  -1 -1;...
   0  0  1  0   0  0 -1  0   0  -1   0 -1  0 -1   0  -1  -1  -1 -1;...
   0  0  0  1   0  0  0 -1   0   0  -1  0 -1 -1  -1  -1   0  -1 -1;...

   1  1  0  0  -1 -1  0  0  -1  -1  -1 -1 -1  0  -1  -1  -1  -1 -1;...
   1  0  1  0  -1  0 -1  0  -1  -1  -1 -1  0 -1  -1  -1  -1  -1 -1;...
   1  0  0  1  -1  0  0 -1  -1  -1  -1  0 -1 -1  -1  -1  -1  -1 -1;...
   0  1  1  0   0 -1 -1  0  -1  -1   0 -1 -1 -1  -1  -1  -1  -1 -1;...
   0  1  0  1   0 -1  0 -1  -1   0  -1 -1 -1 -1  -1  -1  -1  -1 -1;...
   0  0  1  1   0  0 -1 -1   0  -1  -1 -1 -1 -1  -1  -1  -1  -1 -1;...

   1  1  1  0  -1 -1 -1  0  -1  -1  -1 -1 -1 -1  -1  -1  -1  -1 -1;...
   1  1  0  1  -1 -1  0 -1  -1  -1  -1 -1 -1  0  -1  -1  -1  -1 -1;...
   1  0  1  1  -1  0 -1 -1  -1  -1  -1 -1 -1 -1  -1  -1  -1  -1 -1;...
   0  1  1  1   0 -1 -1 -1  -1  -1  -1 -1 -1 -1  -1  -1  -1  -1 -1;...

  -1  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0 -1  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0 -1  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0 -1   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0  -1  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0 -1  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0 -1  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0 -1   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0  -1   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0  -1   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0  -1  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0 -1  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0 -1  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0 -1   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0  -1   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0  -1   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0  -1   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0  -1  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0 -1;...

   1  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  1  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  1  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  1   0  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   1  0  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  1  0  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  1  0   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  1   0   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   1   0   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   1   0  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   1  0  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  1  0  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  1  0   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  1   0   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   1   0   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   1   0   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   1   0  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   1  0;...
   0  0  0  0   0  0  0  0   0   0   0  0  0  0   0   0   0   0  1]; %52x19

b=[zeros(1,14) -eps*ones(1,19) (1-eps)*ones(1,19)]; %1x52

Aeq=[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;...
     0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; %2x19
beq=[1 1]; %1x2

f=zeros(1,19); %1x19

x=linprog(f,A,b,Aeq,beq); 

if ~isempty(x)
    pj=x;
else
    pj=NaN;


end

What I believe (but I don't known how to show it analytically) is that the inequalities/equalities I have put into the linprogr algorithm inside the function solve are such that the solution that Matlab produces will satisfy another inequality as indicated below:

clear
rng default

%solve system 
p1=solve(unifrnd(0,0.05));

%solve system
p2=solve(unifrnd(0,0.05));


%solve system
p3=solve(unifrnd(0,0.05));



if ~isnan(p1) & ~isnan(p2) & ~isnan(p3) %#ok<AND2>

%LHS
lhs=(p1(2)+p1(3)+p1(4))*1*1+...
     p1(1)*(p2(1)+p2(4))*1+...
     p1(1)*p2(2)*(p3(2)+p3(3)+p3(4))+...
     p1(1)*p2(3)*(p3(1)+p3(2)+p3(3)); 


%RHS 
 rhs=(1-(p1(5)+p1(9)+p1(10)+p1(11)+p1(15)+p1(16)+p1(17)+p1(19)))*...
      1*...
      1+... 
      ... +
     (p1(5)+p1(9)+p1(10)+p1(11)+p1(15)+p1(16)+p1(17)+p1(19))*...
     (p2(5)+p2(8)+p2(11))*...
      1+...
      ... +
     (p1(5)+p1(9)+p1(10)+p1(11)+p1(15)+p1(16)+p1(17)+p1(19))*...
     (p2(7)+p2(10)+p2(14)+p2(16))*...
     ((p3(5)+p3(9)+p3(10)+p3(17))+(p3(6)+p3(7)+p3(12)))+...
     ... +
     (p1(5)+p1(9)+p1(10)+p1(11)+p1(15)+p1(16)+p1(17)+p1(19))*...
     (p2(6)+p2(9)+p2(13)+p2(15))*...
     ((p3(8)+p3(13)+p3(14)+p3(18))+(p3(6)+p3(7)+p3(12)))+...
     ... +
     (p1(5)+p1(9)+p1(10)+p1(11)+p1(15)+p1(16)+p1(17)+p1(19))*...
     (p2(12)+p2(17)+p2(19))*...
     (p3(6)+p3(7)+p3(12)); 

check=(lhs>=rhs); %I expect check to be 1
else
end

I believe that the solutions p1,p2,p3 will deliver check=1.

Question: As mentioned above, I don't know how to show this analytically; is there a way to produce a satisfactory numerical argument? Or, can you kill my belief and provide solutions p1,p2,p3 delivering check=0?

Halfback

None can help? The bounty is almost expiring – Halfback

Answer 1

What you can do is extend the problem from 19 parameters to 3x19 = 57 parameters, write a function calculating difference between lhs and rhs and try to minimize this function with given constraints. If you go below zero it means rhs > lhs. Extending the constraints:

ACell = repmat({A}, 1, 3); % diagonal matrix with A repeated 3 times on diagonal
A = blkdiag(ACell{:}); % 156x57
b = [b b b]; % 1x156
AeqCell = repmat({Aeq}, 1, 3);
Aeq = blkdiag(AeqCell{:}); % 6x57
beq = [beq beq beq]; % 1x6

Function to calculate error:

function error = errorFun( p )

p1 = p(1:19);
p2 = p(20:38);
p3 = p(39:57);

if ~isnan(p1) & ~isnan(p2) & ~isnan(p3) %#ok<AND2>

    lhs= (...); 
    rhs= (...);
    error= lhs - rhs;

else
    error= 1;
end
end

Now you can use some optimization tool to find minimum of errorFun, for example fmincon:

x = fmincon(@errorFun, zeros(57, 1), A, b, Aeq, beq, [], [], [], options);

Solution that I found was positive, which means your assumption should be right.

Woodbridge