% 读取地块数据
filename1 = '附件1.xlsx';
landDataTable = readtable(filename1);
landData = table2array(landDataTable(2:end, 1:3));
landTypes = landData(:, 2);
landAreas = landData(:, 3);
% 提取不同类型地块的索引和面积
flatDryLandIndices = contains(landTypes, '平旱地');
terraceIndices = contains(landTypes, '梯田');
hillsideLandIndices = contains(landTypes, '山坡地');
irrigatedLandIndices = contains(landTypes, '水浇地');
ordinaryGreenhouseIndices = contains(landTypes, '普通大棚');
smartGreenhouseIndices = contains(landTypes, '智慧大棚');
flatDryLandAreas = landAreas(flatDryLandIndices);
terraceAreas = landAreas(terraceIndices);
hillsideLandAreas = landAreas(hillsideLandIndices);
irrigatedLandAreas = landAreas(irrigatedLandIndices);
ordinaryGreenhouseAreas = landAreas(ordinaryGreenhouseIndices);
smartGreenhouseAreas = landAreas(smartGreenhouseIndices);
% 读取农作物数据
filename2 = '附件2.xlsx';
cropDataTable = readtable(filename2, 'Sheet', 1, 'Range', 'A2:G');
cropIDs = table2array(cropDataTable(:, 2));
cropYields = table2array(cropDataTable(:, 5));
cropCosts = table2array(cropDataTable(:, 6));
cropPrices = table2array(cropDataTable(:, 7));
cropCategories = table2array(cropDataTable(:, 4));
% 提取不同作物的索引
beanCropIndices = find(contains(cropCategories, '豆类'));
grainCropIndices = find(contains(cropCategories, '粮食') & ~ismember(cropIDs, beanCropIndices));
vegetableCropIndices = find(contains(cropCategories, '蔬菜'));
mushroomCropIndices = find(contains(cropCategories, '食用菌'));
% 定义时间范围、季节范围
t_range = 2024:2030;
j_range = 1:2;
% 设定豆类作物索引,假设黄豆为豆类代表作物
k_bean = find(cropIDs == 1, 1);
% 设定最小种植面积阈值
min_area = 1;
% 设定作物种植分散程度方差上限(简化假设)
var_upper_limit = 10;
% 目标函数:情况(1)超过部分滞销
function obj = objectiveFunction(x, landData, cropData, landAreas, cropYields, cropCosts, cropPrices)
numCrops = length(unique(cropData(:, 2)));
numSeasons = length(j_range);
numLands = length(landAreas);
numYears = length(t_range);
x_ijtk = reshape(x(1:numCrops * numSeasons * numLands * numYears), numCrops, numSeasons, numLands, numYears);
z_ijtk = reshape(x(numCrops * numSeasons * numLands * numYears + 1:end), numCrops, numSeasons, numLands, numYears);
obj = 0;
for t = 1:numYears
for j = 1:numSeasons
for i = 1:numLands
for k = 1:numCrops
yield = cropYields((cropIDs == k) & (j == 1) & (landData(i, 2) == landData(:, 2)));
sale = cropData((cropIDs == k) & (j == 1), 8);
price = cropPrices((cropIDs == k) & (j == 1) & (landData(i, 2) == landData(:, 2)));
cost = cropCosts((cropIDs == k) & (j == 1) & (landData(i, 2) == landData(:, 2)));
obj = obj + price * min(yield * x_ijtk(k, j, i, t), sale) - cost * x_ijtk(k, j, i, t);
end
end
end
end
obj = -obj; % ga函数求最小值,原问题是求最大值,所以取负
end
function [c, ceq] = constraintFunction(x, landData, landAreas, min_area, var_upper_limit)
numCrops = length(unique(cropData(:, 2)));
numSeasons = length(j_range);
numLands = length(landAreas);
numYears = length(t_range);
x_ijtk = reshape(x(1:numCrops * numSeasons * numLands * numYears), numCrops, numSeasons, numLands, numYears);
z_ijtk = reshape(x(numCrops * numSeasons * numLands * numYears + 1:end), numCrops, numSeasons, numLands, numYears);
c = [];
ceq = [];
% 土地资源约束
for t = 1:numYears
for j = 1:numSeasons
c = [c; sum(x_ijtk(:, j, flatDryLandIndices, t)) - sum(flatDryLandAreas);
sum(x_ijtk(:, j, terraceIndices, t)) - sum(terraceAreas);
sum(x_ijtk(:, j, hillsideLandIndices, t)) - sum(hillsideLandAreas);
sum(x_ijtk(:, j, irrigatedLandIndices, t)) - sum(irrigatedLandAreas);
sum(x_ijtk(:, j, ordinaryGreenhouseIndices, t)) - sum(ordinaryGreenhouseAreas);
sum(x_ijtk(:, j, smartGreenhouseIndices, t)) - sum(smartGreenhouseAreas)];
end
end
% 各地块农作物种植面积下限
for t = 1:numYears
for j = 1:numSeasons
for i = 1:numLands
for k = 1:numCrops
c = [c; min_area * z_ijtk(k, j, i, t) - x_ijtk(k, j, i, t)];
end
end
end
end
% 决策变量的约束关系
for t = 1:numYears
for j = 1:numSeasons
for i = 1:numLands
for k = 1:numCrops
c = [c; x_ijtk(k, j, i, t) - landAreas(i) * z_ijtk(k, j, i, t)];
end
end
end
end
% 禁止连续重茬种植
for t = 1:numYears-1
for j = 1:numSeasons
for i = 1:numLands
for k = 1:numCrops
c = [c; x_ijtk(k, j, i, t + 1) - (1 - z_ijtk(k, j, i, t)) * landAreas(i)];
end
end
end
end
% 豆类作物种植要求
for i = 1:numLands
bean_count = 0;
for t = 1:numYears
for j = 1:numSeasons
bean_count = bean_count + z_ijtk(k_bean, j, i, t);
end
c = [c; floor(numYears/3) - bean_count];
end
% 各类地块的种植限制
% 平旱地、梯田、山坡地只能种粮食类作物(水稻除外)
for t = 1:numYears
for j = 1:numSeasons
for i = find(flatDryLandIndices | terraceIndices | hillsideLandIndices)
for k = setdiff(1:numCrops, [grainCropIndices; find(cropIDs == 16)])
c = [c; z_ijtk(k, j, i, t)];
end
end
end
end
% 水浇地种植限制
for t = 1:numYears
for i = find(irrigatedLandIndices)
% 第一季可种多种蔬菜(大白菜、白萝卜和红萝卜除外),第二季只能种大白菜、白萝卜和红萝卜中的一种
c = [c; sum(z_ijtk(setdiff(vegetableCropIndices, [35, 36, 37]), 1, i, t)) - 1;
1 - sum(z_ijtk([35, 36, 37], 2, i, t))];
end
end
% 普通大棚种植限制
for t = 1:numYears
for i = find(ordinaryGreenhouseIndices)
% 第一季可种多种蔬菜(大白菜、白萝卜和红萝卜除外),第二季只能种食用菌
c = [c; sum(z_ijtk(setdiff(vegetableCropIndices, [35, 36, 37]), 1, i, t)) - 1;
1 - sum(z_ijtk(mushroomCropIndices, 2, i, t))];
end
end
% 智慧大棚种植限制
for t = 1:numYears
for i = find(smartGreenhouseIndices)
% 只能种蔬菜(大白菜、白萝卜和红萝卜除外)
c = [c; sum(z_ijtk(setdiff(vegetableCropIndices, [35, 36, 37]), :, i, t)) - 2];
end
end
% 作物种植的分散程度限制(简化处理)
for t = 1:numYears
for j = 1:numSeasons
for k = 1:numCrops
x_k_vec = reshape(x_ijtk(k, j, :, t), 1, []);
var_x_k = var(x_k_vec);
c = [c; var_x_k - var_upper_limit];
end
end
end
end
% 决策变量个数
numVars = numCrops * length(j_range) * numLands * length(t_range) * 2;
% 设定选项
options = optimoptions('ga','PopulationSize',50,'Generations',100,'CrossoverFraction',0.8,'MutationFcn',@mutationuniform);
% 调用ga函数求解
[x, fval, exitflag, output, population, scores] = ga(@(x) objectiveFunction(x, landData, cropData, landAreas, cropYields, cropCosts, cropPrices),...
numVars,[],[],[],[],[],[],...
@(x) constraintFunction(x, landData, landAreas, min_area, var_upper_limit),[],options);
% 解码结果
x_ijtk_opt = reshape(x(1:numCrops * length(j_range) * numLands * length(t_range)), numCrops, length(j_range), numLands, length(t_range));
z_ijtk_opt = reshape(x(numCrops * length(j_range) * numLands * length(t_range) + 1:end), numCrops, length(j_range), numLands, length(t_range));
disp('最优解:');
disp(x_ijtk_opt);
disp(z_ijtk_opt);
对于以上代码,matlab说函数"objectiveFunction"以
'end'结束,但至少有一个其他函数定义没有以此结束。脚本中的所有函数都必须以'end’结束。你帮我改正,并输出改正后的全部代码。
最新发布