vcphub / cspsol

Try column generation techniques on cutting stock problem
GNU General Public License v3.0
0 stars 0 forks source link

Strange optimization result ? #1

Closed GoogleCodeExporter closed 8 years ago

GoogleCodeExporter commented 8 years ago
I just noticed strange results when given quite simple data, shown below.

Reading order data from file o1.txt
Total orders read from file = 2
2000 4
1000 4
Add order demand constraints. Total rows, cols = 2, 0
Added initial patterns. Total rows, cols = 2, 2
Node    1: new patterns =    2 Obj Func Value = 2 Branch.
Node    2: new patterns =    0 Obj Func Value = 2.33333 Branch.
Node    3: new patterns =    0 Obj Func Value = 5.33333 Branch.
Node    4: new patterns =    0 Obj Func Value = 3 INTEGER ***
Node    5: new patterns =    0 Obj Func Value = 3 INTEGER ***
Node    6: new patterns =    0 Obj Func Value = 6 INTEGER ***
Node    7: new patterns =    0 Obj Func Value = 6 INTEGER ***

Branch and bound tree exhausted.

 # Solution Report # 

Best integer obj. func. value = 3
Pattern count =    1:  1000 x  6, 
Pattern count =    2:  2000 x  3, 

# Total runtime = 0 Secs

I would expect it to give results of

Pattern count = 2: 2000 x 2, 1000 x 2

This would be optimal solution, no waste and less material needed. The
lenght of stock in this example is 6000. Attached is o1.txt file.

I'm using 0.3 version of this software. Is there something I could do to
debug this more ?

Original issue reported on code.google.com by henry.pa...@gmail.com on 7 Aug 2008 at 6:54

Attachments:

GoogleCodeExporter commented 8 years ago
Thank you for reporting the problem. I am trying to find out what is causing 
this
problem, will update you soon.

Original comment by vijay.pa...@gmail.com on 11 Aug 2008 at 11:08

GoogleCodeExporter commented 8 years ago
I was able to reproduce the issue on my machine, using version 0.3. I found the 
cause
behind the issue. Basically the pattern (2000 x 2, 1000 x 2) is NOT being 
generated
by the subproblem.

At one stage the subproblem becomes:
--------------------------------------------------------------
\* Problem: SubProb *\

Maximize
 obj: + 0.333333333333333 x_1 + 0.166666666666667 x_2

Subject To
 PatternWidthConstraint: + 1000 x_2 + 2000 x_1 <= 6000

Generals
 x_1
 x_2

End
--------------------------------------------------------------
The objective function values are copies from dual values in the master problem.
Unfortunately the dual values are approximate. The real coefficients are: 

Maximize
 obj: + 1/3 x_1 + 1/6 x_2
 PatternWidthConstraint: + 1000 x_2 + 2000 x_1 <= 6000

The coefficients 1/3 and 1/6 are irrational numbers and cannot be represented
exactly. So values are approximated and integer solution is biased towards 
variable
x_2 (generating pattern 1000 x 6). 

If we could represent the coefficients exactly, then there are multiple optimal
integer solutions to the sub-problem.
s1. x_1 = 3, x_2 = 0, obj. fun. value = 1.0
s2. x_1 = 0, x_2 = 6, obj. fun. value = 1.0
s3. x_1 = 2, x_2 = 2, obj. fun. value = 1.0

In this case, if somehow program chooses solution s3, then pattern (2000 x 2, 
1000 x
2) will be generated and added to the master problem. Branch and bound 
algorithm will
then  select this pattern in the final optimal integer solution.

To summarize, to fix this issue, we need to:

p1. Represent obj. fun. coefficients of subproblem exactly. Since they are 
copied
from dual values of master problems, the dual values themselves need to be
represented exactly.

p2. Devise some heuristic to select one of 3 multiple optimal integer solution. 
This
can be done by checking is pattern generated is duplicate.

Right now I do not know how solve problem p1. I will give it some more time. 
Let me
know if above analysis makes sense.

Original comment by vijay.pa...@gmail.com on 11 Aug 2008 at 1:32

