Mathematica教學：Simplex method in tableau form

myConstrains[constrains_] := 
 Block[{num, s, mycon}, num = Length[constrains];
  mycon = 
   Table[If[Head[constrains[[i]]] === LessEqual, 
     constrains[[i]][[1]] + s[i] == constrains[[i]][[2]], 
     If[Head[constrains[[i]]] === Equal, 
      constrains[[i]][[1]] + R[i] == constrains[[i]][[2]], 
      constrains[[i]][[1]] - s[i] + R[i] == 
       constrains[[i]][[2]]]], {i, num}];
  mycon /. x__[y_] :> ToExpression[ToString[x] <> ToString[y]]]

LPstd[obj_, mycons_] := 
 Block[{myvar, basicvars, nonbasicvars, AB, BM, c, A, B, rowindex1, 
   rowindex2}, myvar = Variables[mycons[[All, 1]]];
  myvar = 
   Flatten[{Sort@
      Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
     Sort@Select[myvar, StringStartsQ[SymbolName[#], "R"] &]}];
  AB = {Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, Not@StringStartsQ[SymbolName[#], "s" | "R"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "s"] &], 
    Coefficient[mycons[[All, 1]], #] & /@ 
     Select[myvar, StringStartsQ[SymbolName[#], "R"] &]};
  AB = Transpose@Flatten[AB, 1];
  AB = Join[{Flatten[
      Table[Select[myvar, 
        StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]]}, 
    AB];
  c = Flatten[{Table[
      Coefficient[obj, #] & /@ 
       Select[myvar, StringStartsQ[SymbolName[#], i] &], {i, {"x", 
        "s"}}], Table[-BM, {Length[
        Select[myvar, StringStartsQ[SymbolName[#], "R"] &]]}], 0}];
  basicvars = 
   Select[AB[[All, #]] & /@ (Flatten /@ 
        Map[Position[myvar, #] &, 
         Tuples[myvar, {Length[mycons]}], {2}]), #[[2 ;; -1]] === 
       IdentityMatrix[Length[mycons]] &][[1, 1]];
  nonbasicvars = 
   Flatten[Table[
     Select[Complement[myvar, basicvars], 
      StringStartsQ[SymbolName[#], i] &], {i, {"x", "s", "R"}}]];
  c = Flatten@{c[[Flatten[Position[myvar, #] & /@ nonbasicvars]]], 
     c[[Flatten[Position[myvar, #] & /@ basicvars]]], 0};
  AB = Join[{Flatten@{nonbasicvars, basicvars}}, 
    Flatten /@ 
     Transpose@{AB[[2 ;; -1]][[All, 
         Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
       IdentityMatrix[Length[mycons]]}];
  rowindex1 = c[[-Length[basicvars] - 1 ;; -2]];
  rowindex2 = Range[Length[mycons]];
  Flatten[{c - 
       Total[rowindex1[[#]]*
           Flatten[{AB[[2 ;; -1]][[rowindex2[[#]]]], 
             mycons[[rowindex2[[#]], 2]]}] & /@ 
         Range[Length@rowindex1]], 
      AB[[2 ;; -1]][[All, 
        Sort[Flatten[Position[myvar, #] & /@ nonbasicvars]]]], 
      IdentityMatrix[Length[mycons]], mycons[[All, 2]], 
      Flatten@{Sort@
         Select[nonbasicvars, StringStartsQ[SymbolName[#], "x"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "s"] &], 
        Sort@Select[nonbasicvars, 
          StringStartsQ[SymbolName[#], "R"] &]}, basicvars} & /@ 
    Range[2], 1]]

simplexmethod[c_, A_, B_, b_, xnb_, xb_, cp_, Ap_, Bp_, bp_, xnbp_, 
  xbp_] := 
 Block[{cB, nA, nB, mysymbol, pivotCol, AB, leaveVarpos, Bnew, xbnew, 
   xnbnew, Anew, cnew, znew, z, BM, ABp},
   (*初始化變數*)
   nA = Length[xnb];
   nB = Length[xb];
   
  (*建立符號變數*)
  mysymbol = Flatten@{xnbp, xbp};

  (*找出進入變數（最負的檢驗數）*)
  pivotCol = PositionSmallest[-c[[1 ;; -2]] /. BM -> 10^6][[1]];
  
  (*合併A和B矩陣*)
  AB = Flatten /@ Transpose[{A, B}];
  Print@MatrixForm@AB;
  ABp = Flatten /@ Transpose[{Ap, Bp}];
  leaveVarpos = 
   PositionSmallest[
     b/(If[# <= 0, 10^-6, #] & /@ AB[[All, pivotCol]] /. BM -> 10^6) /. 
      x_ /; x <= 0 -> 10^6][[1]];
  
  (*更新基變數和非基變數*)
  xbnew = xb /. xb[[leaveVarpos]] -> mysymbol[[pivotCol]];
  xnbnew = xnb /. mysymbol[[pivotCol]] -> xb[[leaveVarpos]];
  Bnew = 
   ABp[[All, Flatten[Position[mysymbol, #] & /@ Flatten@{xbnew}]]];
  cB = cp[[Flatten[Position[mysymbol, #] & /@ xbnew]]];
  
  (*計算新的檢驗數和目標函數值*)
  znew = Flatten[{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
      cB . Inverse[Bnew], cB . Inverse[Bnew] . bp - cp[[-1]]}] // 
    FullSimplify;
  
  (*顯示單形法表格*)
  Print[Grid[
    Insert[Flatten /@ 
      Transpose[{Flatten[{z, xbnew}], 
        Insert[Map[Flatten, 
          Transpose[{Inverse[Bnew] . Ap, Inverse[Bnew], 
            Inverse[Bnew] . bp}]], znew, 1]}], 
     Flatten[{Item["BasicVariable", 
        Background -> Lighter[Yellow, 0.9]], mysymbol, "Solution"}], 
     1], Frame -> All, 
    Background -> {{Lighter[Gray, 0.8]}, {Lighter[Yellow, 0.9], 
       Lighter[Pink, 0.9]}}, ItemStyle -> Directive[FontSize -> 16], 
    Spacings -> {Table[5, {Length[mysymbol] + 3}], 2}]];
  
  (*返回更新後的參數*)
  {-Flatten@{cB . Inverse[Bnew] . Ap - cp[[1 ;; nA]], 
       cB . Inverse[Bnew]}, Inverse[Bnew] . Ap, Inverse[Bnew], 
    Inverse[Bnew] . bp, xnbnew, xbnew, cp, Ap, Bp, bp, xnbp, xbp} // 
   FullSimplify]

Nest[simplexmethod @@ # &, 
 LPstd[-15 x1 - 20 x2, 
  myConstrains[{x1 + 2 x2 >= 10, 2 x1 - 3 x2 <= 6, x1 + x2 >= 6}]], 2]

Tags: Mathematica, Mathematica 教學, Mathematica 教學 程式, Mathematica 教學講義

• 編輯文章
• 此文章網址
• 轉寄此文章
• 目前有0 意見