数模培训第一周——优化模型

易拉罐下料问题

问题: 某公司采用一套冲压设备生产一种罐装饮料的易拉罐,这种易拉罐是用镀锡板冲压制成的.易拉罐为圆柱形,包括罐身、上盖和下底,罐身高10$cm$,上盖和下底的直径均为5$cm$.该公司使用两种不同规格的镀锡板原料:规格1的镀锡板为正方形,边长24$cm$;规格2的镀锡板为长方形,长、宽分别为32$cm$和28$cm$.由于生产设备和生产工艺的限制,对于规格$Ⅰ$的镀锡板原料,只可以按照图1中的模式1、模式⒉或模式3进行冲压;对于规格2的镀锡板原料只能按照模式4进行冲压.使用模式1模式⒉模式3、模式4进行每次冲压所需要的时间分别为1.5 s,2 s,1 s,3 s.
该工厂每周工作40 h,每周可供使用的规格1规格⒉的镀锡板原料分别为5万张和2万张.目前每只易拉罐的利润为0.10元,原料余料损失为0. 001元/cm3(如果周末有罐身、上盖或下底不能配套组装成易拉罐出售,也看做是原料余料损失).
  问工厂应该如何安排每周的生产？

问题分析

  我们已知切割的模式，只要计算各种切割模式下的余料的损失，已知上盖和下底的直径$d$，可以求得其面积$s=\frac{\pi d^2}{4}$，周长$L=\pi d$，罐高$h$，可以求得其面积$S=hL$，于是可以求出来四种模式下的余料的损失。问题要求净利润最大，我们可以列出其约束条件通过$Lingo$来求解最优化问题。

模型建立

  我们可以用$x_i$来表示第$i$种模式下冲压的次数$\left(i=1,2,3,4\right)$，$k$来表示一周所生产的易拉罐的个数，用$L_1$来表示不能配套的罐身的个数，用$L_2$表示不能配套的底、盖的个数，其中$x_i$和$k，L_1，L_2$均为整数变量。
  假设每周的易拉罐全部可以售出，则公司的销售的利润为$0.1\times k$，其中余料的损失包括两部分：

• 4种冲压模式下的余料的损失
• 不配套的罐身和底、盖造成的原料的损失

$$\begin{array}{l}\max z=0.1\times k-0.001\times \left({\sum }_1^4{r}_i{x}_i+157.08L_1+19.64L_2\right).\\ \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}\begin{array}{l}1.5x_1+2x_2+x_3+3x_4\le 144000\\ x_1+x_2+x_3\le 50000\\ x_4\le 20000\\ k\le x_1+2x_2+4x_4\end{array}\\ \begin{array}{l}2k\le 10x_1+4x_2+16x_3+5x_4\\ L_1=x_1+2x_2+4x_4-k\\ L_2=10x_1+4x_2+16x_3+5x_4-2k\\ x_1,x_2,x_3,x_4\text{取整}\\ k,L_1,L_2\text{取整}\end{array}\end{array}\right)\\ \text{其中}{r}_1=222.57,r_2=183.3,r_3=261.84,r_4=169.5\end{array}$$

求解&结果

model :

sets :
row/1..4/ : remain, x, stampingTime, spec1, spec2, body, bottom;

endsets

calc :
pi = @acos(-1); !pi值;
s1 = pi * 5 *10; !罐身面积;
s2 = pi * 5 ^ 2 / 4; !罐底或盖面积;
remain(1) = 24 * 24 - s1 - 10 * s2; !模式1余料损失;
remain(2) = 24 * 24 - 2 * s1 - 4 * s2; !模式2余料损失;
remain(3) = 24 * 24 - 16 * s2; !模式3余料损失;
remain(4) = 32 * 28 - 4 * s1 - 5 * s2; !模式4余料损失;
totalTime = 40 * 3600; !总时间;

endcalc