GoogleCodeExporter commented 8 years ago
"The coefficients 1/3 and 1/6 are irrational numbers and cannot be represented
exactly. So values are approximated and integer solution is biased towards 
variable
x_2 (generating pattern 1000 x 6). "

I should have used term "repeating decimal" instead of "irrational". 
http://mathworld.wolfram.com/RepeatingDecimal.html

Original comment by vijay.pa...@gmail.com on 12 Aug 2008 at 11:39

GoogleCodeExporter commented 8 years ago
[deleted comment]
GoogleCodeExporter commented 8 years ago
Same result (function value = 3) is returned for input
6500
2000 4
1000 4
Here no rounding errors are to be expected.

Original comment by carl.dui...@gmx.de on 12 Aug 2008 at 11:03

GoogleCodeExporter commented 8 years ago
6500
2000 4
1000 4

For above input, I am able to reproduce the reported behavior with version 0.3. 
The
optimal number patterns is 3.0, which is incorrect. The correct value of 
optimal obj.
func. is 2.0

Original comment by vijay.pa...@gmail.com on 13 Aug 2008 at 11:07

GoogleCodeExporter commented 8 years ago
I have come up with an explanation for the issue and a workaround.
Following are the details.

# Background information:
Program cspsol has two MIPs, a master problem and a subproblem. Master problem 
is
solved by custom branch and bound (B&B) algorithm. Subproblem is solved by 
using GLPK
API glp_intopt.

While solving root node LP of master problem's B&B, subproblem MIP is solved 
many
time. Each subproblem solution (optimal integer) produces a pattern which is 
then
added to the master problem as a variable/column. This is done till subproblem
produces a duplicate pattern (i.e a pattern already existing in the master 
problem).
Note that subproblem MIP might have multiple optimal integer solutions, but 
only one
solution is returned to the master problem.

Example:

Subproblem MIP 1 -> Optimal Int. Solution -> Pattern 1 -> Accept Pattern 1 
Subproblem MIP 2 -> Optimal Int. Solution -> Pattern 2 -> Accept Pattern 2
Subproblem MIP 3 -> Optimal Int. Solution -> Pattern 1 -> Pattern 1 is 
Duplicate Stop.

So far so good.

# Issue
For certain input to cspsol, it was observed that some patterns, which are 
necessary
for obtaining optimal solution to master problem are never produced by the 
subproblem.

For example: Consider following subproblem MIP.
--------------------------------------------------------------
Maximize
 obj: + 0.333333333333333 x_1 + 0.166666666666667 x_2

Subject To
 PatternWidthConstraint: + 1000 x_2 + 2000 x_1 <= 6000

Generals
 x_1
 x_2

End
--------------------------------------------------------------
It has following 4 (at least) multiple optimal integer solutions.

s1. x_1 = 3, x_2 = 0, obj. fun. value = 1.0
s2. x_1 = 0, x_2 = 6, obj. fun. value = 1.0
s3. x_1 = 2, x_2 = 2, obj. fun. value = 1.0
s4. x_1 = 1, x_2 = 4, obj. fun. value = 1.0

Unfortunately, we are using GLPK to return only one of above solution, so here 
is
what happens.

Subproblem MIP 1 -> Optimal Int. Solution -> Pattern 1 -> Accept Pattern 1 
Subproblem MIP 2 -> Optimal Int. Solution -> Pattern 2 -> Accept Pattern 2
Subproblem MIP 3 -> Multiple Optimal Int. Solutions -> Pattern 1, 2, 3, 4 -> 
Pattern
1 is Duplicate Stop.

As you can see, patterns 3 and 4 are never used in master problem. If we could 
find a
way to add these patterns to master problem, issue will be resolved.

# Solution 
To resolve the issue, we need to make changes to part the code where pattern is
checked for duplicates. In case of duplicate pattern, we should check is another
optimal integer solution exists. If yes, use that solution to create new 
pattern and
add it to master problem. 

Right now I do not know how to use GLPK API to access multiple optimal integer
solutions, so I came up with following work around.

# Workaround

Workaround Part 1: Somehow get another integer optimal solution from GLPK. 
Assume
atleast once variable branching is done.

