388 lines
12 KiB
Matlab
388 lines
12 KiB
Matlab
%% 数据处理
|
||
% 指定文件名和工作表名
|
||
filename1 = '附件1.xlsx';
|
||
sheetname1 = '乡村的现有耕地';
|
||
filename2 = '附件2.xlsx';
|
||
sheetname2 = '2023年统计的相关数据';
|
||
|
||
% 读取数据
|
||
data_ranges = {'C2:C27', 'C28:C35', 'C36:C51', 'C52:C55'};
|
||
A = readmatrix(filename1, 'Sheet', sheetname1, 'Range', data_ranges{1});
|
||
B = readmatrix(filename1, 'Sheet', sheetname1, 'Range', data_ranges{2});
|
||
C = readmatrix(filename1, 'Sheet', sheetname1, 'Range', data_ranges{3});
|
||
D = readmatrix(filename1, 'Sheet', sheetname1, 'Range', data_ranges{4});
|
||
|
||
% 读取数值数据
|
||
data_ranges = {'F2:F16', 'F17:F31', 'F32:F46', 'F47:F65', 'F84:F86', 'F66:F83', 'F87:F90', 'F91:F108'};
|
||
data_labels = {'MP', 'MT', 'MS', 'MSb', 'MSbb', 'MDc', 'MDcc', 'MZdd'};
|
||
for i = 1:numel(data_ranges)
|
||
eval([data_labels{i} ' = readmatrix(filename2, ''Sheet'', sheetname2, ''Range'', data_ranges{i})'';']);
|
||
end
|
||
MZd = MDc; % 将MDc的值赋给MZd
|
||
|
||
% 读取G列的数据
|
||
data_ranges_G = {'G2:G16', 'G17:G31', 'G32:G46', 'G47:G65', 'G84:G86', 'G66:G83', 'G87:G90', 'G91:G108'};
|
||
data_labels_G = {'CP', 'CT', 'CS', 'CSb', 'CSbb', 'CDc', 'CDcc', 'CZdd'};
|
||
for i = 1:numel(data_ranges_G)
|
||
eval([data_labels_G{i} ' = readmatrix(filename2, ''Sheet'', sheetname2, ''Range'', data_ranges_G{i})'';']);
|
||
end
|
||
CZd = CDc; % 将CDc的值赋给CZd
|
||
|
||
% 读取H列的数据并处理
|
||
data_ranges_H = {'H32:H47', 'H48:H65', 'H84:H108'};
|
||
data_labels_H = {'SDJ', 'SDYJ', 'SDEJ'};
|
||
for k = 1:numel(data_ranges_H)
|
||
data = readcell(filename2, 'Sheet', sheetname2, 'Range', data_ranges_H{k});
|
||
values = zeros(length(data), 1);
|
||
for i = 1:length(data)
|
||
str = data{i};
|
||
if ismissing(str)
|
||
values(i) = NaN;
|
||
else
|
||
nums = str2double(split(str, '-'));
|
||
values(i) = mean(nums);
|
||
end
|
||
end
|
||
eval([data_labels_H{k} ' = values'';']);
|
||
end
|
||
|
||
|
||
% 指定文件名和工作表名
|
||
sheetname = '2023年的农作物种植情况';
|
||
|
||
% 读取整个B列和E列的数据
|
||
B_data = readmatrix(filename, 'Sheet', sheetname, 'Range', 'B2:B88');
|
||
E_data = readmatrix(filename, 'Sheet', sheetname, 'Range', 'E2:E88');
|
||
%% 任务1: 创建并填充a0矩阵
|
||
a0 = zeros(26, 15);
|
||
for i = 1:26
|
||
colIndex = B_data(i);
|
||
if colIndex >= 1 && colIndex <= 15
|
||
a0(i, colIndex) = E_data(i);
|
||
end
|
||
end
|
||
|
||
%% 任务2: 创建并填充b0矩阵
|
||
b0 = zeros(8, 19); % 行数增加到8以容纳最后的设置
|
||
evenRows = [28, 30, 32, 34, 36, 38];
|
||
for i = 1:numel(evenRows)
|
||
colIndex = B_data(evenRows(i) - 1) - 16;
|
||
if colIndex >= 1 && colIndex <= 19
|
||
b0(i, colIndex) = E_data(evenRows(i) - 1);
|
||
end
|
||
end
|
||
b0(7, 1) = 22; % 手动设置的值
|
||
b0(8, 1) = 20;
|
||
|
||
%% 任务3: 创建并填充bb0矩阵
|
||
bb0 = zeros(8, 3);
|
||
oddRows = [29, 31, 33, 35, 37, 39];
|
||
for i = 1:numel(oddRows)
|
||
colIndex = B_data(oddRows(i) - 1) - 34;
|
||
if colIndex >= 1 && colIndex <= 3
|
||
bb0(i, colIndex) = E_data(oddRows(i) - 1);
|
||
end
|
||
end
|
||
|
||
%% 任务4: 创建并填充c0矩阵
|
||
rowsC0 = [42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 71, 73];
|
||
B_data_c0 = readmatrix(filename, 'Sheet', sheetname, 'Range', ['B' num2str(min(rowsC0)) ':B' num2str(max(rowsC0))]);
|
||
E_data_c0 = readmatrix(filename, 'Sheet', sheetname, 'Range', ['E' num2str(min(rowsC0)) ':E' num2str(max(rowsC0))]);
|
||
|
||
% 初始化零矩阵c0，长16, 宽18
|
||
c0 = zeros(16, 18);
|
||
|
||
% 行映射，注意第15行用两次
|
||
rowMapC0 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 15, 16];
|
||
|
||
% 填充c0矩阵
|
||
for i = 1:numel(rowsC0)
|
||
colIndex = B_data_c0(i) - 16; % 计算横坐标
|
||
rowIndex = rowMapC0(i); % 获取映射的行索引
|
||
if colIndex >= 1 && colIndex <= 18
|
||
c0(rowIndex, colIndex) = E_data_c0(i); % 存储数据
|
||
end
|
||
end
|
||
|
||
%% 任务5: 创建并填充cc0矩阵
|
||
cc0 = zeros(16, 4);
|
||
selectedRowsCC0 = [43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 72, 74] - 1;
|
||
for i = 1:numel(selectedRowsCC0)
|
||
colIndex = B_data(selectedRowsCC0(i)) - 37;
|
||
if colIndex >= 1 && colIndex <= 4
|
||
cc0(i, colIndex) = E_data(selectedRowsCC0(i));
|
||
end
|
||
end
|
||
|
||
%% 任务6: 创建并填充d0矩阵
|
||
d0 = zeros(4, 18);
|
||
rowsD0 = [75, 76, 79, 80, 83, 86] - 1;
|
||
rowIndices = [1, 1, 2, 2, 3, 4]; % 映射到相应的行数
|
||
for i = 1:numel(rowsD0)
|
||
colIndex = B_data(rowsD0(i)) - 16;
|
||
if colIndex >= 1 && colIndex <= 18
|
||
d0(rowIndices(i), colIndex) = E_data(rowsD0(i));
|
||
end
|
||
end
|
||
|
||
%% 任务7: 创建并填充dd0矩阵
|
||
dd0 = zeros(4, 18);
|
||
rowsDD0 = [77, 78, 81, 82, 84, 85, 87, 88] - 1;
|
||
rowIndices = [1, 1, 2, 2, 3, 3, 4, 4]; % 映射到相应的行数
|
||
for i = 1:numel(rowsDD0)
|
||
colIndex = B_data(rowsDD0(i)) - 16;
|
||
if colIndex >= 1 && colIndex <= 18
|
||
dd0(rowIndices(i), colIndex) = E_data(rowsDD0(i));
|
||
end
|
||
end
|
||
|
||
%% 计算预期收获
|
||
GDJ = zeros(1, 16);
|
||
GDYJ = zeros(1, 18);
|
||
GDEJ = zeros(1, 25);
|
||
|
||
% 计算GDJ
|
||
GDJ(1:15) = MP .* sum(a0(1:6, :), 1) + MT .* sum(a0(7:20, :), 1) + MS .* sum(a0(21:26, :), 1);
|
||
GDJ(16) = MSb(1) * sum(b0(:, 1), 1);
|
||
|
||
% 计算GDYJ
|
||
GDYJ = MSb(2:19) .* sum(b0(:, 2:19), 1) + MDc .* sum(c0, 1) + MZd .* sum(d0, 1);
|
||
|
||
% 计算GDEJ
|
||
GDEJ(1:3) = MSbb .* sum(bb0, 1);
|
||
GDEJ(4:7) = MDcc .* sum(cc0, 1);
|
||
GDEJ(8:end) = MZdd .* sum(dd0, 1);
|
||
%% 变量设计
|
||
|
||
% 定义变量
|
||
a = sdpvar(26,15,8);
|
||
b = sdpvar(8,19,8);
|
||
bb = sdpvar(8,3,8);
|
||
c = sdpvar(16,18,8);
|
||
cc = sdpvar(16,4,8);
|
||
d = sdpvar(4,18,8);
|
||
dd = sdpvar(4,18,8);
|
||
b_temp = binvar(8,1,8);
|
||
|
||
ar = binvar(26,15,8);
|
||
br = binvar(8,19,8);
|
||
bbr = binvar(8,3,8);
|
||
cr = binvar(16,18,8);
|
||
ccr = binvar(16,4,8);
|
||
dr = binvar(4,18,8);
|
||
ddr = binvar(4,18,8);
|
||
|
||
mDJ = sdpvar(1,16,8);
|
||
mDYJ = sdpvar(1,18,8);
|
||
mDEJ = sdpvar(1,25,8);
|
||
cDJ = sdpvar(1,16,8);
|
||
cDYJ = sdpvar(1,18,8);
|
||
cDEJ = sdpvar(1,25,8);
|
||
|
||
% 约束设置
|
||
Constraints = [];
|
||
% 2023年数据载入
|
||
data = {a0, b0, bb0, c0, cc0, d0, dd0};
|
||
vars = {a, b, bb, c, cc, d, dd};
|
||
for i = 1:numel(data)
|
||
Constraints = [Constraints, vars{i}(:,:,1) == data{i}];
|
||
end
|
||
|
||
% 水浇地约束
|
||
Constraints = [Constraints, sum(bb,2) <= b_temp * 10e5];
|
||
Constraints = [Constraints, 10e-4 - b(:,1,2:8) <= b_temp(:,:,2:8) * 100];
|
||
Constraints = [Constraints, 100 - b(:,1,2:8) >= b_temp(:,:,2:8) * 100];
|
||
|
||
% 种植面积约束
|
||
areas = {a, b, bb, c, cc, d, dd};
|
||
bins = {ar, br, bbr, cr, ccr, dr, ddr};
|
||
limits = {A, B, B, C, C, D, D};
|
||
for i = 1:numel(areas)
|
||
Constraints = [Constraints, sum(bins{i},2) <= 1];
|
||
Constraints = [Constraints, 0.1 * bins{i}(:,:,2:8) <= areas{i}(:,:,2:8)];
|
||
Constraints = [Constraints, repmat(limits{i}, [1, size(bins{i},2), 7]) .* bins{i}(:,:,2:8) >= areas{i}(:,:,2:8)];
|
||
end
|
||
|
||
% 连续两年种植不同
|
||
crops = {ar, br, bbr, cr, ccr, dr, ddr};
|
||
for i = 1:numel(crops)
|
||
Constraints = [Constraints, crops{i}(:,:,1:7) + crops{i}(:,:,2:8) <= 1];
|
||
end
|
||
|
||
% 三年种一次豆类
|
||
for i = 1:6
|
||
Constraints = [Constraints, sum(a(:,1:5,i:i+2),2) >= A * 0.5];
|
||
Constraints = [Constraints, sum(b(:,2:4,i:i+2),2) >= B * 0.5];
|
||
Constraints = [Constraints, sum(c(:,1:3,i:i+2),2) >= C * 0.5];
|
||
Constraints = [Constraints, sum(d(:,1:3,i:i+2) + dd(:,1:3,i:i+2),2) >= D * 0.5];
|
||
end
|
||
|
||
|
||
%% 目标函数
|
||
% 收获(斤)
|
||
mDJ(1,1:15,:) = repmat(MP, [1, 1, 8]) .* sum(a(1:6,:,:),1) + ...
|
||
repmat(MT, [1, 1, 8]) .* sum(a(7:20,:,:),1) + ...
|
||
repmat(MS, [1, 1, 8]) .* sum(a(21:26,:,:),1);
|
||
mDJ(1,16,:) = repmat(MSb(1), [1, 1, 8]) .* sum(b(:,1,:),1);
|
||
mDYJ = repmat(MSb(2:19), [1, 1, 8]) .* sum(b(:,2:19,:),1) + ...
|
||
repmat(MDc, [1, 1, 8]) .* sum(c(:,:,:),1) + ...
|
||
repmat(MZd, [1, 1, 8]) .* sum(d(:,:,:),1);
|
||
mDEJ(1,1:3,:) = repmat(MSbb, [1, 1, 8]) .* sum(bb(:,:,:),1);
|
||
mDEJ(1,4:7,:) = repmat(MDcc, [1, 1, 8]) .* sum(cc(:,:,:),1);
|
||
mDEJ(1,8:end,:) = repmat(MZdd, [1, 1, 8]) .* sum(dd(:,:,:),1);
|
||
|
||
% 总成本
|
||
cDJ(1,1:15,:) = repmat(CP, [1, 1, 8]) .* sum(a(1:6,:,:),1) + ...
|
||
repmat(CT, [1, 1, 8]) .* sum(a(7:20,:,:),1) + ...
|
||
repmat(CS, [1, 1, 8]) .* sum(a(21:26,:,:),1);
|
||
cDJ(1,16,:) = repmat(CSb(1), [1, 1, 8]) .* sum(b(:,1,:),1);
|
||
cDYJ = repmat(CSb(2:19), [1, 1, 8]) .* sum(b(:,2:19,:),1) + ...
|
||
repmat(CDc, [1, 1, 8]) .* sum(c(:,:,:),1) + ...
|
||
repmat(CZd, [1, 1, 8]) .* sum(d(:,:,:),1);
|
||
cDEJ(1,1:3,:) = repmat(CSbb, [1, 1, 8]) .* sum(bb(:,:,:),1);
|
||
cDEJ(1,4:7,:) = repmat(CDcc, [1, 1, 8]) .* sum(cc(:,:,:),1);
|
||
cDEJ(1,8:end,:) = repmat(CZdd, [1, 1, 8]) .* sum(dd(:,:,:),1);
|
||
|
||
|
||
|
||
% 最终收益
|
||
Constraints = [Constraints, mDJ(:,:,2:8)<=repmat(GDJ, [1, 1, 7])*1.5];
|
||
Constraints = [Constraints, mDYJ(:,:,2:8)<=repmat(GDYJ, [1, 1, 7])*1.5];
|
||
Constraints = [Constraints, mDEJ(:,:,2:8)<=repmat(GDEJ, [1, 1, 7])*1.5];
|
||
|
||
|
||
OBJ2 = sum(sum(sum(mDJ(:,:,2:8) .* repmat(SDJ, [1, 1, 7])))) + ...
|
||
sum(sum(sum(mDYJ(:,:,2:8) .* repmat(SDYJ, [1, 1, 7])))) + ...
|
||
sum(sum(sum(mDEJ(:,:,2:8) .* repmat(SDEJ, [1, 1, 7])))) - ...
|
||
sum(sum(sum(cDJ(:,:,2:8)))) - ...
|
||
sum(sum(sum(cDYJ(:,:,2:8)))) - ...
|
||
sum(sum(sum(cDEJ(:,:,2:8))));
|
||
|
||
% 单目标求解
|
||
ops = sdpsettings(...
|
||
'solver', 'gurobi', ...
|
||
'gurobi.timelimit', 60, ...
|
||
'gurobi.MIPFocus', 1, ...
|
||
'gurobi.Heuristics', 0.5, ...
|
||
'gurobi.NodeLimit', 1000, ...
|
||
'allowmilp', 1, ...
|
||
'debug', 1, ...
|
||
'verbose', 0 ...
|
||
);
|
||
|
||
obj = -OBJ2;
|
||
result = optimize(Constraints, obj);
|
||
if result.problem == 0
|
||
disp('找到最优解');
|
||
elseif result.problem == 1
|
||
disp('未找到最优解');
|
||
else
|
||
disp('超时，显示当前最优解');
|
||
end
|
||
value(obj)
|
||
a = value(a);
|
||
b = value(b);
|
||
c = value(c);
|
||
d = value(d);
|
||
bb = value(bb);
|
||
cc = value(cc);
|
||
dd = value(dd);
|
||
%% 数据写入
|
||
|
||
% 定义工作表名称和矩阵数据
|
||
years = 2024:2030;
|
||
sheets = arrayfun(@num2str, years, 'UniformOutput', false);
|
||
|
||
% 创建 Excel 文件
|
||
filename = 'result1_1.xlsx';
|
||
|
||
for i = 1:length(years)
|
||
year = years(i);
|
||
sheetname = sheets{i};
|
||
|
||
% 创建矩阵数据
|
||
a1 = a(:,:,i);
|
||
b1 = b(:,:,i);
|
||
bb1 = bb(:,:,i);
|
||
c1 = c(:,:,i);
|
||
cc1 = cc(:,:,i);
|
||
d1 = d(:,:,i);
|
||
dd1 = dd(:,:,i);
|
||
|
||
% 写入 Excel 文件
|
||
% 1. 首先将 C2 到 AQ83 全部写为 0
|
||
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');
|
||
|
||
% 2. 将数据写入指定区域
|
||
xlswrite(filename, a1, sheetname, 'C2');
|
||
xlswrite(filename, b1, sheetname, 'R28');
|
||
xlswrite(filename, bb1, sheetname, 'AK56');
|
||
xlswrite(filename, c1, sheetname, 'R36');
|
||
xlswrite(filename, cc1, sheetname, 'AN36');
|
||
xlswrite(filename, d1, sheetname, 'R52');
|
||
xlswrite(filename, dd1, sheetname, 'R80');
|
||
end
|
||
|
||
|
||
% 收获 (斤) 和总成本
|
||
mDJ = repmat(MP, [1, 1, 8]).*sum(a(1:6,:,:),1) + ...
|
||
repmat(MT, [1, 1, 8]).*sum(a(7:20,:,:),1) + ...
|
||
repmat(MS, [1, 1, 8]).*sum(a(21:26,:,:),1);
|
||
mDJ(1,16,:) = repmat(MSb(1), [1, 1, 8]).*sum(b(:,1,:),1);
|
||
mDYJ = repmat(MSb(2:19), [1, 1, 8]).*sum(b(:,2:19,:),1) + ...
|
||
repmat(MDc, [1, 1, 8]).*sum(c(:,:,:),1) + ...
|
||
repmat(MZd, [1, 1, 8]).*sum(d(:,:,:),1);
|
||
mDEJ = repmat(MSbb, [1, 1, 8]).*sum(bb(:,:,:),1);
|
||
mDEJ(1,4:7,:) = repmat(MDcc, [1, 1, 8]).*sum(cc(:,:,:),1);
|
||
mDEJ(1,8:end,:) = repmat(MZdd, [1, 1, 8]).*sum(dd(:,:,:),1);
|
||
|
||
cDJ = repmat(CP, [1, 1, 8]).*sum(a(1:6,:,:),1) + ...
|
||
repmat(CT, [1, 1, 8]).*sum(a(7:20,:,:),1) + ...
|
||
repmat(CS, [1, 1, 8]).*sum(a(21:26,:,:),1);
|
||
cDJ(1,16,:) = repmat(CSb(1), [1, 1, 8]).*sum(b(:,1,:),1);
|
||
cDYJ = repmat(CSb(2:19), [1, 1, 8]).*sum(b(:,2:19,:),1) + ...
|
||
repmat(CDc, [1, 1, 8]).*sum(c(:,:,:),1) + ...
|
||
repmat(CZd, [1, 1, 8]).*sum(d(:,:,:),1);
|
||
cDEJ = repmat(CSbb, [1, 1, 8]).*sum(bb(:,:,:),1);
|
||
cDEJ(1,4:7,:) = repmat(CDcc, [1, 1, 8]).*sum(cc(:,:,:),1);
|
||
cDEJ(1,8:end,:) = repmat(CZdd, [1, 1, 8]).*sum(dd(:,:,:),1);
|
||
|
||
% 计算最终收益
|
||
benefits = 0;
|
||
|
||
% 定义收益计算函数
|
||
function b = calculate_benefits(m, G, S)
|
||
[dim1, dim2, dim3] = size(m);
|
||
b = 0;
|
||
for i = 1:dim1
|
||
for j = 1:dim2
|
||
for k = 1:dim3
|
||
if m(i,j,k) <= G(i,j,k)
|
||
b = b + m(i,j,k) * S(i,j,k);
|
||
else
|
||
b = b + G(i,j,k) * S(i,j,k) + (m(i,j,k) - G(i,j,k)) * S(i,j,k) / 2;
|
||
end
|
||
end
|
||
end
|
||
end
|
||
end
|
||
|
||
for i = 1:dim1
|
||
for j = 1:dim2
|
||
for k = 1:dim3
|
||
if mDEJ1(i,j,k) <= GDEJ1(i,j,k)
|
||
% 如果 mDJ 小于或等于 GDJ
|
||
benefits = benefits + mDEJ1(i,j,k) * SDEJ1(i,j,k);
|
||
else
|
||
% 如果 mDJ 大于 GDJ
|
||
benefits = benefits + GDEJ1(i,j,k) * SDEJ1(i,j,k) ...
|
||
+ (mDEJ1(i,j,k) - GDEJ1(i,j,k)) * 0;
|
||
end
|
||
end
|
||
end
|
||
end
|
||
|
||
benefits
|
||
|
||
|