data :
stampingTime = 1.5, 2, 1, 3; !冲压时间;
spec1 = 1, 1, 1, 0; !规格1;
spec2 = 0, 0, 0, 1; !规格2;
totalSpec1 = 50000; !规格1总数;
totalSpec2 = 20000; !规格2总数;
body = 1, 2, 0, 4; !每个模式的罐身个数;
bottom = 10, 4, 16, 5; !每个模式的罐底或者盖子个数;

enddata


max = 0.1 * count - 0.001 * (@sum(row : remain * x) + s1 * L1 + s2 * L2);

@sum(row : stampingTime * x) <= totalTime;
@sum(row : spec1 * x) <= totalSpec1;
@sum(row : spec2 * x) <= totalSpec2;
count <= @sum(row : body * x);
2 * count <= @sum(row : bottom * x);
L1 = @sum(row : body * x) - count;
L2 = @sum(row : bottom * x) - 2 * count;
@for(row : @gin(x));
@gin(count); @gin(L1); @gin(L2);

end

  Global optimal solution found.
  Objective value:                              4298.014
  Objective bound:                              4298.014

                       Variable           Value        Reduced Cost
                          COUNT        160250.0          -0.1000000
                             L1        0.000000           0.1570796
                             L2        0.000000           0.1963495E-01
                          X( 1)        0.000000           0.2225708
                          X( 2)        40125.00           0.1833009
                          X( 3)        3750.000           0.2618407
                          X( 4)        20000.00           0.1695067

即冲压第一种模式的锡板0张，第二种模式的锡板40125张，第三种模式的锡板3750张，第四种模式的锡板20000张，总收益4298.014元

天然肠衣问题

问题:天然肠衣制作加工是我国的一个传统产业。肠衣经过清洗整理后被分割成长度不等的小段(原料)，进入组装工序。传统的生产方式依靠人工，边丈量原料长度边心算，将原材料按指定根数和总长度组装出成品(捆)。原料按长度分档，通常以0.5米为一档，如: 14-14.4米按14米计算，14.5米-14.9米按14.5米计算，其余的依此类推。为了提高生产效率，公司计划改变组装工艺，先丈量所有原料，建立一个原料表。下表某批次原料描述。

  根据以上成品和原料描述，设计一个原料搭配方案，工人根据这个方案“照方抓药”进行生产。每捆标准长度为89米，根数为5根。

  公司对搭配方案有以下具体要求：

1. 对于给定的一批原料，装出的成品捆数越多越好；
2. 对于成品捆数相同的方案，最短长度最长的成品越多，方案越好；
3. 为提高原料使用率，总长度允许有± 0.5米的误差，总根数允许比标准少1根；
4. 为减少组装的复杂性，要求组装的模式尽可能最少。这里一种模式表示各种长度的肠衣构成情况相同。

  请建立上述问题的数学模型，给出求解方法，并对表1给出的实际数据进行
求解，给出搭配方案。

问题分析

最大捆数：

进行数据清洗将材料为0的数据剔除掉，找到20种原材料，不妨设第$i$种的原材料的长度为$l_i$，根数为$a_i，i=1,2,\cdots ,20$，其中$l=14,14.5,15,15.5,16,16.6,17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,22.5,23.5,25.5$

$a=35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,6,1$

20种原材料中长度为$L={\sum }_{i=1}^{20}{a}_i\cdot l_i=12159.5$，每捆长度最少为$88.5$米，因此最多的困数为：$K⩽[12195.5/88.5]=[137.4]=137$，$20$种原材料的总共的根数为$T={\sum }_{i=1}^{20}{a}_i=677$，每捆最少为$4$根，因此最多捆数为$K⩽[677/4]=[169.25]=169$根。取最小值$K⩽137$。

每捆组成模式

每捆成品可以有不同的构成模式，每种模式由一个向量$\left(c_1,c_2,\cdots ,c_{20}\right)$组成，其中$c_i$代表第$i$种材料的根数，取最大值为：

$$M_i=\min \left\{a_i,\left[\frac{89.5}{l_i}\right],5\right\},i=1,2,\cdots ,20$$

计算得到最大取值为：$5,5,5,5,5,5,5,5,4,4,4,4,4,4,4,4,4,2,3,1$

