Simplex Method updated to V20.13 (cas mode)
Update note:
Updated: 2016-12-17
sim20_13 : If there are infinitely many solutions, program will try to choose the smallest value.
Updated: 2016-12-9
sim20_11 : user definable variables name when importing data from Spreadsheet.
Use of argument:
sim20(0) : new input, no previous data.
sim20(1) : Load previous data.
sim20(2),sim20(3).... : Load and save data to "2","3".....
sim20(99) : Import data from spreadsheet.
Features:
1)Solve linear programming problem by simplex method.
2)Provide pure integer solution.
3)Provide Mixed integer, continuous and binary solution. (Read POST#6)
4)Save and reload data (Read POST#6)
5)Set target Z value (Read POST#6)
6)Import data from spreadsheet. (Read POST#10)
7)User definable variables name. (Read POST#11)
Hello,
Simplex program is updated to v20.
First of all, try this example: Type sim20(0) in CAS mode.
Max X1 + 0.64X2
s.t. 50X1 + 31X2 <= 250
3X1 - 2X2 >= -4
X1 and X2 >=0 integers
Generally speaking, B&B should be more reliable. But above B&B solution is not optimal.
Why? It is because the searching level is not deep enough.
Pay attention to Reach Level: 5 ->6, it means the level limit is 5 but the program need to go to level 6 for more searching.
So, try to increase Search Level to 10.
Now, the solution is better and the Reach Level is 7.
It means the program is stopped at level 7 and no need to go deeper.
So, the solution should be optimal.
Example to use this program to find shortest path.
Set up the equations:
The solution is {1,0,1,0,0,1} Z=6, it means X12->X23->X34 is the shortest path.
Note that searching integer algorithm is time consuming, some problem may take one hour to solve. So, be patient.
Troubleshooting:
A) If you find sim20(0) is not working properly.
solution:
1) Delete sim17/18/19 if you have it before. Some function names are same.
2) Check cas variables by pressing [SHIFT] [Mem], then press "CAS Vars" , delete all rubbish function and sim17/18/19 related functions. Then press "Esc".
3) Press [SHIFT][Program],edit sim20_8, press softkey "Check", then "Esc" and "Esc".
4) Type sim20(0) in CAS mode to see whether it is working well. If not, repeat step 2 and 3.
B) Running non recursive evaluator ERROR
solution:
1) Press [SHIFT][CAS setting], go to page 2, set "Recursive Function" to 100.
C) Get Syntax error when pressing "CHECK" in program editor.
solution:
Put Spreadsheet App as the current App of calculator. Calculator title bar will show Spreadsheet App name. Then edit and CHECK program again.
sim20(0) codes:
Code:
EXPORT IMPORTDATA()
BEGIN
LOCAL C,V,P,S1,S2,S3,S4,S5;
LOCAL O,I,M,T,A,B,LEVEL;
LOCAL TAB,E,K,S,Z,FIX;
LOCAL METHOD,METHOD1,METHOD2,METHOD3,SUBLEVEL,INTGR;
LOCAL G,H,METHOD4,IMASK,BMASK;
LOCAL S6,S7,S8,S9;
LOCAL TT,SS,EE,OO,CC,MIXTAB,TEMP;
LOCAL ENDMARK1,ENDMARK2,VVV,CCC,NAMELIST;
LEVEL:=5;SUBLEVEL:=2;METHOD:=4;METHOD1:=0;METHOD2:=0;METHOD3:=0;INTGR:=0;METHOD4:=0;
I:=0;
REPEAT
I:=I+1;
ENDMARK1:=(STRING(Cell(2,I))=="""end""" OR STRING(Cell(2,I))=="""END""" OR STRING(Cell(2,I))=="""MAX""" OR STRING(Cell(2,I))=="""max""" OR STRING(Cell(2,I))=="""MIN""" OR STRING(Cell(2,I))=="""min""");
UNTIL ENDMARK1;
V:=I-1;
I:=2;
REPEAT
I:=I+1;
ENDMARK2:=(STRING(Cell(I,1))=="""end""" OR STRING(Cell(I,1))=="""END""");
UNTIL ENDMARK2 ;
C:=I-3;
VVV:=V;CCC:=C;
INPUT({{VVV,[0],{55,20,0}},{CCC,[0],{55,20,1}},{INTGR,1,{55,20,2}},{LEVEL,[0],{55,20,3}},{METHOD,{"BRANCH & BOUND","CUTTING PLANE","Mixed B&C","ALL"},{55,42,4}},{METHOD4,1,{92,5,5}}},"SIMPLEX",{"NUMBER OF VARIABLES:","NUMBER OF CONSTRINTS:","INTEGER SOLUTION:","INTEGER SEARCH LEVEL:","METHOD:","MIXED INTEGER-CONTINUOUS-BINARY[0,1]:"});
IF CCC==0 OR VVV==0 THEN RETURN(0);END;
IF METHOD==1 AND INTGR==1 THEN METHOD1:=1;END;
IF METHOD==2 AND INTGR==1 THEN METHOD2:=1;END;
IF METHOD==3 AND INTGR==1 THEN METHOD3:=1;END;
IF METHOD==4 AND INTGR==1 THEN METHOD1:=1;METHOD2:=1;METHOD3:=1;END;
////////////// prepare variable //////////////////
O:=MAKELIST(0,X,1,V);
T:=MAKELIST(O,X,1,C);
E:=MAKELIST(0,X,1,C);
S:=MAKELIST(2,X,1,C);
IMASK:=MAKELIST(0,X,1,V);
BMASK:=MAKELIST(0,X,1,V);
////////////// prepare variable //////////////////
////////////// Read Spreadsheet /////////////
H:=1;G:=0;NAMELIST:={};
FOR I FROM 1 TO V DO
IF TYPE(Cell(1,I))==2 THEN
NAMELIST(I):=Cell(1,I);
ELSE
NAMELIST(I):="X"+I;
END;
O(I):=Cell(2,I);
END;//FOR I , Objective, variable name
I:=V+1;
IF TYPE(Cell(1,I))==2 THEN
NAMELIST(I):=Cell(1,I);
ELSE
NAMELIST(I):="Z";
END;
IF STRING(Cell(2,I))=="""MAX""" OR STRING(Cell(2,I))=="""max""" THEN P:=1;END;
IF STRING(Cell(2,I))=="""MIN""" OR STRING(Cell(2,I))=="""min""" THEN P:=2;END;
I:=V+2;
IF STRING(Cell(2,I))==""">=""" OR STRING(Cell(2,I))=="""≥""" THEN H:=2;END;
IF STRING(Cell(2,I))=="""<=""" OR STRING(Cell(2,I))=="""≤""" THEN H:=3;END;
IF STRING(Cell(2,I))=="""=""" THEN H:=4;END;
I:=V+3;
G:=Cell(2,I);
FOR K FROM 1 TO C DO
FOR I FROM 1 TO V DO
T(K,I):=Cell(K+2,I);
END;//FOR I, ==V
I:=V+1;
IF STRING(Cell(K+2,I))==""">=""" OR STRING(Cell(K+2,I))=="""≥""" THEN S(K):=1;END;
IF STRING(Cell(K+2,I))=="""<=""" OR STRING(Cell(K+2,I))=="""≤""" THEN S(K):=2;END;
IF STRING(Cell(K+2,I))=="""=""" THEN S(K):=3;END;
I:=V+2;
E(K):=Cell(K+2,I);
END;//FOR K, ==C
////////////// Read Spreadsheet /////////////
///// MIXED INT CONT BIN /////////////
IF METHOD4==1 THEN
S6:="";
S7:="";
S8:="";
M:=-1;
FOR I FROM 1 TO V DO
IF ((I MOD 3)==1) THEN L:=40;END;
IF ((I MOD 3)==2) THEN L:=40+27;END;
IF ((I MOD 3)==0) THEN L:=40+27*2;END;
IF ((I MOD 3)==1) THEN M:=M+1;END;
S6:=S6+"{IMASK("+I+"),1,{"+L+",5,"+M+"}},";
S7:=S7+""""+LEFT(NAMELIST(I),7)+""",";
END;//FOR I
FOR I FROM 1 TO V DO
IF ((I MOD 3)==1) THEN L:=40;END;
IF ((I MOD 3)==2) THEN L:=40+27;END;
IF ((I MOD 3)==0) THEN L:=40+27*2;END;
IF ((I MOD 3)==1) THEN M:=M+1;END;
S6:=S6+"{BMASK("+I+"),1,{"+L+",5,"+M+"}},";
S8:=S8+""""+LEFT(NAMELIST(I),7)+""",";
END;//FOR I
S6:=LEFT(S6,DIM(S6)-1);
S7:=RIGHT(S7,DIM(S7)-1);
S8:=LEFT(S8,DIM(S8)-1);
S8:=RIGHT(S8,DIM(S8)-1);
S9:="INPUT({"+S6+"},""CHOOSE INT. BIN. VARIABLES"",{""INTEGER: "+S7+"""BINARY: "+S8+"})";
EXPR(S9);
END;//METHOD4==1
///// MIXED INT CONT BIN /////////////
//PRINT EQUATIONS
PRINT();
S1:="";
IF P==1 THEN S1:="MAXIMIZE: "+NAMELIST(V+1)+"= ";
ELSE S1:="MINIMIZE: "+NAMELIST(V+1)+"= ";END;
FOR I FROM 1 TO V DO
IF SIGN(O(I))==−1 THEN S1:=S1+" - ";END;
IF SIGN(O(I))≥0 AND I≠1 THEN S1:=S1+" + ";END;
S1:=S1+STRING(ABS(O(I)),7,-1)+"_"+NAMELIST(I);
END;//FOR I
/////// Target Z string //////////
IF H==2 THEN S1:=S1+" ≥ ";S1:=S1+G;END;
IF H==3 THEN S1:=S1+" ≤ ";S1:=S1+G;END;
IF H==4 THEN S1:=S1+" = ";S1:=S1+G;END;
/////// Target Z string //////////
PRINT(S1);
PRINT("SUBJECT TO:");
FOR K FROM 1 TO C DO
S1:="";
FOR I FROM 1 TO V DO
IF SIGN(T(K,I))==−1 THEN S1:=S1+" - ";END;
IF SIGN(T(K,I))≥0 AND I≠1 THEN S1:=S1+" + ";END;
S1:=S1+STRING(ABS(T(K,I)),7,-1)+"_"+NAMELIST(I);
END;//FOR I
IF S(K)==1 THEN S1:=S1+" ≥ ";END;
IF S(K)==2 THEN S1:=S1+" ≤ ";END;
IF S(K)==3 THEN S1:=S1+" = ";END;
IF SIGN(E(K))==−1 THEN S1:=S1+"-";END;
S1:=S1+STRING(ABS(E(K)),7,-1);
PRINT(S1);
END;//FOR K
PRINT("---------------------");
///////////////////////////////////
////////// Target Z equation /////////
IF H>1 THEN
T:=append(T,O);
E:=append(E,G);
S:=append(S,H-1);
C:=C+1;
END;//IF H>1
////////// Target Z equation /////////
///////////////////////////////////
//Equal
FIX:=SIZE(S);
FOR I FROM 1 TO FIX DO
IF S(I)==3 THEN
S(I):=2;
C:=C+1;
T(C):=T(I);
E(C):=E(I);
S(C):=1;
END;//IF S(I)
END;//FOR I
///////////////////////////////////
///// SAVE T,E,S,O,C /////////////////////
TT:=T;
EE:=E;
SS:=S;
OO:=O;
CC:=C;
///// SAVE T,E,S,O,C /////////////////////
/////////// Construct TAB //////////////////
TAB:=append(T,(−1)^P*O);
TAB:=list2mat(TAB);
A:=MAKEMAT(0,1,C+1);
FOR I FROM 1 TO C DO
A:=A*0;
A(1,I):=(−1)^S(I);
ADDCOL(TAB,A(1),V+I);
END;
Z:=MAKEMAT(0,1,C+1);
Z(1,C+1):=1;
ADDCOL(TAB,Z(1),C+V+1);
E:=append(E,0);
B:=list2mat(E,C+1);
ADDCOL(TAB,B(1),C+V+2);
/////////// Construct TAB //////////////////
/////////// Construct Mix TAB //////////////
IF METHOD4==1 THEN
//// constraint add X? <=1 BIN/////////
FOR I FROM 1 TO V DO
IF BMASK(I)==1 THEN
TEMP:=MAKELIST(0,X,1,V);
TEMP(I):=1;
TT:=append(TT,TEMP);
SS:=append(SS,2);
EE:=append(EE,1);
CC:=CC+1;
END;//IF BMASK==1
END;//FOR I
//// constraint add X? <=1 BIN/////////
MIXTAB:=append(TT,(−1)^P*OO);
MIXTAB:=list2mat(MIXTAB);
A:=MAKEMAT(0,1,CC+1);
FOR I FROM 1 TO CC DO
A:=A*0;
A(1,I):=(−1)^SS(I);
ADDCOL(MIXTAB,A(1),V+I);
END;
Z:=MAKEMAT(0,1,CC+1);
Z(1,CC+1):=1;
ADDCOL(MIXTAB,Z(1),CC+V+1);
EE:=append(EE,0);
B:=list2mat(EE,CC+1);
ADDCOL(MIXTAB,B(1),CC+V+2);
END;//IF method4==1
/////////// Construct Mix TAB //////////////
IMASK:=BITOR(IMASK,BMASK);
RETURN({TAB,C,V,P,LEVEL,SUBLEVEL,METHOD1,METHOD2,METHOD3,METHOD4,IMASK,MIXTAB,CC,NAMELIST});
END;
EXPORT INPUTDATA(REEDIT,REEDITBAG)
BEGIN
LOCAL C,V,P,S1,S2,S3,S4,S5;
LOCAL O,I,M,T,A,B,LEVEL;
LOCAL TAB,E,K,S,Z,FIX;
LOCAL METHOD,METHOD1,METHOD2,METHOD3,SUBLEVEL,INTGR;
LOCAL G,H,METHOD4,IMASK,BMASK;
LOCAL S6,S7,S8,S9,NAMELIST;
LOCAL TT,SS,EE,OO,CC,MIXTAB,TEMP;
LEVEL:=5;SUBLEVEL:=2;METHOD:=4;METHOD1:=0;METHOD2:=0;METHOD3:=0;INTGR:=0;METHOD4:=0;
IF REEDIT>=1 AND SIZE(REEDITBAG)==15 THEN
LEVEL:=REEDITBAG(15);
METHOD4:=REEDITBAG(14);
METHOD:=REEDITBAG(13);
INTGR:=REEDITBAG(12);
C:=REEDITBAG(5);
V:=REEDITBAG(6);
END;
INPUT({{V,[0],{55,20,0}},{C,[0],{55,20,1}},{INTGR,1,{55,20,2}},{LEVEL,[0],{55,20,3}},{METHOD,{"BRANCH & BOUND","CUTTING PLANE","Mixed B&C","ALL"},{55,42,4}},{METHOD4,1,{92,5,5}}},"SIMPLEX",{"NUMBER OF VARIABLES:","NUMBER OF CONSTRINTS:","INTEGER SOLUTION:","INTEGER SEARCH LEVEL:","METHOD:","MIXED INTEGER-CONTINUOUS-BINARY[0,1]:"});
//INPUT({{V,[0],{55,20,0}},{C,[0],{55,20,1}},{LEVEL,[0],{55,20,2}},{METHOD1,1,{55,20,3}},{METHOD2,1,{55,20,4}},{SUBLEVEL,[0],{55,20,5}},{METHOD3,1,{55,20,6}}},"SIMPLEX",{"NUMBER OF VARIABLES:","NUMBER OF CONSTRINTS:","INTEGER SEARCH LEVEL:","CUTTING PLANE:","BRANCH & BOUND:","SUB CUT LEVEL:","MIXED BRANCH & CUT:"});
//INPUT({{V,[0],{55,20,0}},{C,[0],{55,20,1}},{LEVEL,[0],{55,20,2}},{METHOD2,1,{55,20,3}},{METHOD1,1,{55,20,4}},{SUBLEVEL,[0],{55,20,5}},{METHOD3,1,{55,20,6}}},"SIMPLEX",{"NUMBER OF VARIABLES:","NUMBER OF CONSTRINTS:","INTEGER SEARCH LEVEL:","BRANCH & BOUND:","CUTTING PLANE:","SUB CUT LEVEL:","MIXED BRANCH & CUT:"});
IF C==0 OR V==0 THEN RETURN(0);END;
IF METHOD==1 AND INTGR==1 THEN METHOD1:=1;END;
IF METHOD==2 AND INTGR==1 THEN METHOD2:=1;END;
IF METHOD==3 AND INTGR==1 THEN METHOD3:=1;END;
IF METHOD==4 AND INTGR==1 THEN METHOD1:=1;METHOD2:=1;METHOD3:=1;END;
////////////// prepare variable //////////////////
O:=MAKELIST(0,X,1,V);
T:=MAKELIST(O,X,1,C);
E:=MAKELIST(0,X,1,C);
S:=MAKELIST(2,X,1,C);
IMASK:=MAKELIST(0,X,1,V);
BMASK:=MAKELIST(0,X,1,V);
IF REEDIT>=1 AND SIZE(REEDITBAG)==15 THEN
IF V>REEDITBAG(6) THEN
FOR I FROM REEDITBAG(6)+1 TO V DO
REEDITBAG(4):=append(REEDITBAG(4),0);
REEDITBAG(10):=append(REEDITBAG(10),0);
REEDITBAG(11):=append(REEDITBAG(11),0);
FOR K FROM 1 TO REEDITBAG(5) DO
REEDITBAG(1,K):=append(REEDITBAG(1,K),0);
END;//FOR K
END;//FOR I
END;//IF larger V
IF V<REEDITBAG(6) THEN
REEDITBAG(4):=SUB(REEDITBAG(4),1,V);
REEDITBAG(10):=SUB(REEDITBAG(10),1,V);
REEDITBAG(11):=SUB(REEDITBAG(11),1,V);
FOR K FROM 1 TO REEDITBAG(5) DO
REEDITBAG(1,K):=SUB(REEDITBAG(1,K),1,V);
END;//FOR K
END;//IF smaller V
IF C>REEDITBAG(5) THEN
FOR I FROM REEDITBAG(5)+1 TO C DO
REEDITBAG(2):=append(REEDITBAG(2),2);
REEDITBAG(3):=append(REEDITBAG(3),0);
REEDITBAG(1):=append(REEDITBAG(1),O);
END;//FOR I
END;//IF larger C
IF C<REEDITBAG(5) THEN
REEDITBAG(2):=SUB(REEDITBAG(2),1,C);
REEDITBAG(3):=SUB(REEDITBAG(3),1,C);
REEDITBAG(1):=SUB(REEDITBAG(1),1,C);
END;//IF larger C
////// RELOAD DATA ///////
T:=REEDITBAG(1);
S:=REEDITBAG(2);
E:=REEDITBAG(3);
O:=REEDITBAG(4);
P:=REEDITBAG(7);
G:=REEDITBAG(8);
H:=REEDITBAG(9);
IMASK:=REEDITBAG(10);
BMASK:=REEDITBAG(11);
////// RELOAD DATA ///////
END;//IF reedit==1
////////////// prepare variable //////////////////
//PRINT("XXXXXXXX");
//PRINT(REEDITBAG(1));
//PRINT(REEDITBAG(4));
//PRINT(REEDITBAG(10));
//PRINT(REEDITBAG(11));
//PRINT(REEDITBAG(2));
//PRINT(REEDITBAG(3));
//PRINT("XXXXXXXX");
///// MIXED INT CONT BIN /////////////
IF METHOD4==1 THEN
S6:="";
S7:="";
S8:="";
M:=-1;
FOR I FROM 1 TO V DO
IF ((I MOD 4)==1) THEN L:=25;END;
IF ((I MOD 4)==2) THEN L:=25+15;END;
IF ((I MOD 4)==3) THEN L:=25+15*2;END;
IF ((I MOD 4)==0) THEN L:=25+15*3;END;
IF ((I MOD 4)==1) THEN M:=M+1;END;
S6:=S6+"{IMASK("+I+"),1,{"+L+",5,"+M+"}},";
S7:=S7+"""X"+I+""",";
END;//FOR I
FOR I FROM 1 TO V DO
IF ((I MOD 4)==1) THEN L:=25;END;
IF ((I MOD 4)==2) THEN L:=25+15;END;
IF ((I MOD 4)==3) THEN L:=25+15*2;END;
IF ((I MOD 4)==0) THEN L:=25+15*3;END;
IF ((I MOD 4)==1) THEN M:=M+1;END;
S6:=S6+"{BMASK("+I+"),1,{"+L+",5,"+M+"}},";
S8:=S8+"""X"+I+""",";
END;//FOR I
S6:=LEFT(S6,DIM(S6)-1);
S7:=RIGHT(S7,DIM(S7)-1);
S8:=LEFT(S8,DIM(S8)-1);
S8:=RIGHT(S8,DIM(S8)-1);
S9:="INPUT({"+S6+"},""CHOOSE INT. BIN. VARIABLES"",{""INTEGER: "+S7+"""BINARY: "+S8+"})";
EXPR(S9);
END;//METHOD4==1
///// MIXED INT CONT BIN /////////////
S1:="{P,{""Max"",""Min""},{1,17,0}}";
S3:="";
S4:="""Z"",";
M:=−1;
FOR I FROM 1 TO V DO
IF ((I MOD 3)==1) THEN L:=27;END;
IF ((I MOD 3)==2) THEN L:=27+27;END;
IF ((I MOD 3)==0) THEN L:=27+27*2;END;
IF (I MOD 3)==1 THEN M:=M+1;END;
S3:=S3+"{O("+I+"),[0],{"+L+",18,"+M+"}},";
S4:=S4+"""+X"+I+""",";
END;
//////// Target Z /////////////
I:=V+1;
IF ((I MOD 3)==1) THEN L:=27;END;
IF ((I MOD 3)==2) THEN L:=27+27;END;
IF ((I MOD 3)==0) THEN L:=27+27*2;END;
IF (I MOD 3)==1 THEN M:=M+1;END;
S3:=S3+"{H,{""Z"",""≥"",""≤"",""=""},{"+(L-8)+",9,"+(M)+"}},{G,[0],{"+(L+2)+",16,"+M+"}},";
S4:=S4+""""","""",";
//////// Target Z /////////////
S3:=LEFT(S3,DIM(S3)-1);
S5:="";
M:=M+1;
FOR K FROM 1 TO C DO
FOR I FROM 1 TO V DO
IF ((I MOD 3)==1) THEN L:=6;END;
IF ((I MOD 3)==2) THEN L:=6+24;END;
IF ((I MOD 3)==0) THEN L:=6+24*2;END;
IF ((I MOD 3)==1) THEN M:=M+1;END;
S5:=S5+"{T("+K+","+I+"),[0],{"+L+",17,"+(M)+"}},";
S4:=S4+"""X"+I+""",";
END;
S5:=S5+"{S("+K+"),{""≥"",""≤"",""=""},{72,9,"+(M)+"}},{E("+K+"),[0],{82,17,"+(M)+"}},";
S4:=S4+""""","""",";
END;
S5:=LEFT(S5,DIM(S5)-1);
S4:=LEFT(S4,DIM(S4)-1);
S2:="INPUT({"+S1+","+S3+","+S5+"},""SIMPLEX"",{"+S4+"})";
EXPR(S2);
//////////////////////////////
/////////////// SAVE VARIABLES FOR REEDIT ///////////////
REEDITBAG:=MAKELIST(0,X,1,15);
REEDITBAG(1):=T;
REEDITBAG(2):=S;
REEDITBAG(3):=E;
REEDITBAG(4):=O;
REEDITBAG(5):=C;
REEDITBAG(6):=V;
REEDITBAG(7):=P;
REEDITBAG(8):=G;
REEDITBAG(9):=H;
REEDITBAG(10):=IMASK;
REEDITBAG(11):=BMASK;
REEDITBAG(12):=INTGR;
REEDITBAG(13):=METHOD;
REEDITBAG(14):=METHOD4;
REEDITBAG(15):=LEVEL;
/////////////// SAVE VARIABLES FOR REEDIT ///////////////
//PRINT EQUATIONS
PRINT();
S1:="";
IF P==1 THEN S1:="MAXIMIZE: Z= ";
ELSE S1:="MINIMIZE: Z= ";END;
FOR I FROM 1 TO V DO
IF SIGN(O(I))==−1 THEN S1:=S1+" - ";END;
IF SIGN(O(I))≥0 AND I≠1 THEN S1:=S1+" + ";END;
S1:=S1+STRING(ABS(O(I)),7,-1)+"_X"+I;
END;//FOR I
/////// Target Z string //////////
IF H==2 THEN S1:=S1+" ≥ ";S1:=S1+G;END;
IF H==3 THEN S1:=S1+" ≤ ";S1:=S1+G;END;
IF H==4 THEN S1:=S1+" = ";S1:=S1+G;END;
/////// Target Z string //////////
PRINT(S1);
PRINT("SUBJECT TO:");
FOR K FROM 1 TO C DO
S1:="";
FOR I FROM 1 TO V DO
IF SIGN(T(K,I))==−1 THEN S1:=S1+" - ";END;
IF SIGN(T(K,I))≥0 AND I≠1 THEN S1:=S1+" + ";END;
S1:=S1+STRING(ABS(T(K,I)),7,-1)+"_X"+I;
END;//FOR I
IF S(K)==1 THEN S1:=S1+" ≥ ";END;
IF S(K)==2 THEN S1:=S1+" ≤ ";END;
IF S(K)==3 THEN S1:=S1+" = ";END;
IF SIGN(E(K))==−1 THEN S1:=S1+"-";END;
S1:=S1+STRING(ABS(E(K)),7,-1);
PRINT(S1);
END;//FOR K
PRINT("---------------------");
///////////////////////////////////
////////// Target Z equation /////////
IF H>1 THEN
T:=append(T,O);
E:=append(E,G);
S:=append(S,H-1);
C:=C+1;
END;//IF H>1
////////// Target Z equation /////////
///////////////////////////////////
//Equal
FIX:=SIZE(S);
FOR I FROM 1 TO FIX DO
IF S(I)==3 THEN
S(I):=2;
C:=C+1;
T(C):=T(I);
E(C):=E(I);
S(C):=1;
END;//IF S(I)
END;//FOR I
///////////////////////////////////
///// SAVE T,E,S,O,C /////////////////////
TT:=T;
EE:=E;
SS:=S;
OO:=O;
CC:=C;
///// SAVE T,E,S,O,C /////////////////////
/////////// Construct TAB //////////////////
TAB:=append(T,(−1)^P*O);
TAB:=list2mat(TAB);
A:=MAKEMAT(0,1,C+1);
FOR I FROM 1 TO C DO
A:=A*0;
A(1,I):=(−1)^S(I);
ADDCOL(TAB,A(1),V+I);
END;
Z:=MAKEMAT(0,1,C+1);
Z(1,C+1):=1;
ADDCOL(TAB,Z(1),C+V+1);
E:=append(E,0);
B:=list2mat(E,C+1);
ADDCOL(TAB,B(1),C+V+2);
/////////// Construct TAB //////////////////
/////////// Construct Mix TAB //////////////
IF METHOD4==1 THEN
//// constraint add X? <=1 BIN/////////
FOR I FROM 1 TO V DO
IF BMASK(I)==1 THEN
TEMP:=MAKELIST(0,X,1,V);
TEMP(I):=1;
TT:=append(TT,TEMP);
SS:=append(SS,2);
EE:=append(EE,1);
CC:=CC+1;
END;//IF BMASK==1
END;//FOR I
//// constraint add X? <=1 BIN/////////
MIXTAB:=append(TT,(−1)^P*OO);
MIXTAB:=list2mat(MIXTAB);
A:=MAKEMAT(0,1,CC+1);
FOR I FROM 1 TO CC DO
A:=A*0;
A(1,I):=(−1)^SS(I);
ADDCOL(MIXTAB,A(1),V+I);
END;
Z:=MAKEMAT(0,1,CC+1);
Z(1,CC+1):=1;
ADDCOL(MIXTAB,Z(1),CC+V+1);
EE:=append(EE,0);
B:=list2mat(EE,CC+1);
ADDCOL(MIXTAB,B(1),CC+V+2);
END;//IF method4==1
/////////// Construct Mix TAB //////////////
NAMELIST:={};
FOR I FROM 1 TO V DO
NAMELIST(I):="X"+I;
END;//FOR I variable Name
I:=V+1;
NAMELIST(I):="Z";
IMASK:=BITOR(IMASK,BMASK);
RETURN({TAB,C,V,P,LEVEL,SUBLEVEL,METHOD1,METHOD2,METHOD3,METHOD4,IMASK,MIXTAB,CC,REEDITBAG,NAMELIST});
END;
#cas
caspivot(m,r,c):=BEGIN
local t,k,si;
t:=m;
si:=SIZE(m);
t[r]:=m[r]*(1/m[r,c]);
FOR k FROM 1 to si(1) DO
IF k≠r THEN
t[k]:=m[k]-m[k,c]*t[r];
END;//IF K≠RO
END;//FOR K
return(t);
END;
#end
#cas
casphaseone(tab,c,v,stp1,intgr)
BEGIN
//PHASE ONE
local la,lb,alln,onezc;
local ill,lll,nse,nzn,nzr,k,i;
local rrr,mcl,w,maxi;
local ii,kk,ia,ib,nzb,nzbr;
onezc:=0;
REPEAT
FOR la FROM 1 TO c+1 DO
alln:=1;ill:=0;
FOR lb FROM 1 TO c+v+2 DO
IF tab[la,lb]>0 THEN alln:=0;END;
END;//FOR-LB
IF alln THEN tab[la]:=−1*tab[la];END;
END;//FOR-LA
nse:=0;lll:=0;
FOR i FROM 1 TO c+v DO
nzn:=0;nzr:=0;
FOR k FROM 1 TO c+1 DO
IF tab[k,i]≠0 THEN nzn:=nzn+1;nzr:=k;END;
END;//FOR-K
IF nzn≠0 THEN
IF (nzn==1) AND (nzr!=c+1) AND SIGN(tab[nzr,c+v+2])<0 AND SIGN(tab[nzr,i])<0 THEN tab[nzr]:=-1*tab[nzr];END;
IF (nzn==1) AND (nzr!=c+1) AND tab[nzr,c+v+2]==0 AND SIGN(tab[nzr,i])<0 THEN tab[nzr]:=-1*tab[nzr];END;
IF (nzn==1) AND (nzr!=c+1) AND (SIGN(tab[nzr,c+v+2])*SIGN(tab[nzr,i])<0 OR (tab[nzr,c+v+2]==0 AND SIGN(tab[nzr,i])<0)) AND (onezc==0) THEN
alln:=1;
FOR lb FROM 1 TO c+v DO
IF SIGN(tab[nzr,c+v+2])*SIGN(tab[nzr,lb])>0 OR (tab[nzr,c+v+2]==0 AND SIGN(tab[nzr,lb])>0) THEN alln:=0;END;
END;//FOR-LB
IF alln==0 THEN
nse:=1;rrr:=nzr;lll:=i;
/////////////////// ADD same row check
FOR ia FROM 1 to c+v DO
IF ia!=lll THEN
nzb:=0;nzbr:=0;
FOR ib FROM 1 to c+1 DO
IF tab[ib,ia]!=0 THEN nzb:=nzb+1;nzbr:=ib;END;
END;//FOR ib (row)
END;//IF ia!=lll
IF nzb==1 AND nzbr==rrr AND tab[nzbr,ia]>0 THEN nse:=0;BREAK;END;
END;//FOR ia (COL)
////////////////// ADD same row check
ELSE
ill:=1;
END;//IF-ALLN
IF nse==1 THEN BREAK;END;
END;//IF-(NZN==1)
END;//IF-NZN≠0
//IF (NZN==1) AND (TAB(NZR,C+V+2)==0) AND (ONEZC==0)THEN NSE:=1;RRR:=NZR;LLL:=I;ONEZC:=1;END;
END;//FOR-I
IF nse==1 THEN
maxi:=0;mcl:=0;
FOR i FROM 1 TO c+v DO
IF (SIGN(tab[rrr,c+v+2])*tab[rrr,i]>maxi OR (tab[rrr,c+v+2]==0 AND tab[rrr,i]>maxi)) AND i≠lll THEN mcl:=i;maxi:=MAX(SIGN(tab[rrr,c+v+2])*tab[rrr,i],tab[rrr,i]);END;
END;
//TAB(RRR):=TAB(RRR)/TAB(RRR,MCL);
tab:=caspivot(tab,rrr,mcl);
////////////////////////
///////////////////////
IF intgr==0 THEN PRINT("PHASE ONE: ROW="+rrr+" COL="+mcl);END;
END;
stp1:=stp1+1;
UNTIL nse==0 OR stp1>=100;
return({tab,c,v,stp1});
END;
#end
#cas
casphasetwo(tab,c,v,stp2,intgr):=BEGIN
local mini,neg,la,lb,alln,ook;
local i,k,j,r,q1,q2,pack,stp1;
local ookcol,ii,kk,nonzero,nonzerorow,nonzerorowlist;
stp2:=0;stp1:=0;
REPEAT
mini:=inf;neg:=0;
IF tab[c+1,c+v+1]<0 THEN tab[c+1]:=-1*tab[c+1];END;
FOR la FROM 1 TO c+1 DO
alln:=1;
FOR lb FROM 1 TO c+v+2 DO
IF tab[la,lb]>0 THEN alln:=0;END;
END;//FOR-LB
IF alln THEN tab[la]:=−1*tab[la];END;
END;//FOR-LA
ook:=2;ookcol:=0;
FOR i FROM 1 TO v+c DO
IF (tab[c+1,i]<0) AND (tab[c+1,i]<mini) THEN
ook:=0;ookcol:=i;
FOR k FROM 1 TO c DO
IF tab[k,i]>0 THEN
ook:=1;
END;//IF2
END;//FOR-K
IF ook==1 THEN
mini:=tab[c+1,i];j:=i;neg:=1;
END;//IF3
END;//IF (tab[c+1,i]<0) AND (tab[c+1,i]<mini)
END;//FOR-I
/////////// TESTING ///////////////
IF ook==0 THEN
nonzerorowlist:=MAKELIST(X,X,1,c);
FOR ii FROM 1 TO c+v DO
nonzero:=0;nonzerorow:=0;
FOR kk FROM 1 TO c+1 DO
IF tab[kk,ii]!=0 THEN nonzero:=nonzero+1;nonzerorow:=kk;END;
END;//FOR kk
IF nonzero==1 THEN nonzerorowlist:=remove(nonzerorow,nonzerorowlist);END;
END;//FOR ii
IF SIZE(nonzerorowlist)!=0 THEN
IF tab[nonzerorowlist[1],ookcol]!=0 THEN
tab[nonzerorowlist[1]]:=-1*tab[nonzerorowlist[1]];
tab:=caspivot(tab,nonzerorowlist[1],ookcol);
IF intgr==0 THEN PRINT("TESTING: ROW="+nonzerorowlist[1]+" COL="+ookcol);END;
IF tab[c+1,c+v+1]<0 THEN tab[c+1]:=-1*tab[c+1];END;
//PRINT("phaseone");
pack:=casphaseone(tab,c,v,stp1,intgr);
tab:=pack[1];
stp1:=pack[4];
neg:=2;
END;//IF tab[nonzerorowlist[1],ookcol]!=0
END;//IF SIZE(nonzerorowlist)!=0
END;//IF ook==0
/////////// TESTING ///////////////
IF neg==1 THEN
mini:=inf;r:=0;
FOR k FROM 1 TO c DO
IF tab[k,j]>0 THEN
IF RANDINT(1,2)==1 THEN
q1:=(tab[k,c+v+2]/tab[k,j])<mini;
ELSE
q1:=(tab[k,c+v+2]/tab[k,j])<=mini;
END;
q2:=(tab[k,c+v+2]/tab[k,j])≥0;
IF q1 THEN
mini:=tab[k,c+v+2]*(1/tab[k,j]);
r:=k;
END;
END;
END;
//
tab:=caspivot(tab,r,j);
//////////////////////////
IF intgr==0 THEN PRINT("SIMPLEX: ROW="+r+" COL="+j);END;
IF tab[c+1,c+v+1]<0 THEN tab[c+1]:=-1*tab[c+1];END;
//////////
//PRINT("phaseone");
pack:=casphaseone(tab,c,v,stp1,intgr);
tab:=pack[1];
stp1:=pack[4];
//////////
//RETURN(TAB);
END;
//PRINT("OOK="+OOK);
stp2:=stp2+1;
//PRINT("stp2="+stp2);
UNTIL neg==0 OR stp2>=100;
return({tab,c,v,stp2,stp1,ook});
END;
#end
#cas
cassimplex(tab,c,v,p,intgr,namelist):=BEGIN
local org,tempsimdata,retry;
local stp1,pack,pack2;
local s1,s2,s3;
local ABC,INFSOL,NSOL,SIG,temptab;
local stp2;
local nzr,nzn,r,t,mini;
local q1,q2,ook;
local ill2;
local mts,mtab,lmx;
local latab,ill3;
local i,k,simdata;
local nitem;
local solutionrow,m,samerow,firstpos,solutionlist;
local equ,posi,comment;
tab:=exact(tab);
org:=tab;
IF SIZE(namelist)!=(v+1) OR TYPE(namelist)!=6 THEN
namelist:={};
FOR i FROM 1 TO v DO
namelist[i]:="X"+i;
END;//FOR
namelist[v+1]:="Z";
END;//IF
///////////////////////////////////////////////////////////////////
//PHASE ONE
retry:=0;
REPEAT
stp1:=0;
pack:={};
pack:=casphaseone(tab,c,v,stp1,intgr);
tab:=pack[1];
//c:=pack(2);
//v:=pack(3);
stp1:=pack[4];
///////////////////////////////////////////////////////////////////
//SIMPLEX (Phase two)
stp2:=0;
pack:=casphasetwo(tab,c,v,stp2,intgr);
tab:=pack[1];
stp2:=pack[4];
stp1:=pack[5];
ook:=pack[6];
//////////////////////////////////////////////////////////////////
IF tab==org AND retry==0 THEN
retry:=1;
//PRINT("TEST");
FOR i FROM 1 to c+v DO
IF tab[c+1,i]==0 THEN
nitem:=0;
FOR k FROM 1 to c DO
IF tab[k,i]!=0 THEN nitem:=nitem+1;END;
END;//for k
IF nitem>1 THEN
mini:=inf;r:=0;
FOR t FROM 1 TO c DO
IF tab[t,i]>0 THEN
q1:=(tab[t,c+v+2]*(1/tab[t,i]))<mini;
q2:=(tab[t,c+v+2]*(1/tab[t,i]))≥0;
IF q1 THEN
mini:=tab[t,c+v+2]*(1/tab[t,i]);
r:=t;
END;//IF q1 and q2
END;//IF SIGN
END;//FOR t
IF r!=0 THEN
tab:=caspivot(tab,r,i);
IF intgr==0 THEN PRINT("PIVOT: ROW="+r+" COL="+i);END;
BREAK;
END;//IF r!=0
END;//if nitem>1 (zero col)
END;//if last row==0
END;//for i
ELSE //END;//if tab==org
retry:=0;
END;
UNTIL retry==0;
////////////////////////////////////////////////
//orgtab:=tab;
////////////////////////////////////////////////
//FIND VALID SOLUTION
IF intgr==0 THEN PRINT("---------------------");END;
latab:={};
latab[1]:=org;
/////////// Find multiple solutions list/////////////////////////
mts:={};
FOR i FROM 1 to c+v DO
nzn:=0;nzr:=0;
FOR k FROM 1 to c+1 DO
IF tab[k,i]!=0 THEN nzn:=nzn+1;nzr:=k;END;
END;//FOR k
IF nzn>1 AND tab[c+1,i]==0 THEN mts:=append(mts,i);END;
END;//FOR i
/////////// Find multiple solutions list/////////////////////////
lmx:={};lmx[1]:=tab;
/////////// PIVOT multiple col //////////////////////////////////
IF SIZE(mts)!=0 THEN
FOR i FROM 1 to SIZE(mts) DO
mini:=inf;r:=0;
FOR k FROM 1 TO c DO
IF tab[k,mts[i]]>0 THEN
q1:=(tab[k,c+v+2]*(1/tab[k,mts[i]]))<mini;
q2:=(tab[k,c+v+2]*(1/tab[k,mts[i]]))≥0;
IF q1 THEN
mini:=tab[k,c+v+2]*(1/tab[k,mts[i]]);
r:=k;
END;//if q1 and q2
END;//if SIGN
END;//for k
IF r!=0 THEN
stp1:=0;
stp2:=0;
mtab:=caspivot(tab,r,mts[i]);
pack:=casphaseone(mtab,c,v,stp1,1);
pack2:=casphasetwo(pack[1],c,v,stp2,1);
lmx:=append(lmx,pack2[1]);
stp1:=pack2[5];
stp2:=pack2[4];
END;//IF r!=0
END;//FOR i (size(mts))
END;//IF size(mts)!=0
/////////// PIVOT multiple col //////////////////////////////////
/////////// Find Solution ///////////////////////////////////////
solutionlist:={};simdata:={};
FOR m FROM 1 TO SIZE(lmx)(1) DO
// find solution row //////////////
solutionrow:={};mtab:=lmx[m];
FOR i FROM 1 TO c+v DO
nzn:=0;nzr:=0;
FOR k FROM 1 TO c+1 DO
IF lmx[m,k,i]!=0 THEN nzn:=nzn+1;nzr:=k;END;
END;//FOR k
IF nzn==1 AND nzr!=c+1 THEN solutionrow[i]:=nzr;
ELSE solutionrow[i]:=0;END;
END;//FOR i
// find solution row //////////////
// Infinitely many solutions //////
NSOL:=0;INFSOL:={};INFSOL[1]:="";
FOR i FROM c+v DOWNTO 1 DO
equ:=" = ";
IF solutionrow[i]!=0 THEN
IF lmx[m,solutionrow[i],i]>0 AND (lmx[m,solutionrow[i],c+v+2]>0 OR lmx[m,solutionrow[i],c+v+2]==0)THEN
INFSOL[NSOL+1]:="";
firstpos:=0;
IF i<=v THEN INFSOL[NSOL+1]:=" "+lmx[m,solutionrow[i],i]+"_"+namelist[i];firstpos:=1;END;
IF i>v THEN equ:=" ≤ " ;END;
samerow:=0;
FOR k FROM i-1 DOWNTO 1 DO
IF solutionrow[k]==solutionrow[i] THEN
IF k<=v AND lmx[m,solutionrow[i],k]>0 THEN
samerow:=1;
firstpos:=1;
SIG:=" - ";
IF lmx[m,solutionrow[i],k]>= 0 THEN SIG:=" + ";END;
INFSOL[NSOL+1]:=INFSOL[NSOL+1]+SIG+ABS(lmx[m,solutionrow[i],k])+"_"+namelist[k];
END;//if k<=v
IF k<=v AND lmx[m,solutionrow[i],k]<0 THEN
IF firstpos==1 THEN samerow:=1;END;
SIG:=" - ";
IF lmx[m,solutionrow[i],k]>= 0 THEN SIG:=" + ";END;
INFSOL[NSOL+1]:=INFSOL[NSOL+1]+SIG+ABS(lmx[m,solutionrow[i],k])+"_"+namelist[k];
END;
solutionrow[k]:=0;
END;//if same row
END;//for k
IF samerow==1 THEN INFSOL[NSOL+1]:=INFSOL[NSOL+1]+equ+lmx[m,solutionrow[i],c+v+2];NSOL:=NSOL+1;INFSOL[NSOL+1]:="";END;
END;//if tab>0
IF lmx[m,solutionrow[i],i]<0 AND (lmx[m,solutionrow[i],c+v+2]>0 OR lmx[m,solutionrow[i],c+v+2]==0) THEN
INFSOL[NSOL+1]:="";
IF i<=v THEN INFSOL[NSOL+1]:=" "+lmx[m,solutionrow[i],i]+"_"+namelist[i];END;
IF i>v THEN equ:=" ≥ ";END;
samerow:=0;
firstpos:=0;
posi:=i;
FOR k FROM i-1 DOWNTO 1 DO
IF solutionrow[k]==solutionrow[posi] THEN
IF k>v AND lmx[m,solutionrow[posi],k]>0 AND firstpos!=0 THEN
solutionrow[k]:=0;
END;
IF k>v AND lmx[m,solutionrow[posi],k]>0 AND firstpos==0 THEN
equ:=" ≤ ";
firstpos:=1;
solutionrow[posi]:=0;
posi:=k;
END;
IF k>v AND lmx[m,solutionrow[posi],k]<0 THEN
solutionrow[k]:=0;
END;
IF k<=v AND lmx[m,solutionrow[posi],k]>0 AND firstpos!=0 THEN
solutionrow[k]:=0;
samerow:=1;
SIG:=" - ";
IF lmx[m,solutionrow[posi],k]>= 0 THEN SIG:=" + ";END;
INFSOL[NSOL+1]:=INFSOL[NSOL+1]+SIG+ABS(lmx[m,solutionrow[posi],k])+"_"+namelist[k];
END;
IF k<=v AND lmx[m,solutionrow[posi],k]>0 AND firstpos==0 THEN
firstpos:=1;
solutionrow[posi]:=0;
posi:=k;
samerow:=1;
SIG:=" - ";
IF lmx[m,solutionrow[posi],k]>= 0 THEN SIG:=" + ";END;
INFSOL[NSOL+1]:=INFSOL[NSOL+1]+SIG+ABS(lmx[m,solutionrow[posi],k])+"_"+namelist[k];
END;
IF k<=v AND lmx[m,solutionrow[posi],k]<0 THEN
solutionrow[k]:=0;
SIG:=" - ";
IF lmx[m,solutionrow[posi],k]>= 0 THEN SIG:=" + ";END;
INFSOL[NSOL+1]:=INFSOL[NSOL+1]+SIG+ABS(lmx[m,solutionrow[posi],k])+"_"+namelist[k];
END;
//IF k<=v AND lmx[m,solutionrow[posi],k]<0 AND firstpos==0 THEN
// solutionrow[k]:=0;
//END;
END;//if solutionrow==solutionrow and tab<0
END;//for k
IF samerow==1 THEN INFSOL[NSOL+1]:=INFSOL[NSOL+1]+equ+lmx[m,solutionrow[posi],c+v+2];NSOL:=NSOL+1;INFSOL[NSOL+1]:="";END;
//IF firstpos==1 THEN solutionrow[i]:=0;END;
END;//if tab<0
END;//if solutionrow!=0;
END;//for i
IF NSOL>0 THEN ill3:=1; ELSE ill3:=0;END;
// Infinitely many solutions //////
///////////// check infeasible solution /////////////////////////
ill2:=0;
FOR i FROM 1 to c+v DO
IF solutionrow[i]!=0 THEN
IF lmx[m,solutionrow[i],i]*SIGN(lmx[m,solutionrow[i],c+v+2])<0 THEN ill2:=1;END;
END;//IF solutionrow[i]!=0
END;//FOR i
///////////// check infeasible solution /////////////////////////
s1:="";s2:="";s3:="";tempsimdata:={};ABC:=0;
///////////// prepare solution string ///////////////////////////
FOR i FROM 1 TO v DO
IF solutionrow[i]!=0 THEN
temptab:=lmx[m];
temptab[solutionrow[i]]:=temptab[solutionrow[i]]*(1/temptab[solutionrow[i],i]);
s1:=s1+namelist[i]+"="+STRING(temptab[solutionrow[i],c+v+2],7,-1)+" , ";
s2:=s2+namelist[i]+"="+STRING(temptab[solutionrow[i],c+v+2],14,-1)+" , ";
s3:=s3+namelist[i]+"="+STRING(temptab[solutionrow[i],c+v+2],2,4)+" , ";
tempsimdata:=append(tempsimdata,temptab[solutionrow[i],c+v+2]);
IF ROUND(temptab[solutionrow[i],c+v+2],0)!=temptab[solutionrow[i],c+v+2] THEN ABC:=1;END;
ELSE //ELSE solutionrow !=0
s1:=s1+namelist[i]+"=0 , ";
s2:=s2+namelist[i]+"=0 , ";
s3:=s3+namelist[i]+"=0 , ";
tempsimdata:=append(tempsimdata,0);
END;//IF solutionrow !=0
END;//FOR i
//////////////// Z /////
IF lmx[m,c+1,c+v+1]!=0 THEN
temptab:=lmx[m];
temptab[c+1]:=temptab[c+1]*(1/temptab[c+1,c+v+1]);
s1:=s1+namelist[v+1]+"="+STRING(((-1)^(p+1)*temptab[c+1,c+v+2]),7,-1);
s2:=s2+namelist[v+1]+"="+STRING(((-1)^(p+1)*temptab[c+1,c+v+2]),14,-1);
s3:=s3+namelist[v+1]+"="+STRING(((-1)^(p+1)*temptab[c+1,c+v+2]),2,4);
tempsimdata:=append(tempsimdata,(-1)^(p+1)*temptab[c+1,c+v+2]);
IF ROUND(temptab[c+1,c+v+2],0)!=temptab[c+1,c+v+2] THEN ABC:=1;END;
ELSE
s1:=s1+namelist[v+1]+"=NA"
s2:=s2+namelist[v+1]+"=NA"
s3:=s3+namelist[v+1]+"=NA"
tempsimdata:=append(tempsimdata,(-1)^p*inf);
END;
//////////////// P /////
///////////// prepare solution string ///////////////////////////
///////////// check NEW solution ///////////////////////////////
IF contains(solutionlist,tempsimdata)==0 THEN
solutionlist:=append(solutionlist,tempsimdata);
latab:=append(latab,lmx[m]);
///////////// PRINT solution ////////////
IF intgr==0 THEN
PRINT(s1);
IF s1!=s2 THEN
PRINT(s2);
END;
IF ABC==1 THEN
PRINT(s3);
END;
PRINT("---------------------");
comment:=0;
IF ill3 AND (NOT ill2) THEN PRINT("*Infinitely many solutions*");
PRINT("*Solution may satisfy below equation:*");
FOR i FROM 1 TO NSOL DO
PRINT(INFSOL[i]);
END;//FOR NSOL
comment:=1;
END;//IF ill3
IF ook==0 AND ill2==0 THEN PRINT("*Unbounded*");comment:=1;END;
IF ill2 THEN PRINT("*Infeasible solution*");comment:=1;END;
IF comment==1 THEN PRINT("---------------------");END;
END;//if intgr==0
///////////// PRINT solution ////////////
//////////// prepare simdata ////////////
//IF ill4 THEN simdata[m*2-1]:=2;END;//multi
simdata:=append(simdata,0);
IF ook==0 AND ill2==0 THEN simdata[SIZE(simdata)]:=3;END;//unbounded
IF ill2==1 THEN simdata[SIZE(simdata)]:=1;END;//infeasible
simdata:=append(simdata,tempsimdata);
//PRINT("simdata="+simdata);
//////////// prepare simdata ////////////
END;//if contains (solutionlist,tempsimdata)==0
///////////// check NEW solution ///////////////////////////////
END;//FOR m (solution tableau)
/////////// Find Solution ///////////////////////////////////////
IF SIZE(simdata)>2 and intgr==0 THEN PRINT("*Multiple optimal solutions*");END;
IF intgr==0 THEN return(latab); END;
IF intgr==1 THEN return(simdata); END;
IF intgr==2 THEN return({latab,simdata});END;
END;
#end
#cas
casintegerBB(tab,c,v,p,level,st,maxlevel,imask,bbormix,namelist):=BEGIN
local i,simdata,allinteger,newvector;
local newtab1,m,whichlevel,reachlevel;
local s1,s2,ABC;
whichlevel:=level;reachlevel:=level;
IF level==1 THEN
IF p==1 THEN best:={-inf};END;
IF p==2 THEN best:={inf};END;
END;
allinteger:=0;
//PRINT("before simplex");
simdata:=cassimplex(tab,c,v,p,1,{});
//PRINT(level+":"+(simdata));
//PRINT("Level="+level+" best="+best+" new="+simdata);
//*****
//IF contains(simdata,{4,0,1,0,5})!=0 and b==0 THEN a:=tab;b:=1;END;
//*****
m:=0;
REPEAT
m:=m+1;
//if feasible solution
IF simdata[m*2-1]!=1 THEN
IF ((simdata[m*2,v+1]>=best[1]) AND p==1) OR ((simdata[m*2,v+1]<=best[1]) AND p==2) THEN
allinteger:=1;
FOR i FROM 1 to v DO
IF imask[i]==1 THEN
IF ROUND(simdata[m*2,i],0)!=simdata[m*2,i] THEN
allinteger:=0;
reachlevel:=level+1;
IF level<maxlevel THEN
//PRINT("in level="+level+" simdata="+simdata+" m="+m+" i="+i);
//FLOOR ////////////////////////////////////////
newvector:=MAKELIST(0,X,1,c+v+2);
newvector[i]:=1;
newvector[c+v+2]:=FLOOR(simdata[m*2,i]);
newtab1:=tab;
ADDROW(newtab1,newvector,c+1);
newvector:=MAKELIST(0,X,1,c+1+1);
newvector[c+1]:=1;
ADDCOL(newtab1,newvector,v+c+1);
///////////////////////////////////////////////
//PRINT("X"+i+" <= "+FLOOR(simdata[m,i]));
//st+CHAR(10)+"X"+i+" <= "+FLOOR(simdata[m,i])
whichlevel:=casintegerBB(newtab1,c+1,v,p,level+1,"",maxlevel,imask,bbormix,{});
IF whichlevel>reachlevel THEN reachlevel:=whichlevel;END;
//PRINT("pass FLOOR");
//CEILING//////////////////////////////////////
newvector:=MAKELIST(0,X,1,c+v+2);
newvector[i]:=1;
newvector[c+v+2]:=CEILING(simdata[m*2,i]);
newtab1:=tab;
ADDROW(newtab1,newvector,c+1);
newvector:=MAKELIST(0,X,1,c+1+1);
newvector[c+1]:=-1;
ADDCOL(newtab1,newvector,v+c+1);
///////////////////////////////////////////////
//PRINT("X"+i+" >= "+CEILING(simdata[m,i]));
//st+CHAR(10)+"X"+i+" >= "+CEILING(simdata[m,i])
whichlevel:=casintegerBB(newtab1,c+1,v,p,level+1,"",maxlevel,imask,bbormix,{});
IF whichlevel>reachlevel THEN reachlevel:=whichlevel;END;
///////////////////////////////////////////////
END;//IF level<maxlevel
//return(newtab1);
END;//IF not integer
END;//IF imask==1
END;//FOR i
END;//if not better than best
END;//IF feasible solution
IF allinteger==1 THEN
IF ((simdata[m*2,v+1]>=best[1]) AND p==1) OR ((simdata[m*2,v+1]<=best[1]) AND p==2) THEN
IF simdata[m*2,v+1]==best[1] THEN
IF contains(best,simdata[m*2])==0 THEN
best:=append(best,simdata[m*2]);
//best:=INTERSECT(best,best);
END;//contains
ELSE
best:={0};
best[1]:=simdata[m*2,v+1];
best[2]:=simdata[m*2];
END;//equal
END;//if best P
END;//IF allinteger
//PRINT("m="+m+" simdata="+SIZE(simdata)+" "+simdata);
UNTIL m*2==SIZE(simdata);
IF level==1 THEN
//PRINT(best);
IF SIZE(best)>1 THEN
PRINT("---------------------");
IF bbormix==0 THEN PRINT("Integer Solution [BRANCH & BOUND]");
ELSE PRINT("Mixed Solution [Mixed INT.CONT.BIN.]");END;
IF simdata[1]==3 THEN PRINT("*a possible solution*");END;
//IF reachlevel >= maxlevel THEN PRINT("*Try to increase search level for more solutions*");END;
//best:=INTERSECT(best,best);
//best:=mat2list(best);
m:=1;
REPEAT
m:=m+1;s1:="";s2:="";ABC:=0;
FOR i FROM 1 to v DO
s1:=s1+namelist[i]+"="+best[m,i]+" , ";
s2:=s2+namelist[i]+"="+STRING(best[m,i],2,4)+" , ";
IF ROUND(best[m,i],0)!=best[m,i] THEN ABC:=1;END;
END;//FOR i
s1:=s1+namelist[v+1]+"="+best[m,v+1];
s2:=s2+namelist[v+1]+"="+STRING(best[m,v+1],2,4);
IF ROUND(best[m,v+1],0)!=best[m,v+1] THEN ABC:=1;END;
PRINT(s1);
IF bbormix==1 AND ABC==1 THEN PRINT(s2);END;
UNTIL m==SIZE(best);
ELSE
PRINT("---------------------");
IF bbormix==0 THEN PRINT("*No solution [BRANCH & BOUND]*");
ELSE PRINT("*No Solution [Mixed INT.CONT.BIN.]*");END;
IF reachlevel>maxlevel THEN PRINT("*Try to increase search level*");END;
END;//if size(best)
//PRINT("-END-");
return(reachlevel);
END;//if level 1
return(reachlevel);
END;
#end
#cas
casintegercut(tab,c,v,p,level,st,maxlevel,namelist):=BEGIN
local i,simdata,allinteger,newvector,pack,simtab;
local newtab1,m,k,cutrow,whichlevel,reachlevel;
local s1,noplaneerror,newdata,s2,ABC;
whichlevel:={};reachlevel:=level;noplaneerror:=1;newdata:={};
IF level==1 THEN
IF p==1 THEN best:={-inf};END;
IF p==2 THEN best:={inf};END;
END;
allinteger:=0;
//PRINT("Simplex");
//a:=tab;
pack:={};
pack:=cassimplex(tab,c,v,p,2,{});
simtab:=pack[1];
simdata:=pack[2];
//PRINT("pass simplex");
//a:=simtab;
//b:=simdata;
//PRINT(level+":"+(simdata));
//PRINT("Level="+level+" best="+best+" new="+simdata);
//IF level==1 THEN a:=simtab; END;
//PRINT(simdata);
//PRINT(exact(simtab));
m:=0;
REPEAT
m:=m+1;
//if feasible solution
IF simdata[m*2-1]!=1 THEN
//PRINT("feasbile="+m);
IF ((simdata[m*2,v+1]>=best[1]) AND p==1) OR ((simdata[m*2,v+1]<=best[1]) AND p==2) THEN
allinteger:=1;
//PRINT("Integer="+m);
FOR i FROM 1 to v DO
IF ROUND(simdata[m*2,i],0)!=simdata[m*2,i] THEN
allinteger:=0;
//PRINT("non integer="+m);
reachlevel:=level+1;
IF level<maxlevel THEN
/////
/////
//FIND ROW /////////////////////////////////////
//PRINT("find row="+m);
//a:=simtab;
//b:=simdata;
cutrow:=0;
FOR k FROM 1 to c DO
IF simtab[m+1,k,i]==1 THEN cutrow:=k;END;
END;//for k
IF cutrow!=0 THEN
newvector:=MAKELIST(0,X,1,c+v+2);
FOR k FROM 1 to c+v+2 DO
newvector[k]:=simtab[m+1,cutrow,k]-FLOOR(simtab[m+1,cutrow,k]);
END;//FOR k
///no Plane ERROR
noplaneerror:=0;
FOR k FROM 1 to c+v+1 DO
IF newvector[k]!=0 THEN noplaneerror:=noplaneerror+1;END;
END;//FOR error K
////////
IF noplaneerror!=0 THEN
///////////////
newtab1:=simtab[m+1];
ADDROW(newtab1,newvector,c+1);
newvector:=MAKELIST(0,X,1,c+1+1);
newvector[c+1]:=-1;
ADDCOL(newtab1,newvector,v+c+1);
//return(newtab1);
//casintegercut ////////////////////////////////////////
////////////////////////////////////////////////////////
//PRINT("enter next level");
//PRINT("simdata="+simdata);
whichlevel:=casintegercut(newtab1,c+1,v,p,level+1,"",maxlevel,{});
IF whichlevel[1]>reachlevel THEN reachlevel:=whichlevel[1]; END;
IF noplaneerror!=0 AND whichlevel[2]==0 THEN noplaneerror:=0;END;
newdata:=whichlevel[3];
////////////////////////////////////////////////////////
END;//IF noplaneerror!=0
END;//IF cutrow!=0
END;//IF level<maxlevel
//return(newtab1);
END;//IF not integer
END;//FOR i
END;//if not better than best
END;//IF feasible solution
IF allinteger==1 THEN
//PRINT("integer==1");
IF ((simdata[m*2,v+1]>=best[1]) AND p==1) OR ((simdata[m*2,v+1]<=best[1]) AND p==2) THEN
IF simdata[m*2,v+1]==best[1] THEN
IF contains(best,simdata[m*2])==0 THEN
best:=append(best,simdata[m*2]);
//best:=INTERSECT(best,best);
END;//contains
ELSE
best:={0};
best[1]:=simdata[m*2,v+1];
best[2]:=simdata[m*2];
END;//equal
END;//if best P
END;//IF allinteger
//PRINT("repeat");
UNTIL m*2==SIZE(simdata);
//PRINT("XXXXXXXXXXXXXX");
//a:=simdata;
//b:=newdata;
IF newdata=={} THEN
newdata:=simdata;
END;
//PRINT("YYYYYYYYYYYYYY");
IF level==1 THEN
//PRINT(best);
IF SIZE(best)>1 THEN
PRINT("---------------------");
PRINT("Integer Solution [CUTTING PLANE]");
IF simdata[1]==3 THEN PRINT("*a possible solution*");END;
//IF reachlevel >= maxlevel THEN PRINT("*Try to increase search level for more solutions*");END;
m:=1;
REPEAT
m:=m+1;s1:="";s2:="";ABC:=0
FOR i FROM 1 to v DO
s1:=s1+namelist[i]+"="+best[m,i]+" , ";
s2:=s2+namelist[i]+"="+STRING(best[m,i],2,4)+" , ";
IF ROUND(best[m,i],0)!=best[m,i] THEN ABC:=1;END;
END;//FOR i
s1:=s1+namelist[v+1]+"="+best[m,v+1];
s2:=s2+namelist[v+1]+"="+STRING(best[m,v+1],2,4);
IF ROUND(best[m,v+1],0)!=best[m,v+1] THEN ABC:=1;END;
PRINT(s1);
//IF ABC==1 THEN PRINT(s2);END;
UNTIL m==SIZE(best);
ELSE
PRINT("---------------------");
PRINT("*No integer solution [CUTTING PLANE]*");
IF reachlevel > maxlevel THEN PRINT("*Try to increase search level*");END;
IF noplaneerror==0 THEN PRINT("*Try to use B&B*");END;
END;//if size(best)
//PRINT("-END-");
return({reachlevel,noplaneerror,newdata});
END;//if level 1
return({reachlevel,noplaneerror,newdata});
END;
#end
#cas
casintegerBcut(tab,c,v,p,level,st,maxlevel,sublevel,namelist):=BEGIN
local i,simdata,allinteger,newvector;
local newtab1,m,whichlevel,reachlevel;
local s1,s2,ABC;
whichlevel:=level;reachlevel:=level;
IF level==1 THEN
IF p==1 THEN best:={-inf};END;
IF p==2 THEN best:={inf};END;
END;
allinteger:=0;
//PRINT("integer cut");
//simdata:=cassimplex(tab,c,v,p,1);
//casintegercut(newtab1,c+1,v,p,level+1,"",maxlevel);
whichlevel:=casintegercut(tab,c,v,p,2,"",sublevel,{});
simdata:=whichlevel[3];
IF (whichlevel[1]-1)>reachlevel THEN reachlevel:=(whichlevel[1]-1); END;
//b:=simdata;
//PRINT(level+":"+(simdata));
//PRINT("Level="+level+" best="+best+" new="+simdata);
//a:=tab;b:=simdata;
//*****
//IF contains(simdata,{4,0,1,0,5})!=0 and b==0 THEN a:=tab;b:=1;END;
//*****
m:=0;
REPEAT
m:=m+1;
//if feasible solution
IF simdata[m*2-1]!=1 AND simdata!={} THEN
IF ((simdata[m*2,v+1]>=best[1]) AND p==1) OR ((simdata[m*2,v+1]<=best[1]) AND p==2) THEN
allinteger:=1;
FOR i FROM 1 to v DO
IF ROUND(simdata[m*2,i],0)!=simdata[m*2,i] THEN
allinteger:=0;
reachlevel:=level+1;
IF level<maxlevel THEN
//FLOOR ////////////////////////////////////////
newvector:=MAKELIST(0,X,1,c+v+2);
newvector[i]:=1;
//y:=simdata;
//z:=simdata[m,i];
newvector[c+v+2]:=FLOOR(simdata[m*2,i]);
newtab1:=tab;
ADDROW(newtab1,newvector,c+1);
newvector:=MAKELIST(0,X,1,c+1+1);
newvector[c+1]:=1;
ADDCOL(newtab1,newvector,v+c+1);
///////////////////////////////////////////////
//PRINT("X"+i+" <= "+FLOOR(simdata[m,i]));
//st+CHAR(10)+"X"+i+" <= "+FLOOR(simdata[m,i])
whichlevel:=casintegerBcut(newtab1,c+1,v,p,level+1,"",maxlevel,sublevel,{});
IF whichlevel>reachlevel THEN reachlevel:=whichlevel;END;
//CEILING//////////////////////////////////////
newvector:=MAKELIST(0,X,1,c+v+2);
newvector[i]:=1;
newvector[c+v+2]:=CEILING(simdata[m*2,i]);
newtab1:=tab;
ADDROW(newtab1,newvector,c+1);
newvector:=MAKELIST(0,X,1,c+1+1);
newvector[c+1]:=-1;
ADDCOL(newtab1,newvector,v+c+1);
///////////////////////////////////////////////
//PRINT("X"+i+" >= "+CEILING(simdata[m,i]));
//st+CHAR(10)+"X"+i+" >= "+CEILING(simdata[m,i])
whichlevel:=casintegerBcut(newtab1,c+1,v,p,level+1,"",maxlevel,sublevel,{});
IF whichlevel>reachlevel THEN reachlevel:=whichlevel;END;
///////////////////////////////////////////////
END;//IF level<maxlevel
//return(newtab1);
END;//IF not integer
END;//FOR i
END;//if not better than best
END;//IF feasible solution
IF allinteger==1 THEN
IF ((simdata[m*2,v+1]>=best[1]) AND p==1) OR ((simdata[m*2,v+1]<=best[1]) AND p==2) THEN
IF simdata[m*2,v+1]==best[1] THEN
IF contains(best,simdata[m*2])==0 THEN
best:=append(best,simdata[m*2]);
//best:=INTERSECT(best,best);
END;//contains
ELSE
best:={0};
best[1]:=simdata[m*2,v+1];
best[2]:=simdata[m*2];
END;//equal
END;//if best P
END;//IF allinteger
UNTIL m*2==SIZE(simdata);
IF level==1 THEN
//PRINT(best);
IF SIZE(best)>1 THEN
PRINT("---------------------");
PRINT("Integer Solution [MIXED BRANCH & CUT]");
IF simdata[1]==3 THEN PRINT("*a possible solution*");END;
//IF reachlevel >= maxlevel THEN PRINT("*Try to increase search level for more solutions*");END;
//best:=INTERSECT(best,best);
//best:=mat2list(best);
m:=1;
REPEAT
m:=m+1;s1:="";s2:="";ABC:=0;
FOR i FROM 1 to v DO
s1:=s1+namelist[i]+"="+best[m,i]+" , ";
s2:=s2+namelist[i]+"="+STRING(best[m,i],2,4)+" , ";
IF ROUND(best[m,i],0)!=best[m,i] THEN ABC:=1;END;
END;//FOR i
s1:=s1+namelist[v+1]+"="+best[m,v+1];
s2:=s2+namelist[v+1]+"="+STRING(best[m,v+1],2,4);
IF ROUND(best[m,v+1],0)!=best[m,v+1] THEN ABC:=1;END;
PRINT(s1);
//IF ABC==1 THEN PRINT(s2);END;
UNTIL m==SIZE(best);
ELSE
PRINT("---------------------");
PRINT("*No integer solution [MIXED BRANCH & CUT]*");
IF reachlevel>maxlevel THEN PRINT("*Try to increase search level*");END;
END;//if size(best)
//PRINT("-END-");
return(reachlevel);
END;//if level 1
return(reachlevel);
END;
#end
#cas
sim20(reedit):=BEGIN
local inputdata,tab,c,v,p,level;
local result,time1,time2,reachlevel;
local method1,method2,method3,sublevel;
local imask,method4,mixtab,cc,tempmask;
local namelist;
IF SIZE(reeditbag)<2 OR TYPE(reeditbag)!=6 THEN reeditbag:={{},{}};END;
IF reedit>=1 AND SIZE(reeditbag[reedit])!=15 THEN reeditbag[reedit]:={};END;
IF reedit==99 THEN
inputdata:=IMPORTDATA();
IF inputdata==0 THEN return(0);END;
tab:=exact(inputdata[1]);
c:=inputdata[2];
v:=inputdata[3];
p:=inputdata[4];
level:=inputdata[5];
sublevel:=inputdata[6];
method1:=inputdata[7];
method2:=inputdata[8];
method3:=inputdata[9];
method4:=inputdata[10];
imask:=inputdata[11];
mixtab:=inputdata[12];
cc:=inputdata[13];
namelist:=inputdata[14];
ELSE //IMPORTDATA
inputdata:=INPUTDATA(reedit,reeditbag[reedit]);
IF inputdata==0 THEN return(0);END;
tab:=exact(inputdata[1]);
c:=inputdata[2];
v:=inputdata[3];
p:=inputdata[4];
level:=inputdata[5];
sublevel:=inputdata[6];
method1:=inputdata[7];
method2:=inputdata[8];
method3:=inputdata[9];
method4:=inputdata[10];
imask:=inputdata[11];
mixtab:=inputdata[12];
cc:=inputdata[13];
reeditbag[1]:=inputdata[14];
namelist:=inputdata[15];
IF reedit>=2 THEN reeditbag[reedit]:=inputdata[14];END;
END;//INPUTDATA
IF method1==1 OR method2==1 OR method3==1 THEN intgr==1;
ELSE intgr==0;
END;
result:=cassimplex(tab,c,v,p,0,namelist);
//IF intgr==0 THEN PRINT("-END-");END;
IF method1==1 THEN
tempmask:=MAKELIST(1,X,1,v);
time1:=TICKS;
reachlevel:=casintegerBB(tab,c,v,p,1,"",level,tempmask,0,namelist);
time2:=TICKS;
IF reachlevel<=level THEN PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+reachlevel);
ELSE PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+level+" -> "+reachlevel);
END;//IF reachlevel
END;//if Branch
IF method2==1 THEN
time1:=TICKS;
reachlevel:=casintegercut(tab,c,v,p,1,"",level,namelist);
time2:=TICKS;
IF reachlevel[1]<=level THEN PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+reachlevel[1]);
ELSE PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+level+" -> "+reachlevel[1]);
END;//if reachlevel
END;//if cut
IF method3==1 THEN
time1:=TICKS;
reachlevel:=casintegerBcut(tab,c,v,p,1,"",level,sublevel+1,namelist);
time2:=TICKS;
IF reachlevel<=level THEN PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+reachlevel);
ELSE PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+level+" -> "+reachlevel);
END;//IF reachlevel
END;//if Branch
IF method4==1 THEN
time1:=TICKS;
reachlevel:=casintegerBB(mixtab,cc,v,p,1,"",level,imask,1,namelist);
time2:=TICKS;
IF reachlevel<=level THEN PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+reachlevel);
ELSE PRINT("TIME="+STRING((time2-time1)/1000,2,3)+"sec Reach Level: "+level+" -> "+reachlevel);
END;//IF reachlevel
END;//if Branch
PRINT("-END-");
return(result);
END;
#end
 sim20_13.hpprgm (Size: 87.94 KB / Downloads: 167)
|
Search
Member List
Calendar
Help
/
Simplex Method updated to V20.13 (cas mode)
simA.hpappdir.zip (Size: 1.95 KB / Downloads: 30)
