二维剪板机下料问题（2-D Guillotine Cutting Stock Problem) 的混合整数规划精确求解——数学规划的计算智能特征
二维剪板机下料问题（2-D Guillotine Cutting Stock Problem) 的混合整数规划精确求解——数学规划的计算智能特征
二维剪板机下料（2D-GCSP） 的混合整数规划是最优美的整数规划模型之一。以往很多人认为像2D-GCSP这样的问题由于本质上的递归性，很难建立成混合整数规划模型，但是持这种观点的人忽略了数学规划的智能特征。所谓计算智能是如果将问题描述给机器，则机器可以自行得到问题的结果。例如很多元启发算法都被称为计算智能算法，就是因为他们好像能够通用性地解决一些问题。这也给我们一个启示，就是不怕问题刁钻，只要我们把问题的特征用建模逻辑描述出来，那么它就是合理的数学规划模型，就能为我们工作。
2D-GCSP问题
首先2D-GCSP中的G来自西文Guillotine，Sire Guillotine 是历史上一个著名法国的医生，他的医学知识包括人被割下脑袋后就不会活，进而发明了很有视觉冲击力下落式铡刀。为了纪念这位了不起的法国医生，以后凡是瞬间切断东西的方式都被命名为Guillotine。在航母、飞机、建筑等制造工业上经常会使用一种剪板机对金属板材进行切割下料，每次必须横向或纵向切断整个板材，下料图案十分独特，无论多少次切割，就是你总能找到一刀切位置。下图是剪板机。
问题数据
下表是BWMV bin packing数据库中第40个例子（需展开查看）。第一行表示例子编号，第二行30，30表示原始板材大小，80表示总计下料需求板材的种类数。第三行开始是需求定义，例如第三行表示需要5X8的小板材1件，第四行表示需要8X8的小板材1件。所有板材必须从WXH=30X30的大板材上G方式剪裁下来，目标是规划下料方式以极小化所需要的大板材数量。对原始数据，用以下符号表示：W-原始板材宽度，H-原始板材高度，n-需求小板材数量，a[i], b[i]-第i类小板材宽度和高度，m[i]-第i类小板材的需求数量。
30 80 8 1 8 1 7 1 6 1 5 1 2 1 4 1 4 1 7 1 2 1 3 1 3 1 2 1 2 1 9 1 2 1 2 1 6 1 2 1 3 1 3 1 3 1 4 1 2 1 5 1 6 1 10 1 7 1 8 1 10 1 10 1 1 1 9 1 3 1 5 1 10 1 7 1 9 1 1 1 2 1 4 1 2 1 7 1 9 1 5 1 3 1 5 1 7 1 9 1 1 1 1 1 6 1 9 1 5 1 3 1 6 1 9 1 1 1 4 1 4 1 10 1 3 1 10 1 9 1 3 1 6 1 9 1 3 1 3 1 9 1 3 1 1 1 10 1 4 1 9 1 2 1 6 1 7 1 4 1 8 1
问题数据（请展开查看）
整数规划模型
Silva等[1]首次给出G-下料的整数规划模型，此模型描述真实下料时的逻辑。如下图所示，工人在下料时候要依照给出的需求下料，那么下料后的余料被放起来准备之后的小板材下料。于是除了原始板材外，随着下料的进展，待切割板材会堆积成板材库。那么逻辑是要保证某类待切割的板材数量大于从该类下料小板材的数量。
用$x_{ij}$表示第i种需求小板材从第j种余料（包括原始板材）中下料的数量。用0-1数值矩阵$a_{ijk}$表示从j类余料中下料i类小板擦会产生一件余料k，于是根据逻辑（即工人始终能够找到需要的余料，或者说余料的数量大于将要在该类余料进行下料的数量）：
$\sum_{j=0,…,M;i=1,…,n}(a_{ij}x_{ij})\geq\sum_{i=1,…,n}(x_{ik}), k=1,…,M$ // (1)
这里的$M$是所有的余料种类，$j=0$表示原始板材（即$W\times H$板材）。
另外，下料需求需要被满足：
$\sum_{j=0,…,M}(x_{ij})\geq m_i, i=1,…,n$ // (2)
目标当然是原始板材使用数量最少：
$\min \sum_{i=1,…,n}(x_{i0})$ // (3)
下图是根本逻辑：从余料j（j=0表示原始板材）上下料i会产生两块 新生的余料备用。生成何种余料用三维0-1常数矩阵$a_{ijk}$表示。约束（1）表示一定要保证有足够的余料可用于需求小板材的下料。
至于如何生成$a_{ijk}$, Silva等[1]建议首先做枚举。本文以下直接将这个枚举写进模型里！
模型实现
上面的模型（1）-（3）比较抽象，现在把它具体化，并把对$a_{ijk}$的枚举写进模型。
首先，对余料从0到$M$的编码抽象性太强，这里把他们具体化为板材的长度和宽度，即使用x[w][h][i]表示从大小为wXh余料上切割第i类小板材的数量。于是目标是极小化从WxH板材（原始板材）上切割出所有小板材的数量，即：
min sum{i=1,…,n}x[W][H][i] //（4）
也就是使用（4）代替上面的（3）。
其次必须满足所有需求：
sum{w=1,…,W;h=1,…,H;w>=a[i];h>=b[i]}x[w][h][i] >= m[i] | i=1,…,n //(5)
也就是用（5）取代上面的（2）。
最后用下面的（6）取代上面的（1）：
sum{w=1,…,W;h=1,…,H;i=1,…,n;w>=a[i];h>=b[i];w==(a[i]+w1);b[i]==h1}x[w][h][i]+–>
sum{w=1,…,W;h=1,…,H;i=1,…,n;w>=a[i];h>=b[i];h==(b[i]+h1);w1==w}x[w][h][i]–>
>= sum{i=1,…,n;w1>=a[i];h1>=b[i]}x[w1][h1][i]–>
| w1=1,…,W;h1=1,…H;(h1+w1)<(H+W) //(6)
这个看起来复杂，仔细分析会发现跟上面的（1）是一丘之貉，不过好处是扔掉了$a_{ijk}$。
整体模型看起来是这个样子的（数据区采用小案例，节省空间）：
min sum{i=1,...,n}x[W][H][i] //(4)
subject to
sum{w=1,...,W;h=1,...,H;w>=a[i];h>=b[i]}x[w][h][i] >=m[i]|i=1,...,n //(5)
sum{w=1,...,W;h=1,...,H;i=1,...,n;w>=a[i];h>=b[i];w==(a[i]+w1);b[i]==h1}x[w][h][i]+-->
sum{w=1,...,W;h=1,...,H;i=1,...,n;w>=a[i];h>=b[i];h==(b[i]+h1);w1==w}x[w][h][i]-->
>=sum{i=1,...,n;w1>=a[i];h1>=b[i]}x[w1][h1][i]-->
| w1=1,...,W;h1=1,...H;(h1+w1)<(H+W) //(6)
where
n is an integer
W,H are numbers
dt is a set of number
a[i],b[i],m[i] are numbers|i=1,...,n
x[w][h][i] is a variable of nonnegative integer|w=1,...,W;h=1,...,H;i=1,...,n
data_relation
n=_$(dt)/3
a[i]=dt[3(i-1)+1] | i=1,...,n
b[i]=dt[3(i-1)+2] | i=1,...,n
m[i]=dt[3i] | i=1,...,n
data
W=10
H=10
dt={
2 2 1
8 6 1
2 10 1
3 1 1
4 8 1
10 3 1
9 1 1
5 1 1
3 6 1
1 1 1
2 4 1
2 9 1
9 1 1
5 9 1
7 4 1
2 2 1
4 3 1
7 9 1
1 4 1
8 9 1
}
上面的模型用+Leapms解析并用cplex命令调用cplex求解即可秒解。
下料排样图的生成
上面的模型只给出x[w][h][i]的最优值，即从wxh余料中切割第i类小板材的具体数量，由此信息得到排样图是另外一个挑战。首先看结果样子（以上面80小板材数据为例）：
而后看技术。这里技术包括用+Leampms编程接口，还有细致的逻辑思考，这里不再赘述，只贴出c++代码，此源码能够一次求出文件中500个案例：
1 #include<iostream> 2 #include<fstream> 3 using namespace std; 4 5 extern int W,H,n; 6 extern int w[1000],h[1000],m[1000]; 7 extern ofstream ofcad; 8 extern int nInstance; 9 10 class c_sheet{ 11 public: 12 int w,h; //尺寸 13 int n,nl; //数量，已定位的数量 14 int x[1000],y[1000]; //坐标 15 int item[1000]; //下料 16 bool good; 17 c_sheet(){ 18 w=W; 19 h=H; 20 n=nl=0; 21 } 22 c_sheet(int w,int h){ 23 this->w=w; 24 this->h=h; 25 n=0; 26 good=false; 27 } 28 void reset(){ 29 w=h=n=nl=0; 30 good=false; 31 } 32 void addItem(int item); 33 void print(); 34 void draw(ofstream &off); 35 bool cut(); 36 bool cut(int x, int y); 37 }; 38 39 int prepare(int instance=0);//准备模型 40 int findSheet(int w, int h); 41 extern int nSheet; 42 extern c_sheet sheet[1000];
#include<iostream> #include<fstream> #include "leap.h" #include "Guillotine.h" using namespace std; int W,H,n; int w[1000],h[1000],m[1000]; ofstream ofcad; void addSheet(int w, int h, int item, int ns); int dealInstance();
//主函数，循环计算所有500案例 int main(){ while(1){ if(dealInstance()<0)break; } return 0; } int dealInstance(){ //设置时间 char time[]="600"; //准备模型 if(prepare()<0)return -1; //if(nInstance<41)return 0; //求解模型 c_leap lp; lp.loadLeap("current.leap"); lp.cplex("",time,""); //生成板材 int ns; for(int w=W;w>0;w--){ for(int h=H;h>0;h--){ for(int i=0;i<n;i++){ ns=(int)lp.getVar("x",3,w,h,i+1); if(ns>0){ addSheet(w,h,i,ns); } } } } //剪切板材 if(nSheet>1)sheet[0].cut(); //输出板材图形 char cadfname[64]; if(lp.getState()!=IP_OPTIMAL)sprintf_s(cadfname,64,"graphic(%d)(%d_%d)!.scr",nInstance,(int)((nInstance-1)/50)+1,nInstance-(int)((nInstance-1)/50)*50); else sprintf_s(cadfname,64,"graphic(%d)(%d_%d).scr",nInstance,(int)((nInstance-1)/50)+1,nInstance-(int)((nInstance-1)/50)*50); ofcad.open(cadfname); ofcad<<"erase all "<<endl; for(int i=0;i<nSheet;i++){ sheet[i].draw(ofcad); } ofcad<<"zoom e"<<endl; ofcad.close(); //清空sheet for(int i=0;i<nSheet;i++){ sheet[i].reset(); } return 0; } void addSheet(int w, int h, int item, int ns){ int sht=findSheet(w,h); if(sht<0){ sht=nSheet++; sheet[sht].w=w; sheet[sht].h=h; } for(int k=0;k<ns;k++){ sheet[sht].addItem(item); } }
#include<iostream> #include<fstream> #include "leap.h" #include "Guillotine.h" using namespace std; ifstream iffd;//原始数据文件句柄 bool iffdOpen=false; int nInstance=0;//案例编号 int prepare(int instance){ //准备模型 if(nInstance>500)return -1; //如果没有打开原始数据文件，打开之====================== const char iffdname[]="BWMV bin packing instances.txt"; if(!iffdOpen){ int ntmp; iffd.open(iffdname); iffd>>ntmp; cout<<"ntmp="<<ntmp<<endl; iffdOpen=true; } if(!iffd){ cout<<"Can not open the file"<<iffdname<<endl; exit(0); } //打开模板文件并复制==================================== ifstream iffm; ofstream off; iffm.open("Guillotine.leap"); if(!iffm){ cout<<"\tCannot open file 'Guillotine.leap'"<<endl; return 0; } off.open("current.leap"); if(!off){ cout<<"\tCannot create file 'current.txt'"<<endl; return 0; } char line[4096]; while(!iffm.eof()){ iffm.getline(line,4096); off<<line<<endl; } //生成Leap模型数据区===================================== nInstance++; int dummy; iffd>>dummy; iffd>>W>>H>>n; off<<"\tW="<<W<<endl; off<<"\tH="<<H<<endl; off<<"\tdt={"<<endl; for(int i=0;i<n;i++){ iffd>>w[i]>>h[i]>>m[i]; off<<"\t\t"<<w[i]<<"\t"<<h[i]<<"\t"<<m[i]<<endl; } off<<"\t}"<<endl; iffm.close();off.close(); //===================================================== nSheet=0; return nInstance; }
#include<iostream> #include<fstream> #include "Guillotine.h" int nSheet=0; c_sheet sheet[1000]; int findSheet(int w,int h){ for(int i=0;i<nSheet;i++){ if(w==sheet[i].w&&h==sheet[i].h)return i; } return -1; } void c_sheet::addItem(int item){ this->item[n]=item; n++; } void c_sheet::print(){ cout<<"Sheet"<<w<<"X"<<h; cout<<"\titems:"; for(int i=0;i<n;i++){ cout<<item[i]+1<<","; } cout<<endl; } bool c_sheet::cut(){ bool rt=true; int w1,h1,w2,h2,s1,s2; if(w==W&&h==H){ nl=n; for(int i=0;i<=n;i++){ x[i]=(int)(i*1.2*W+.5); y[i]=0; w1=w-::w[item[i]]; h1=::h[item[i]]; s1=findSheet(w1,h1); if(s1>=0)sheet[s1].cut(x[i]+::w[item[i]],y[i]); w2=w; h2=h-::h[item[i]]; s2=findSheet(w2,h2); if(s2>=0)sheet[s2].cut(x[i],y[i]+::h[item[i]]); } } return rt; } bool c_sheet::cut(int xx, int yy){ if(nl==n)return false; int w1,h1,w2,h2,s1,s2; x[nl]=xx;y[nl]=yy; w1=w-::w[item[nl]]; h1=::h[item[nl]]; s1=findSheet(w1,h1); if(s1>=0)sheet[s1].cut(xx+::w[item[nl]],yy); w2=w; h2=h-::h[item[nl]]; s2=findSheet(w2,h2); if(s2>=0)sheet[s2].cut(xx,yy+::h[item[nl]]); nl++; return true; } void c_sheet::draw(ofstream &off){ bool hangle=true; for(int i=0;i<nl;i++){ if(w==W&&h==H)off<<"rectangle "<<x[i]<<","<<y[i]<<" @"<<w<<","<<h<<endl; off<<"rectangle "<<x[i]<<","<<y[i]<<" @"<<::w[item[i]]<<","<<::h[item[i]]<<endl; off<<"-hatch p ansi31 0.25 "<<hangle*90<<" "<<endl; hangle=!hangle; off<<"-hatch "<<x[i]+::w[item[i]]/2.0<<","<<y[i]+::h[item[i]]/2.0<<" "<<endl; off<<"line "<<x[i]<<","<<y[i]+::h[item[i]]<<" @"<<w<<",0"<<endl<<endl; } }
参考文献
[1] Silva E , Alvelos F , J.M. Valério de Carvalho. An integer programming model for two- and three-stage two-dimensional cutting stock problems[J]. European Journal of Operational Research, 2010, 205(3):