模型建立

求解&结果

我们可以通过DFS来求得所有的模式数，也可以通过枚举剪枝来减少无效枚举的次数，通过枚举我们求得所有符合的模式数$2783$种。我们将其以矩阵的形式存放入$1.txt$，其中$b_{i}j$表示第$i$种模式中第$j$种长度的肠衣的根数，其中$i=1,2,\cdots ,2783，j=1,2,\cdots 20$，不妨设第$i$种模式$x_i$捆，则成捆最多的目标函数为：
$$\max Z_1=\sum _{i=1}^{2783}{x}_i$$

通过求得的所有符合的模式中，前三种模式为最短长度最长的模式，则有第二目标函数求最短长度最长的捆数最大：
$$\max Z_2=x_1+x_2+x_3$$

满足约束为各种长度的原料的数量有

$$\sum _{i=1}^{2783}{x}_i{b}_{ij}⩽a_j，j=1,2,\cdots 20$$

其中$a_j$为第$j$种长度的原料的根数

为使搭配的模式尽量少，我们在前面两个目标情况下考虑使搭配模式最少的方案。建立决策变量：

$$y_i\{\begin{array}{l}1\text{第}i\text{种模式选中}\\ 0\text{第}i\text{种模式未选中}\end{array}\left(i=1,2,\cdots ,2783\right)$$

设第三目标函数为总模式数最小，则有

$$\min Z_3=\sum _{i=1}^{2783}{y}_i$$

满足的约束为：

$$y_i⩽x_i⩽M{y}_i,\left(i=1,2,\cdots ,20\right)$$

其中$M$为一个足够大的整数，可以满足上述的是否选中。

综上总模型为：

$$\begin{gathered} \max Z_1=\sum _{i=1}^{2783}{x}_i \\ \max Z_2=x_1+x_2+x_3 \\ \min Z_3=\sum _{i=1}^{2783}{y}_i \\ s.t\{\begin{array}{l}{\sum }_{i=1}^{2783}{x}_i{b}_{ij}⩽a_j，j=1,2,\cdots 20\\ y_i⩽x_i⩽M{y}_i,i=1,2,\cdots ,20\\ x_i\text{为整数},y=0\text{或}1,i=1,2,\cdots ,2783\end{array} \end{gathered}$$

所有模式

%%
clear, close, clc
tic

%%
l = [14 : 0.5 : 22.5, 23.5, 25.5]; % 20种长度
a = [35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,6,1]; %20种原料的根数
cnt1 = floor(sum(a) / 4); % 由总根数求得最多捆数
cnt2 = floor(l * a' / 88.5); % 由总米数求得最多捆数
fprintf('最大捆数为：%5d\n', min(cnt1, cnt2))
n = size(l, 2);
aMin = zeros(1, n); % 初始化各种长度的类型肠衣的最多根数
for i = 1 : n
    aMin(i) = min([5, floor(89.5 / l(i)), a(i)]); % 各种长度的类型肠衣的最多根数
end
count = 0;
res = zeros(1, 20);
result = zeros(1, 20);
for i1 = 0 : aMin(1)
    result(1) = i1;
    Len = result * [l(1 : 1), zeros(1,n - 1)]';
    Gen = sum(result(1 : 1));
    if Len > 89.5 || Gen > 5
        break;
    end
    for i2 = 0 : aMin(2)
        result(2) = i2;
        Len =  result * [l(1 : 2), zeros(1,n - 2)]';
        Gen = sum(result(1 : 2));
        if Len > 89.5 || Gen > 5
            break;
        end
        for i3 = 0 : aMin(3)
            result(3) = i3;
            Len =  result * [l(1 : 3), zeros(1,n - 3)]';
            Gen = sum(result(1 : 3));
            if Len > 89.5 || Gen > 5
                break;
            end
            for i4 = 0 : aMin(4)
                result(4) = i4;
                Len =  result * [l(1 : 4), zeros(1,n - 4)]';
                Gen = sum(result(1 : 4));
                if Len > 89.5 || Gen > 5
                    break;
                end
                for i5 = 0 : aMin(5)
                    result(5) = i5;
                    Len =  result * [l(1 : 5), zeros(1,n - 5)]';
                    Gen = sum(result(1 : 5));
                    if Len > 89.5 || Gen > 5
                        break;
                    end
                    for i6 = 0 : aMin(6)
                        result(6) = i6;
                        Len =  result * [l(1 : 6), zeros(1,n - 6)]';
                        Gen = sum(result(1 : 6));
                        if Len > 89.5 || Gen > 5
                            break;
                        end
                        for i7 = 0 : aMin(7)
                            result(7) = i7;
                            Len =  result * [l(1 : 7), zeros(1,n - 7)]';
                            Gen = sum(result(1 : 7));
                            if Len > 89.5 || Gen > 5
                                break;
                            end
                            for i8 = 0 : aMin(8)
                                result(8) = i8;
                                Len =  result * [l(1 : 8), zeros(1,n - 8)]';
                                Gen = sum(result(1 : 8));
                                if Len > 89.5 || Gen > 5
                                    break;
                                end
                                for i9 = 0 : aMin(9)
                                    result(9) = i9;
                                    Len =  result * [l(1 : 9), zeros(1,n - 9)]';
                                    Gen = sum(result(1 : 9));
                                    if Len > 89.5 || Gen > 5
                                        break;
                                    end
                                    for i10 = 0 : aMin(10)
                                        result(10) = i10;
                                        Len =  result * [l(1 : 10), zeros(1,n - 10)]';
                                        Gen = sum(result(1 : 10));
                                        if Len > 89.5 || Gen > 5
                                            break;
                                        end
                                        for i11 = 0 : aMin(11)
                                            result(11) = i11;
                                            Len =  result * [l(1 : 11), zeros(1,n - 11)]';
                                            Gen = sum(result(1 : 11));
                                            if Len > 89.5 || Gen > 5
                                                break;
                                            end
                                            for i12 = 0 : aMin(12)
                                                result(12) = i12;
                                                Len =  result * [l(1 : 12), zeros(1,n - 12)]';
                                                Gen = sum(result(1 : 12));
                                                if Len > 89.5 || Gen > 5
                                                    break;
                                                end
                                                for i13 = 0 : aMin(13)
                                                    result(13) = i13;
                                                    Len =  result * [l(1 : 13), zeros(1,n - 13)]';
                                                    Gen = sum(result(1 : 13));
                                                    if Len > 89.5 || Gen > 5
                                                        break;
                                                    end
                                                    for i14 = 0 : aMin(14)
                                                        result(14) = i14;
                                                        Len =  result * [l(1 : 14), zeros(1,n - 14)]';
                                                        Gen = sum(result(1 : 14));
                                                        if Len > 89.5 || Gen > 5
                                                            break;
                                                        end
                                                        for i15 = 0 : aMin(15)
                                                            result(15) = i15;
                                                            Len =  result * [l(1 : 15), zeros(1,n - 15)]';
                                                            Gen = sum(result(1 : 15));
                                                            if Len > 89.5 || Gen > 5
                                                                break;
                                                            end
                                                            for i16 = 0 : aMin(16)
                                                                result(16) = i16;
                                                                Len =  result * [l(1 : 16), zeros(1,n - 16)]';
                                                                Gen = sum(result(1 : 16));
                                                                if Len > 89.5 || Gen > 5
                                                                    break;
                                                                end
                                                                for i17 = 0 : aMin(17)
                                                                    result(17) = i17;
                                                                    Len =  result * [l(1 : 17), zeros(1,n - 17)]';
                                                                    Gen = sum(result(1 : 17));
                                                                    if Len > 89.5 || Gen > 5
                                                                        break;
                                                                    end
                                                                    for i18 = 0 : aMin(18)
                                                                        result(18) = i18;
                                                                        Len =  result * [l(1 : 18), zeros(1,n - 18)]';
                                                                        Gen = sum(result(1 : 18));
                                                                        if Len > 89.5 || Gen > 5
                                                                            break;
                                                                        end
                                                                        for i19 = 0 : aMin(19)
                                                                            result(19) = i19;
                                                                            Len =  result * [l(1 : 19), zeros(1,n - 19)]';
                                                                            Gen = sum(result(1 : 19));
                                                                            if Len > 89.5 || Gen > 5
                                                                                break;
                                                                            end
                                                                            for i20 = 0 : aMin(20)
                                                                                result(20) = i20;
                                                                                Len =  result * [l([1 : 20]), zeros(1,n - 20)]';
                                                                                Gen = sum(result);
                                                                                if Len >= 88.5 && Len <= 89.5 && Gen >= 4 && Gen <= 5
                                                                                    count = count + 1;
                                                                                    res(count, :) = result;
                                                                                    if i20 == 1
                                                                                        result(20) = 0;
                                                                                        break;
                                                                                    end
                                                                                end
                                                                            end
                                                                        end
                                                                    end
                                                                end
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
%%
toc

%% 文件写入
fid = fopen('changyi.txt', 'w');
for i = 1 : size(res, 1);
   for j = 1 : 20
       fprintf(fid, "%3d", res(i, j));
   end
   fprintf(fid, "\n");
end
fclose(fid);

结果

所有符合的模式1.txt

求解最优模式

Model :

Sets :
row/1..2783/ : x, y;
col/1..20/ : l, a;
link(row, col) : b;

Endsets

Data :
l=14,14.5,15,15.5,16,16.5,17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,22.5,23.5,25.5;
a=35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,6,1;
b = @file('1.txt');
@text('30.txt') = @writefor(row(i) | x(i) #gt# 0 : @name(x), ' ', x(i), @newline(1));

Enddata

Min = z3;
137 = @sum(row : x); !目标函数1; 
3 = x(1) + x(2) + x(3); !目标函数2;
z3 = @sum(row : y); !目标函数3;

@for(col(j) : @sum(row(i) : x(i) * b(i, j)) <= a(j)); !原料约束;
@for(row : x <= 10000 * y);
@for(row : x >= y);
@for(row : @gin(x));
@for(row : @bin(y));

End

  Feasible solution found.
  Objective value:                              14.00000
  Objective bound:                              8.000000
  Infeasibilities:                              0.000000
  Extended solver steps:                          126503
  Total solver iterations:                       4318795

打印

%%
clear, close, clc

%% 读入数据
load 1.txt
totalModel = X1;
[~, modelRow, ~, modelNum] = textread('10.txt','%s %d %s %d');
x = [modelRow, modelNum];
%%
l = [14 : 0.5 : 22.5, 23.5, 25.5]; % 20种长度
a = [35,29,30,42,28,42,45,49,50,64,52,63,49,35,27,16,12,2,6,1]; %20种原料的根数
totalGen = 0; % 总根数
totalKun = 0; % 总捆数
totalLen = 0; % 总长度
for i = 1 : length(modelRow)
    preLen = 0; % 当前长度
    preGen = 0; % 当前根数
    i1 = modelRow(i);
    for k = 1 : length(l)
       preLen = preLen + totalModel(i1, k) * l(k);
       preGen = preGen + totalModel(i1, k);
    end
    fprintf("第%02d行，模式%04d为：", i, x(i, 1));
    for j = 1 : length(l)
       if totalModel(i1, j) > 0
           fprintf('(%3.1f米, %1d根 )', l(j), totalModel(i1, j))
       end
    end
    fprintf('长度%3.1f 米，%1d根，%2d捆 \n', preLen, preGen, x(i, 2))
    totalLen = totalLen + preLen * x(i, 2);
    totalKun = totalKun + x(i, 2);
    totalGen = totalGen + preGen * x(i, 2);
end
fprintf('总长度%6.1f， 总根数%6d， 总捆数%5d\n',totalLen, totalGen, totalKun);

结果

---

如有错误以及可以改进的地方欢迎在下方评论区留言！