%% 分布式农产品生产系统优化
clear; clc; close all;
%% 1. 系统参数设置
fprintf('初始化分布式生产系统...\n');
% 系统结构
num_branch_factories = 3; % 3个分工厂
machines_per_branch = 2; % 每个分工厂2台机器
machines_central = 4; % 总工厂4台机器
% 生产参数
proc_time_S1 = 60; % S1生产时间(分钟)
proc_time_S2 = 30; % S2生产时间(分钟)
proc_time_S3 = 10; % S3生产时间(分钟)
% 变质时间
spoil_time_S1 = 240; % 4小时
spoil_time_S2 = 180; % 3小时
spoil_time_S3 = 120; % 2小时
% 运输参数(分工厂到总工厂)
transport_time_branch_to_central = [15, 20, 25]; % 三个分工厂到总工厂的运输时间
% 目标权重(所有产品同等重要)
weights = [1, 1, 1]; % S1, S2, S3
%% 2. 订单数据(分配到不同分工厂)
fprintf('生成订单数据...\n');
% 为每个分工厂生成订单
orders_per_factory = 6; % 进一步减少订单数以加快求解
total_orders = num_branch_factories * orders_per_factory;
orders = [];
order_id = 1;
for factory = 1:num_branch_factories
for i = 1:orders_per_factory
% 随机生成订单参数
order_time = 8 + rand() * 4; % 8:00-12:00之间
product_type = randi([1, 3]); % 随机产品类型
location_id = factory; % 所在分工厂
orders = [orders; order_id, factory, order_time, location_id, product_type];
order_id = order_id + 1;
end
end
orders_table = array2table(orders, 'VariableNames', ...
{'OrderID', 'FactoryID', 'OrderTime', 'LocationID', 'ProductType'});
%% 3. 创建优化问题
fprintf('创建分布式生产优化问题...\n');
prob = optimproblem('Description', '分布式生产优化', 'ObjectiveSense', 'max');
%% 4. 定义决策变量
fprintf('定义决策变量...\n');
% 分工厂生产决策:订单 × 分工厂机器
y_branch = optimvar('y_branch', total_orders, num_branch_factories * machines_per_branch, ...
'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
% 总工厂生产决策:订单 × 总工厂机器
y_central = optimvar('y_central', total_orders, machines_central, ...
'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
% 原料转移决策:订单 × 分工厂 → 总工厂
x_transfer = optimvar('x_transfer', total_orders, num_branch_factories, ...
'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
% 完成标志
z = optimvar('z', total_orders, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
%% 5. 设置目标函数:最大化总产量
fprintf('设置目标函数...\n');
objective = sum(z);
prob.Objective = objective;
%% 6. 添加约束条件
fprintf('添加约束条件...\n');
% 约束1: 每个订单只能在一个地方生产(分工厂或总工厂)
for i = 1:total_orders
branch_production = sum(y_branch(i, :));
central_production = sum(y_central(i, :));
prob.Constraints.(['single_production_', num2str(i)]) = ...
branch_production + central_production <= 1;
end
% 约束2: 分工厂机器容量约束
for factory = 1:num_branch_factories
machine_start = (factory-1) * machines_per_branch + 1;
machine_end = factory * machines_per_branch;
for machine = machine_start:machine_end
machine_load = sum(y_branch(:, machine));
prob.Constraints.(['branch_machine_cap_f', num2str(factory), '_m', num2str(machine)]) = ...
machine_load <= 1;
end
end
% 约束3: 总工厂机器容量约束
for machine = 1:machines_central
machine_load = sum(y_central(:, machine));
prob.Constraints.(['central_machine_cap_', num2str(machine)]) = ...
machine_load <= 1;
end
% 约束4: 原料转移逻辑约束
for i = 1:total_orders
factory_id = orders_table.FactoryID(i);
% 如果订单在总工厂生产,必须从对应分工厂转移原料
central_production = sum(y_central(i, :));
transfer_from_factory = x_transfer(i, factory_id);
prob.Constraints.(['transfer_required_', num2str(i)]) = ...
central_production <= transfer_from_factory;
end
% 约束5: 分工厂原料平衡约束
for factory = 1:num_branch_factories
% 该分工厂的订单
factory_orders = find(orders_table.FactoryID == factory);
% 分工厂本地生产的订单
machine_start = (factory-1) * machines_per_branch + 1;
machine_end = factory * machines_per_branch;
local_production = sum(sum(y_branch(factory_orders, machine_start:machine_end), 2));
% 转移到总工厂的订单
transfer_out = sum(x_transfer(factory_orders, factory));
% 约束:本地生产 + 转移出去 ≤ 总订单数
prob.Constraints.(['material_balance_f', num2str(factory)]) = ...
local_production + transfer_out <= length(factory_orders);
end
% 约束6: 变质时间约束 - 使用预计算的方式
fprintf('预计算变质时间约束...\n');
% 预计算哪些订单可能超时
time_constrained_orders = [];
for i = 1:total_orders
product_type = orders_table.ProductType(i);
order_time = orders_table.OrderTime(i) * 60; % 转换为分钟
factory_id = orders_table.FactoryID(i);
% 根据产品类型确定生产时间
if product_type == 1
production_time = proc_time_S1;
max_spoil_time = spoil_time_S1;
elseif product_type == 2
production_time = proc_time_S2;
max_spoil_time = spoil_time_S2;
else
production_time = proc_time_S3;
max_spoil_time = spoil_time_S3;
end
% 计算最坏情况(总工厂生产 + 运输时间)
worst_case_time = production_time + transport_time_branch_to_central(factory_id);
% 如果最坏情况都超时,强制该订单不能完成
if order_time + worst_case_time >= max_spoil_time
time_constrained_orders = [time_constrained_orders; i];
prob.Constraints.(['spoil_constraint_', num2str(i)]) = z(i) == 0;
end
end
fprintf('有 %d 个订单因变质时间限制无法在总工厂生产\n', length(time_constrained_orders));
% 约束7: 完成订单必须生产
for i = 1:total_orders
total_production = sum(y_branch(i, :)) + sum(y_central(i, :));
prob.Constraints.(['completion_requires_production_', num2str(i)]) = ...
z(i) <= total_production;
end
% 约束8: 转移原料必须对应总工厂生产
for i = 1:total_orders
factory_id = orders_table.FactoryID(i);
% 如果转移原料,必须在总工厂生产
prob.Constraints.(['transfer_implies_production_', num2str(i)]) = ...
x_transfer(i, factory_id) <= sum(y_central(i, :));
end
%% 7. 使用启发式方法生成好的初始解
fprintf('使用启发式方法生成初始分配方案...\n');
% 运行启发式算法获取好的初始解
heuristic_sol = distributed_heuristic(orders_table, num_branch_factories, ...
machines_per_branch, machines_central, transport_time_branch_to_central, ...
proc_time_S1, proc_time_S2, proc_time_S3, spoil_time_S1, spoil_time_S2, spoil_time_S3);
fprintf('启发式解完成订单: %d/%d\n', sum(heuristic_sol.z), total_orders);
%% 8. 快速求解
fprintf('开始求解分布式生产问题...\n');
% 修复优化选项参数
options = optimoptions('intlinprog', ...
'Display', 'iter', ... % 显示迭代过程
'MaxTime', 30, ... % 减少到30秒
'MaxNodes', 500, ... % 减少节点数
'Heuristics', 'advanced', ...
'CutGeneration', 'none', ... % 关闭割平面
'IntegerTolerance', 1e-4, ... % 使用有效范围内的值
'RelativeGapTolerance', 0.05); % 允许5%的间隙
tic;
try
[sol, fval, exitflag, output] = solve(prob, 'Options', options);
solve_time = toc;
% 检查求解结果
if exitflag <= 0
fprintf('优化求解失败,使用启发式解\n');
sol = heuristic_sol;
fval = sum(heuristic_sol.z);
exitflag = 99;
else
fprintf('优化求解成功\n');
end
catch ME
fprintf('求解过程中出现错误: %s\n', ME.message);
fprintf('使用启发式解\n');
sol = heuristic_sol;
fval = sum(heuristic_sol.z);
exitflag = 98;
solve_time = toc;
end
fprintf('求解时间: %.2f 秒\n', solve_time);
fprintf('总产量: %d 个订单\n', fval);
fprintf('退出标志: %d\n', exitflag);
%% 9. 结果分析
fprintf('\n=== 分布式生产结果分析 ===\n');
% 统计各工厂产量
[branch_production, central_production, transfer_count, total_production] = ...
analyze_production_results(sol, orders_table, num_branch_factories, machines_per_branch);
fprintf('\n各工厂生产统计:\n');
fprintf('工厂\t分厂产量\t总厂产量\t原料转移\t总产量\n');
fprintf('----\t--------\t--------\t--------\t------\n');
for factory = 1:num_branch_factories
factory_total = branch_production(factory) + central_production(factory);
fprintf('F%d\t%d\t\t%d\t\t%d\t\t%d\n', ...
factory, branch_production(factory), central_production(factory), ...
transfer_count(factory), factory_total);
end
fprintf('总计\t%d\t\t%d\t\t%d\t\t%d\n', ...
sum(branch_production), sum(central_production), sum(transfer_count), total_production);
%% 10. 详细分配结果
fprintf('\n=== 详细分配结果 ===\n');
fprintf('显示前8个完成订单的分配情况:\n');
fprintf('订单\t工厂\t产品\t生产地点\t机器\n');
fprintf('----\t----\t----\t--------\t----\n');
display_count = 0;
max_display = 8;
for i = 1:total_orders
if display_count >= max_display
break;
end
if sol.z(i) > 0.5 % 只显示完成的订单
factory_id = orders_table.FactoryID(i);
product_type = orders_table.ProductType(i);
% 确定生产地点和机器
[production_location, machine_id] = get_production_location(sol, i, factory_id, ...
num_branch_factories, machines_per_branch, machines_central);
product_str = get_product_string(product_type);
fprintf('O%d\tF%d\t%s\t%s\t\tM%d\n', ...
i, factory_id, product_str, production_location, machine_id);
display_count = display_count + 1;
end
end
%% 11. 性能分析
fprintf('\n=== 性能分析 ===\n');
completion_rate = total_production / total_orders * 100;
central_production_ratio = sum(central_production) / max(total_production, 1) * 100;
fprintf('系统总订单数: %d\n', total_orders);
fprintf('系统总产量: %d\n', total_production);
fprintf('系统完成率: %.1f%%\n', completion_rate);
fprintf('总工厂产量占比: %.1f%%\n', central_production_ratio);
fprintf('原料转移订单数: %d\n', sum(transfer_count));
% 机器利用率分析
branch_utilization = calculate_branch_utilization(sol, num_branch_factories, machines_per_branch);
central_utilization = calculate_central_utilization(sol, machines_central);
fprintf('\n机器利用率分析:\n');
fprintf('分工厂机器利用率: %.1f%%\n', branch_utilization);
fprintf('总工厂机器利用率: %.1f%%\n', central_utilization);
% 产品类型分析
fprintf('\n产品类型完成情况:\n');
for product_type = 1:3
type_orders = find(orders_table.ProductType == product_type);
type_completed = sum(sol.z(type_orders) > 0.5);
type_total = length(type_orders);
product_name = get_product_string(product_type);
fprintf('%s: %d/%d (%.1f%%)\n', product_name, type_completed, type_total, ...
type_completed/max(type_total,1)*100);
end
%% 12. 优化建议
fprintf('\n=== 优化建议 ===\n');
if completion_rate < 70
fprintf('⚠ 系统完成率较低,建议:\n');
fprintf(' - 检查变质时间约束是否过严\n');
fprintf(' - 考虑增加生产能力\n');
fprintf(' - 优化运输路线减少运输时间\n');
elseif completion_rate >= 90
fprintf('✓ 系统运行良好,完成率很高\n');
else
fprintf('✓ 系统运行正常,有改进空间\n');
end
if central_utilization < 50 && branch_utilization > 80
fprintf('💡 分工厂负载较高,建议增加向总工厂的原料转移\n');
elseif central_utilization > 80 && branch_utilization < 50
fprintf('💡 总工厂负载较高,建议充分利用分工厂生产能力\n');
end
if sum(transfer_count) == 0 && central_utilization < 30
fprintf('💡 未充分利用总工厂产能,建议适当转移原料\n');
end
%% 13. 简单可视化
fprintf('\n生成简单可视化...\n');
figure('Position', [100, 100, 1000, 600]);
% 子图1: 产量分布
subplot(2, 2, 1);
production_data = [branch_production, central_production];
bar(1:num_branch_factories, production_data, 'stacked');
xlabel('分工厂');
ylabel('产量(订单数)');
title('各工厂产量分布');
legend({'分工厂生产', '总工厂生产'}, 'Location', 'best');
grid on;
% 子图2: 完成率分析
subplot(2, 2, 2);
completion_data = [total_production, total_orders - total_production];
pie(completion_data, {'已完成', '未完成'});
title(sprintf('总完成率: %.1f%%', completion_rate));
% 子图3: 机器利用率
subplot(2, 2, 3);
utilization_data = [branch_utilization, central_utilization];
bar(1:2, utilization_data);
set(gca, 'XTickLabel', {'分工厂', '总工厂'});
ylabel('利用率 (%)');
title('机器利用率对比');
ylim([0, 100]);
grid on;
% 子图4: 生产地点选择
subplot(2, 2, 4);
location_data = [sum(branch_production), sum(central_production)];
pie(location_data, {'分工厂生产', '总工厂生产'});
title('生产地点分布');
fprintf('\n分布式生产优化完成!\n');
%% 辅助函数定义
function sol = distributed_heuristic(orders_table, num_branch_factories, ...
machines_per_branch, machines_central, transport_times, ...
proc_time_S1, proc_time_S2, proc_time_S3, spoil_time_S1, spoil_time_S2, spoil_time_S3)
total_orders = height(orders_table);
% 初始化解
sol.y_branch = zeros(total_orders, num_branch_factories * machines_per_branch);
sol.y_central = zeros(total_orders, machines_central);
sol.x_transfer = zeros(total_orders, num_branch_factories);
sol.z = zeros(total_orders, 1);
% 按工厂和产品类型排序(S1优先)
[~, order_idx] = sort(orders_table.ProductType);
% 资源状态
branch_machines_used = zeros(num_branch_factories * machines_per_branch, 1);
central_machines_used = zeros(machines_central, 1);
for idx = 1:total_orders
i = order_idx(idx);
factory_id = orders_table.FactoryID(i);
product_type = orders_table.ProductType(i);
order_time = orders_table.OrderTime(i) * 60;
% 确定生产时间
if product_type == 1
production_time = proc_time_S1;
spoil_time = spoil_time_S1;
elseif product_type == 2
production_time = proc_time_S2;
spoil_time = spoil_time_S2;
else
production_time = proc_time_S3;
spoil_time = spoil_time_S3;
end
% 先尝试在分工厂生产
assigned = false;
machine_start = (factory_id-1) * machines_per_branch + 1;
machine_end = factory_id * machines_per_branch;
for machine = machine_start:machine_end
if branch_machines_used(machine) == 0
% 检查变质时间
total_time = production_time;
if order_time + total_time <= spoil_time
sol.y_branch(i, machine) = 1;
branch_machines_used(machine) = 1;
sol.z(i) = 1;
assigned = true;
break;
end
end
end
% 如果分工厂无法生产,尝试总工厂
if ~assigned
for machine = 1:machines_central
if central_machines_used(machine) == 0
% 检查变质时间(包括运输时间)
total_time = production_time + transport_times(factory_id);
if order_time + total_time <= spoil_time
sol.y_central(i, machine) = 1;
sol.x_transfer(i, factory_id) = 1;
central_machines_used(machine) = 1;
sol.z(i) = 1;
assigned = true;
break;
end
end
end
end
end
end
function [branch_production, central_production, transfer_count, total_production] = ...
analyze_production_results(sol, orders_table, num_branch_factories, machines_per_branch)
total_orders = height(orders_table);
branch_production = zeros(num_branch_factories, 1);
central_production = zeros(num_branch_factories, 1);
transfer_count = zeros(num_branch_factories, 1);
for factory = 1:num_branch_factories
factory_orders = find(orders_table.FactoryID == factory);
% 分工厂本地生产
machine_start = (factory-1) * machines_per_branch + 1;
machine_end = factory * machines_per_branch;
for i = factory_orders
for machine = machine_start:machine_end
if sol.y_branch(i, machine) > 0.5
branch_production(factory) = branch_production(factory) + 1;
break;
end
end
end
% 总工厂生产
for i = factory_orders
for machine = 1:size(sol.y_central, 2)
if sol.y_central(i, machine) > 0.5
central_production(factory) = central_production(factory) + 1;
break;
end
end
end
% 原料转移数量
for i = factory_orders
if sol.x_transfer(i, factory) > 0.5
transfer_count(factory) = transfer_count(factory) + 1;
end
end
end
total_production = sum(branch_production) + sum(central_production);
end
function [production_location, machine_id] = get_production_location(sol, order_id, factory_id, ...
num_branch_factories, machines_per_branch, machines_central)
production_location = '无';
machine_id = 0;
% 检查分工厂生产
machine_start = (factory_id-1) * machines_per_branch + 1;
machine_end = factory_id * machines_per_branch;
for machine = machine_start:machine_end
if sol.y_branch(order_id, machine) > 0.5
production_location = '分工厂';
machine_id = machine - machine_start + 1;
return;
end
end
% 检查总工厂生产
for machine = 1:machines_central
if sol.y_central(order_id, machine) > 0.5
production_location = '总工厂';
machine_id = machine;
return;
end
end
end
function product_str = get_product_string(product_type)
if product_type == 1
product_str = 'S1';
elseif product_type == 2
product_str = 'S2';
else
product_str = 'S3';
end
end
function utilization = calculate_branch_utilization(sol, num_branch_factories, machines_per_branch)
total_machines = num_branch_factories * machines_per_branch;
used_machines = 0;
for machine = 1:total_machines
if any(sol.y_branch(:, machine) > 0.5)
used_machines = used_machines + 1;
end
end
utilization = used_machines / total_machines * 100;
end
function utilization = calculate_central_utilization(sol, machines_central)
used_machines = 0;
for machine = 1:machines_central
if any(sol.y_central(:, machine) > 0.5)
used_machines = used_machines + 1;
end
end
utilization = used_machines / machines_central * 100;
end
代码纠错
最新发布