简单单纯形法

Posted

技术标签:

【中文标题】简单单纯形法【英文标题】:Simple Simplex method 【发布时间】:2015-01-29 12:34:28 【问题描述】:

我编写了一个求解单纯形法的程序,但它仅适用于约束数量等于或小于目标函数中的变量数量的方程,如果有任何其他方程存在 OutOfBoundsException 而我没有如何解决这个问题呢。如果有人知道,请告诉我或分享工作算法的链接。

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;
            
        
    
    

我输入了一些方程来求解。解决的就对了。 尝试对此进行更改

    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 ;

【问题讨论】:

如果您准确指出OutOfBoundsException 出现的位置,那么您可以更快地得到响应。 temp = calctemp(temp, constraintsLeft, targetFunc, basic); 那个 calctemp 方法在哪里? 同时发布完整的堆栈跟踪。 【参考方案1】:

这是我的 scala 版本。我在退化的情况下尝试过,我认为它支持“品牌规则”。

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();
    


这个用例有效。我试过用一个循环的,它也可以工作......

    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)
    

【讨论】:

以上是关于简单单纯形法的主要内容,如果未能解决你的问题,请参考以下文章

单纯形法剖析,一句话描述单纯形法

matlab 线性规划 单纯形法

matlab 线性规划 单纯形法

线性规划中的单纯形法与内点法(原理步骤以及matlab实现)

GLOP - Google OR 工具 - 选择单纯形法

单纯形法 -- 求解线性规划