In GLPK branch & bound algorithm (glp_intopt), a LP node is NOT explored 
further if
it's optimal obj. function value (bound) is less than best integer solution 
(obtained
so far) + eps. This is when obj function dir is maximization (which is the case 
for
subproblem). Following is relevant code.

----------
Function: static int is_branch_hopeful(glp_tree *tree, int p)
File: glpios03.c

eps = tree->parm->tol_obj * (1.0 + fabs(tree->mip->mip_obj));
if (bound <= tree->mip->mip_obj + eps) ret = 0;
----------

The logic behind this code is that NOT to waste time finding an integer solution
which is better than the existing integer solution by only very small amount 
(eps).
This is fine as a default behavior of GLPK. Besides parameter 
"tree->parm->tol_obj"
can be used defined.

Now remember, to solve the issue, we want GLPK to explore solution which is as 
good
as existing integer solution. In fact we even want solution which is worse by 
small
amount, because our aim is to find alternate optimal integer solution. By 
alteranate
I mean what GLPK (glp_intopt) would have found by default. This can be achieved 
by
setting "tree->parm->tol_obj" to a small negative value.

    Call glp_intopt. 
    Use default optimal integer solution.
    Set "tree->parm->tol_obj" to a small negative value.
    Call glp_intopt. 
    Use alternate optimal integer solution.

Workaround Part 2: 
Above workaround alone will NOT work in case there is no branching. This can 
happen
if root node lp of BB tree has a integer solution. So we need another idea. As
suggested by Xypron (on GLPK mailing list) we could ensure branching by 
adjusting
maximum pattern width. 

For example: If user input for max. pattern width is 6000, and all demand 
widths are
integer then we could replace 6000 with 6000.5. More thought needed.

# Code Changes

By default GLPK does not allow 0 or -ve value for "tree->parm->tol_obj". We 
need to
comment out following code.

Function: int glp_intopt(glp_prob *mip, const glp_iocp *parm)
File: glpapi08.c
if (!(0.0 < parm->tol_obj && parm->tol_obj < 1.0))
    xfault("glp_intopt: tol_obj = %g; invalid parameter\n",
            parm->tol_obj);

Next release of cspsol will contain code related the workaround. Please note 
however
that we need to make change to GLPK source code also. I have done some testing 
and
able to obtain desired result. Also note that this is a workaround to fix this 
issue,
ideal and general solution is to access all optimal integer solutions from GLPK.

Comments are welcome. Attached file is a C program to demonstrate the 
workaround and
subprob lp file.

Original comment by vijay.pa...@gmail.com on 13 Aug 2008 at 12:24

Attachments:

GoogleCodeExporter commented 8 years ago

Released version 0.4, which contains the workaround to temporarily fix the 
issue. To
download source code tarball visit download section. For more details please see
release notes.

http://code.google.com/p/cspsol/source/browse/trunk/ReleaseNotes.txt

Original comment by vijay.pa...@gmail.com on 14 Aug 2008 at 5:57

GoogleCodeExporter commented 8 years ago
Segment of communication with Xypron.

"What if we use DP to solve this 0-N (not 0-1) knapsack problem? Would be 
slower than
solving MIP using GLPK?. I am writing a prototype, will let you know how it 
goes. It
would be easier to access multiple optimal solution in DP. If this works, then 
we
will have two different/separate options to solve the subproblem.
1. Using GLPK glp_intopt.
2. Using DP."

I completed prototype implementation of dynamic programming solution to knapsack
problem, which could be used in cspsol, as long term solution to original issue.
See attached source code file. Hopefully it is not very inefficient.

$ g++ -o dp -Wall subprob.cpp 
$ dp test.in 
6000
2000 0.333333
1000 0.166667

# Optimal integer solutions by DP : 

1000 x 6,  obj. val = 1

1000 x 4, 2000 x 1,  obj. val = 1

1000 x 2, 2000 x 2,  obj. val = 1

2000 x 3,  obj. val = 1

-----------------

Original comment by vijay.pa...@gmail.com on 14 Aug 2008 at 6:34

