function [xmin, fval, flag] = LinProg(f, AInEq, bInEq, AEq, bEq, debug)
if nargin < 6
debug = false;
end
if nargin < 5
bEq = [];
end
if nargin < 4
AEq = [];
end
InEq = ~isempty(AInEq);
Eq = ~isempty(AEq);
if (InEq)
validate_input(AInEq,bInEq,f);
A = [AInEq, eye(size(AInEq, 1))];
c = [f, zeros(1, size(AInEq, 1))];
b = bInEq;
end
if (Eq)
if (InEq)
validate_input(AEq,bEq,f);
p = size(A, 2) - size(AEq, 2);
A = [A; [AEq, zeros(size(AEq, 1), p)]];
b = [b; bEq];
else
A = AEq;
c = f;
b = bEq;
end
end
for j=1:length(b)
if b(j)<0
A(j,:)=-A(j,:);
b(j)=-b(j);
end
end
[A,b,c] = make_phase_one(A,b,c,debug);
try
[throw_away, x] = SimplexPhase2(A,b,c,debug,'phase two');
xmin = x(1:numel(f));
fval = f*xmin;
flag = 'Solution was found';
catch
xmin = x(1:numel(f));
fval = f*xmin;
flag = 'Solution was not found: primal is unbounded';
end
end
function [A,b,c] = make_phase_one(A,b,c,debug)
[m,n] = size(A);
tab = zeros(m+2,n+m+1);
tab(1:m,1:n) = A;
tab(end,n+1:end-1) = 1;
tab(1:m,end) = b(:);
tab(m+1,1:n) = c;
tab(1:m,n+1:n+m) = eye(m);
if debug
fprintf('Current tableau [phase one]\n');
disp(tab);
end
for i = 1:m
tab(end,:) = tab(end,:)-tab(i,:);
end
A = tab(1:m+1,1:n+m);
b = tab(1:m+1,end);
c = tab(end,:);
tab = SimplexPhase1(A,b,c,debug,'phase one');
A = tab(1:m,1:n);
b = tab(1:m+1,end);
c = tab(m+1,1:n);
end
function tab = SimplexPhase1(A,b,c,debug,phase_name)
[m,n] = size(A);
tab = zeros(m+1,n+1);
tab(1:m,1:n) = A;
tab(m+1,:) = c;
tab(1:m,end) = b(:);
keep_running = true;
while keep_running
if debug
fprintf('***********************\n');
fprintf('Current tableau [%s] \n',phase_name);
disp(tab);
end
if any(tab(end,1:n)<0)
[throw_away,J] = min(tab(end,1:n));
if all(tab(1:m-1,J)<=0)
error('problem unbounded. All entries <= 0 in column %d',J);
else
pivot_row = 0;
min_found = inf;
for i = 1:m-1
if tab(i,J)>0
tmp = tab(i,end)/tab(i,J);
if tmp < min_found
min_found = tmp;
pivot_row = i;
end
end
end
if debug
fprintf('pivot row is %d\n',pivot_row);
fprintf('pivot col is %d\n',J);
end
tab(pivot_row,:) = tab(pivot_row,:)/tab(pivot_row,J);
for i=1:m+1
if i ~= pivot_row
tab(i,:)=tab(i,:)-tab(i,J)*tab(pivot_row,:);
end
end
end
if debug
fprintf('current basic feasible solution is\n');
disp(get_current_x());
end
else
keep_running=false;
end
end
function current_x = get_current_x()
current_x = zeros(n,1);
for j=1:n
if length(find(tab(:,j)==0))==m
idx= tab(:,j)==1;
current_x(j)=tab(idx,end);
end
end
end
end
function [tab, x] = SimplexPhase2(A,b,c,debug,phase_name)
[m,n] = size(A);
tab = zeros(m+1,n+1);
tab(1:m,1:n) = A;
tab(m+1,1:n) = c;
tab(1:m+1,end) = b(:);
inneriter = 0;
keep_running = true;
while keep_running
if debug
fprintf('***********************\n');
fprintf('Current tableau [%s] \n',phase_name);
disp(tab);
end
if any(tab(end,1:n)<0)
[throw_away,J] = min(tab(end,1:n));
if all(tab(1:m,J)<=0)
error('problem unbounded. All entries <= 0 in column %d',J);
else
pivot_row = 0;
min_found = inf;
for i = 1:m
if tab(i,J)>0
tmp = tab(i,end)/tab(i,J);
if tmp < min_found
min_found = tmp;
pivot_row = i;
end
end
end
if debug
fprintf('pivot row is %d\n',pivot_row);
fprintf('pivot col is %d\n',J);
end
tab(pivot_row,:) = tab(pivot_row,:)/tab(pivot_row,J);
for i=1:m+1
if i ~= pivot_row
tab(i,:)=tab(i,:)-tab(i,J)*tab(pivot_row,:);
end
end
end
x = get_current_x();
if debug
fprintf('current basic feasible solution is\n');
disp(x);
end
elseif(any(sum(abs(tab(1:m,1:n-m)))) > 1 && inneriter < m*5)
colmax = zeros(1,n-m);
for i = 1:n-m
colmax(i) = max(tab(:,i));
end
[throw_away,J] = max(colmax);
pivot_row = 0;
min_found = inf;
for i = 1:m
if tab(i,J)>0
tmp = tab(i,end)/tab(i,J);
if tmp < min_found
min_found = tmp;
pivot_row = i;
end
end
end
if debug
fprintf('pivot row is %d\n',pivot_row);
fprintf('pivot col is %d\n',J);
end
tab(pivot_row,:) = tab(pivot_row,:)/tab(pivot_row,J);
for i=1:m+1
if i ~= pivot_row
tab(i,:)=tab(i,:)-tab(i,J)*tab(pivot_row,:);
end
end
else
x = get_current_x();
keep_running=false;
end
end
function current_x = get_current_x()
current_x = zeros(n,1);
for j=1:n
if length(find(tab(:,j)==0))==m
idx= tab(:,j)==1;
current_x(j)=tab(idx,end);
end
end
end
end
function validate_input(A,b,c)
if ~isvector(b)
error('b must be vector');
end
if ~isvector(c)
error('c must be vector');
end
[m,n]=size(A);
if length(b) ~= m
error('b must have same size as number of rows of A');
end
if length(c) ~= n
error('c must have same size as number of columns of A');
end
end