Hello, Please forgive me if this is not in the correct format. If there is a better way to submit this report please advise me.
A research group I work with has been using GLPK and CPLEX for the past four years to solve queueing network problems. Most of our work occurs during the summer. Recently I upgraded the version of GLPK we were using from 4.29 to 4.38 but noticed that some programs written in MathProg that used to work correctly (i.e. GLPK gave the same results as CPLEX) no longer did so. After a bit of work, I figured out: Version 4.29 and Verison 4.30 work okay, Version 4.31 and Version 4.38 do not. I assume that the versions between 4.31 and 4.38 also fail, but I've not tested them. Attached is a MathProg program that illustrates the problem. I run with "glpsol --math seriesexp3bad.mod" The output from Version 4.30 isReading model section from seriesexp3bad.mod... Reading data section from seriesexp3bad.mod... 128 lines were read Display statement at line 31 lambda = 1 Lambda = 3.75 B = 0.833333333333333 Display statement at line 32 i = 1 c[1] = 1 mu[1] = 1.5 i = 2 c[2] = 2 mu[2] = 1.25 Generating cost_per_hour... Generating Type_1... Generating Type10... Generating Type00... Generating Type1a... Generating Type1b... Generating Type2a1... Generating Type2a2... Generating Type2b... Generating Type4a... Generating Quadratic... Generating EXP2... Model has been successfully generated glp_simplex: original LP has 131 rows, 9 columns, 1074 non-zeros glp_simplex: presolved LP has 128 rows, 7 columns, 826 non-zeros lpx_adv_basis: size of triangular part = 128 * 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0) * 17: obj = 6.612454212e+00 infeas = 3.264e-14 (0) OPTIMAL SOLUTION FOUND Time used: 0.0 secs Memory used: 0.4 Mb (378322 bytes) The output from Version 4.38 isGLPSOL: GLPK LP/MIP Solver 4.38 Reading model section from seriesexp3bad.mod... Reading data section from seriesexp3bad.mod... 128 lines were read Display statement at line 31 lambda = 1 Lambda = 3.75 B = 0.833333333333333 Display statement at line 32 i = 1 c[1] = 1 mu[1] = 1.5 i = 2 c[2] = 2 mu[2] = 1.25 Generating cost_per_hour... Generating Type_1... Generating Type10... Generating Type00... Generating Type1a... Generating Type1b... Generating Type2a1... Generating Type2a2... Generating Type2b... Generating Type4a... Generating Quadratic... Generating EXP2... Model has been successfully generated glp_simplex: original LP has 131 rows, 9 columns, 1074 non-zeros glp_simplex: presolved LP has 128 rows, 7 columns, 826 non-zeros Scaling... A: min|aij| = 1.355e-20 max|aij| = 7.438e+01 ratio = 5.488e+21 GM: min|aij| = 6.755e-06 max|aij| = 1.480e+05 ratio = 2.192e+10 EQ: min|aij| = 4.563e-11 max|aij| = 1.000e+00 ratio = 2.192e+10 Crashing... Size of triangular part = 128 * 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0) Warning: numerical instability (primal simplex, phase II) 6: obj = 5.142018849e+00 infeas = 7.140e-02 (0) * 30: obj = 5.172971112e+00 infeas = 4.714e-09 (0) Warning: numerical instability (primal simplex, phase II) 31: obj = 5.085248263e+00 infeas = 2.024e-01 (0) * 60: obj = 5.138363589e+00 infeas = 3.042e-09 (0) * 62: obj = 1.308182055e+01 infeas = 9.367e-08 (0) OPTIMAL SOLUTION FOUND Time used: 0.0 secs Memory used: 0.4 Mb (459402 bytes) Jonathan R. Senning Chair, Department of Mathematics and Computer Science Gordon College, 255 Grapevine Road, Wenham MA 01984 Voice: 978-867-4376, FAX: 978-867-4666 Email: [email protected] Web: http://www.math-cs.gordon.edu/~senningHello,
Please forgive me if this is not in the correct format. If there is a better way to submit this report please advise me.
A research group I work with has been using GLPK and CPLEX for the past four years to solve queueing network problems. Most of our work occurs during the summer. Recently I upgraded the version of GLPK we were using from 4.29 to 4.38 but noticed that some programs written in MathProg that used to work correctly (i.e. GLPK gave the same results as CPLEX) no longer did so.
After a bit of work, I figured out:
Version 4.29 and Verison 4.30 work okay, Version 4.31 and Version 4.38 do not. I assume that the versions between 4.31 and 4.38 also fail, but I've not tested them.
Attached is a MathProg program that illustrates the problem. I run with "glpsol --math seriesexp3bad.mod"
The output from Version 4.30 is
Reading model section from seriesexp3bad.mod... Reading data section from seriesexp3bad.mod... 128 lines were read Display statement at line 31 lambda = 1 Lambda = 3.75 B = 0.833333333333333 Display statement at line 32 i = 1 c[1] = 1 mu[1] = 1.5 i = 2 c[2] = 2 mu[2] = 1.25 Generating cost_per_hour... Generating Type_1... Generating Type10... Generating Type00... Generating Type1a... Generating Type1b... Generating Type2a1... Generating Type2a2... Generating Type2b... Generating Type4a... Generating Quadratic... Generating EXP2... Model has been successfully generated glp_simplex: original LP has 131 rows, 9 columns, 1074 non-zeros glp_simplex: presolved LP has 128 rows, 7 columns, 826 non-zeros lpx_adv_basis: size of triangular part = 128 * 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0) * 17: obj = 6.612454212e+00 infeas = 3.264e-14 (0) OPTIMAL SOLUTION FOUND Time used: 0.0 secs Memory used: 0.4 Mb (378322 bytes)
The output from Version 4.38 is
GLPSOL: GLPK LP/MIP Solver 4.38 Reading model section from seriesexp3bad.mod... Reading data section from seriesexp3bad.mod... 128 lines were read Display statement at line 31 lambda = 1 Lambda = 3.75 B = 0.833333333333333 Display statement at line 32 i = 1 c[1] = 1 mu[1] = 1.5 i = 2 c[2] = 2 mu[2] = 1.25 Generating cost_per_hour... Generating Type_1... Generating Type10... Generating Type00... Generating Type1a... Generating Type1b... Generating Type2a1... Generating Type2a2... Generating Type2b... Generating Type4a... Generating Quadratic... Generating EXP2... Model has been successfully generated glp_simplex: original LP has 131 rows, 9 columns, 1074 non-zeros glp_simplex: presolved LP has 128 rows, 7 columns, 826 non-zeros Scaling... A: min|aij| = 1.355e-20 max|aij| = 7.438e+01 ratio = 5.488e+21 GM: min|aij| = 6.755e-06 max|aij| = 1.480e+05 ratio = 2.192e+10 EQ: min|aij| = 4.563e-11 max|aij| = 1.000e+00 ratio = 2.192e+10 Crashing... Size of triangular part = 128 * 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0) Warning: numerical instability (primal simplex, phase II) 6: obj = 5.142018849e+00 infeas = 7.140e-02 (0) * 30: obj = 5.172971112e+00 infeas = 4.714e-09 (0) Warning: numerical instability (primal simplex, phase II) 31: obj = 5.085248263e+00 infeas = 2.024e-01 (0) * 60: obj = 5.138363589e+00 infeas = 3.042e-09 (0) * 62: obj = 1.308182055e+01 infeas = 9.367e-08 (0) OPTIMAL SOLUTION FOUND Time used: 0.0 secs Memory used: 0.4 Mb (459402 bytes)
|
# ALP for Series Queue, EXP3(N) h approximation
# Adapted from M. LeClair's "EXP1" ALP
# by M. Veatch 10/18/04
#
# h(x)=(1/2)x'Qx+p'x+r1*B^x2+r2*x1*B^x2+r3*x2*B^x2
# where B = mu2/mu1
# Reduced over x1, truncated x2 <= N
param k := 2; # of classes
param lambda; #the arrival rate into the system
param mu{i in 1..k}; #the departure rate for the different machines
param c{i in 1..k}; #cost of being in the different states
param N; #truncation x2 <= N
param Lambda:= lambda + mu[1] + mu[2]; #total Event Rate
param B:= mu[2]/mu[1];
#
var J; #Average cost per stage
var Q{i in 1..k, j in 1..k}; #varibles in the upper triangular matrix
var p{i in 1..k}; # the linear x term variable
var r{i in 1..3}; # the exponential x term variables
display lambda, Lambda, B;
display {i in 1..2}:i, c[i], mu[i];
#---------------------------------------------------------------------
#Objective
#---------------------------------------------------------------------
maximize cost_per_hour: Lambda*J;
#---------------------------------------------------------------------
#Constraints
#---------------------------------------------------------------------
s.t. Type_1 #u=(1,1) (x1=1) and u=(0,1) (x1=0)
{x1 in 0..1, x2 in 1..N}: #u1=x1, u2=1
Lambda*J<= c[1]*x1+c[2]*x2
+lambda*(
.5*(x1+1)^2*Q[1,1]+(x1+1)*x2*Q[1,2]+.5*x2^2*Q[2,2]+p[1]*(x1+1)+p[2]*x2
+r[1]*B^x2 + r[2]*(x1+1)*B^x2 +r[3]*x2*B^x2
)
#the above is the lamda*h(x+e1) term
+mu[1]*(
.5*0^2*Q[1,1]+0*(x2+x1)*Q[1,2]+.5*(x2+x1)^2*Q[2,2]
+p[1]*0+p[2]*(x2+x1)
+r[1]*B^(x2+x1) + r[2]*0*B^(x2+x1) + r[3]*(x2+x1)*B^(x2+x1)
)
#the above is the mu1*h(x-e1+e2) term, or h(x) if x1=0
+mu[2]*(
.5*x1^2*Q[1,1]+x1*(x2-1)*Q[1,2]+.5*(x2-1)^2*Q[2,2]+p[1]*x1+p[2]*(x2-1)
+r[1]*B^(x2-1) + r[2]*x1*B^(x2-1) + r[3]*(x2-1)*B^(x2-1)
)
#the above term is the mu2*h(x-e2) term
-Lambda*(
.5*x1^2*Q[1,1]+x1*x2*Q[1,2]+.5*x2^2*Q[2,2]+p[1]*x1+p[2]*x2
+r[1]*B^x2 + r[2]*x1*B^x2 + r[3]*x2*B^x2
);
# the above is the -Lambda*h(x) term
s.t. Type10 #u=(1,0)
{x1 in 1..1, x2 in 0..0}: #u1=x1, u2=1
Lambda*J<= c[1]*x1+c[2]*x2
+lambda*(
.5*(x1+1)^2*Q[1,1]+(x1+1)*x2*Q[1,2]+.5*x2^2*Q[2,2]+p[1]*(x1+1)+p[2]*x2
+r[1]*B^x2 + r[2]*(x1+1)*B^x2 +r[3]*x2*B^x2
)
#the above is the lamda*h(x+e1) term
+mu[1]*(
.5*0^2*Q[1,1]+0*(x2+x1)*Q[1,2]+.5*(x2+x1)^2*Q[2,2]
+p[1]*0+p[2]*(x2+x1)
+r[1]*B^(x2+x1) + r[2]*0*B^(x2+x1) + r[3]*(x2+x1)*B^(x2+x1)
)
#the above is the mu1*h(x-e1+e2) term, or h(x) if x1=0
-(lambda + mu[1])*(
.5*x1^2*Q[1,1]+x1*x2*Q[1,2]+.5*x2^2*Q[2,2]+p[1]*x1+p[2]*x2
+r[1]*B^x2 + r[2]*x1*B^x2 + r[3]*x2*B^x2
);
# the above is the -(lambda + mu[1])*h(x) term
s.t. Type00: Lambda*J <= lambda*(.5*Q[1,1]+p[1]+r[2]);
s.t. Type1a: c[1]+(lambda-mu[1])*Q[1,1]+(mu[1]-mu[2])*Q[1,2] >= 0;
s.t. Type1b: c[2]+(lambda-mu[1])*Q[1,2]+(mu[1]-mu[2])*Q[2,2] >= 0;
s.t. Type2a1: c[1]+lambda*Q[1,1]-mu[2]*Q[1,2]+r[2]*(mu[1]-mu[2])*B >= 0;
s.t. Type2a2: c[1]+lambda*Q[1,1]-mu[2]*Q[1,2] >= 0;
s.t. Type2b: c[2]+lambda*Q[1,2]-mu[2]*Q[2,2] >= 0;
s.t. Type4a: c[1]+(lambda-mu[1])*Q[1,1]+mu[1]*Q[1,2]-r[2]*(mu[1]-mu[2])>=0;
#s.t. Qbound: Q[2,2] >= Q[1,2]; #enforces switching curve#
#s.t. EXP1: r[1] = 0;
s.t. Quadratic: r[2] = 0;
s.t. EXP2: r[3] = 0;
#Quadratic: set all r = 0
#-------------------------------------------------------------------------------
# Data
#-------------------------------------------------------------------------------
data;
param lambda := 1.0;
param mu:= 1 1.5 2 1.25;
param c := 1 1 2 2;
param N := 60;
end;
_______________________________________________ Bug-glpk mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-glpk
