【数学建模】天然肠衣搭配问题

2011高教社杯全国大学生数学建模竞赛D题
天然肠衣（以下简称肠衣）制作加工是我国的一个传统产业，出口量占世界首位。肠衣经过清洗整理后被分割成长度不等的小段（原料），进入组装工序。传统的生产方式依靠人工，边丈量原料长度边心算，将原材料按指定根数和总长度组装出成品（捆）。
原料按长度分档，通常以0.5米为一档，如：3-3.4米按3米计算，3.5米-3.9米按3.5米计算，其余的依此类推。表1是几种常见成品的规格，长度单位为米，∞表示没有上限，但实际长度小于26米。
\space   \space   \space   \space   \space   \space   \space   表1 成品规格表

最短长度	最大长度	根数	总长度
3	6.5	20	89
7	13.5	8	89
14	$\infty$	5	89

为了提高生产效率，公司计划改变组装工艺，先丈量所有原料，建立一个原料表。表2为某批次原料描述。

根据以上成品和原料描述，设计一个原料搭配方案，工人根据这个方案“照方抓药”进行生产。
公司对搭配方案有以下具体要求：
(1) 对于给定的一批原料，装出的成品捆数越多越好；
(2) 对于成品捆数相同的方案，最短长度最长的成品越多，方案越好；
(3) 为提高原料使用率，总长度允许有

提出假设

假设降级使用仅可降一级，不可多级降，如如长度为14米的原料不可以和长度介于3-6.5米的进行捆扎

问题简单复述:

用以0.5米为一档给出的原料按指定根数和总长度组装出成品，为此设计一个原料搭配方案。
方案好坏比较首先是装出的成品捆数越多越好；其次捆数相同的方案，最短长度为14的成品越多，方案越好
使用附带原则有总长度允许有$\pm 0.5$米的误差，总根数允许比标准少$1$根；某种规格对应原料如果出现剩余，可以降级使用。

建立数学模型:
首先考虑装出的成品捆数，设三种不同最短长度(从小到大)的成品个数依次为$x_1,x_2,x_3$
装出的成品捆数越多越好有:
$max={\sum }_{1\le i\le 3}{x}_i$

为了简化问题，解决最短长度为14的成品越多，方案越好的方案为:
设已求得一组可行解$x_1^{\prime },x_2^{\prime },x_3^{\prime }$
加入限制条件
$x_3>x_3^{\prime }$

考虑原料如果出现剩余，可以降级使用。那么解决顺序应该是成品三，成品二，成品一。
因为每个成品最短长度到最长长度的区间没有重叠，所以三个成品可以分开讨论。

对应每一个任意种类的成品建立第一种模型

设原料按长度由小到大的使用个数依次为 $y_1,y_2,...,y_n$,总长度为$a$,根数为$b_{i^{\prime }}$,最小大长度$c_{i^{\prime }}$,最大长度$c_{j^{\prime }}$
都要满足
成品总长度在标准的$\pm 0.5$米误差范围内:
$a-0.5\le (3+0.5(i-1))y_i+(3+0.5i)y_{i+1}+...(3+0.5(j-1))y_j\le a+0.5$
总根数和标准一致或总根数比标准少$1$根:
$y_i+y_{i+1}+...+y_j=b_{i^{\prime }}-z,z=0/1$
使用的原料长度在成品规格表的标准长度范围内
$2(c_{i^{\prime }}-3)-1\le y_k\le 2(c_{j^{\prime }}-3)-1,i\le k\le j$
考虑原料如果出现剩余，可以降级使用。我们设上一级(如果有的话)的的最大长度为$c_{j^{\prime }+1}$
$2(c_{i^{\prime }}-3)-1\le y_k\le 2(c_{j^{\prime }+1}-3)-1,i\le k\le j$

因为做出一个成品，对应使用的原料数量会降低，设原料按长度由小到大的减少总个数依次为 $u_1,u_2,...,u_n$,原料按长度由小到大的总个数依次为 $d_1{d}_2,...,d_n$

我们要满足每一次原料都不能凭空产生
$y_i\le d_i-u_i,1\le i\le n$

