24高教社杯数学建模国赛C题成品论文
一、问题一模型建立与求解
1.1模型建立
(1)决策变量设计
表示一个26×15×8的矩阵,其中26是平旱地梯田和山坡地的总数,15是在这几类土地上可以种植的农作物数量,8则表示从2023到2030这8年时间。上述决策变量可以反映平旱地、梯田和山坡地在这8年时间里的种植情况。
表示水浇地在8年时间里第一季的种植情况,19种农作物种包含了水稻,且其排在第一位。
表示水浇地在8年时间里第二季的种植情况,仅包含三种农作物。
分别表示普通大棚在8年时间里第一季和第二季的种植情况。
别表示智慧大棚在8年时间里第一季和第二季的种植情况。
(2)约束变量设计
本题的约束条件主要有一般的边界条件以及题目给出的下述约束,逐个建立数学表达式即可:
1. 平旱地、梯田和山坡地每年适宜单季种植粮食类作物(水稻除外)。由于设置了变量,可以保证平旱地、梯田和山坡仅单季种植粮食类作物。
2. 水浇地每年可以单季种植水稻或两季种植蔬菜作物。针对这一条件,首先引入中间0-1变量,然后进行如下表达:
上述公式可以保证当中存在某数字大于0,则bb全为0。其物理意义即种植水稻后不再种植其他植物。
3. 若在某块水浇地种植两季蔬菜,第一季可种植多种蔬菜(大白菜、白萝卜和红萝卜除外+);第二季只能种植大白菜、白萝卜和红萝卜中的一种(便于管理)。变量b和bb可以满足上述要求。
4. 根据季节性要求,大白菜、白萝卜和红萝卜只能在水浇地的第二季种植。变量b和bb可以满足上述要求。
5. 普通大棚每年种植两季作物,第一季可种植多种蔬菜(大白菜、白萝卜和红萝卜除外),第二季只能种植食用菌。变量c和cc可以满足上述要求。
6. 因食用菌类适应在较低且适宜的温度和湿度环境中生长,所以只能在秋冬季的普通大棚里种植。变量c和cc可以满足上述要求.
7. 智慧大棚每年都可种植两季蔬菜(大白菜、白萝卜和红萝卜除外)。变量d和dd可以满足上述要求。
8. 普通大棚适宜每年种植一季蔬菜和一季食用菌,智慧大棚适宜每年种植两季蔬菜。变量c,cc,d,dd可以满足上述要求.
9. 同一地块(含大棚) 每季可以合种不同的作物。详见附件 1。
10. 种植方案应考虑到方便耕种作业和田间管理, 譬如: 每种作物每季的种植地不能太分散, 每种作物在单个地块(含大棚) 种植的面积不宜太小。
注意,由附件1中可以得到,一般一块田地不会种植超过两种植物,且由于要求作物不能连续种植以及三年要种植一次豆类。因此在同意土地上种植多种植物极其不方便管理且容易导致减产。因此本模型假设一般田地均只种植同类型植物。在此基础上,引入下面中间0-1变量:
并设置如下约束:
即可以保证每块土地在某一时间段内仅种植一种植物。如果想种植多种,仅需将公式中的1修改为其他整数。
11. 每种作物在同一地块(含大棚)都不能连续重茬种植,否则会减产。
上式可保证连续两年内同一块土地不会种植相同的植物。
12. 所有土地三年内至少种植一次豆类作物。
在此,首先假设种植的豆类植物应当占据该土地面积的一半以上,则有如下约束
其中ABCD分别表示四种土地的最大种植面积,A包含了平旱地梯田和山坡地三种土地。
13. 最后,是种植面积的基本上下界,假设至少种植0.1亩地,则有如下约束:
分别表示了各种植物在不同时间不同土地的种植上下界约束。
(3)目标函数设计
首先,引入下图9个常量
分别表示平旱地、梯田、山坡地、水浇地第一季、水浇地第二季、普通大棚第一季、普通大棚第二季、智慧大棚第一季、智慧大棚第二季的亩产量。
上述常量则表示平旱地、梯田、山坡地、水浇地第一季、水浇地第二季、普通大棚第一季、普通大棚第二季、智慧大棚第一季、智慧大棚第二季的种植成本。
上述则表示平旱地、梯田、山坡地、水浇地第一季、水浇地第二季、普通大棚第一季、普通大棚第二季、智慧大棚第一季、智慧大棚第二季的销售单价。注意,本文中的销售单价均取平均值。
接着,计算模型的整体产量为。目标函数1为:
其中为单季各作物产量,为第一季产量,为第二季产量。目标函数2为:
同理,单季、第一季和第二季的成本可按照如上计算。
由于附件未直接提供预期销售量,因此以2023年产量作为预期销售量,并分别标记为
将单季、第一季、第二季的价格标记为
引入中间变量,使其满足
即可保证用于计算的产量值取预期销售额和最终产量的最小值,以此计算最终收益,即目标函数如下:
1.2模型求解与分析
(1)数据处理
针对题目提供的数据,首先做如下处理,以方便matlab计算:
- 读取表格“附件1.xlsx”中“乡村的现有耕地”这一页,将C2到C27的数据存储在变量A,C28到C35的数据存储在变量B,C36到C51的数据存储在变量C,C52到C55的数据存储在变量D。
- 读取表格“附件2.xlsx”中”2023年统计的相关数据"这一页,将F2到F16的数据存储在MP,F17到F31的数据存储在MT,F32到F45的数据存储在MS,F46到F65的数据存储在MSb,F84到F86的数据存储在MSbb,F66到F83的数据存储在MDc,F87到F90的数据存储在MDcc,F91到F108的数据存储在MZdd,最后,令MZd=MDc。
- 读取表格“附件2.xlsx”中”2023年统计的相关数据"这一页,将F2到F16的数据存储在MP,F17到F31的数据存储在MT,F32到F45的数据存储在MS,F46到F65的数据存储在MSb,F84到F86的数据存储在MSbb,F66到F83的数据存储在MDc,F87到F90的数据存储在MDcc,F91到F108的数据存储在MZdd,最后,令MZd=MDc。
- 读取表格“附件2.xlsx”中”2023年统计的相关数据"这一页,将H32到H47的数据处理后存储在SDJ,将H48到H65的数据处理后存储在SDYJ,将H84到H108的数据处理后存储在SDEJ。处理方式为,将单元格内的字符串中‘-’之前的字符串和‘-’之后的字符串均转化为浮点数并取均值。
- 读取表格“附件2.xlsx”中”2023年的农作物种植情况"这一页,进行如下操作:
- 针对第2行到第27行,设置一个零矩阵a0,长26,宽15,以1到26为纵坐标,B2到B27中的数值为横坐标,将E2到E27的数值存储进a0。
- 针对第28,30,32,34,36,38行,首先设置一个零矩阵b0,长8,宽19,以1到6为纵坐标,B列中的数值减去16为横坐标,将对应的E列的数值存储进b0。最后,设置b0(7,1)=22, b0(8,1)=20.
- 针对第29,31,33,35,37,39行,首先设置一个零矩阵bb0,长8,宽3,以1到8为纵坐标,B列中奇数行的数值减去34为横坐标,将对应的E列的数值存储进bb0。
- 针对第42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,71,73行,首先设置一个零矩阵c0,长16,宽18,以1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,15,16分别为纵坐标,B列中的数值减去16为横坐标,将对应的E列的数值存储进c0。
- 针对第43,45,47,49,51,53,55,57,59,61,63,65,67,69,72,74行,首先设置一个零矩阵cc0,长16,宽4,以1到16为纵坐标,B列中的数值减去37为横坐标,将对应的E列的数值存储进cc0。
- 针对75,76,79,80,83,86行,首先设置一个零矩阵d0,长4,宽18,以1,1,2,2,3,4,4为纵坐标,B列中的数值减去16为横坐标,将对应的E列的数值存储进d0。
- 针对77,78,81,82,84,85,87,88行,首先设置一个零矩阵dd0,长4,宽18,以1,1,2,2,3,3,4,4为纵坐标,B列中的数值减去16为横坐标,将对应的E列的数值存储进dd0。
(2)代码计算
以matlab结合gurubi直接求得最终解。注意,由于目标函数过于复杂,因此在求解过程中将其松弛为
同时,引入如下约束以确保种植农作物数量不会超出预期销售量太多,α为松弛系数,针对第一小问和第二小问分别设置为1和1.5。
(3)结果展示分析
当超出部分没有收入时的总收入为:4.4269e+07。当超出部分有百分之五十收入时为:4.6788e+07。具体原因在于由于超出部分具有收入,因此会选针对一些经济价值更高的农作物进行种植,以提高最终的收入。具体的结果在表格中可见。也可以考虑增设一些图片。
1.3代码演示
%% 数据处理
% 指定文件名和工作表名
filename1 = '附件1.xlsx';
sheetname1 = '乡村的现有耕地';% 读取数据
A = readmatrix(filename1, 'Sheet', sheetname1, 'Range', 'C2:C27');
B = readmatrix(filename1, 'Sheet', sheetname1, 'Range', 'C28:C35');
C = readmatrix(filename1, 'Sheet', sheetname1, 'Range', 'C36:C51');
D = readmatrix(filename1, 'Sheet', sheetname1, 'Range', 'C52:C55');% 指定文件名和工作表名
filename2 = '附件2.xlsx';
sheetname2 = '2023年统计的相关数据';% 读取数据
MP = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F2:F16')';
MT = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F17:F31')';
MS = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F32:F46')';
MSb = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F47:F65')';
MSbb = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F84:F86')';
MDc = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F66:F83')';
MDcc = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F87:F90')';
MZdd = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'F91:F108')';
MZd = MDc; % 将MDc的值赋给MZd% 读取G列的数据
CP = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G2:G16')';
CT = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G17:G31')';
CS = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G32:G46')';
CSb = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G47:G65')';
CSbb = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G84:G86')';
CDc = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G66:G83')';
CDcc = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G87:G90')';
CZdd = readmatrix(filename2, 'Sheet', sheetname2, 'Range', 'G91:G108')';
CZd = CDc; % 将CDc的值赋给CZd% 读取H列的数据
data_SDJ = readcell(filename2, 'Sheet', sheetname2, 'Range', 'H32:H47');
data_SDYJ = readcell(filename2, 'Sheet', sheetname2, 'Range', 'H48:H65');
data_SDEJ = readcell(filename2, 'Sheet', sheetname2, 'Range', 'H84:H108');% 初始化存储平均值的数组
SDJ = zeros(length(data_SDJ), 1);
SDYJ = zeros(length(data_SDYJ), 1);
SDEJ = zeros(length(data_SDEJ), 1);% 处理SDJ数据
for i = 1:length(data_SDJ)str = data_SDJ{i}; % 获取单元格数据if ismissing(str)SDJ(i) = NaN; % 处理缺失数据elseparts = split(str, '-'); % 按'-'分割字符串nums = str2double(parts); % 将分割后的字符串转换为数字SDJ(i) = mean(nums); % 计算平均值end
end% 处理SDYJ数据
for i = 1:length(data_SDYJ)str = data_SDYJ{i};if ismissing(str)SDYJ(i) = NaN;elseparts = split(str, '-');nums = str2double(parts);SDYJ(i) = mean(nums);end
end% 处理SDEJ数据
for i = 1:length(data_SDEJ)str = data_SDEJ{i};if ismissing(str)SDEJ(i) = NaN;elseparts = split(str, '-');nums = str2double(parts);SDEJ(i) = mean(nums);end
endSDJ = SDJ';
SDYJ = SDYJ';
SDEJ = SDEJ';% 指定文件名和工作表名
filename = '附件2.xlsx';
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:26colIndex = B_data(i);if colIndex >= 1 && colIndex <= 15a0(i, colIndex) = E_data(i);end
end% 任务2: 创建并填充b0矩阵
b0 = zeros(8, 19); % 行数增加到8以容纳最后的设置
evenRows = [28, 30, 32, 34, 36, 38];
for i = 1:6colIndex = B_data(evenRows(i) - 1) - 16;if colIndex >= 1 && colIndex <= 19b0(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:6colIndex = B_data(oddRows(i) - 1) - 34;if colIndex >= 1 && colIndex <= 3bb0(i, colIndex) = E_data(oddRows(i) - 1);end
end% 任务4: 创建并填充c0矩阵
filename = '附件2.xlsx';
sheetname = '2023年的农作物种植情况';% 读取B列和E列的指定行数据
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:length(rowsC0)colIndex = B_data_c0(i) - 16; % 计算横坐标rowIndex = rowMapC0(i); % 获取映射的行索引if colIndex >= 1 && colIndex <= 18c0(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:length(selectedRowsCC0)colIndex = B_data(selectedRowsCC0(i)) - 37;if colIndex >= 1 && colIndex <= 4cc0(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:length(rowsD0)colIndex = B_data(rowsD0(i)) - 16;if colIndex >= 1 && colIndex <= 18d0(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:length(rowsDD0)colIndex = B_data(rowsDD0(i)) - 16;if colIndex >= 1 && colIndex <= 18dd0(rowIndices(i), colIndex) = E_data(rowsDD0(i));end
endGDJ =zeros(1,16);
GDYJ =zeros(1,18);
GDEJ =zeros(1,25);
% 预期收获(斤)
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 = MSb(2:19).*sum(b0(:,2:19),1) + MDc.*sum(c0(:,:),1) + MZd.*sum(d0(:,:),1);
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);
QDJ= sdpvar(1,16,8);
QDYJ= sdpvar(1,18,8);
QDEJ= sdpvar(1,25,8);% = sdpvar(,,);
% = sdpvar(,,);% 约束设置
Constraints = [];% 2023年数据载入
Constraints = [Constraints,a(:,:,1) == a0];
Constraints = [Constraints,b(:,:,1) == b0];
Constraints = [Constraints,bb(:,:,1) == bb0];
Constraints = [Constraints,c(:,:,1) == c0];
Constraints = [Constraints,cc(:,:,1) == cc0];
Constraints = [Constraints,d(:,:,1) == d0];
Constraints = [Constraints,dd(:,:,1) == dd0];% % 水浇地每年可以单季种植水稻或两季种植蔬菜作物。
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];% 种植面积约束
Constraints = [Constraints, sum(ar,2) <= 1];
Constraints = [Constraints, 0.1*ar(:,:,2:8)<=a(:,:,2:8)];
Constraints = [Constraints, repmat(A, [1, 15, 7]).*ar(:,:,2:8)>=a(:,:,2:8)];Constraints = [Constraints, sum(br,2) <= 1];
Constraints = [Constraints, 0.1*br(:,:,2:8)<=b(:,:,2:8)];
Constraints = [Constraints, repmat(B, [1, 19, 7]).*br(:,:,2:8)>=b(:,:,2:8)];Constraints = [Constraints, sum(bbr,2) <= 1];
Constraints = [Constraints, 0.1*bbr(:,:,2:8)<=bb(:,:,2:8)];
Constraints = [Constraints, repmat(B, [1, 3, 7]).*bbr(:,:,2:8)>=bb(:,:,2:8)];
Constraints = [Constraints, sum(cr,2) <= 1];
Constraints = [Constraints, 0.1*cr(:,:,2:8)<=c(:,:,2:8)];
Constraints = [Constraints, repmat(C, [1, 18, 7]).*cr(:,:,2:8)>=c(:,:,2:8)];
Constraints = [Constraints, sum(ccr,2) <= 1];
Constraints = [Constraints, 0.1*ccr(:,:,2:8)<=cc(:,:,2:8)];
Constraints = [Constraints, repmat(C, [1, 4, 7]).*ccr(:,:,2:8)>=cc(:,:,2:8)];Constraints = [Constraints, sum(dr,2) <= 1];
Constraints = [Constraints, 0.1*dr(:,:,2:8)<=d(:,:,2:8)];
Constraints = [Constraints, repmat(D, [1, 18, 7]).*dr(:,:,2:8)>=d(:,:,2:8)];Constraints = [Constraints, sum(ddr,2) <= 1];
Constraints = [Constraints, 0.1*ddr(:,:,2:8)<=dd(:,:,2:8)];
Constraints = [Constraints, repmat(D, [1, 18, 7]).*ddr(:,:,2:8)>=dd(:,:,2:8)];% % 初始种植面积约束
% Constraints = [Constraints, ar(:,:,1)>=a(:,:,1)/100+0.99999];
% Constraints = [Constraints, ar(:,:,1)<=a(:,:,1)/100+0.99999];
% Constraints = [Constraints, br(:,:,1)>=b(:,:,1)/100+0.99999];
% Constraints = [Constraints, br(:,:,1)<=b(:,:,1)/100+0.99999];
% Constraints = [Constraints, bbr(:,:,1)>=bb(:,:,1)/100+0.99999];
% Constraints = [Constraints, bbr(:,:,1)<=bb(:,:,1)/100+0.99999];
% Constraints = [Constraints, cr(:,:,1)>=c(:,:,1)/100+0.99999];
% Constraints = [Constraints, cr(:,:,1)<=c(:,:,1)/100+0.99999];
% Constraints = [Constraints, ccr(:,:,1)>=cc(:,:,1)/100+0.99999];
% Constraints = [Constraints, ccr(:,:,1)<=cc(:,:,1)/100+0.99999];
% Constraints = [Constraints, dr(:,:,1)>=d(:,:,1)/100+0.99999];
% Constraints = [Constraints, dr(:,:,1)<=d(:,:,1)/100+0.99999];
% Constraints = [Constraints, ddr(:,:,1)>=dd(:,:,1)/100+0.99999];
% Constraints = [Constraints, ddr(:,:,1)<=dd(:,:,1)/100+0.99999];% 连续两年种植不同
Constraints = [Constraints, ar(:,:,1:7)+ar(:,:,2:8)<=1];
Constraints = [Constraints, br(:,:,1:7)+br(:,:,2:8)<=1];
Constraints = [Constraints, bbr(:,:,1:7)+bbr(:,:,2:8)<=1];
Constraints = [Constraints, cr(:,:,1:7)+cr(:,:,2:8)<=1];
Constraints = [Constraints, ccr(:,:,1:7)+ccr(:,:,2:8)<=1];
Constraints = [Constraints, dr(:,:,1:7)+dr(:,:,2:8)<=1];
Constraints = [Constraints, ddr(:,:,1:7)+ddr(:,:,2:8)<=1];% 三年种一次豆类for i =1:6Constraints = [Constraints, sum(a(:,1:5,i)+a(:,1:5,i+1)+a(:,1:5,i+2),2)>=A*0.5];Constraints = [Constraints, sum(b(:,2:4,i)+b(:,2:4,i+1)+b(:,2:4,i+2),2)>=B*0.5];Constraints = [Constraints, sum(c(:,1:3,i)+c(:,1:3,i+1)+c(:,1:3,i+2),2)>=C*0.5];Constraints = [Constraints, sum(d(:,1:3,i)+d(:,1:3,i+1)+d(:,1:3,i+2)+dd(:,1:3,i)+dd(:,1:3,i+1)+dd(:,1:3,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, QDJ(:,:,2:8)<=mDJ(:,:,2:8)];
% Constraints = [Constraints, QDJ(:,:,2:8)<=repmat(GDJ, [1, 1, 7])];
% Constraints = [Constraints, QDYJ(:,:,2:8)<=mDYJ(:,:,2:8)];
% Constraints = [Constraints, QDYJ(:,:,2:8)<=repmat(GDYJ, [1, 1, 7])];
% Constraints = [Constraints, QDEJ(:,:,2:8)<=mDEJ(:,:,2:8)];
% Constraints = [Constraints, QDEJ(:,:,2:8)<=repmat(GDEJ, [1, 1, 7])];Constraints = [Constraints, mDJ(:,:,2:8)<=repmat(GDJ, [1, 1, 7])*1.1];
Constraints = [Constraints, mDYJ(:,:,2:8)<=repmat(GDYJ, [1, 1, 7])*1.1];
Constraints = [Constraints, mDEJ(:,:,2:8)<=repmat(GDEJ, [1, 1, 7])*1.1];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)))) ;
% OBJ = sum(sum(sum(a)))+sum(sum(sum(b)))+sum(sum(sum(bb)))+sum(sum(sum(c)))+sum(sum(sum(cc)))+sum(sum(sum(d)))+sum(sum(sum(dd)));% 单目标求解
%ops = sdpsettings('solver', 'gurobi', 'gurobi.IterationLimit', 100000, 'guro-bi.TimeLimit', 3600);
%ops = sdpsettings('solver', 'gurobi', 'gurobi.FeasibilityTol', 1e-9, 'guro-bi.OptimalityTol', 1e-9);
ops = sdpsettings(...'solver', 'gurobi', ... % 使用 Gurobi 求解器'gurobi.timelimit', 60, ... % 设置时间限制为60秒'gurobi.MIPFocus', 1, ... % 专注于找到可行解(MIPFocus = 1)'gurobi.Heuristics', 0.5, ... % 启用启发式求解以提高找到可行解的概率'gurobi.NodeLimit', 1000, ... % 设置节点限制(可选,根据需要调整)'allowmilp', 1, ... % 允许混合整数规划'debug', 1, ... % 启用调试信息'verbose', 0 ... % 关闭详细输出(根据需要调整)
);
obj=-OBJ2;
result=optimize(Constraints,obj);
if result.problem==0disp('!!!!!!!!!!!!Solver thinks it is feasible')
elseif result.problem == 1disp('!!!!!!!!!!!!Solver thinks it is infeasible')
elsedisp('!!!!!!!!!!!!Timeout, Display the current optimal solution')
end
value(obj) a = value(a);
b = value(b);
c = value(c);
d = value(d);
bb = value(bb);
cc = value(cc);
dd = value(dd);%% 数据写入% 创建矩阵数据2024
a1 = a(:,:,2);
b1 = b(:,:,2);
bb1 = bb(:,:,2);
c1 = c(:,:,2);
cc1 = cc(:,:,2);
d1 = d(:,:,2);
dd1 = dd(:,:,2); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2024';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 创建矩阵数据2025
a1 = a(:,:,3);
b1 = b(:,:,3);
bb1 = bb(:,:,3);
c1 = c(:,:,3);
cc1 = cc(:,:,3);
d1 = d(:,:,3);
dd1 = dd(:,:,3); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2025';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 创建矩阵数据2026
a1 = a(:,:,4);
b1 = b(:,:,4);
bb1 = bb(:,:,4);
c1 = c(:,:,4);
cc1 = cc(:,:,4);
d1 = d(:,:,4);
dd1 = dd(:,:,4); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2026';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 创建矩阵数据2027
a1 = a(:,:,5);
b1 = b(:,:,5);
bb1 = bb(:,:,5);
c1 = c(:,:,5);
cc1 = cc(:,:,5);
d1 = d(:,:,5);
dd1 = dd(:,:,5); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2027';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 创建矩阵数据2028
a1 = a(:,:,6);
b1 = b(:,:,6);
bb1 = bb(:,:,6);
c1 = c(:,:,6);
cc1 = cc(:,:,6);
d1 = d(:,:,6);
dd1 = dd(:,:,6); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2028';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 创建矩阵数据2029
a1 = a(:,:,7);
b1 = b(:,:,7);
bb1 = bb(:,:,7);
c1 = c(:,:,7);
cc1 = cc(:,:,7);
d1 = d(:,:,7);
dd1 = dd(:,:,7); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2029';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 创建矩阵数据2030
a1 = a(:,:,8);
b1 = b(:,:,8);
bb1 = bb(:,:,8);
c1 = c(:,:,8);
cc1 = cc(:,:,8);
d1 = d(:,:,8);
dd1 = dd(:,:,8); % Excel 文件和工作表名称
filename = 'result1_2.xlsx';
sheetname = '2030';% 写入 Excel 文件
% 1. 首先将 C2 到 AQ83 全部写为 0
xlswrite(filename, zeros(82, 43), sheetname, 'C2:AQ83');% 2. 将 a1 写在 C2 到 Q27
xlswrite(filename, a1, sheetname, 'C2');% 3. 将 b1 写在 R28 到 AJ35
xlswrite(filename, b1, sheetname, 'R28');% 4. 将 bb1 写在 AK56 到 AM63
xlswrite(filename, bb1, sheetname, 'AK56');% 5. 将 c1 写在 R36 到 AJ51
xlswrite(filename, c1, sheetname, 'R36');% 6. 将 cc1 写在 AN36 到 AQ51
xlswrite(filename, cc1, sheetname, 'AN36');% 7. 将 d1 写在 R52 到 AJ55
xlswrite(filename, d1, sheetname, 'R52');% 8. 将 dd1 写在 R80 到 AJ83
xlswrite(filename, dd1, sheetname, 'R80');% 收获(斤)
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);% 最终收益% 初始化 benefits 为 0
benefits = 0;mDJ1 = mDJ(:,:,2:8);
SDJ1 = repmat(SDJ, [1, 1, 7]);
GDJ1 = repmat(GDJ, [1, 1, 7]);% 获取 mDJ1, GDJ1 和 SDJ1 的维度
[dim1, dim2, dim3] = size(mDJ1);% 遍历所有元素并更新 benefits
for i = 1:dim1for j = 1:dim2for k = 1:dim3if mDJ1(i,j,k) <= GDJ1(i,j,k)% 如果 mDJ 小于或等于 GDJbenefits = benefits + mDJ1(i,j,k) * SDJ1(i,j,k);else% 如果 mDJ 大于 GDJbenefits = benefits + GDJ1(i,j,k) * SDJ1(i,j,k) ...+ (mDJ1(i,j,k) - GDJ1(i,j,k)) * SDJ1(i,j,k) / 2;endendend
endmDYJ1 = mDYJ(:,:,2:8);
SDYJ1 = repmat(SDYJ, [1, 1, 7]);
GDYJ1 = repmat(GDYJ, [1, 1, 7]);% 获取 mDJ1, GDJ1 和 SDJ1 的维度
[dim1, dim2, dim3] = size(mDYJ1);% 遍历所有元素并更新 benefits
for i = 1:dim1for j = 1:dim2for k = 1:dim3if mDYJ1(i,j,k) <= GDYJ1(i,j,k)% 如果 mDJ 小于或等于 GDJbenefits = benefits + mDYJ1(i,j,k) * SDYJ1(i,j,k);else% 如果 mDJ 大于 GDJbenefits = benefits + GDYJ1(i,j,k) * SDYJ1(i,j,k) ...+ (mDYJ1(i,j,k) - GDYJ1(i,j,k)) * SDYJ1(i,j,k) / 2;endendend
endmDEJ1 = mDEJ(:,:,2:8);
SDEJ1 = repmat(SDEJ, [1, 1, 7]);
GDEJ1 = repmat(GDEJ, [1, 1, 7]);% 获取 mDJ1, GDJ1 和 SDJ1 的维度
[dim1, dim2, dim3] = size(mDEJ1);% 遍历所有元素并更新 benefits
for i = 1:dim1for j = 1:dim2for k = 1:dim3if mDEJ1(i,j,k) <= GDEJ1(i,j,k)% 如果 mDJ 小于或等于 GDJbenefits = benefits + mDEJ1(i,j,k) * SDEJ1(i,j,k);else% 如果 mDJ 大于 GDJbenefits = benefits + GDEJ1(i,j,k) * SDEJ1(i,j,k) ...+ (mDEJ1(i,j,k) - GDEJ1(i,j,k)) * SDEJ1(i,j,k) / 2;endendend
endbenefits