Follow up on http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/cutting-stock-problem-example.html.
>I decided to have a go at the challenge set in this paper - hope you don't mind.
>We receive bars of length 1200 and must cut them into the following
> [length, quantity] pairs:[100,100],[200,2],[300,2],[493,250],[590,2],
> [630,2],[780,2],[930,2]
....
>So the solution equates to:12 * [100,100,100,100,100,100,100,493]2 *
> [100,100,100,100,100,630]2 * [100,100,930]2 * [100,300,780]2 *
>[100,493,590]118 * [200,493,493]
>
>Total bars used = 138
>
>Theoretical min:
>(100*100+200*2+300*2+493*250+590*2+630*2+780*2+930*2)/1200=116.76
>
>...and, finally, the numbers of each length of bar cut:100: 12 * 7 + 2 * 5 + 2 * 2 +
>2 * 1 + 2 * 1 = 102200: 118 * 1 = 118300: 2 * 1 = 2493: 118 * 2 + 2 * 1 + 12 * 1 =
>250590: 2 * 1 = 2630: 2 * 1 = 2780: 2 * 1 = 2930: 2 * 1 = 2
See also: http://people.sabanciuniv.edu/ertekg/papers/2006/alp_et_al_ims2006.pdf.
I can find a better solution using the column generation model. The following input was used:
*-----------------------------------------------------
* Data
*-----------------------------------------------------
scalar r 'raw width' /1200/;
table demanddata(i,*)
width demand
width1 100 100
width2 200 2
width3 300 2
width4 493 250
width5 590 2
width6 630 2
width7 780 2
width8 930 2
;
The result is a solution (not proven optimal) of 131:
---- 174 PARAMETER pat pattern usage
p4 p5 p9 p10 p11 p15 Total
width1 1 2 100
width2 1 2
width3 1 2
width4 2 1 2 250
width5 2 2
width6 1 2
width7 1 2
width8 1 2
Count 75 1 2 2 2 49 131
**** SOLVER STATUS 8 USER INTERRUPT
**** MODEL STATUS 14 NO SOLUTION RETURNED
**** OBJECTIVE VALUE 3885.0000
RESOURCE USAGE, LIMIT 6.618 3600.000
ITERATION COUNT, LIMIT 0 10000000
MOSEK Link May 1, 2008 22.7.1 WIN 3927.4700 VIS x86/MS Windows
M O S E K version 5.0.0.79 (Build date: Jan 29 2008 13:41:45)
Copyright (C) MOSEK ApS, Fruebjergvej 3, Box 16
DK-2100 Copenhagen, Denmark
http://www.mosek.com
No solution returned
set si /i1,i2/;
parameter duedateadj(si) /i1 10, i2 10/;
set k /k1,k2/;
parameter duration(si,k) / i1.k1 1 /;
parameters
minduration(si)
mintime(si) "calculate in one step"
mintime2(si) "calculate in two steps"
;
mintime(si) = duedateadj(si)-smin(k$duration(si,k),duration(si,k));
minduration(si) = smin(k$duration(si,k),duration(si,k));
mintime2(si) = duedateadj(si)-minduration(si);
display mintime, mintime2;---- 17 PARAMETER mintime calculate in one step
i1 9.000, i2 -INF
---- 17 PARAMETER mintime2 calculate in two steps
i1 9.000, i2 10.000
set
j /x,y/
i /i1,i2,i3/
;
parameter data(i,j);
* may 30, may 31, june 1
data(i,'x') = jdate(2008,5,30)+ord(i)-1;
data(i,'y') = ord(i);
* check dates
parameter check(i,*);
check(i,'year') = gyear(data(i,'x'));
check(i,'month') = gmonth(data(i,'x'));
check(i,'day') = gday(data(i,'x'));
option check:0:1:1;
display data,check;
execute_unload 'data',data;
---- 20 PARAMETER data
x y
i1 39597.000 1.000
i2 39598.000 2.000
i3 39599.000 3.000
---- 20 PARAMETER check
year month day
i1 2008 5 30
i2 2008 5 31
i3 2008 6 1
Sub color()
ActiveSheet.UsedRange.Style = "Normal"
Call colorsub(xlNumbers, "Accent1")
Call colorsub(xlTextValues, "Accent2")
Call colorsub(xlErrors, "Accent3")
Call colorsub(xlLogical, "Accent4")
End Sub
Sub colorsub(SpecialCells As Integer, styletype As String)
Dim r As Range
On Error GoTo err
Set r = ActiveSheet.UsedRange.SpecialCells(xlCellTypeFormulas, SpecialCells)
r.Style = styletype
err:
End Sub
S O L V E S U M M A R Y
MODEL TIME OBJECTIVE LL
TYPE NLP DIRECTION MAXIMIZE
SOLVER PATHNLP FROM LINE 160
**** SOLVER STATUS 4 TERMINATED BY SOLVER
**** MODEL STATUS 7 INTERMEDIATE NONOPTIMAL
**** OBJECTIVE VALUE 47.2365
RESOURCE USAGE, LIMIT 35.697 1000.000
ITERATION COUNT, LIMIT 0 100000
EVALUATION ERRORS 0 0
PATH-NLP May 1, 2008 22.7.2 WEX 3906.4799 WEI x86_64/MS Windows
NLP size: 829 rows, 881 cols, 10145 non-zeros, 1.39% dense.
MCP size: 1708 rows/cols, 26626 non-zeros, 0.91% dense.
S O L V E S U M M A R Y
MODEL m OBJECTIVE cost
TYPE MIP DIRECTION MINIMIZE
SOLVER CPLEX FROM LINE 168
**** SOLVER STATUS 1 NORMAL COMPLETION
**** MODEL STATUS 1 OPTIMAL
**** OBJECTIVE VALUE 3712.0000
RESOURCE USAGE, LIMIT 630.895 3600.000
ITERATION COUNT, LIMIT 48023779 10000000
Example output file:
* miptrace file trace.csv: ID = Cplex 11.0.1
* fields are lineNum, seriesID, node, seconds, bestFound, bestBound
1, S, 0, 0.075, na, 4053.06
2, N, 0, 0.084, na, 4053.06
3, N, 0, 0.09, 1740.37, 3852.17
4, N, 0, 0.095, 1740.37, 3804.84
5, N, 0, 0.101, 1740.37, 3798.37
6, N, 0, 0.107, 1740.37, 3781.77
7, N, 0, 0.113, 1740.37, 3781.53
8, N, 0, 0.149, 1826.58, 3779.84
9, N, 100, 0.229, 2769.62, 3327.77
10, N, 200, 0.304, 2891.26, 3135.02
11, N, 300, 0.391, 2891.26, 2910.57
12, E, 320, 0.397, 2891.26, 2891.26
* miptrace file trace.csv closed
See also http://www.ampl.com/BOOKLETS/amplcplex110userguide.pdf (the AMPL link has a number of these features implemented) and http://www.gams.com/dd/docs/solvers/cplex.pdf.
I sent it to GAMS for further investigation. This problem seems only to happen with the 64 bit version. The work around would be to use the 32 bit version.
MOSEK Link May 1, 2008 22.7.2 WEX 3927.4799 WEI x86_64/MS Windows
M O S E K version 5.0.0.79 (Build date: Jan 29 2008 13:47:58)
semicont variable x;
variable y;
equation e;
y.lo=0;
x.fx=1;
e.. y =e= x;
model m /e/;
solve m minimizing y using mip;
ROWS
N obj
E c1
COLUMNS
sc1 c1 -1
x2 obj 1
x2 c1 1
RHS
BOUNDS
FX bnd sc1 1
ENDATA
ROWS
N obj
E c1
COLUMNS
sc1 c1 -1
x2 obj 1
x2 c1 1
RHS
BOUNDS
LO bnd sc1 1
SC bnd sc1 1
ENDATA
S O L V E S U M M A R Y
MODEL thainavy OBJECTIVE obj
TYPE MIP DIRECTION MINIMIZE
SOLVER CPLEX FROM LINE 84
**** SOLVER STATUS 4 TERMINATED BY SOLVER
**** MODEL STATUS 14 NO SOLUTION RETURNED
**** OBJECTIVE VALUE 0.0000
RESOURCE USAGE, LIMIT 0.037 1000.000
ITERATION COUNT, LIMIT 0 5
The SOLVER STATUS is also wrong, this should be 2 ITERATION INTERRUPT. Turns out I have reported this already in january 2006, but it does not seem to get fixed.
Correct values for these settings are often important when running Cplex from a more complex algorithm formulated in GAMS.
One way of achieving this is to add the constraints:
a(i) ≥ 0.001·p(i)
b(i) ≥ 0.001·q(i)
c(i) ≥ 0.001·r(i)
sum(i,p(i)) ≥ 2
sum(i,q(i)) ≥ 2
sum(i,r(i)) ≥ 2
p(i),q(i),r(i) binary variables
p(i)=0 => a(i)=0
x ≤ (1-w)·M where M=x.up
x ≥ 0.001·(1-w)
| Method | Time | Accuracy | Time | Accuracy | Time | Accuracy |
|---|---|---|---|---|---|---|
| n=50 | n=100 | n=200 | ||||
| Nonlinear Equations | 0.148 | 9.6e-8 | 0.917 | 1.2e-4 | 8.991 | 0.621 |
| External Solver | 0.024 | 1.6e-13 | 0.025 | 1.9e-10 | 0.056 | 1.4e-5 |
| GAMS Algorithm | 0.034 | 1.9e-13 | 0.482 | 6.3e-11 | 7.571 | 7.7e-6 |
I am surprised that the GAMS algorithm is competitive. Notice the accuracy issues when using nonlinear equations due to ill-conditioning of the system of non-linear equations (small changes in a(i,j) lead to large changes in L(i,j)).
Here is some info about the external solver:
c:\>cholesky
CHOLESKY: matrix decomposition A=LL^T
Usage
> cholesky gdxin i a gdxout L
where
gdxin : name of gdxfile with matrix
i : name of set used in matrix
a : name of 2 dimensional parameter inside gdxin
gdxout : name of gdxfile for results (factor L)
L : name of 2 dimensional parameter inside gdxout
Calculates the Cholesky decomposition A=LL^T of a symmetric
positive definite matrix A=a(i,j) where i and j are aliased sets.
L will contain the Cholesky factor L(i,j)
C:\>
if(1,
$set name "hello"
else
$set name "world"
);
display "%name%";
scalar s /NA/;
s$(s>1) = 3.14;
display s;
> touch f1.txt f2.txt f3.txtHowever it fails miserably:
> ren f*.txt gg*.txt
A ("somewhat" complicated) solution is presented here.C:\projects\test\test1>touch f1.txt f2.txt f3.txt
C:\projects\test\test1>ren f*.txt gg*.txt
A duplicate file name exists, or the file cannot be found.
A duplicate file name exists, or the file cannot be found.
C:\projects\test\test1>dir
Volume in drive C has no label.
Volume Serial Number is 7563-3993
Directory of C:\projects\test\test1
05/16/2008 08:42 PM <DIR> .
05/16/2008 08:42 PM <DIR> ..
05/16/2008 08:42 PM 0 f2.txt
05/16/2008 08:42 PM 0 f3.txt
05/16/2008 08:42 PM 0 gg.txt
3 File(s) 0 bytes
2 Dir(s) 159,707,815,936 bytes free
C:\projects\test\test1>
>I would like to write the results of a GAMS program into a text file
>whose name should depend on the date and time of the run. How could
>that be done?
>
>Ex: output-05.15.08-09.39.txt
>The above is just a sample and it does not have to be exactly the same.
The suggestion to use myfile /filename_%system.date%_%system.time%/ does not work for several reasons: quotes are needed to handle embedded blanks and the ‘/’ characters in the system date cause problems when used inside a file name. Here is a different approach that actually works (I ran it before suggesting it):
The generated include file will look like:$onecho > createfile.awk
BEGIN { print "file f / \"f_" strftime("%Y-%m-%d_%H.%M.%S") ".txt\"/" }
$offecho
$call "awk -f createfile.awk > ftime.inc"
$include ftime.inc
put f "hello"/;
For more info on AWK see gawk.htmlfile f / "f_2008-05-16_14.28.59.txt"/
Trouble polishing MIP start.
> There is a worked example of how to solve a basic Cutting Stock Problem at
> http://docs.google.com/View?docid=dfkkh8qj_64rrpz86db
> It uses simple methods, and the open source program lp_solve. Briefly,
> the solution shown is:
>
> * create a list of cut combinations (not permutations - there would be
> too many of those to work with in most cases)
>
> * from this list, generate a set of linear programming (simplex) expressions,
> and output them to a file
>
> * this file becomes the input to lp_solve, which calculates an optimal solution
>
> * most of the cut combinations in the solution will have a score of zero.
> For those that have a score greater than zero, look them up on the list
> created in step 1
A clever way to solve these problems is using a Delayed Column Generation approach. Such a model is not too difficult to code in GAMS. See http://www.amsterdamoptimization.com/pdf/colgen.pdf. A nice write up on the cutting stock problem is available from http://en.wikipedia.org/wiki/Cutting_stock_problem.
The last restriction is essentially a big-M constraint:
\[\boxed{\begin{split}\min\{0,x^{lo}\} &\le z \le \max\{0,x^{up}\} \\
x^{lo} \cdot \delta &\le z \le x^{up} \cdot \delta \\
x − x^{up} \cdot (1 − \delta) &\le z \le x − x^{lo} \cdot (1 − \delta)\end{split}}\]
where \(M_1\), \(M_2\) are chosen as tightly as possible while not excluding \(z=0\) when \(\delta=0\).
\[\boxed{x − M_1 \cdot (1 − \delta) \le z \le x + M_2 \cdot (1 − \delta)}\]
The construct \(z=x· \delta\) can be used to model an OR condition: "\(z = 0\) OR \(z = x\)".
\[\boxed{\begin{split}& z \in [0,x^{up}]\\ & z ≤ x^{up} · \delta\\ & z ≤ x\\ & z ≥ x − x^{up} · (1 − \delta)\end{split}}\]
I used GNUPLOT to plot results, e.g.:
Notice that we relaxed \(z\) to be continuous between 0 and 1: \(z\) is automatically integer; a property we can use to reduce the number of discrete variable. Note that modern high performance solvers may reintroduce these variables as binary as they can exploit this in their algorithms.
\[\boxed{\begin{split}0 &≤ z ≤ 1\\ z &≤ x_1\\ z &≤ x_2\\ z &≥ x_1 + x_2 − 1\end{split}}\]
This formulation has fewer equations but has two major disadvantages:
\[\boxed{\begin{split}z &\in \{0,1\}\\ z & ≤ \frac{x_1+x_2}{2}\\z &≥ x_1 + x_2 − 1\end{split}}\]
An application of this reformulation is a QP (Quadratic Program) with binary variables:
\[\boxed{\begin{split} \min\> & x'Qx\\ & Ax=b\\ & x \in \{0,1\} \end{split}}\]
This can be modeled as:
variable y(i,j);
set ij(i,j);
ij(i,j)$( ord(i)<ord(j) ) = yes;
y.lo(ij)=0;
y.up(ij)=1;
obj.. z =e= sum(i, q(i,i)*x(i)) + sum(ij(i,j), (q(i,j)+q(j,i))*y(i,j) );
ymul1(ij(i,j)).. y(i,j) =L= x(i);
ymul2(ij(i,j)).. y(i,j) =L= x(j);
ymul3(ij(i,j)).. y(i,j) =G= x(i)+x(j)-1;
The more general product \(z = \prod_i x_i\) with \(x_i\in \{0,1\}\) can be linearized in the same way:
where \(n\) is the cardinality of the set \(I\). Again if we want we can relax \(z\) to \(z \in [0,1]\).
\[\boxed{\begin{split} &z ≤ x_i\\ & z ≥ \sum_i x_i – (n-1)\\& z \in \{0,1\}\end{split}}\]
$set stochastic 0
$set erslivestockmodel 1
$ifthen %stochastic%==1
$ ifthen %erslivestockmodel%==1
b(i,"Exports",k) = tempbx(k);
b(i,"Season Average Price",k) = tempbp(k);
b(i,"Feed Use",k) = tempbdf(k);
b(i,"Food",k) = tempbdi(k);
cfixed(i,"Exports",k) = 0;
cfixed(i,"Season Average Price",k) = 0;
cfixed(i,"Feed Use",k) = 0;
cfixed(i,"Food",k) = 0;
$ else
$ include grains.inc
$ endif
$endif
This does not work as nesting is not allowed. You don't get an error message either indicating nesting is not allowed. In some cases you may get a possibly misleading error message, in the worst case it just behaves incorrectly. The documentation in the McCarl User Guide does not mention this limitation. Looks like $ifthen will be less useful in practice than I had hoped.
The workaround is to use a run-time if statement:
$set stochastic 0
$set erslivestockmodel 1
scalars
stochastic /%stochastic%/
erslivestockmodel /%erslivestockmodel%/
;
if(stochastic=1,
if(erslivestockmodel=1,
b(i,"Exports",k) = tempbx(k);
b(i,"Season Average Price",k) = tempbp(k);
b(i,"Feed Use",k) = tempbdf(k);
b(i,"Food",k) = tempbdi(k);
cfixed(i,"Exports",k) = 0;
cfixed(i,"Season Average Price",k) = 0;
cfixed(i,"Feed Use",k) = 0;
cfixed(i,"Food",k) = 0;
else
$ include grains.inc
);
);
This has additional advantages: all branches are syntax-checked even if they are not executed, and we check the type of the parameters using the scalar initializations. Of course we lose some of the compile time features using this approach.
Another workaround would be to place the inner $ifthen in an include file:
$set stochastic 0Conclusion: GAMS users be aware when using nested $ifthen constructs. They are not implemented (correctly) and no error message is produced to warn you about this.
$set erslivestockmodel 1
$ifthen %stochastic%==1
$include gamsbugworkaround.inc
$endif