并且当满足$\{\begin{array}{l}a-0.5\le (3+0.5(i-1))y_i+(3+0.5i)y_{i+1}+...(3+0.5(j-1))y_j\le a+0.5\\ y_i+y_{i+1}+...+y_j=b_i-z,z=0/1\\ 2(c_{i^{\prime }}-3)-1\le y_k\le 2(c_{j^{\prime }+1}-3)-1,i\le k\le j\\ y_i\le d_i-u_i,1\le i\le n\end{array}$的时候对应种类的成品数量加一，即$x_i=x_i+1$,并且原材料使用后对应的原材料数目减少，即$u_i=u_i+y_i,1\le i\le n$

C++

#include <iostream>
#include <vector>
#include <array>
using namespace std;
int cnt[] ={43,59,39,41,27,28,34,21,24,24,20,25,21,23,21,18,31,23,22,59,18,25,35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,0,6,0,0,0,1}; // 3+0.5i
int x[4];
int x_3 = 0;//已知x3的最大解
void solve();double min_size = 3,max_size =6.5,total_size =89;
int roots_Number = 20;
int f_xi = 1;double now_total_size=0;
int now_roots_Number = 0;
void f(  int k ,int u){if(now_total_size <= 0.5 + total_size && now_total_size >= total_size - 0.5){if(now_roots_Number == roots_Number || now_roots_Number == roots_Number-1){x[f_xi] = max(x[f_xi],u+1);double kt = now_total_size;int kr = now_roots_Number;now_total_size = now_roots_Number = 0;f(((int)max_size - 3)*2,u+1);now_total_size = kt , now_roots_Number = kr;return ;}else return ;}//if(now_total_size>0.5 + total_size)return ; // cut down
#if 0 //small to bigif(3+0.5*k > max_size)return ; // can't choseint ct = min(cnt[k] ,(int)((total_size -  now_total_size)/(3+0.5*k)));for(int i=0;i<=ct;i++){cnt[k]-=i;f(now_total_size+i*(3+0.5*k) , now_roots_Number+i ,k+1,u );cnt[k]+=i;}
#endif
#if 1 // big to smallif(3+0.5*k < min_size)return ; // can't choseint ct = min(cnt[k] ,(int)((total_size -  now_total_size)/(3+0.5*k)));ct = min(ct,roots_Number - now_roots_Number);//if(now_roots_Number > roots_Number)return ;for(int i=0;i<=ct;i++){cnt[k]-=i;now_total_size += i*(3+0.5*k);now_roots_Number+=i;f(k-1,u);now_total_size -= i*(3+0.5*k);now_roots_Number-=i;cnt[k]+=i;}
#endif
}void solve(){//排除要求四/*min_size = 3,max_size =6.5,total_size =89;roots_Number = 20;f_xi = 1;now_total_size=0;now_roots_Number = 0;f(((int)max_size - 3)*2,0);//3+0.5*k = min_size//k = (min_size - 3) / 0.5min_size = 7,max_size =13.5,total_size =89;roots_Number = 8;f_xi = 2;now_total_size=0;now_roots_Number = 0;f(((int)max_size - 3)*2,0);*/min_size = 14,max_size =25.5,total_size =89;roots_Number = 5;f_xi = 3;now_total_size=0;now_roots_Number = 0;f((max_size - 3)*2,0);cout << x[3];
}int main() {solve();return 0;
}

基于成品分配方案的第二种模型

考虑到第一种模型跑程序时候，就算剪枝也会跑许多多余的分支，考虑到本题是分配问题，并且每一种成品规格限制比较死，所以我们可以列出该成品分配的所有方案。

可以自行运行一下程序查看所有方案,对应成品三有1823种
C++

#include <iostream>
#include <vector>
#include <array>
using namespace std;
int cnt[46] ={43,59,39,41,27,28,34,21,24,24,20,25,21,23,21,18,31,23,22,59,18,25,35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,0,6,0,0,0,1}; // 3+0.5i
int x[4];
int x_3 = 0;//已知x3的最大解
void solve();double min_size = 3,max_size =6.5,total_size =89;
int roots_Number = 20;
int f_xi = 1;double now_total_size=0;
int now_roots_Number = 0;
int bt = 1;
int a[46];
void f(  int k ,int u){if(now_total_size <= 0.5 + total_size && now_total_size >= total_size - 0.5){if(now_roots_Number == roots_Number || now_roots_Number == roots_Number-1){cout << "\n\n方案" << bt << ":\n|";int kk = 0;for(int i=0;i<=45;i++){if(a[i]){cout << 3+0.5*i << '|';kk++;}}cout << "\n|";for(int i=0;i<kk;i++){cout << "--|";}cout << "\n|";for(int i=0;i<=45;i++){if(a[i])cout << a[i] << '|';}bt++;return ;}else return ;}//if(now_total_size>0.5 + total_size)return ; // cut down
#if 0 //small to bigif(3+0.5*k > max_size)return ; // can't choseint ct = min(cnt[k] ,(int)((total_size -  now_total_size)/(3+0.5*k)));for(int i=0;i<=ct;i++){cnt[k]-=i;f(now_total_size+i*(3+0.5*k) , now_roots_Number+i ,k+1,u );cnt[k]+=i;}
#endif
#if 1 // big to smallif(3+0.5*k < min_size)return ; // can't choseint ct = min(cnt[k] ,(int)((total_size -  now_total_size)/(3+0.5*k)));ct = min(ct,roots_Number - now_roots_Number);//if(now_roots_Number > roots_Number)return ;for(int i=0;i<=ct;i++){cnt[k]-=i;now_total_size += i*(3+0.5*k);now_roots_Number+=i;a[k]+=i;f(k-1,u);now_total_size -= i*(3+0.5*k);now_roots_Number-=i;cnt[k]+=i;a[k]-=i;}
#endif
}void solve(){//排除要求四/*min_size = 3,max_size =6.5,total_size =89;roots_Number = 20;f_xi = 1;now_total_size=0;now_roots_Number = 0;f(((int)max_size - 3)*2,0);//3+0.5*k = min_size//k = (min_size - 3) / 0.5min_size = 7,max_size =13.5,total_size =89;roots_Number = 8;f_xi = 2;now_total_size=0;now_roots_Number = 0;f(((int)max_size - 3)*2,0);*/min_size = 14,max_size =25.5,total_size =89;roots_Number = 5;f_xi = 3;now_total_size=0;now_roots_Number = 0;f((max_size - 3)*2,0);}int main() {solve();return 0;
}

下面列出其中几个方案
方案1:

17.5	18
3	2

方案2:

17	17.5	18
1	1	3

方案3:

17.5	18
2	3

方案1467:

14	14.5	16	21.5	22.5
1	1	1	1	1

方案1653:

14	14.5	15	21.5	23.5
1	1	1	1	1

方案1800:

20	20.5	22.5	25.5
1	1	1	1

方案1823:

16.5	23.5	25.5
1	2	1

列出该成品分配的所有方案后，设计第二种模型:
设使用方案$i$的次数为$y_i$，原料按长度由小到大的总个数依次为 $d_1{d}_2,...,d_n$
方案总数为$a$;方案$i$使用第$j$种原料数量为$z_{ij}$
对于每一个原料我们都得
$d_i\ge {\sum }_{1\le j\le a}{y}_j\ast z_{ji},1\le i\le n$
求 $\max {\sum }_{1\le i\le a}{y}_i$

对于成品三:
Lingo

sets:aa/1..46/:d;bb/1..1823/:y;cc(bb,aa):z;
endsets
data:d = 43,59,39,41,27,28,34,21,24,24,20,25,21,23,21,18,31,23,22,59,18,25,35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,0
,6,0,0,0,1;enddata
max = @sum(bb(i):y(i));
@for(bb(i):@gin(y(i)));
@for(bb(i):y(i)>0);
@for(aa(i):d(i)>@sum(bb(j):y(j)*z(j,i)));

将下面C++程序运行结果带入上面Lingo中 , 在data: enddata 里面加 z = …

#include <iostream>
#include <vector>
#include <array>
using namespace std;
int cnt[46] ={43,59,39,41,27,28,34,21,24,24,20,25,21,23,21,18,31,23,22,59,18,25,35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,0,6,0,0,0,1}; // 3+0.5i
void solve();double min_size = 3,max_size =6.5,total_size =89;
int roots_Number = 20;
double now_total_size=0;
int now_roots_Number = 0;
int bt = 1;
int a[46];
void f(  int k ,int u){if(now_total_size <= 0.5 + total_size && now_total_size >= total_size - 0.5){if(now_roots_Number == roots_Number || now_roots_Number == roots_Number-1){//方案数在后for(int i=0;i<=45;i++){//cout << "z(" << i  << "," << bt << ")=" << a[i] << '\n';cout << a[i] << ",";}bt++;return ;}else return ;}if(3+0.5*k < min_size)return ; // can't choseint ct = min(cnt[k] ,(int)((total_size -  now_total_size)/(3+0.5*k)));ct = min(ct,roots_Number - now_roots_Number);//if(now_roots_Number > roots_Number)return ;for(int i=0;i<=ct;i++){cnt[k]-=i;now_total_size += i*(3+0.5*k);now_roots_Number+=i;a[k]+=i;f(k-1,u);now_total_size -= i*(3+0.5*k);now_roots_Number-=i;cnt[k]+=i;a[k]-=i;}
}void solve(){min_size = 14,max_size =25.5,total_size =89;roots_Number = 5;f_xi = 3;now_total_size=0;now_roots_Number = 0;f((max_size - 3)*2,0);
}int main() {solve();return 0;
}

Lingo求解得

Objective value:                              137.0000
Objective bound:                              137.0000

修改一下Lingo,添加下面这一行可以查看剩余原料数量

@for(aa(i):u(i) = d(i)-@sum(bb(j):y(j)*z(j,i)));

同上处理成品一二得到最终程序
Lingo

sets:aa/1..46/:d,u,v;bb/1..1823/:y;cc(bb,aa):z;dd/1..466598/:x;ee(dd,aa):e;ff/1..2861814/:fgg(ff,aa):h;
endsets
data:d = 43,59,39,41,27,28,34,21,24,24,20,25,21,23,21,18,31,23,22,59,18,25,35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,0
,6,0,0,0,1;enddata
!成品三
!max = @sum(bb(i):y(i));
@sum(bb(i):y(i)) > 137;
@for(bb(i):@gin(y(i)));
@for(bb(i):y(i)>0);
@for(aa(i):d(i)>@sum(bb(j):y(j)*z(j,i)));
@for(aa(i):u(i) = d(i)-@sum(bb(j):y(j)*z(j,i)));
!成品二
@sum(dd(i):x(i)) > ...;!求成品三后成品二最大组成数量+@for(dd(i):@gin(x(i)));
@for(dd(i):x(i)>0);
@for(aa(i):u(i)>@sum(dd(j):x(j)*e(j,i)));
@for(aa(i):v(i) = d(i)-@sum(dd(j):x(j)*e(j,i)));
!成品一
@for(ff(i):@gin(f(i)));
@for(ff(i):f(i)>0);
@for(aa(i):v(i)>@sum(ff(j):f(j)*h(j,i)));
max = @sum(bb(i):y(i)) + @sum(dd(j):x(j)) + @sum(dd(j):x(j));

需要在C++/其他程序中跑出成品一的2861814种方案的列表和成品二的466598种方案的列表。加在上述代码中。或者用其他方式给e和f赋值

e = ... !成品二的方案
h = ... !成品一的方案

模型二总结

思路总结：配完成品三，得出剩余的原材料数量，再配成品二，得出剩余原材料数量，最后配成品一。最后得到最优配比

目前的问题：无法有效的解决剩余这个关键，剩余条件为剩下的无法组成任何一种成品方案。
仅靠代码@sum(bb(i):y(i)) > 137;去约束剩余这个条件，只是局部最优，而不一定是全局最优

代码优化思路：例如成品三对于前面20多个都没用，直接优化掉;成品二一同上

难点：目前已知成品一方案数多达百万量级

手动实现局部最优解

虽然得出的答案不一定是全局最优，但可以先试试手动模拟成品三，二，一的步骤，先得出一种(多种)局部最优。如果有规律或者新的上限模型，就可以推出某一个局部最优解即是全局最优
按照模型二跑一次成品三

  Objective value:                              137.0000Objective bound:                              137.0000U( 30)        1.000000            0.000000

长度范围为$17.5-17.9$剩余1根
将剩余的长度的带入，得出成品二的分配方案有12729种，相比于带入无处理成品三后数据的直接计算的分配方案少了一个量级。

将成品二的分配方案带入模型和程序中得:

Objective value:                              37.00000Objective bound:                              37.00000U( 1)        43.00000            0.000000U( 2)        59.00000            0.000000U( 3)        39.00000            0.000000U( 4)        41.00000            0.000000U( 5)        27.00000            0.000000U( 6)        28.00000            0.000000U( 7)        34.00000            0.000000U( 8)        21.00000            0.000000U( 9)        23.00000            0.000000U( 10)        23.00000            0.000000U( 11)        7.000000            0.000000U( 12)        1.000000            0.000000U( 14)        2.000000            0.000000U( 15)        2.000000            0.000000U( 17)        1.000000            0.000000

将剩余的长度的带入，得出成品一的分配方案有1213327种，相比于带入无处理成品二后数据的直接计算的分配方案量级没变。

但如果不算降级处理的东西分配方案仅有59942种

单独解决成品一:

  Global optimal solution found.Objective value:                              14.00000Objective bound:                              14.00000Infeasibilities:                              0.000000Extended solver steps:                               0Total solver iterations:                          8821Variable           Value        Reduced CostU( 2)        5.000000            0.000000U( 4)        6.000000            0.000000U( 5)        1.000000            0.000000U( 6)        3.000000            0.000000U( 9)        23.00000            0.000000U( 10)        23.00000            0.000000U( 11)        7.000000            0.000000U( 12)        1.000000            0.000000U( 14)        2.000000            0.000000U( 15)        2.000000            0.000000U( 17)        1.000000            0.000000

局部最优解

成品三137个，成品二37个，成品三14个，总188个