### 问题1的整数规划模型框架与步骤
#### 1. 模型框架
我们建立混合整数线性规划（MILP）模型，核心框架如下：
- **目标函数**：最大化总利润 = 销售收入 - 种植成本
- **决策变量**：二元种植决策 + 连续面积分配
- **核心约束**：
- 耕地类型匹配
- 重茬约束（禁止连续种植相同作物）
- 三年豆类轮作
- 最小种植面积
- 产量与销售平衡
#### 2. 建模步骤
##### 步骤1：定义集合与参数
```matlab
% 基础集合
years = 2024:2030; % 规划年份
plots = land_data.地块名称; % 地块集合
crops = stats_data.作物名称; % 作物集合
seasons = ["第一季", "第二季", "单季"]; % 季节集合
legumes = ["黄豆","黑豆","红豆","绿豆","爬豆","豇豆","刀豆","芸豆"]; % 豆类作物
% 关键参数
min_area = 0.5; % 最小种植面积阈值(亩)
yield = containers.Map(); % 亩产量字典
cost = containers.Map(); % 种植成本字典
price = containers.Map(); % 销售价格字典
demand = containers.Map();% 预期销售量
```
##### 步骤2：参数计算（预处理）
```matlab
% 计算2023年各作物预期销售量（按季）
for c = crops
for s = seasons
mask = (planting_data.作物名称 == c) & (planting_data.种植季次 == s);
if any(mask)
area_planted = planting_data.种植面积_亩(mask);
plot_types = land_data.地块类型(ismember(land_data.地块名称, planting_data.种植地块(mask)));
% 获取对应亩产量
total_yield = 0;
for i = 1:sum(mask)
plot_type = plot_types(i);
y = stats_data.亩产量_斤(...
stats_data.作物名称 == c & stats_data.地块类型 == plot_type);
total_yield = total_yield + y * area_planted(i);
end
demand(strjoin([c,s], '_')) = total_yield;
else
demand(strjoin([c,s], '_')) = 0;
end
end
end
% 填充亩产量、成本、价格字典（使用中值）
for i = 1:height(stats_data)
crop = stats_data.作物名称{i};
plot_type = stats_data.地块类型{i};
key = strjoin([crop, plot_type], '_');
yield(key) = stats_data.亩产量_斤(i);
cost(key) = stats_data.种植成本_元_亩_(i);
% 处理价格区间
if ischar(stats_data.销售单价_元_斤_{i})
prices = sscanf(stats_data.销售单价_元_斤_{i}, '%f-%f');
price(key) = mean(prices);
else
price(key) = stats_data.销售单价_元_斤_{i};
end
end
```
##### 步骤3：定义决策变量
```matlab
% 二元决策变量：是否种植
z = optimvar('z', length(years), length(plots), length(seasons),...
length(crops), 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
% 连续决策变量：种植面积
a = optimvar('a', length(years), length(plots), length(seasons),...
length(crops), 'LowerBound', 0);
% 销售变量（情况1和2共用）
y_sale = optimvar('y_sale', length(years), length(crops),...
length(seasons), 'LowerBound', 0);
```
##### 步骤4：建立目标函数
```matlab
% 情况1：滞销模型
revenue1 = optimexpr;
for t = 1:length(years)
for c = 1:length(crops)
for s = 1:length(seasons)
revenue1 = revenue1 + price(crops{c}) * y_sale(t, c, s);
end
end
end
% 情况2：降价模型
revenue2 = optimexpr;
for t = 1:length(years)
for c = 1:length(crops)
for s = 1:length(seasons)
% 总产量计算
total_yield = optimexpr;
for i = 1:length(plots)
plot_type = land_data.地块类型(land_data.地块名称 == plots{i});
key = strjoin([crops{c}, plot_type], '_');
total_yield = total_yield + yield(key) * a(t, i, s, c);
end
revenue2 = revenue2 + 0.5 * price(crops{c}) * (total_yield + y_sale(t, c, s));
end
end
end
% 种植成本（共用）
cultivation_cost = optimexpr;
for t = 1:length(years)
for i = 1:length(plots)
for s = 1:length(seasons)
for c = 1:length(crops)
plot_type = land_data.地块类型(land_data.地块名称 == plots{i});
key = strjoin([crops{c}, plot_type], '_');
cultivation_cost = cultivation_cost + cost(key) * a(t, i, s, c);
end
end
end
end
% 最终目标函数
model1 = optimproblem('ObjectiveSense', 'maximize');
model1.Objective = revenue1 - cultivation_cost;
model2 = optimproblem('ObjectiveSense', 'maximize');
model2.Objective = revenue2 - cultivation_cost;
```
##### 步骤5：添加约束条件
```matlab
% 约束1：地块面积约束
area_constr = optimconstr(length(years), length(plots), length(seasons));
for t = 1:length(years)
for i = 1:length(plots)
for s = 1:length(seasons)
total_area = sum(a(t, i, s, :));
area_constr(t, i, s) = total_area <= land_data.地块面积_亩(i);
end
end
end
model1.Constraints.area_constr = area_constr;
model2.Constraints.area_constr = area_constr;
% 约束2：最小种植面积
min_area_constr = optimconstr(length(years), length(plots), length(seasons), length(crops));
for t = 1:length(years)
for i = 1:length(plots)
for s = 1:length(seasons)
for c = 1:length(crops)
min_area_constr(t, i, s, c) = a(t, i, s, c) >= min_area * z(t, i, s, c);
min_area_constr(t, i, s, c) = a(t, i, s, c) <= land_data.地块面积_亩(i) * z(t, i, s, c);
end
end
end
end
model1.Constraints.min_area_constr = min_area_constr;
model2.Constraints.min_area_constr = min_area_constr;
% 约束3：重茬约束
replant_constr = optimconstr(length(years)-1, length(plots), length(seasons), length(crops));
for t = 1:length(years)-1
for i = 1:length(plots)
for s = 1:length(seasons)
for c = 1:length(crops)
replant_constr(t, i, s, c) =...
z(t, i, s, c) + z(t+1, i, s, c) <= 1;
end
end
end
end
model1.Constraints.replant_constr = replant_constr;
model2.Constraints.replant_constr = replant_constr;
% 约束4：豆类轮作
legume_constr = optimconstr(length(years)-2, length(plots));
for t = 1:length(years)-2
for i = 1:length(plots)
legume_sum = optimexpr;
for k = t:t+2
for s = 1:length(seasons)
for c = 1:length(crops)
if ismember(crops{c}, legumes)
legume_sum = legume_sum + z(k, i, s, c);
end
end
end
end
legume_constr(t, i) = legume_sum >= 1;
end
end
model1.Constraints.legume_constr = legume_constr;
model2.Constraints.legume_constr = legume_constr;
% 约束5：耕地类型匹配
land_type_constr = optimconstr(length(years), length(plots), length(seasons), length(crops));
for t = 1:length(years)
for i = 1:length(plots)
plot_type = land_data.地块类型(land_data.地块名称 == plots{i});
for s = 1:length(seasons)
for c = 1:length(crops)
crop_type = stats_data.作物类型(stats_data.作物名称 == crops{c});
% 根据地块类型和作物类型设置约束
if plot_type == "平旱地" || plot_type == "梯田" || plot_type == "山坡地"
% 仅允许粮食作物，单季
if ~contains(crop_type, '粮食') || s ~= find(seasons == "单季")
land_type_constr(t, i, s, c) = z(t, i, s, c) == 0;
end
elseif plot_type == "水浇地"
% 特殊处理：水稻或蔬菜轮作
if crops{c} == "水稻"
% 水稻只能单季，且第二季不能种植
if s == find(seasons == "单季")
% 允许种植
else
land_type_constr(t, i, s, c) = z(t, i, s, c) == 0;
end
else
% 蔬菜作物：第一季任意蔬菜，第二季特定蔬菜
if s == find(seasons == "第二季") &&...
~ismember(crops{c}, ["大白菜","白萝卜","红萝卜"])
land_type_constr(t, i, s, c) = z(t, i, s, c) == 0;
end
end
elseif plot_type == "普通大棚"
% 第一季蔬菜，第二季食用菌
if s == find(seasons == "第一季") &&...
~contains(crop_type, '蔬菜')
land_type_constr(t, i, s, c) = z(t, i, s, c) == 0;
elseif s == find(seasons == "第二季") &&...
~contains(crop_type, '食用菌')
land_type_constr(t, i, s, c) = z(t, i, s, c) == 0;
end
elseif plot_type == "智慧大棚"
% 两季均蔬菜，禁种特定作物
if ~contains(crop_type, '蔬菜') ||...
ismember(crops{c}, ["大白菜","白萝卜","红萝卜"])
land_type_constr(t, i, s, c) = z(t, i, s, c) == 0;
end
end
end
end
end
end
model1.Constraints.land_type_constr = land_type_constr;
model2.Constraints.land_type_constr = land_type_constr;
% 约束6：销售平衡
sale_constr1 = optimconstr(length(years), length(crops), length(seasons));
sale_constr2 = optimconstr(length(years), length(crops), length(seasons));
for t = 1:length(years)
for c = 1:length(crops)
for s = 1:length(seasons)
key = strjoin([crops{c}, seasons{s}], '_');
% 产量计算
total_yield = optimexpr;
for i = 1:length(plots)
plot_type = land_data.地块类型(land_data.地块名称 == plots{i});
y_key = strjoin([crops{c}, plot_type], '_');
total_yield = total_yield + yield(y_key) * a(t, i, s, c);
end
% 销售量不能超过产量和预期需求
sale_constr1(t, c, s) = y_sale(t, c, s) <= total_yield;
sale_constr2(t, c, s) = y_sale(t, c, s) <= demand(key);
end
end
end
model1.Constraints.sale_constr1 = sale_constr1;
model1.Constraints.sale_constr2 = sale_constr2;
model2.Constraints.sale_constr1 = sale_constr1;
model2.Constraints.sale_constr2 = sale_constr2;
```
##### 步骤6：模型求解与结果输出
```matlab
% 求解情况1（滞销）
opts = optimoptions('intlinprog',...
'Display', 'iter',...
'MaxTime', 3600,...
'Heuristics', 'advanced',...
'CutGeneration', 'advanced');
[sol1, fval1] = solve(model1, 'Options', opts);
% 求解情况2（降价）
[sol2, fval2] = solve(model2, 'Options', opts);
% 结果提取与保存
results_table = table();
for t = 1:length(years)
for i = 1:length(plots)
for s = 1:length(seasons)
for c = 1:length(crops)
if sol1.a(t, i, s, c) > 0.01 % 忽略极小面积
new_row = {years(t), plots{i}, seasons{s},...
crops{c}, sol1.a(t, i, s, c)};
results_table = [results_table; new_row];
end
end
end
end
end
results_table.Properties.VariableNames =...
{'年份', '地块', '季节', '作物', '种植面积'};
% 保存到Excel
writetable(results_table, 'result_1.xlsx');
将这段代码优化一下，基于MATLAB（2024A）
最新发布
本文介绍了MATLAB中if-else-end指令的基本用法及其三种不同的结构形式，包括单分支、双分支及多分支结构，适用于不同场景下的条件判断。
扫一扫
258
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
