%% 原料数据初始化
lengths = [3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, ...
7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, ...
11, 11.5, 12, 12.5, 13, 13.5, 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, 23.5, 24, 24.5, 25, 25.5];
stocks = [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];
% 原料分组索引
spec1_direct = 1:8; % 规格1直接使用 (3-6.5m)
spec2_direct = 9:22; % 规格2直接使用 (7-13.5m)
spec3_direct = 23:46; % 规格3直接使用 (14m+)
% 可降级原料索引
downgrade1_idx = []; % 可降级至规格1的原料 (6.5 < len <= 8.5)
downgrade2_idx = []; % 可降级至规格2的原料 (13.5 < len <= 22.5)
for j = spec2_direct
if lengths(j) > 6.5 && lengths(j) <= 8.5
downgrade1_idx = [downgrade1_idx, j];
end
end
for j = spec3_direct
if lengths(j) > 13.5 && lengths(j) <= 22.5
downgrade2_idx = [downgrade2_idx, j];
end
end
%% 估算最大捆数
total_stocks = sum(stocks);
n_max = min([floor(total_stocks/4), ... % 规格3最小根数
floor(total_stocks/7), ... % 规格2最小根数
floor(total_stocks/19)]); % 规格1最小根数
fprintf('最大潜在捆数估算: %d\n', n_max);
%% 建立模型基础结构
n = n_max; % 最大潜在捆数
num_vars = 3*n + 46*n + length(downgrade1_idx)*n + length(downgrade2_idx)*n;
% 变量索引映射
x_start = 1; % X_ik (3n个变量)
c_start = 3*n + 1; % C_ij (46n个变量)
d1_start = c_start + 46*n;% D_ij^{(1)}
d2_start = d1_start + length(downgrade1_idx)*n; % D_ij^{(2)}
%% 第一阶段:最大化总捆数
fprintf('\n=== 第一阶段:最大化总捆数 ===\n');
% 目标函数:最大化总捆数
f = zeros(num_vars, 1);
f(x_start:x_start+3*n-1) = 1; % 所有X_ik系数为1
% 约束计数器
constr_count = 0;
% 约束1:每个捆选择恰好一种规格 (n个等式约束)
A1 = zeros(n, num_vars);
for i = 1:n
idx = x_start + 3*(i-1);
A1(i, idx:idx+2) = 1;
end
b1 = ones(n, 1);
constr_count = constr_count + n;
% 约束2:规格1长度约束 (2n个不等式)
A2 = zeros(2*n, num_vars);
b2 = zeros(2*n, 1);
for i = 1:n
% 下界约束: 88.5*X_i1 <= sum
idx = x_start + 3*(i-1);
A2(2*i-1, idx) = -88.5; % -88.5*X_i1
% 上界约束: sum <= 89.5*X_i1
A2(2*i, idx) = 89.5; % 89.5*X_i1
% 原料求和项
for j = spec1_direct
var_idx = c_start + 46*(i-1) + j - 1;
A2(2*i-1, var_idx) = lengths(j); % 正系数
A2(2*i, var_idx) = -lengths(j); % 负系数
end
for j = downgrade1_idx
d_idx = find(downgrade1_idx == j);
var_idx = d1_start + length(downgrade1_idx)*(i-1) + d_idx - 1;
A2(2*i-1, var_idx) = lengths(j); % 正系数
A2(2*i, var_idx) = -lengths(j); % 负系数
end
end
b2(1:2:end) = 0; % 下界约束 <= 0
b2(2:2:end) = 0; % 上界约束 <= 0
constr_count = constr_count + 2*n;
% 约束3:规格2长度约束 (2n个不等式)
A3 = zeros(2*n, num_vars);
b3 = zeros(2*n, 1);
for i = 1:n
idx = x_start + 3*(i-1) + 1; % X_i2
A3(2*i-1, idx) = -88.5; % -88.5*X_i2
A3(2*i, idx) = 89.5; % 89.5*X_i2
for j = spec2_direct
var_idx = c_start + 46*(i-1) + j - 1;
A3(2*i-1, var_idx) = lengths(j);
A3(2*i, var_idx) = -lengths(j);
end
for j = downgrade2_idx
d_idx = find(downgrade2_idx == j);
var_idx = d2_start + length(downgrade2_idx)*(i-1) + d_idx - 1;
A3(2*i-1, var_idx) = lengths(j);
A3(2*i, var_idx) = -lengths(j);
end
end
b3(1:2:end) = 0;
b3(2:2:end) = 0;
constr_count = constr_count + 2*n;
% 约束4:规格3长度约束 (2n个不等式)
A4 = zeros(2*n, num_vars);
b4 = zeros(2*n, 1);
for i = 1:n
idx = x_start + 3*(i-1) + 2; % X_i3
A4(2*i-1, idx) = -88.5; % -88.5*X_i3
A4(2*i, idx) = 89.5; % 89.5*X_i3
for j = spec3_direct
var_idx = c_start + 46*(i-1) + j - 1;
A4(2*i-1, var_idx) = lengths(j);
A4(2*i, var_idx) = -lengths(j);
end
end
b4(1:2:end) = 0;
b4(2:2:end) = 0;
constr_count = constr_count + 2*n;
% 约束5:库存限制 (46个不等式)
A5 = zeros(46, num_vars);
b5 = stocks';
for j = 1:46
for i = 1:n
% C_ij
var_idx = c_start + 46*(i-1) + j - 1;
A5(j, var_idx) = 1;
% D_ij^{(1)}
if ismember(j, downgrade1_idx)
d_idx = find(downgrade1_idx == j);
var_idx = d1_start + length(downgrade1_idx)*(i-1) + d_idx - 1;
A5(j, var_idx) = 1;
end
% D_ij^{(2)}
if ismember(j, downgrade2_idx)
d_idx = find(downgrade2_idx == j);
var_idx = d2_start + length(downgrade2_idx)*(i-1) + d_idx - 1;
A5(j, var_idx) = 1;
end
end
end
constr_count = constr_count + 46;
% 合并所有约束
A = [A1; A2; A3; A4; A5];
b = [b1; b2; b3; b4; b5];
% 变量边界
lb = zeros(num_vars, 1);
ub = Inf(num_vars, 1);
% X_ik为0-1变量
ub(x_start:x_start+3*n-1) = 1;
% 整数约束:所有变量为整数
intcon = 1:num_vars;
% 求解第一阶段
options = optimoptions('intlinprog', 'Display', 'iter', 'MaxTime', 600);
[x1, fval1, exitflag1] = intlinprog(-f, intcon, A, b, [], [], lb, ub, options);
if exitflag1 <= 0
error('第一阶段求解失败');
end
total_bundles = round(fval1);
fprintf('第一阶段结果:最大总捆数 = %d\n', total_bundles);
%% 第二阶段:最大化第3类捆数
fprintf('\n=== 第二阶段:最大化第3类捆数 ===\n');
% 添加约束:总捆数等于第一阶段结果
A_eq1 = zeros(1, num_vars);
A_eq1(x_start:x_start+3*n-1) = 1;
b_eq1 = total_bundles;
% 目标函数:最大化第3类捆数
f2 = zeros(num_vars, 1);
for i = 1:n
idx = x_start + 3*(i-1) + 2; % X_i3
f2(idx) = 1;
end
% 求解第二阶段
[x2, fval2, exitflag2] = intlinprog(-f2, intcon, A, b, A_eq1, b_eq1, lb, ub, options);
if exitflag2 <= 0
error('第二阶段求解失败');
end
spec3_bundles = round(-fval2);
fprintf('第二阶段结果:第3类捆数 = %d\n', spec3_bundles);
%% 第三阶段:最大化第2类捆数
fprintf('\n=== 第三阶段:最大化第2类捆数 ===\n');
% 添加约束:第3类捆数等于第二阶段结果
A_eq2 = zeros(1, num_vars);
for i = 1:n
idx = x_start + 3*(i-1) + 2; % X_i3
A_eq2(idx) = 1;
end
b_eq2 = spec3_bundles;
% 目标函数:最大化第2类捆数
f3 = zeros(num_vars, 1);
for i = 1:n
idx = x_start + 3*(i-1) + 1; % X_i2
f3(idx) = 1;
end
% 求解第三阶段
[x3, fval3, exitflag3] = intlinprog(-f3, intcon, A, b, [A_eq1; A_eq2], [b_eq1; b_eq2], lb, ub, options);
if exitflag3 <= 0
error('第三阶段求解失败');
end
spec2_bundles = round(-fval3);
spec1_bundles = total_bundles - spec3_bundles - spec2_bundles;
%% 输出最终结果
fprintf('\n=== 最终优化结果 ===\n');
fprintf('总捆数: %d\n', total_bundles);
fprintf('规格1捆数: %d\n', spec1_bundles);
fprintf('规格2捆数: %d\n', spec2_bundles);
fprintf('规格3捆数: %d\n', spec3_bundles);
% 计算原料使用率
used = zeros(46, 1);
for j = 1:46
for i = 1:n
% C_ij
idx = c_start + 46*(i-1) + j - 1;
used(j) = used(j) + x3(idx);
% D_ij^{(1)}
if ismember(j, downgrade1_idx)
d_idx = find(downgrade1_idx == j);
idx = d1_start + length(downgrade1_idx)*(i-1) + d_idx - 1;
used(j) = used(j) + x3(idx);
end
% D_ij^{(2)}
if ismember(j, downgrade2_idx)
d_idx = find(downgrade2_idx == j);
idx = d2_start + length(downgrade2_idx)*(i-1) + d_idx - 1;
used(j) = used(j) + x3(idx);
end
end
end
utilization = used ./ stocks';
fprintf('\n原料使用情况:\n');
disp(table((1:46)', lengths', stocks', used, utilization, ...
'VariableNames', {'原料ID', '长度', '库存', '使用量', '使用率'}));
第一阶段:最大化总捆数 ===
LP: Optimal objective value is 0.000000.
找到最优解。
Intlinprog 在根节点处停止,因为
目标值在最优值的间隙容差范围内,options.AbsoluteGapTolerance = 0。intcon 变量是
容差范围内的整数,options.IntegerTolerance = 1e-05。
第一阶段结果:最大总捆数 = 0
=== 第二阶段:最大化第3类捆数 ===
LP: Optimal objective value is 0.000000.
找到最优解。
Intlinprog 在根节点处停止,因为
目标值在最优值的间隙容差范围内,options.AbsoluteGapTolerance = 0。intcon 变量是
容差范围内的整数,options.IntegerTolerance = 1e-05。
第二阶段结果:第3类捆数 = 0
=== 第三阶段:最大化第2类捆数 ===
LP: Optimal objective value is 0.000000.
找到最优解。
Intlinprog 在根节点处停止,因为
目标值在最优值的间隙容差范围内,options.AbsoluteGapTolerance = 0。intcon 变量是
容差范围内的整数,options.IntegerTolerance = 1e-05。
=== 最终优化结果 ===
总捆数: 0
规格1捆数: 0
规格2捆数: 0
规格3捆数: 0
为什么都是0?
最新发布