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。
以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; % 处理缺失数据
else
parts = 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;
else
parts = 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;
else
parts = split(str, '-');
nums = str2double(parts);
SDEJ(i) = mean(nums);
end
end
SDJ = 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: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:6
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:6
colIndex = B_data(oddRows(i) - 1) - 34;
if colIndex >= 1 && colIndex <= 3
bb0(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 <= 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:length(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:length(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:length(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(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:6
Constraints = [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==0
disp('!!!!!!!!!!!!Solver thinks it is feasible')
elseif result.problem == 1
disp('!!!!!!!!!!!!Solver thinks it is infeasible')
else
disp('!!!!!!!!!!!!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:dim1
for j = 1:dim2
for k = 1:dim3
if mDJ1(i,j,k) <= GDJ1(i,j,k)
% 如果 mDJ 小于或等于 GDJ
benefits = benefits + mDJ1(i,j,k) * SDJ1(i,j,k);
else
% 如果 mDJ 大于 GDJ
benefits = benefits + GDJ1(i,j,k) * SDJ1(i,j,k) ...
+ (mDJ1(i,j,k) - GDJ1(i,j,k)) * SDJ1(i,j,k) / 2;
end
end
end
end
mDYJ1 = 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:dim1
for j = 1:dim2
for k = 1:dim3
if mDYJ1(i,j,k) <= GDYJ1(i,j,k)
% 如果 mDJ 小于或等于 GDJ
benefits = benefits + mDYJ1(i,j,k) * SDYJ1(i,j,k);
else
% 如果 mDJ 大于 GDJ
benefits = benefits + GDYJ1(i,j,k) * SDYJ1(i,j,k) ...
+ (mDYJ1(i,j,k) - GDYJ1(i,j,k)) * SDYJ1(i,j,k) / 2;
end
end
end
end
mDEJ1 = 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: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)) * SDEJ1(i,j,k) / 2;
end
end
end
end
benefits
4366
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
