MODEL:
! Find alternative primal optimal extreme point solutions (EnumrXtrmC.lng)
to a Linear Program, and report the average of all of them.
All the solutions found are written to the file:
\temp\AltOpt.txt
! Example: The following model has
3 alternative primal extreme point optima
(X1 = 1 or X2 = 1 or X3 =1) and
3 alternative dual extreme point optima,
(DP1 = 1 or DP2 = 1/2 or DP3 = 1/3) for a total of
9 alternative optimal extreme point solutions.
We consider only the primal extreme point solutions as distinct.
MIN= X1 + X2 + X3,
X1 + X2 + X3 >= 1,
2* X1 + 2* X2 + 2* X3 >= 2,
3* X1 + 3* X2 + 3* X3 >= 3,
It reports
The average interior solution is:
X1 = 0.3333333333333333
X2 = 0.3333333333333333
X3 = 0.3333333333333333
Method: It is computationally expensive to enumerate every
distinct alternative optimal solution.
The approach used here is to randomly choose a large number of
random objectives/directions and then
minimize this objective subject to:
the solution is feasible and is
optimal for the original objective.
If enough random directions are generated, then
all alternative primal optimal solutions will be generated.;
! Keywords: Alternative optima, Average solution, Basic solution,
Corner point, Enumeration, Extreme point, K-Best solutions;
SETS:
ROW: RHS, TYPER;
COL: OBJ, X, SUB, OBJPTRB;
NONZ( ROW, COL): COEF;
SLIST;
SXC( SLIST, COL): XSV;
ENDSETS
DATA:
! Here are a number of cases;
! Case Murty: Example based on
"A Problem in Enumerating Extreme Points, and an efficient Algorithm,"
by Katta G. Murty, Dept Ind. and Opns Eng.
U. Michigan, March 2007.
! Number of desired extreme points. (There appear to be 6 distinct primal solutions);
!Cmurty NDXP = 12;
!Cmurty SLIST = 1..NDXP;
!Find all extreme points;
! Names of rows;
!Cmurty ROW= R01 R02 R03 R04 R05 R06 R07 R08 R09 R10 R11 R12 R13 R14 R15 ;
!Cmurty RHS= 1;
!Cmurty TYPER = 0; !Cmurty COL =
X01 X02 X03 X04 S01 S02 S03 S04 S05 S06 S07 S08 S09 S10 S11 S12 S13 S14 S15 ;
!Cmurty OBJ= 0;
! Upper bounds;
!Cmurty SUB=
1 1 1 1 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5;
! Constraint coefficients, with slacks explicitly added;
!Cmurty COEF =
1 1 1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1, 1, -1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
1, 1, -1, -1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
1, -1, 1, 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
1, -1, 1, -1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
1, -1, -1, 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
1, -1, -1 -1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
-1, 1, 1, 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
-1, 1, 1, -1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
-1, 1, -1, 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
-1, 1, -1, -1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
-1, -1, 1, 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
-1, -1, 1, -1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
-1 -1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
-1, -1, -1, -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ;
! Case OML6: Example from chapter 6 of Optimization Modeling with LINGO.
! It has two alternative optima;
!Number of desired extreme points.;
!COLM6 NDXP = 3;
!COLM6 SLIST = 1..3;
!Find all extreme points with profit contribution = optimal value;
! Names of rows;
!COLM6 ROW= MACH CAP HOURS SHIP PROD1 PROD2 PROD3;
!COLM6 RHS= 35 35 36 600 218 114 111 ;
!COLM6 TYPER= 1 1 1 1 -1 -1 -1; !COLM6 COL =
B34 B38 B48 B58 B42 B52 ;
!COLM6 OBJ=
-15.89 -17.89 -16.5 -15.22 -17.5 -16.22 ;
! Upper bounds;
!COLM6 SUB=
999 999 999 999 999 999 ;
! Constraint coefficients triples, ( Row Col Coefficient);
!COLM6 NONZ COEF =
MACH B34 0.11111 MACH B38 0.11111
CAP B48 0.16667 CAP B42 0.16667
HOURS B58 0.22222 HOURS B52 0.22222
SHIP B34 1 SHIP B38 1 SHIP B48 1 SHIP B58 1 SHIP B42 1 SHIP B52 1
PROD1 B34 1
PROD2 B38 1 PROD2 B48 1 PROD2 B58 1
PROD3 B42 1 PROD3 B52 1;
! Case Dyera: Example from Dyer and Proll, Math. Programming, vol. 12, 1977
in <= constraint form;
! It has 10 primal extreme points;
! Number of desired extreme point tries;
!CDyera NDXP = 60;
!Number of desired extreme points.;
!CDyera SLIST = 1..NDXP;
!CDyera ROW= ROW1 ROW2 ROW3 ROW4 ROW5;
!CDyera RHS= 5 16 3 17 10;
!CDyera TYPER = 1; !CDyera COL= X1 X2 X3 ;
!CDyera OBJ= 0 0 0;
!CDyera SUB= 99 99 99;
! Constraint coefficients;
!CDyera COEF =
3 2 -1
3 2 4
3 0 -4
2.25 4 3
1 2 1
;
! Case Xueyu 1;
! It has two alternative optima;
! Number of desired extreme points to try.;
!CX1 NDXP = 5;
!CX1 SLIST = 1..NDXP;
! Names of rows;
!CX1 ROW = CON1 CON2 CON3 CON4;
!CX1 RHS = -6 4 8 6;
!CX1 TYPER = 0; !CX1 COL =
X1 X2 SLK1 SLK2 SLK3 SLK4;
! Upper bounds;
!CX1 SUB=
999 999 999 999 999 999 ;
!CX1 OBJ =
-2 -1 0 0 0 0;
! Constraint coefficients;
!CX1 COEF =
-3 -2 1 0 0 0
2 -1 0 1 0 0
2 1 0 0 1 0
-1 2 0 0 0 1
;
! Case Xueyu 2;
! It has 7 extreme points;
! Number of desired extreme points tries.;
!CX2 NDXP = 10;
!CX2 SLIST = 1..NDXP;
! Names of rows;
!CX2 ROW = CON1 CON2 CON3 CON4 CON5 CON6 CON7;
!CX2 RHS = -1 -4 -12 -36 -26 -3 14;
!CX2 TYPER = 0; !CX2 COL =
X1 X2 SLK1 SLK2 SLK3 SLK4 SLK5 SLK6 SLK7;
! Upper bounds;
!CX2 SUB=
999 999 999 999 999 999 999 999 999;
!CX2 OBJ =
0 0 0 0 0 0 0 0 0;
! Constraint coefficients;
!CX2 COEF =
2 -1 -1 0 0 0 0 0 0
3 -2 0 -1 0 0 0 0 0
1 -2 0 0 -1 0 0 0 0
-3 -2 0 0 0 -1 0 0 0
-4 1 0 0 0 0 -1 0 0
-1 2 0 0 0 0 0 -1 0
2 4 0 0 0 0 0 0 -1
;
! Case: Variant of Astro-Cosmo problem.
! It has two alternative optima;
! Max = 15 A + 30 C
A <= 60,
C <= 50,
A + 2 C <= 120,
! Number of desired extreme points tries.;
!CAC1 NDXP = 3;
!CAC1 SLIST = 1..NDXP;
! The last three variables in the tableau,
cols 3-5, are the slack variables;
!CAC1 ROW = ALIM CLIM LABOR;
!CAC1 RHS = 60 50 120;
!CAC1 TYPER = 0; ! Names of columns;
!CAC1 COL = A__ C__ S_1 S_2 S_3;
!CAC1 OBJ = -15 -30 0 0 0;
! Upper bounds;
!CAC1 SUB = 60 50 999 999 999;
! Names of rows;
! Matrix coefficients, including Obj and RHS;
!CAC1 COEF =
1 0 1 0 0
0 1 0 1 0
1 2 0 0 1 ;
! Case 2: Variant of extended Astro-Cosmo problem.
! It has two alternative optima;
! Max = 15 A + 30 C + 44 D,
A + D <= 60,
C + D <= 50,
A + 2 C + 3 D <= 120,
! Case PD3;
! It has 3 alternative primal extreme point optima, and
3 alternative dual extreme point optima, for a total
of 9 alternative extreme point optima;
! MIN= X1 + X2 + X3,
X1 + X2 + X3 >= 1,
2* X1 + 2* X2 + 2* X3 >= 2,
3* X1 + 3* X2 + 3* X3 >= 3,
! Number of desired extreme points tries.;
!CPD3 NDXP = 8;
!CPD3 SLIST = 1..NDXP;
! Names of rows;
!CPD3 ROW = CON1 CON2 CON3 ;
!CPD3 RHS = 1 2 3;
!CPD3 TYPER = -1; !CPD3 COL =
X1 X2 X3 ;
! Upper bounds;
!CPD3 SUB=
999 999 999 ;
!CPD3 OBJ =
1 1 1 ;
! Constraint coefficients in sparse form;
!CPD3 NONZ, COEF =
CON1 X1 1
CON1 X2 1
CON1 X3 1
CON2 X1 2
CON2 X2 2
CON2 X3 2
CON3 X1 3
CON3 X2 3
CON3 X3 3
;
! Number of desired extreme points tries.;
!CAC2 NDXP = 3;
!CAC2 SLIST = 1..NDXP;
! Names of rows;
!CAC2 ROW = ALIM CLIM LABOR;
!CAC2 RHS = 60 50 120;
!CAC2 TYPER = 0; ! The last three variables in the tableau,
cols 4-6, are the slack variables;
! Names of columns;
!CAC2 COL= A__ C__ D__ S_1 S_2 S_3 ;
!CAC2 OBJ= -15 -30 -44 0 0 0 ;
! Upper bounds;
!CAC2 SUB = 60 50 40 999 999 999;
! Constraint coefficients in sparse form;
!CAC2 NONZ, COEF =
ALIM A__ 1
ALIM D__ 1
ALIM S_1 1
CLIM C__ 1
CLIM D__ 1
CLIM S_2 1
LABOR A__ 1
LABOR C__ 2
LABOR D__ 3
LABOR S_3 1;
! Case 3: Example from Dyer and Proll, Math. Programming, vol. 12, 1977
in equality constraint form;
! It has 10 primal extreme points;
! Number of desired extreme point tries;
!CDyer; NDXP = 60!Number of desired extreme points.;
!CDyer; SLIST = 1..NDXP!CDyer; ROW= ROW1 ROW2 ROW3 ROW4 ROW5!CDyer; RHS= 5 16 3 17 10!CDyer; TYPER = 0!CDyer; COL= X1 X2 X3 X4 X5 X6 X7 X8 !CDyer; OBJ= 0 0 0 0 0 0 0 0!CDyer; SUB= 99 99 99 99 99 99 99 99! Constraint coefficients;
!CDyer; COEF =
3 2 -1 1 0 0 0 0
3 2 4 0 1 0 0 0
3 0 -4 0 0 1 0 0
2.25 4 3 0 0 0 1 0
1 2 1 0 0 0 0 1
;
! We can read the data from a currently open( only one)
spreadsheet with range names matching LINGO names;
!CExcel NDXP = @OLE() ;
!Number of desired extreme points.;
!CExcel SLIST = 1..NDXP;
! Names of rows;
!CExcel ROW = @OLE();
!CExcel RHS = @OLE();
!CExcel TYPER = @OLE(); ! Names of columns;
!CExcel COL= @OLE();
! Objective coefficients;
!CExcel OBJ= @OLE();
! Upper bounds;
!CExcel SUB = @OLE();
! Constraint coefficients in sparse form;
!CExcel NONZ, COEF = @OLE();
ENDDATA
SUBMODEL Simple ! A weighted objective;
MIN = ALPHA* OBVTRUE + (1 - ALPHA)* OBVPRTRB;
@FREE( OBVTRUE); ! The true objective;
OBVTRUE = @SUM( COL( j) : OBJ( j)* X( j));
@FREE( OBVPRTRB); ! The perturbed objective;
OBVPRTRB = @SUM( COL( j) : OBJPTRB( j)* X( j));
! Allow an upper bound constraint on the true Objective value;
OBVTRUE <= OBVUL;
@FREE( OBVUL); ! It could be negative;
@FOR( ROW( I) | TYPER( I) #GT# 0: ! <= rows;
@SUM( NONZ( I, J): COEF( I, J) * X( J)) <= RHS( I)
);
@FOR( ROW( I) | TYPER( I) #EQ# 0:! = rows;
@SUM( NONZ( I, J): COEF( I, J) * X( J)) = RHS( I)
);
@FOR( ROW( I) | TYPER( I) #LT# 0: ! >= rows;
@SUM( NONZ( I, J): COEF( I, J) * X( J)) >= RHS( I)
);
@FOR( COL( J): ! All vars are bounded: 0 <= X( J) <= SUB( J);
@BND( 0, X( J), SUB( J))
);
ENDSUBMODEL
CALC:
@SET( 'OROUTE', 1); ! Buffer size for routing output to window;
@SET( 'TERSEO', 3); ! Output level (0:verb, 1:terse, 2:only errors, 3:none);
@SET( 'SOLVEL', 3); ! Linear solver (0:LINGO decides,1:primal,2:dual,3:barrier);
TOFILE = 1; ! 1 means send solutions to file \temp\AltOpt.txt;
NCOL = @SIZE( COL);
TOLDIF = 0.00000001; ! Tolerance for considering two values different;
ALPHA = 1; ! Use original true objective;
@FOR( COL( j): OBJPTRB( j) = 0); ! Initially, no perturbations;
@SOLVE( SimpleMODEL);
@WRITE( 'The Objective value= ', OBVTRUE, @NEWLINE( 1));
ISTAT = @STATUS();
@WRITE( ' Solve Status= ', ISTAT, @NEWLINE( 1));
ITER = 1;
NSFND = 1; ! Number of solutions found;
! Save this solution;
@FOR( COL( j): XSV( NSFND, j) = X( j)); ! Save current solution;
@WRITE(' Solution# ', NSFND, ' (Iter: ', ITER,')', @NEWLINE( 1));
@DIVERT( '\temp\AltOpt.txt'); ! Send output to a file;
@WRITE(' Solution# ', NSFND, ' (Iter: ', ITER,')', @NEWLINE( 1));
@FOR( COL( j) | X( j) #GT# 0:
TEMP = X( j);
@WRITE(' ', COL( j),' ', TEMP, @NEWLINE( 1));
);
@WRITE(@NEWLINE( 1));
@DIVERT(); ! Revert to output to the screen;
! Restrict to just optimal solutions;
OBVUL = OBVTRUE;
! Switch to perturb objective;
ALPHA = 0;
! How to choose next direction;
MODE = 1; ! 1: random, 2: flip, 3: mod;
! Generate a perturbed objective;
! A arbitrary random starting seed;
URANDV = 0.13245768;
@WHILE( ITER #LT# NDXP:
ITER = ITER + 1;
! Generate a direction, try to make different from previous directions;
@IFC( MODE #EQ# 1: !1) Choose a random direction;
@FOR( COL( j):
URANDV = @RAND( URANDV); ! Get next random uniform;
OBJPTRB( j) = URANDV - 0.5; ! Uniform in ( -0.5, +0.5);
);
);
@IFC( MODE #EQ# 2 #OR# MODE #EQ# 4: !2) Choose opposite direction;
@FOR( COL( j):
OBJPTRB( j) = - OBJPTRB( j);
);
);
@IFC( MODE #EQ# 3: !3) Choose approximately orthogonal direction,
e.g. ( -0.1, +0.4) becomes (+0.4, -0.1);
@FOR( COL( j):
TEMP = OBJPTRB( j) + 0.5;
@IFC( TEMP #GT# 0.5: TEMP = TEMP - 1.0;
OBJPTRB( j) = TEMP;
);
);
);
MODE = MODE + 1; ! Use different direction for next iteration;
@IFC( MODE #GT# 4: MODE = 1);
@SOLVE( SimpleMODEL);
ISTAT = @STATUS();
@WRITE( ' Solve Status= ', ISTAT, @NEWLINE( 1));
! Check if we have already seen this primal solution;
NOTDUP = 1; ! Is it a duplicate? ;
IA = 0; ! Loop over already found solutions;
@WHILE( IA #LT# NSFND #AND# NOTDUP:
IA = IA + 1;
ISAME = 1; ! Duplicate until found otherwise;
IB = 0; ! Loop over variables in solution IA;
@WHILE( IB #LT# NCOL #AND# ISAME:
IB = IB + 1;
@IFC( @ABS( X( IB) - XSV( IA, IB)) #GT# TOLDIF:
ISAME = 0; ! Current soln and soln IA are different;
);
);
@IFC( ISAME: NOTDUP = 0); ! Current soln and soln IA are same;
);
@IFC( NOTDUP:
! We found a new primal alternative optima;
NSFND = NSFND + 1;
@WRITE(' Solution# ', NSFND, ' (Iter: ', ITER,')', @NEWLINE( 1));
@DIVERT( '\temp\AltOpt.txt', 'A'); ! Send output to file, Appending to existing contents;
@WRITE(' Solution# ', NSFND, ' (Iter: ', ITER,')', @NEWLINE( 1));
@FOR( COL( j) :
XSV( NSFND, j) = X( j); ! Save it;
@IFC( X( j) #GT# 0: ! Write out the nonzeroes;
TEMP = X( j);
@WRITE(' ', COL( j),' ', TEMP, @NEWLINE( 1));
);
);
@WRITE(@NEWLINE( 1));
@DIVERT(); ! Revert to output to the screen;
);
);
@WRITE(' Completed ', ITER, ' iterations', @NEWLINE( 1));
@DIVERT( '\temp\AltOpt.txt', 'A'); ! Send output to file, Appending to existing contents;
@WRITE( @NEWLINE( 1),' The average interior solution is:', @NEWLINE( 1));
@FOR( COL( j) :
TEMP = @SUM( SLIST( i) | i #LE# NSFND: XSV( i, j))/ NSFND; ! Save it;
@WRITE(' ', COL( j),' ', TEMP, @NEWLINE( 1));
);
@DIVERT(); ! Revert to output to the screen;
ENDCALC
|
