Simple Simplex method

Question

I wrote a program which is solving the Simplex Method but it works only on equations where the number of constraints is equal or less then the number of variables in the target function. If there is any other equation there is an OutOfBoundsException and I don't know how to solve this problem. If someone knows please tell me or share the link to the working algorithm.

private static int ROW;

private static int COL;

private static Scanner scanner = new Scanner(System.in);

private static double[] calctemp(double[] temp, double[][] constLeft,
        double[] targetFunc, int[] basic) {
    double[] calcTemp = new double[temp.length];
    for (int i = 0; i < COL; i++) {
        calcTemp[i] = 0;
        for (int j = 0; j < ROW; j++) {
            calcTemp[i] += targetFunc[basic[j]] * constLeft[j][i];
        }
        calcTemp[i] -= targetFunc[i];
    }
    return calcTemp;
}

private static int minimum(double[] arr) {
    double arrmin = arr[0];
    int minPos = 0;
    for (int i = 0; i < arr.length; i++) {
        if (arr[i] < arrmin) {
            arrmin = arr[i];
            minPos = i;
        }
    }
    return minPos;
}

private static void printFrame(double[] targetFunc) {
    StringBuilder sb = new StringBuilder();
    sb.append("Cj\t\t\t");
    for (int i = 0; i < targetFunc.length; i++) {
        sb.append(targetFunc[i] + "\t");
    }
    sb.append("\ncB\txB\tb\t");
    for (int i = 0; i < targetFunc.length; i++) {
        sb.append("a" + (i + 1) + "\t");
    }
    System.out.print(sb);
}

private static void printAll(double[] targetFunc, double[] constraintRight,
        double[][] constraintLeft, int[] basic) {
    printFrame(targetFunc);
    StringBuilder sb = new StringBuilder();
    for (int i = 0; i < ROW; i++) {
        sb.append("\n" + targetFunc[basic[i]] + "\tx" + (basic[i] + 1)
                + "\t" + constraintRight[i] + "\t");
        for (int j = 0; j < COL; j++) {
            sb.append(constraintLeft[i][j] + "\t");
        }
        sb.append("\n");
    }
    System.out.println(sb);
}

public static void main(String[] args) {
    double[] targetFunc = { 6, -5, 0, 0};
    ROW = 2;
    COL = 2 + ROW;
    double[][] constraintsLeft = { { 2, 5, 1, 0 },
            { 5, 2, 0, 1 }};
    double[] constraintsRight = { 10, 10 };

    double[] temp = new double[COL];

    int tempMinPos;
    double[] miniRatio = new double[ROW];
    int miniRatioMinPos = 0;
    double key;
    int goOutCol = 0;
    double z;
    double[] x = new double[COL];
    int[] basic = new int[ROW];
    int[] nonBasic = new int[ROW];
    boolean flag = false;

    for (int i = 0; i < ROW; i++) {
        basic[i] = (i + ROW);
        nonBasic[i] = i;
    }
    System.out.println("------------Calculating------------");
    while (!flag) {
        z = 0;
        temp = calctemp(temp, constraintsLeft, targetFunc, basic);

        tempMinPos = minimum(temp);
        printAll(targetFunc, constraintsRight, constraintsLeft, basic);
        System.out.print("Zj-Cj\t\t\t");
        for (int i = 0; i < COL; i++) {
            System.out.print(temp[i] + "\t");
        }
        System.out
                .println("\n--------------------------------------------------");
        System.out.println("Basic variables : ");
        for (int i = 0; i < ROW; i++) {
            x[basic[i]] = constraintsRight[i];
            x[nonBasic[i]] = 0;
            System.out.println("x" + (basic[i] + 1) + " = "
                    + constraintsRight[i]);
        }
        for (int i = 0; i < ROW; i++) {
            z = z + targetFunc[i] * x[i];
        }
        System.out.println("Max(z) = " + z);

        for (int i = 0; i < ROW; i++) {
            if (constraintsLeft[i][tempMinPos] <= 0) {
                miniRatio[i] = 999;
                continue;
            }
            miniRatio[i] = constraintsRight[i]
                    / constraintsLeft[i][tempMinPos];
        }
        miniRatioMinPos = minimum(miniRatio);

        for (int i = 0; i < ROW; i++) {
            if (miniRatioMinPos == i) {
                goOutCol = basic[i];
            }
        }
        System.out.println("Outgoing variable : x" + (goOutCol + 1));
        System.out.println("Incoming variable : x" + (tempMinPos + 1));

        basic[miniRatioMinPos] = tempMinPos;
        nonBasic[tempMinPos] = goOutCol;

        key = constraintsLeft[miniRatioMinPos][tempMinPos];
        constraintsRight[miniRatioMinPos] /= key;
        for (int i = 0; i < COL; i++) {
            constraintsLeft[miniRatioMinPos][i] /= key;
        }
        for (int i = 0; i < ROW; i++) {
            if (miniRatioMinPos == i) {
                continue;
            }
            key = constraintsLeft[i][tempMinPos];
            for (int j = 0; j < COL; j++) {
                constraintsLeft[i][j] -= constraintsLeft[miniRatioMinPos][j]
                        * key;
            }
            constraintsRight[i] -= constraintsRight[miniRatioMinPos] * key;
        }

        for (int i = 0; i < COL; i++) {
            flag = true;
            if (temp[i] < 0) {
                flag = false;
                break;
            }
        }
    }
}    

