-
Notifications
You must be signed in to change notification settings - Fork 0
/
phase1.m~
103 lines (91 loc) · 2.47 KB
/
phase1.m~
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
function [ Bx, xB,message, iter ] = phase1( A,b,c)
%function [ Bx,xB,message, iter ] = phase1( A,b,c)
%
% Phase 1 des Simplexverfahrens
%
% Input: A, b, c - Daten für LP in primaler Standardform
% min c'x s.t. Ax=b, x>=0
% Output: Bx, xB - primal zulaessige Basis, zugehörige Basislösung
% message - Information über Optimallösung oder Unbeschraenktheit
% iter - Anzahl der Iterationen
%
% Autoren: Vladyslav Yushchenko, Jonas Molina Ramirez, Florian Beck
%% Test data: remove before submission
% clear;
% A = [1 1 -1 0;
% 3 -1 0 0;
% 4 2 0 1];
% b = [1 20 100]';
% c = [-2 -1 0 0]';
%% Check compatibility of dimensions
size_b = size(b);
size_A = size(A);
size_c = size(c);
%%
if(size_b(1) ~= size_A(1))
error('Dimensions of A and b are incompatible.')
end
if(size_c(1) ~= size_A(2))
error('Dimensions of A and c are incompatible.')
end
% Check full rank of matrix A
rank_A = rank(A);
if(size_A(1) ~= rank_A)
error('A does not have full row rank.');
end
%% check for non-negative b
negative_b_indizes = find(b<0);
A(negative_b_indizes,:) = A(negative_b_indizes,:)*-1;
b(negative_b_indizes) = b(negative_b_indizes)*-1;
%% prepare function parameters for simplex
n = size_c(1);
m = size_b(1);
input_A = [A eye(size_A(1))];
B = n+1:n+m;
input_c = [zeros(1,n) ones(1,m)];
%% call simplex
[xopt,result_B,messageSimplex,iterSimplex] = primalSimplex(input_A,b,input_c,B);
%% Set iter = 0 for phase 1 algorithm
iter = 0;
%% compute x, y, result_N
x = xopt(1:n);
y = xopt(n+1:n+m);
xB = x(result_B);
if(y ~= 0)
message = 'LP inadmissible';
return;
end
result_N = setdiff(1:n+m,result_B);
%%
Bx = intersect(result_B,1:n,'stable');
By = setdiff(result_B,Bx,'stable');
Nx = intersect(1:n,result_N);
Ny = setdiff(result_N,Nx);
%% Solange :-)
while(length(Bx) < m)
iter = iter + 1; % increment iter
indices_leq_i = [];
j= [];
i = [];
for temp_j = Nx
i = m - length(Bx);
w = [ones(By-n) A(:,Bx)]\A(:,j)
indices_of_w_unequal_zero = find(w ~= 0);
indices_leq_i = indices_of_w_unequal_zero(indices_of_w_unequal_zero <= i);
if(~empty(indices_leq_i))
j = temp_j;
break;
end
end
if(~empty(j))
Bx = union(Bx,j);
By = setdiff(By,By(i));
Nx = setdiff(Nx,j);
continue;
else
message = sprintf('Rank(A) < %d',m);
return;
end
end
message = 'Found base Bx and base solution xB';
return;