Wsrtjtyk

I am trying to implement a simplex algorithm following the rules I was given at my optimization course. The problem is

min c'*x    s.t.

    Ax = b

    x >= 0

All vectors are assumes to be columns, ' denotes the transpose. The algorithm should also return the solution to dual LP. The rules to follow are:

Here, A_J denotes columns from A with indices in J and x_J, x_K denotes elements of vector x with indices in J or K respectively. Vector a_s is column s of matrix A.

Now I do not understand how this algorithm takes care of condition x >= 0, but I decided to give it a try and follow it step by step. I used Matlab for this and got the following code.

X = zeros(n, 1);

Y = zeros(m, 1);



% i. Choose starting basis J and K = {1,2,...,n}  J

J = [4 5 6]  % for our problem

K = setdiff(1:n, J)



% this while is for goto

while 1

    % ii. Solve system A_J*bar{x}_J = b.

    xbar = A(:,J)  b



    % iii. Calculate value of criterion function with respect to current x_J.

    fval = c(J)' * xbar



    % iv. Calculate dual solution y from A_J^T*y = c_J.

    y = A(:,J)'  c(J)



    % v. Calculate bar{c}^T = c_K^T - u^T A_K. If bar{c}^T >= 0, we have

    % found the optimal solution. If not, select the smallest s in K, such

    % that c_s < 0. Variable x_s enters basis.

    cbar = c(K)' - c(J)' * inv(A(:,J)) * A(:,K)

    cbar = cbar'



    tmp = findnegative(cbar)

    if tmp == -1  % we have found the optimal solution since cbar >= 0

        X(J) = xbar;

        Y = y;

        FVAL = fval;

        return

    end



    s = findnegative(c, K)  %x_s enters basis



    % vi. Solve system A_J*bar{a} = a_s. If bar{a} <= 0, then the problem is

    % unbounded.

    abar = A(:,J)  A(:,s)



    if findpositive(abar) == -1  % we failed to find positive number

        disp('The problem is unbounded.')

        return;

    end



    % vii. Calculate v = bar{x}_J / bar{a} and find the smallest rho in J,

    % such that v_rho > 0. Variable x_rho exits basis.

    v = xbar ./ abar

    rho = J(findpositive(v))



    % viii. Update J and K and goto ii.

    J = setdiff(J, rho)

    J = union(J, s)

    K = setdiff(K, s)

    K = union(K, rho)

end

Functions findpositive(x) and findnegative(x, S) return the first index of positive or negative value in x. S is the set of indices, over which we look at. If S is omitted, whole vector is checked. Semicolons are omitted for debugging purposes.

The problem I tested this code on is

c = [-3 -1 -3 zeros(1,3)];

A = [2 1 1; 1 2 3; 2 2 1];

A = [A eye(3)];

b = [2; 5; 6];

The reason for zeros(1,3) and eye(3) is that the problem is inequalities and we need slack variables. I have set starting basis to [4 5 6] because the notes say that starting basis should be set to slack variables.

Now, what happens during execution is that on first run of while, variable with index 1 enters basis (in Matlab, indices go from 1 on) and 4 exits it and that is reasonable. On the second run, 2 enters the basis (since it is the smallest index such that c(idx) < 0 and 1 leaves it. But now on the next iteration, 1 enters basis again and I understand why it enters, because it is the smallest index, such that c(idx) < 0. But here the looping starts. I assume that should not have happened, but following the rules I cannot see how to prevent this.

I guess that there has to be something wrong with my interpretation of the notes but I just cannot see where I am wrong. I also remember that when we solved LP on the paper, we were updating our subjective function on each go, since when a variable entered basis, we removed it from the subjective function and expressed that variable in subj. function with the expression from one of the equalities, but I assume that is different algorithm.

Any remarks or help will be highly appreciated.

The problem has been solved. Turned out that the point 7 in the notes was wrong. Instead, point 7 should be