Attachments:

GoogleCodeExporter commented 8 years ago
I just verified that release 0.4 works ok with both of above cases. Using glpk 
4.25,
with modified sources as above Mr. Vijay stated in comment #7. Next I'll run 
more
extensive testing with actual production data.

Original comment by henry.pa...@gmail.com on 15 Aug 2008 at 9:59

GoogleCodeExporter commented 8 years ago

Second attempt to write dynamic programming (DP) solution. First attempt
(subprob.cpp) is inefficient code. 

Attached is new code which should be faster. Program should print multiple 
optimal
integer solutions. TODO: Need to add logic for avoiding duplicate solutions. I 
was
lacking concentration towards the end, so there may be bugs.

Once ready, DP related can be used to solve subproblem in cspsol, which is the 
long
term solution for this issue. However this approach requires a restriction that 
all
order widths and max. order width should be INTEGER values.

vijay@vijay-hq:~/algo/dp$ g++ -o multdp -g -Wall multdp.cpp 

vijay@vijay-hq:~/algo/dp$ ./multdp test.in 
6000
2000 0.333333
1000 0.166667

# Optimal integer solutions by DP : 

# Solution set. obj. value = 1
S0 ow0:3, 

S1 ow0:2, ow1:2, 

S2 ow0:2, ow1:2, 

S3 ow0:2, ow1:2, 

S4 ow0:1, ow1:4, 

S5 ow0:2, ow1:2, 

S6 ow0:2, ow1:2, 

S7 ow0:1, ow1:4, 

S8 ow0:2, ow1:2, 

S9 ow0:1, ow1:4, 

S10 ow0:1, ow1:4, 

S11 ow0:1, ow1:4, 

S12 ow1:6, 

For test case provided by Xypron, relevant portion of the output is: 

# Optimal integer solutions by DP : 

# Solution set. obj. value = 3.40638
S0 ow14:1, ow78:3, 

S1 ow14:1, ow78:3, 

S2 ow14:1, ow78:3, 

S3 ow14:1, ow78:3, 

Original comment by vijay.pa...@gmail.com on 16 Aug 2008 at 6:36

Attachments:

GoogleCodeExporter commented 8 years ago

Just released version 0.5. Please see release notes.
http://cspsol.googlecode.com/svn/trunk/ReleaseNotes.txt

You can download the source code tarball from Downloads section.
Soon I will close this issue with status "Fixed".

Original comment by vijay.pa...@gmail.com on 18 Aug 2008 at 6:00

GoogleCodeExporter commented 8 years ago
I have certain doubt on the solutions obtained with cspsol V0.4 and cspsol V0.5.

The configuration I used is the following
    CygWin
    GLPK V-4.30
    cspsol V0.4 and cspsol V0.5

/*******************************************************************************
*************/
1) Test with cspsol V0.4
I test the program with a modified version of the example orders1.txt you 
proposed.

jacquenot@PC-GJ-10 % cat ./data/orders1.txt
100
10 4
20 7
30 8
40 6
50 4
60 6
70 7
80 3
90 10
100 2
---- /cygdrive/d/cspsol-0.4 -----

The modification I brought is the following:  I double the length of the stock 
(L=200).
As a consequence, if we know the solution of the problem with stocks of length 
100,
we have an upper bound to the new problem.
This upper bound of the new problem is the upper nearest integer of half the 
solution
of the 1st problem.

Example
The solution of problem orders1.txt where stocks are of length 100 is 32. 
Therefore,
the upper bound of the problem with stocks of length 200 is 16.
However, the solution obtained with cspsol V0.4 is 17!!!
A quick analysis reveals that all demands are satisfied (i.e. all constraints 
are
satisfied) but that more than the demands is provided by the cutting.
Maybe not all patterns were generated?
Do you report the same solution?

See the log below

