Algorithms/LinearProgramming/dual_simplex.m at main · fdistasio/Algorithms · GitHub

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
function [B, x, y, status] = dual_simplex(A, b, c, B)

    % Assume that B is dual feasible
    n = length(b);
    m = length(c);

    while(true)

        teta = +inf;
        y = zeros(1, n);

        A_B = A(B, :);
        b_B = b(B, :);
        A_N = A;
        A_N(B, :) = [];
        b_N = b;
        b_N(B, :) = [];

        B = sort(B);

        % Primal solution
        A_B_inv = inv(A_B);
        x = A_B_inv * b_B;

        % Dual solution
        for i=1:n
            for j=1:m
                if(i == B(j))
                    y(i) = c * A_B_inv(:, j);
                end
            end
        end

        % End if x is primal feasible
        A_N_x = A_N*x;
        if(all(A_N_x <= b_N))
            status = "optimal";
            break;
        end

        % Incoming index k
        for i=1:n-m
            if(A_N_x(i) > b_N(i))
                k = i;
                break;
            end
        end

        % eta_B
        eta_B = A(k, :) * A_B_inv;

        % End if Primal is empty
        if(all(eta_B <= 0))
            status = "Primal empty";
            break;
        end

        % Outgoing index h
        for i=1:n
            if(y(i) == 0)
                continue;
            end

            if (eta_B(i) > 0)
                teta_i = y(i)/eta_B(i);

                if(teta_i < teta)
                    teta = teta_i;
                    h = i;
                end
            end
        end

        % Change Basis
        B(find(B == h)) = k;

    end

end