Mathematica 教學 Linear Proggramming Simplex Method

LinearProggramming這個指令在Mathematica裡面無法像NMinimize可以輸出演算法的迭代過程，所以自己寫了一下。目標函數在下面simplex這個函數是以極大化為主，所以若是極小化的問題必須改變一下。

simplex[varstemp_List,bvartemp_List,nonbvartemp_List,
simplexmatrixtemp_List,objcoeftemp_List,soltemp_List,
objsoltemp_]:=
Block[{vars=varstemp,bvar=bvartemp,nonbvar=nonbvartemp,
simplexmatrix=simplexmatrixtemp,objcoef=objcoeftemp,
sol=soltemp,objsol=objsoltemp,bmatrix,invar,outvar,
pivotcol1,pivotcol2,rowtrans,rowcomp,simplexobj,bvarlength,
outtemp,mytab,z,Solution},
bvarlength=Length@bvar;
simplexobj=-objcoef;

(*選擇進入變數*)
invar=Select[
Transpose@{simplexobj,vars},#[[1]]<0||Coefficient[#[[1]],M]<0&][[1,2]]//Quiet;
pivotcol1=Position[vars,invar][[1,1]];

(*選擇離去變數*)
outvar=SortBy[
Select[Transpose@{bvar,sol/simplexmatrix[[All,pivotcol1]]},
#[[2]]>0&],Last][[1,1]]//Quiet;
pivotcol2=Position[bvar,outvar][[1,1]];

(*對基本變數進行列運算*)
rowtrans=(-simplexmatrix[[#,pivotcol2]]/
simplexmatrix[[pivotcol1,pivotcol2]]&/@
Complement[Range[bvarlength],{pivotcol1}]);

rowcomp=Complement[Range[bvarlength],{pivotcol1}];
objsol=objsol-simplexobj[[pivotcol1]]/simplexmatrix[[pivotcol1,pivotcol2]]*
sol[[pivotcol1]];

simplexobj=simplexobj-simplexobj[[pivotcol1]]/simplexmatrix[[pivotcol1,pivotcol2]]*
simplexmatrix[[pivotcol1]];

sol=SortBy[
Insert[{rowcomp[[#]],sol[[pivotcol1]]*rowtrans[[#]]+sol[[rowcomp[[#]]]]}&/@
Range[Length@rowcomp],{pivotcol2,
sol[[pivotcol1]]/simplexmatrix[[pivotcol1,pivotcol2]]},1],First][[All,2]];

simplexmatrix=
SortBy[Insert[{rowcomp[[#]],
simplexmatrix[[pivotcol1]]*rowtrans[[#]]+
simplexmatrix[[rowcomp[[#]]]]}&/@
Range[Length@rowcomp],{pivotcol2,simplexmatrix[[pivotcol1]]/
simplexmatrix[[pivotcol1,pivotcol2]]},1],First][[All,2]];
bvar=bvar/.outvar->invar;nonbvar=nonbvar/.invar->outvar;

(*將變數名稱加入SimplexMethod表格輸出*)
outtemp=
Insert[Insert[simplexmatrix[[#]],sol[[#]],-1]&/@
Range[Length@bvar],Flatten[{simplexobj,objsol}],1];

(*SimplexMethod表格輸出*)
mytab=Grid[
Prepend[Partition[
Flatten[Transpose@{Flatten[{z,bvar}],outtemp}],
2+Length@vars],Flatten[{,vars,Solution}]],
Spacings->{2,2},Frame->All,
Background->{None,{Lighter[Yellow,.9],
Lighter[Red,.9],{Lighter[Blend[{Blue,Green}],.8],White}}},
FrameStyle->Directive[Thickness[1]]];
Simplify@{vars,bvar,Complement[vars,bvar],
simplexmatrix,-simplexobj,sol,objsol,mytab}
]

(* Example 1*)

bvar={s1,s2,s3,s4};
nonbvar={x1,x2};
vars=Flatten[{nonbvar,bvar}];
objcoef={5,4,0,0,0,0};
nonbmatrix={{6,4},{1,2},{-1,1},{0,1}};
bmatrix=IdentityMatrix[Length@bvar];
sol={24,6,1,2};
objsol=0;

simplexmatrix=MapThread[Flatten[Append[#1,#2],1]&,{nonbmatrix,bmatrix}];
objsol=objsol+sol.(1*objcoef[[#]]&/@Range[Length@nonbvar+1,Length@vars]);
objcoef=objcoef+Plus@@(simplexmatrix*(-1*objcoef[[#]]&/@
Range[Length@nonbvar+1,Length@vars]));

step1=simplex[vars,bvar,nonbvar,simplexmatrix,objcoef,sol,objsol]
step2 = simplex @@ step1[[1 ;; -2]]

(* Example 2*)

vars={x1,x2,x3,R1,R2,x4};
bvar={R1,R2,x4};
nonbvar={x1,x2,x3};
objcoef={-4,-1,0,-M,-M,0};
nonbmatrix={{3,1,0},{4,3,-1},{1,2,0}};
sol={3,6,4};
objsol=0;
bmatrix=IdentityMatrix[Length@bvar];
simplexmatrix=Flatten[Insert[nonbmatrix[[#]],bmatrix[[#]],-1]]&/@Range[Length@bvar];
objsol=objsol+sol.(1*objcoef[[#]]&/@Range[Length@nonbvar+1,Length@vars]);
objcoef=objcoef+Plus@@(simplexmatrix*(-1*objcoef[[#]]&/@
Range[Length@nonbvar+1,Length@vars]));


step1=simplex[vars,bvar,nonbvar,simplexmatrix,objcoef,sol,objsol]
step2=simplex@@step1[[1;;-2]]
step3=simplex@@step2[[1;;-2]]

Tags: Mathematica 教學

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