jacquenot@PC-GJ-49 % cat ./data/orders1_200.txt
200
10 4
20 7
30 8
40 6
50 4
60 6
70 7
80 3
90 10
100 2
---- /cygdrive/d/cspsol-0.4 -----
jacquenot@PC-GJ-50 % ./src/cspsol.exe --wa --data ./data/orders1_200.txt
Reading order data from file ./data/orders1_200.txt
Total orders read from file = 10
10 4
20 7
30 8
40 6
50 4
60 6
70 7
80 3
90 10
100 2
Add order demand constraints. Total rows, cols = 10, 0
Added initial patterns. Total rows, cols = 10, 10
Node    1: new patterns =   12 Obj Func Value = 15.25 Branch.
Node    2: new patterns =    0 Obj Func Value = 16.0667 Branch.
Node    3: new patterns =    0 Obj Func Value = 15.6889 Branch.
Node    4: new patterns =    0 Obj Func Value = 16.7333 Branch.
Node    5: new patterns =    0 Obj Func Value = 16.0833 Branch.
Node    6: new patterns =    0 Obj Func Value = 15.88 Branch.
Node    7: new patterns =    0 Obj Func Value = 17.1667 Branch.
Node    8: new patterns =    0 Obj Func Value = 16.8333 Branch.
.
.
.
Node  156: new patterns =    0 Obj Func Value = 18 INTEGER ***
Node  157: new patterns =    0 Obj Func Value = 18 INTEGER ***
Node  158: new patterns =    0 Obj Func Value = 18 INTEGER ***
Node  159: new patterns =    0 Obj Func Value = 18 INTEGER ***
Node  160: new patterns =    0 Obj Func Value = 18 INTEGER ***
Node  161: new patterns =    0 LP worse than integer incumbent 17.16 >= 17. 
Fathom node. 

Branch and bound tree exhausted.

 # Solution Report # 

Best integer obj. func. value = 17
Pattern count =    1:    40 x  5, 
Pattern count =    1:    50 x  4, 
Pattern count =    1:    20 x  1,    60 x  3, 
Pattern count =    3:    60 x  1,    70 x  2, 
Pattern count =    2:    40 x  1,    80 x  2, 
Pattern count =    5:    20 x  1,    90 x  2, 
Pattern count =    1:   100 x  2, 
Pattern count =    2:    10 x  2,    30 x  6, 
Pattern count =    1:    10 x  1,    20 x  4,    40 x  1,    70 x  1, 

# Total runtime = 0 Secs
---- /cygdrive/d/cspsol-0.4 -----

All demand is satisfied (i.e. all constraints are satisfied) and even more
-  5 components of length 10 are cut whereas only 4 were needed
- 10 components of length 20 are cut whereas only 7 were needed
- 12 components of length 30 are cut whereas only 8 were needed
-  8 components of length 40 are cut whereas only 6 were needed
-  4 components of length 80 are cut whereas only 3 were needed

Maybe it works fine with the new version (0.5), but I didn't test it.

/*******************************************************************************
*************/
2) Test with cspsol V0.4
With version 0.5, I couldn't get the optimal results for the o1 problem
It seems to me that the program didn't generate all possible patterns or return 
the
upper bound of the problem
I obtain the same results with DP and GLPK

jacquenot@PC-GJ-42 % cat ./data/o1.txt
6000
2000 4
1000 4

With the dynamic programming resolution modell

jacquenot@PC-GJ-42 % ./src/cspsol.exe -d ./data/o1.txt 
Reading order data from file ./data/o1.txt
Total orders read from file = 2
2000 4
1000 4
Add order demand constraints. Total rows, cols = 2, 0
Added initial patterns. Total rows, cols = 2, 2
Node    1: new patterns =    1 Obj Func Value = 8 INTEGER ***

Branch and bound tree exhausted.

 # Solution Report # 

Best integer obj. func. value = 8
Pattern count =    4:  2000 x  1, 
Pattern count =    4:  1000 x  1, 

# Total runtime = 0 Secs
---- /cygdrive/d/cspsol-0.5 -----

With GLPK

jacquenot@PC-GJ-43 % ./src/cspsol.exe --wa --si -d ./data/o1.txt
IMP: Must use patched GLPK lib to allow -ve tol_obj.