I entered some equation to solve. It's is solved right. Try to change on this

    double[] targetFunc = { 8, 2, 0, 0, 0};
    
    ROW = 3;
    COL = 2 + ROW;
    double[][] constraintsLeft = { { 1, -4, 1, 0, 0 },
            { -4, 1, 0, 1, 0 },
            { 1, 1, 0, 0, 1}};
    double[] constraintsRight = { 4, 4, 6 };

Answer 1

Here is my scala version. I tried it on degenerated case and I think it supports the "Brand Rule".

object Simplex {

    sealed trait Pivot {};
    case class Next(row: Int, col: Int) extends Pivot;
    object NoSolution extends Pivot;
    object NoMore extends Pivot;


    def minSuch[T,U](array: Array[T])(fn: (T,Int)=>Option[U])(implicit order: scala.math.Ordering[U]): Option[(Int, T, U)] = {
        @scala.annotation.tailrec
        def compute(idx: Int, res: Option[(Int, T, U)]): Option[(Int, T, U)] = if(idx>=array.length) res else (res, fn(array(idx), idx)) match {
            case (r , None) => compute(idx+1,r)
            case (r @ Some((_, _, u1)), Some(u2)) if order.lt(u1, u2) => compute(idx+1, r)
            case (_ , Some(u)) => compute(idx+1, Some((idx, array(idx), u)))
        }
        return compute(0, Option.empty[(Int, T, U)])
    }

    def solve[T](A: Array[Array[T]], Y: Array[T], C: Array[T])(implicit frac:scala.math.Fractional[T], classtag: scala.reflect.ClassTag[T]) : Option[(T,Array[T])] = {
        import scala.math.Fractional.Implicits._
        import scala.math.Ordering.Implicits._

        val N = (0 to  (C.length-1) by +1).toArray
        val B = (1 to -(Y.length  ) by -1).toArray
        val Z = C.map(-_)
        var z = frac.zero

        def pivot(): Pivot = minSuch(Z) { case (z,_) => 
            if( z<frac.zero ) Some(z) else None
        }.map { case (col, _, _) => 
            minSuch(A) { case(cells,row) =>
                if( cells(col)>frac.zero ) Some(Y(row)/cells(col)) else None
            }.map { case (row, _, _) =>
                new Next(row, col)
            }.getOrElse(NoSolution)
        }.getOrElse(NoMore)

        @scala.annotation.tailrec
        def resolve(): Option[(T, Array[T])] = pivot() match {
            case NoSolution => None
            case NoMore => {
                Some((z, Y.zip(B).foldLeft(Array.fill(C.length)(frac.zero))( (result, yb)=> 
                    if( yb._2 >= 0 ) result.updated(yb._2, yb._1) else result
                )))
            }
            case Next(row, col) => {
                val coef = A(row)(col)
                val tmp = B(row)
                B(row) = N(col)
                N(col) = tmp

                Z(col) = -Z(col) / coef
                z = z + Z(col) * Y(row)
                for(c <- 0 to Z.length-1 if(c!=col)) Z(c) = Z(c) + A(row)(c) * Z(col)

                Y(row) =  Y(row) / coef
                for(r <- 0 to Y.length-1 if(r!=row)) Y(r) = Y(r) - A(r)(col) * Y(row)

                A(row)(col) = frac.one / coef
                for(c <- 0 to A(row).length-1 if(col!=c) ) A(row)(c)=A(row)(c)/coef
                for(r <- 0 to A.length-1 if(row!=r); c <- 0 to A(r).length-1 if(col!=c)) A(r)(c)=A(r)(c) - A(r)(col) * A(row)(c)
                for(r <- 0 to A.length-1 if(row!=r)) A(r)(col) = -A(r)(col) / coef

                return resolve()
            }
        }

        return resolve();
    }

}

This use case works. I've tried with a cyclic one and it works too...

    Simplex.solve(
        Array(
            Array(Rational(1), Rational(1)),
            Array(Rational(1), Rational(-2)),
            Array(Rational(-1), Rational(4))
        ),
        Array(
            Rational(2),
            Rational(0),
            Rational(1)
        ),
        Array(Rational(5), Rational(8))
    ).foreach { result=>
        println(result._1)
        result._2.foreach(println)
    }