Reading order data from file ./data/o1.txt
Total orders read from file = 2
2000 4
1000 4
Add order demand constraints. Total rows, cols = 2, 0
Added initial patterns. Total rows, cols = 2, 2
Node    1: new patterns =    1 Obj Func Value = 8 INTEGER ***

Branch and bound tree exhausted.

 # Solution Report # 

Best integer obj. func. value = 8
Pattern count =    4:  2000 x  1, 
Pattern count =    4:  1000 x  1, 

# Total runtime = 0 Secs
---- /cygdrive/d/cspsol-0.5 -----

However, my thanks to the author
guillaume dot jacquenot at gmail dot com

Original comment by guillaum...@gmail.com on 19 Aug 2008 at 10:13

GoogleCodeExporter commented 8 years ago
Thank you for detailed issue report.

Part 1) Test with cspsol V0.4

"The solution of problem orders1.txt where stocks are of length 100 is 32. 
Therefore,
the upper bound of the problem with stocks of length 200 is 16.
However, the solution obtained with cspsol V0.4 is 17!!!"

Indeed version 0.4 (even with workaround) will produce sub-optimal solutions. 
This is
because it is incapable of generating ALL optimal solutions to the subproblem.

Version 0.5 should give optimal solution 16. I just tried it, however it was 
taking
too long, which got be thinking about how to reduce the runtime. I have a good 
idea
to reduce the runtime.

File: main.cpp
Line 101:
if(node->get_opt_obj_val() >= BBNode::get_best_int_obj_val()) 

Replace above line with:
if(node->get_opt_obj_val() >= BBNode::get_best_int_obj_val()-1.0) 

This is based on the observation that objective function value of master 
problem, for
any integer solution will be always a integer value. So if best_int_obj_value 
is 16,
then we are not interested in exploring any nodes whose lp opt obj. value is 
greater
than 15. I made a test run and got optimal solution 16 in 80 seconds on my 
machine.

Summary: You could make above change to version 0.5 source code and it should 
solve
your issue. This change will be available in next version 0.6.

Part 2) 

"With version 0.5, I couldn't get the optimal results for the o1 problem"

This is unlikely. Please download the source code version 0.5 again and build 
the
exe. It should give correct solution for input o1.txt.

Original comment by vijay.pa...@gmail.com on 19 Aug 2008 at 11:34

GoogleCodeExporter commented 8 years ago
[deleted comment]
GoogleCodeExporter commented 8 years ago
I am sorry to insist, but I cannot find the correct solution of problem o1 with 
the
version 0.5.

o1.txt 

6000
2000 4
1000 4

If I compare the lp files generated during the execution of Version 0.4 and 
Version
0.5, I do not obtain the same files.
With version 0.4, I obtain the correct solution, which is 2 but not the version 
0.5

Log:

    Version 0.4
    ./src/cspsol.exe --wa --data ./data/o1.txt
    best.lp
    \* Problem: MasterCutStock *\

    Minimize
     obj: + x_1 + x_2 + x_3 + x_4 + x_5

    Subject To
     Width1: + 2 x_5 + 3 x_4 + x_1 >= 4
     Width2: + 2 x_5 + 6 x_3 + x_2 >= 4

    Bounds
     x_3 = 0

    Generals
     x_1
     x_2
     x_3
     x_4
     x_5

    End

    Version 0.5
    ./src/cspsol.exe -d ./data/o1.txt
    or
    ./src/cspsol.exe --wa --si -d ./data/o1.txt
    (Give the same result)
    best.lp
    \* Problem: MasterCutStock *\

    Minimize
     obj: + x_1 + x_2 + x_3

    Subject To
     Width1: + x_1 >= 4
     Width2: + 6 x_3 + x_2 >= 4

    Generals
     x_1
     x_2
     x_3

    End

This explains why the solutions are different from version 0.4 to version 0.5
However I do not know where is the error

Original comment by guillaum...@gmail.com on 19 Aug 2008 at 1:03

GoogleCodeExporter commented 8 years ago

All problems mentioned in this issue are fixed in version 0.7. Closing this 
issue.

Original comment by vijay.pa...@gmail.com on 22 Aug 2008 at 3:29