2021数学建模国赛思路 已修订
Posted 您好啊数模君
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了2021数学建模国赛思路 已修订相关的知识，希望对你有一定的参考价值。
目录
已出：
国赛A题第一版思路 FAST主动反射面的形状调节（已修订）
国赛B题第一版思路 乙醇偶合制备C4烯烃（已修订）
国赛C题第一版思路 生产企业原材料的订购与运输（已修订）
国赛D题第一版思路 连铸切割的在线优化
国赛E题第一版思路 中药材的鉴别
国赛C题第一问程序
国赛C题第二问程序
国赛C题第三问程序
思路文件下载链接：
链接：https://pan.baidu.com/s/13aSy2-hlkLa7Ps8wknYY1Q
提取码：gap4
这里有百种算法出处整理，国赛题算法可从上面找取：
A “FAST”主动反射面的形状调节
附件1画的图，为每块反射面链接主锁节点的位置
X=xlsread('附件1.csv');
figure
plot3(X(:,1),X(:,2),X(:,3),'*')
附件2画的图，每个主锁有上下两端点，为什么要给呢，一是主锁倾斜方向就是促动器的作用力方向，就算是之后的工作状态主锁的倾斜方向也不会有太大变化，可忽略；二是结合附件1可以得到每块反射面基准态时的倾斜方向角度，因为每块反射面的角点与其对应的主锁上端点的距离是固定的，需要先计算出来
X=xlsread('附件2.csv');
figure
plot3(X(:,1),X(:,2),X(:,3),'b*')
hold on
plot3(X(:,4),X(:,5),X(:,6),'r*')有些同学会发现，附件1和附件2的上端点对不上，确实是的，但是每个上端点与反射面角点的间距是固定的，所以后面改变下拉索伸缩时，反射面角点的位置有了与主锁上端点距离，以及主锁方向上的拉伸可以重新计算出反射面角点的位置
附件3结合附件1画的图
[~,~,X1]=xlsread('附件1.csv');
[~,~,X2]=xlsread('附件3.csv');
X1=string(X1);X2=string(X2);
X1(1,:)=[];X2(1,:)=[];
figure
hold on
a=[];b=[];c=[];
for i=1:size(X2,1)
a=find(X1(:,1)==X2(i,1));
b=find(X1(:,1)==X2(i,2));
c=find(X1(:,1)==X2(i,3));
plot3(double(X1([a,b,c],2)),double(X1([a,b,c],3)),double(X1([a,b,c],4)),'b-*')
end
中国天眼图，这道题可以不用考虑每个主锁节点和反射面的重量，题目也没给，实际中就算是有重量，基准态时促进器也会给予一定作用力来保持球面形貌，在这道题中我们就假设基准球面是一个球形面，并且工作态时的照射面会形成一个椭球形面
题目中的α和β为方位角和仰角，注意α是x与y轴的夹角，β是CS线与xy平面的夹角
来看下附录：第7点需要理解下
基准状态下，促动器顶端径向伸缩量为0，其径向伸缩范围为-0.6～+0.6 米，这个条件就是刚刚说到的基准态时促进器也会给予一定作用力来保持圆球形貌，所以会有一个伸缩控制范围，但是需要注意调整工作抛物面需要促动器联动，那么如何分析呢，其实SC直线与基础球形面的垂线的交点就是照射面的中心，那么我们可以认为该中心的下拉索收缩幅度是最大的，而照射面边界点的下拉索会有一定伸长，其伸长是所有主索最大的，其实我们结合实际想想，最后会形成一个照射面截面为圆的椭球体面，那么空间的椭球体公式为(x-x1)²/a²+(y-y1)²/b²+(z-z1)²/c²=1,其中(x1,y1,z1)为椭球体的中心点，我们可以取一些参考点来推算该公式，为了多将光线反射至馈源舱，肯定照射面中心会有合理的凹陷，那么对于边界点的下拉索也有合理的伸长
下面理想抛物线
上面还不理解就看下图，其实理想抛物面就是椭球体一部分，照射区应当为椭球体形貌，那么本文为了好做分析，可以直接用椭球体方程，只求照射面没什么影响，照射面就是下图红色实线部分。接下来我们找一下参考点去推算公式，第一个照射面中心点就是SC垂线与基础球面的的交点，然后计算工作时下拉索收缩后该点的位置，然而具体收缩多少是我们后面要寻优的；然后找出基础球形面上到达SC直线近似150m距离的点，作为照射区的边界点，对于该类参考点，计算工作时下拉索伸长后该点的位置，同样的对于伸长多少，这里也是需要寻优的，不管怎么样即使第二类参考点也要收缩，那也要比第一个参考点收缩更小一点；有了这两类点，我们还缺椭球体公式中的中心点(x1,y1,z1)，但是大家想想我们知道了边界点和中心点，但是未知其他位置的工作时的形貌，只知道是一个理想的抛物面，也就是的椭球体面，然而椭球体的中心点位置是会影响到其余照射区的形貌的，但是我们能知道椭球体的中心点在SC直线上，这里我们可以对椭球体中心点到照射中心点的距离进行寻优，从而推出椭球体公式；以上三类点足够推出该椭球体公式了，至于照射区域，参考题目给的图4，就考虑最大照射区就行，然后寻优得到他的形貌公式，有了该椭球体公式，椭球主轴又与SC直线重合，接下来分别求基础圆面和椭球体与SC的交点，然后以基础圆面上各节点与SC直线的交点做弧长距离，按该距离取对应椭球体面上找到相应的位置，通过计算照射区内前后位置变化就可以求得各节点对应的主索伸缩情况
这道题本是一个优化问题，可以观察这幅图可以看一下工作区的形貌变化
第一问自变量有照射中心的主索伸缩量、照射区边界点主索伸缩量(照射区边界点主索伸缩量应当相同)、椭球体中心到照射中心的距离(先有照射中心主索伸缩量在生成该参数)，以上参数会得到一个椭球体面，基础圆面和椭球体面与SC直线会得到两个交点，根据弧长距离对应取得照射区主索的更新位置，然后计算主索的伸缩量，要找到最优的抛物面，肯定是要求更多的反射片能够反射至馈源舱内，馈源舱为直径为1m的圆，反射片为三角形，如下图，只要反射的光线能涵盖馈源舱部分都可以算是反射成功，目标函数可以定位能够将目标反射至馈源舱的反射片数，第三问在考虑反射的面积，这里还要计算出各三角形反射片的角点，以确定其倾斜方向，然后去算反射线能达到馈源舱，对于光线反射问题，反射片的垂线作为主轴去算就行，可能程序有点繁琐
第一问就考虑正上方的天体，第二问考虑其他位置的，大家也不用去纠结前面讲的照射区，按α=36.795°，β=78.169°方向，注意方向是以基础圆面中心垂直线来参考的角度的，按SC线算的的照射区口径差不多就是300m了，第二问没有给天体高度距离来算，或者如果想讨论不同天体高度距离也可以，天体所处位置方向不变就行。第二问同样按照第一问模型来做寻优
第三问则是要考虑上每块反射片，反射光线后能照射到馈源舱圆形上的面积之和，目标函数改为这个其他步骤同理，最后就基础球面对比
B 乙醇偶合制备C4烯烃
首先来看下这道题的逻辑，第一问找关系式，第二问比较催化剂组合，第三问寻优，第四问增设实验
可能会有小伙伴去考虑除乙醇转化率、C4烯烃选择性以外的指标，
其实这个表的意思是，在这种催化剂组合下，以及恒定温度及添加的乙醇速率条件下，达到反应平衡时，有2.07的乙醇转化为了右边几类化学物质，题目要求分析C4，那就只分析C4就行，如果看觉得内容少还可以自行增加一些分析
第一问首先分别找到乙醇转化率、C4烯烃的选择性与温度的关系式，求关系式算法可以参考下面几篇推文，也可以用1stopt软件进行拟合，21种组合都要求下温度与乙醇转化率、C4烯烃的选择性的关系式
最小二乘与非线性最小二乘(Matlab)
https://mp.weixin.qq.com/s/Jl4Zkx-TK6fCNN4tgd47jQ
逐步回归补充篇——确定经验公式https://mp.weixin.qq.com/s/zmN7QrbusZLNRjw-8gJcBg
接下来分析附件2结果，附件2中273分钟时，乙醇转化率已稳定，C4烯烃选择性变化都很小，这个时候我们可以认为273分钟已经达到了平衡，接下来就用350摄氏度、乙醇转化率29.9%、C4烯烃选择性39.04%与附件1中的实验进行匹配，选出差别最小的作为该给定的催化剂组合，这里匹配简单做法直接算距离，也可以水论文用机器学习算法做，对附件2不用做时间与各项指标的关系曲线，后问也用不上，想多做一点分析也可以。
第二问，其实就是要让你们找哪一个对照试验比较好，对照试验一般使用的是方差分析，可以分别求得对照组内某因素是否对乙醇转化率、C4烯烃选择性结果显著影响，如果显著影响那么影响程度怎么样，p值是多少，p值月低越好，此外还要看下这组对照试验中出现最好的结果（乙醇转化率以及C4烯烃选择性）分别是多少，这样可以得到四个指标，然后采用评价算法对不同对照试验组进行评价，评价可以用多种算法例如投影寻踪、熵权法、Tosis等评价，然后在通过组合评价方法例如Borda求一个综合性的评价值，然后进行排序，但是注意这里对数据进行标准化时有一个问题，客观算法可能会存在偏离主观的现象，可以这么来做，对于你们认为的重要性小的指标可以通过ln(β+x)函数映射，这样能够减少该指标对最后结果的影响。
方差分析（单因素、双因素、不均匀样本）
方差分析(Matlab)https://mp.weixin.qq.com/s/rdvmOjClKNUhD9_S-x_rqQ
本问主要用方差分析去找到比较好的对照实验方式，题目说探讨不同催化剂组合及温度对乙醇转化率以及C4烯烃选择性大小的影响，其实就是说怎样的对照试验中温度对结果的影响是最大的，并且实验结果中能产生不错的结果，那么第三问就基于此进行研究
接下来找一找对照试验组，题目说的1wt%Co/SiO2没有表述清楚，/SiO2就说明SiO2为主要基体，Co在这个基体中的质量占总质量的1%，也就是如果有100mg1wt%Co/SiO2，那么就有1mgCo，99mgSiO2，石英砂就是SiO2，我们对附件1第二列信息归纳下
我们来找一下对照试验组
①按Co-SiO2-HAP质量比例为1-100-100且三者质量相等的实验组，但乙醇添加速率不同，共找到6组，去掉重复实验实则为4组
然后我们按不均匀样本取数，以不同温度下乙醇转化率结果为例
然后通过方差分析求该对照试验组中的参数条件是否对结果构成显著影响
X=[19.7 29 40 58.6 76 6.3 8.8 13.2 31.7 56.1 1.4 3.5 6.9 19.9 44.5 1.4 3.4 6.7 19.3 43.6 2.1 3 4.7 13.4 40.8 2.1 3.8 5.8 9.8 15.9 45];
group=[ones(1,5),2.*ones(1,5),3.*ones(1,10),4.*ones(1,11)];
[p,anovatab,stats]=anova1(X,group);%单因素方差分析
fa=finv(0.95,anovatab{2,3},anovatab{3,3});%计算fa
F=anovatab{2,5};%F值
if p<=0.01 && F>fa
disp('非常显著')
elseif p<=0.05 && F>fa
disp('显著')
else
disp('不显著')
end
②。。。
其他实验组自己筛选，同①找三个参数相同一个参数不同的对照组实验
在比较出较好对照试验组后，然后同问题一中的方法分析该对照试验组下温度对结果的影响，如果是Co或SiO2或HAP或乙醇速率对照实验组，那还需要分析哪一个参数在同一温度下对结果（对乙醇转化率、C4烯烃选择性大小）的影响，可以拟合出关系式，第二问最后的分析，需要看你们得到的最优对照试验组结果。
第三问，基于第二问结果以C4烯烃收率为目标函数（C4烯烃收率=乙醇转化率*C4烯烃的选择性），基于第三问的关系是，并设立，实验参数及温度范围对目标函数进行寻优，优化算法系列推文见下链接
https://mp.weixin.qq.com/mp/homepage?__biz=MzU0OTEyNjAyMQ==&hid=25&sn=4dd5f9b697c058388b9efba1372f7e34&scene=18&uin=&key=&devicetype=Windows+10+x64&version=63030532&lang=zh_CN&ascene=7&fontgear=2第四问，这里的实验肯定是第二问定下来的实验，一般做实验，都会等间距的设置参数，是找出最佳实验条件所处大致区间，这里可以采用黄金分割法这类方法去增设对照实验，例如第二问确定以Co、SiO2、HAP比例为1:100:100且不同乙醇添加速率的对照实验，参考第三问结果，我会去看最佳乙醇速率设置多少较好，就会在周围设立不同精度的参数，再次得到结果看是否会出现更好的结果，一般设置这样的参数实验，采用的是黄金分割法或者1/2法去设置
C 生产企业原材料的订购与运输
这道题是优化问题，先看看题目，以下为题目告知的条件和目标函数
条件：
①生产48周，每周产能2.82wm³，折合每周需要A类原材料1.692wm³或B原材料1.8612wm³或C原材料2.0304wm³；
②24周的原材料订购和转运计划，那么相当于平均每周需要送达的材料能够满足5.64wm³，这里是包含了订购及转运时间加起来共计24周以内要完成，如果考虑8家转运商满负荷运转以及平均损耗，那么一天送达的就有4.75wm³,就算是按最低的材料置换产能的比例也能大于5.64wm³，所以运输这块没啥问题，接下来看供应商供应能力，如果也按取平均，那么肯定是不够的，所以可以取最大供应量的周数据作为该供应商最大供应能力的体现
③A类和B类原材料的采购单价分别比C类原材料高20%和10%，那么就假设C单位立方米成本为1，A成本为1.2，B成本为1.1；
④三类原材料运输和储存的单位费用相同，一般运输一批货物只算单次运输的费用，这里的储存费用就不用按照每周多少来考虑，如果要考虑，那么就要涉及到储存的材料使用情况，毕竟边供应边生产；
⑤我们可以就假设当周订货当周送达
目标函数：
一家供应商每周供应的原材料尽量由一家转运商运输，也就是说供应商数要少，供应商数可作为一目标函数
接下来看问题
第一问确定50家最重要的供应商，一般的都是从供应情况去看，这里可以确定一些指标，然后进行评价，最后评分最高Top50作为本问结果，可以考虑以下指标：
①供应商都是只供应一种原料，那么就会有材料得性价比，算一下哪一个原材料经济效益最高，用1/[0.6 0.66 0.72]=[1.67 1.52 1.39]，然后再去除下单价/[1.2 1.1 1]=[1.3889 1.3774 1.3889]，然后各供应商按供应的原材料对应，将这个结果作为第一个指标，指标越大越好
②供应商最大供应能力，最大供应量并换算为产能，指标越大越好
③供应商平均供应能力，平均供应量并换算为产能，指标越大越好
④供应商供应得稳定性，对历史供应数据，做下方差，指标越小越好
⑤历史供应量/历史订单量比，指标越大越好
。。。
除了以上指标大家还可以多增加一些，然后用评价算法评价后排序即可，第一问评价可以用多种算法例如投影寻踪、熵权法、Tosis等评价，然后在通过组合评价方法例如Borda求一个综合性的评价值，然后进行排序，但是注意这里对数据进行标准化时有一个问题，客观算法可能会存在偏离主观的现象，可以这么来做，对于你们认为的重要性小的指标可以通过ln(β+x)函数映射，这样能够减少该指标对最后结果的影响
第二问，A、C类原材料经济效益最高，我们发现C类原材料运输费用和储存费用要比AB较高，因为置换为产能需要更多的供应量，并且供应C类材料的供应商是最少的，在本问，考虑供应商的最大供应能力即可，但是为了方便运算，最好是再增加一个产能计算，最大供应能力就以历史5年中供应量最大的一周作为起最大供应能力，对于转运商的运输损耗，可以直接取平均，也可以预测未来24周的每周的损耗都可以，前者更简单，对于本问建模部分，采用优化算法进行多目标求解，目标函数有：1、供应商数最少，2、总成本=运输成本+储存成本=2*运输成本最低；约束条件有：1、24周从供应到转运结束的原材料转化为产能要达标，2、只有24周转运时间，3、每周至多送达8*6000*损耗，4、每周至多订购8*6000的供货量。为了降低寻优维度，可以增加这样一个优先条件，C原材料置换为产能置换率最低，那么C原材料就优先用损耗较大转运商转运，这样对整体的产能损耗较小，其次为B
第二问设计启发式算法如下：
第三问，尽量多地采购 A 类和尽量少地采购 C 类原材料，以减少转运及仓储的成本，同时希望转运商的转运损耗率尽量少，本问在第二问基础上增加一个目标函数A类原材料-C材料采购量为目标函数，同样按第二问方法来做
第四问，只考虑产能，转运商运能拉满情况下去选择供应商，每周最多订6000m³*8=4.8wm³，考虑所有Top50供应商日均ABC类最大供应量，在第二问基础上，第四问增加一个目标函数为产能的最大化，同样按第二问做法来做，但题目说的很明确了，A类材料才是最合适的材料，但是第二问算法中涉及到k个供应商的选取，那么第四问实则就只是选择最佳供应商而已，对于供应量，按A、B、C类原材料的供应商依次按最大供应量选取，直到供应量达到4.8wm³为止，对于选择转运商还是按第二问默认的规则来，C类材料给损耗率最大的转运商
本题目优与供应商数会发生变化，建议通过模拟退火算法框架进行寻优
D 连铸切割的在线优化
切割机切断一块钢坯的时间为3分钟切割后，返回到工作起点的时间为1分钟，连铸拉坯的速度为1.0米/分钟，也就是说切割机最高效率切割间隔为4m。切割后的钢坯长度必须在 4.8~12.6 米之间，否则无法运走，报废段有两种定义，一个是多切割的部分，如果切割下来的长度未达标则全部报废，一个是结晶器出现异常时，突然产生的报废段，此时的报废段长度为0.8m
来看第一问，注意尾坯是除开了从结晶器中心到切割机工作起点处的钢坯长度，有60m，所以第一问钢坯长度并不是问题中的，还需要加上60，第一问不考虑结晶器异常，看题目对m的精确到最后一位小数，那么在第一问中就考虑9:0.1:10共11种切割方式，但是不能确定会切割多少段，这里可以自己设计一个启发式算法，可用模拟退火算法等算法进行外包装
第二问，在第一问的基础上增加判断条件：
第三问，同样考虑不同的目标值区间，按0.1m等间隔取值作为不同的切割方式
E 中药材的鉴别
这里先给大家普及下红外光谱仪，药品会一次性加入到红外光谱机里，然后内部会水循环，水循环速度是不变的，不同药材随水的方向移动是不同的，毕竟这涉及到颗粒大小，重量，扩散性等因素影响，并且红外辐射会激发分子振动从而产生红外吸收光谱，不同药物分子情况不同，因此局造成了每次采集数据时，检测到的吸光度不同，其实简单理解就是说红外光谱仪能根据分子的在水中的扩散性、化学稳定性、颗粒分子大小等等对药物检测出不同的波信号，可能会有波信号相近的，这是仪器避免不了的误差问题
这道题千万别拿着就开始先降维数据了，图1和图3中的药材大概经历了一个水循环周期，图2大概经历了4个水循环周期，所以说这道题直接用降维算法处理是不严谨的，这道题可以直接用全量数据带入机器学习算法中训练然后识别类别，也可以通过时间序列中EMA、AMA做下曲线平滑会对训练效果有一点提升
第一问为无监督式分类，可以这么来做，以某个k聚类算法作为目标函数，即分好类别后计算下样本与所在聚类中心的差值之和，k和聚类中心数作为自变量，外循环配一个优化算法进行单目标寻优
第二问有监督分类，使用机器学习方法比如说dbn、随机森林、Xgboost、支持向量机都可以用用，带入已分类数据进去训练，然后带为止类别进去识别出是哪种类别；这个问也可以通过聚类算法和相关性/距离等方式去做，例如，先对已知类别数据，针对每类算一个聚类中心，然后用聚类中心的数据向量对未知类别的样本数据做相关性/计算距离，哪一个相关信号/距离近就是哪一个类别；当然也可以用协同过滤推荐算法去做
第三问，两种红外有一定差别，一般的近红外要好点，其实这个问就是先用近红外数据识别，然后用中红外数据去矫正，反过来也可以，看你们使用的算法的效果来定，这个问有两种方式去解决，一个是将近红外和中红外检测的数据进行点乘得到一个新的红外数据，直接带入算法中训练然后识别，第二种就是分别对两个红外数据集中已知类别样本进行训练，得出两种结果，然后将这两种结合和输出类别重新构建一个训练集，带入算法再次训练并输出类别
第四问不用考虑同一产地同一类别作为一个类型，一般都是先检测药物类别，在针对每个类别去识别产地，算法同上，也可以每个问换一下算法都行
国赛C第一问程序
评价方法自己选
指标提取
| clear X1=xlsread('附件1 近5年402家供应商的相关数据.xlsx','企业的订货量（m?）'); X2=xlsread('附件1 近5年402家供应商的相关数据.xlsx','供应商的供货量（m?）'); [~,X3,~]=xlsread('附件1 近5年402家供应商的相关数据.xlsx','企业的订货量（m?）'); X3=string(X3); X3=X3(2:end,2); X3(find(X3=="A"))=1.3889; X3(find(X3=="B"))=1.3774; X3(find(X3=="C"))=1.3889; %从供应角度看供应商的供应情况 Y(:,1)=max(X2')';%供应商最大供应能力，正向指标 Y(:,2)=mean(X2,2);%供应商平均供应能力，正向指标 Y(:,3)=1./std(X2')';%供应商供应的稳定性,对历史供应数据，做方差，负向指标，可以做个倒数 Y(:,4)=double(X3);%供应商材料经济效益，正向指标 Y(:,5)=sum(X2,2)./sum(X1,2);%历史供应量/历史订单量比，正向指标 %从企业角度看对供应商的认可情况 Y(:,6)=max(X1')';%供应商最大订单量，正向指标 Y(:,7)=mean(X1,2);%供应商平均订单量，正向指标 Y(:,8)=1./std(X1')';%供应商订单量的稳定性,对历史供应数据，做方差，负向指标，可以做个倒数 %还有其他指标可自加 Y=log(1+Y);%对指标ln函数映射，以减少数据本身差距，将评分重心转移至权重 Y=mapminmax(Y',0,1)'; |
层次分析评价
| a=[1 5/4 5/4 5/9 5/8 5/6 1 1 4/5 1 1 4/9 4/8 4/6 4/5 4/5 4/5 1 1 4/9 4/8 4/6 4/5 4/5 9/5 9/4 9/4 1 9/8 9/6 9/5 9/5 8/5 2 2 8/9 1 8/6 8/5 8/5 6/5 6/4 6/4 6/9 6/8 1 6/5 6/5 1 5/4 5/4 5/9 5/8 5/6 1 1 1 5/4 5/4 5/9 5/8 5/6 1 1]; %构建判断矩阵 n1=length(a); %a矩阵的长度 [x,y]=eig(a); %求特征值 lamda=max(diag(y)) %最大特征向量 num=find(diag(y)==lamda); w0=x(:,num)/sum(x(:,num)) %准则层权向量 ri=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45]; %平均随机一致性指标 CR1=(lamda-n1)/(n1-1); %一致性指标 CR=CR1/ri(n1); %一致性比例 if CR>=0.10 disp('矩阵没通过一致性检验，请重新调整矩阵'); else disp('矩阵通过一致性检验'); end f=Y*w0; [f,paixu]=sort(f,'descend'); |
熵权法评价
| [n,m]=size(Y); for i=1:n for j=1:m p(i,j)=Y(i,j)/sum(Y(:,j)); end end %% 计算第 j 个指标的熵值 e(j) k=1/log(n); for j=1:m e(j)=-k*sum(p(:,j).*log(p(:,j))); end d=ones(1,m)-e; % 计算信息熵冗余度 w=d./sum(d) % 求权值 w f=Y*w'; [f,paixu]=sort(f,'descend'); |
秩和比
| w=[0.0532 0.3138 0.0865 0.0000 0.1076 0.0318 0.1549 0.2523];%通过其他方法求得的权重 ra=tiedrank(Y); [n,m]=size(ra);%求矩阵维度 RSR=mean(ra,2)/n; W=repmat(w,[n,1]);%用于复制平铺矩阵，相当于讲w矩阵看成一个元素，形成n×1维的矩阵并用W矩阵记录 WRSR=sum(ra.*W,2)/n;%加权秩和比 [sWRSR,ind]=sort(WRSR,'descend');%sort为排序函数 p=[1:n]/n; p(end)=0.99999; Probit=norminv(p,0,1); X=[ones(n,1),Probit']; [ab,abint,r,rint,stats]=regress(sWRSR,X); WRSRfit=ab(1)+ab(2)*Probit; f=WRSRfit'; paixu=ind; |
Topsis
| [n,m]=size(Y); %n为A矩阵的行数，m为A矩阵的列数 c=sqrt(sum(Y.*Y)); %规范化决策矩阵 d=Y./c; w=[0.0532 0.3138 0.0865 0.0000 0.1076 0.0318 0.1549 0.2523];%通过其他方法求得的权重 c=w.*d; cmax=max(c); cmin=min(c); for i=1:n c1=c(i,:)-cmax; s1(i)=norm(c1); c2=c(i,:)-cmin; s2(i)=norm(c2); T(i)=s2(i)/(s1(i)+s2(i)); end %排名 [~,pm]=sort(T,'descend'); disp('评分结果，评分区间[0,1]') disp(T) disp('方案排名') disp(pm) |
因子分析
| X=Y; mapping.mean = mean(X, 1); X = X - repmat(mapping.mean, [size(X, 1) 1]);%去均值 C = cov(X);%协方差矩阵 C(isnan(C)) = 0; C(isinf(C)) = 0; [M, lambda] = eig(C);%求C矩阵特征向量及特征值 [lambda, ind] = sort(diag(lambda), 'descend');%排序 lambda=lambda./sum(lambda); L=lambda; lambda=cumsum(lambda); mapping.lambda = lambda; k=find(lambda>0.95); M = M(:,ind(1:k(1)));%%取前k列 mappedX = X * M;%降维后的X mappedX=mapminmax(mappedX',0,1)'; L=L(1:k(1))./sum(L(1:k(1)));%贡献率归一化 f=mappedX*L; [f,paixu]=sort(f,'descend'); |
投影寻踪-聪明的鲨鱼(自己改名)
| w1=[0.0532 0.3138 0.0865 0.0000 0.1076 0.0318 0.1549 0.2523];%通过其他方法求得的权重 [bestx,f]=youhua(w1,Y); |
| function [bestx,f1]=youhua(w1,p) %% 启发式算法+投影寻踪 %可调参数 sub=w1.*0.5; %自变量下限 up=w1.*1.5; %自变量上限 opt=1;%-1为最小化，1为最大化 % 程序为最小化寻优，如果是最大化寻优更改排序部分程序即可 n=length(sub); %自变量个数 num=100*n^2; %种群大小 if num>10000 num=10000; end det=10+10*n;%迭代次数 k=1.5+0.1*n;%k为最大环绕圈数 R=0.01*n;%当鲨鱼进入猎物该范围，则直接对猎物位置进行逼近 Mc=(up-sub).*0.1;%猎物行动最大范围 x=[]; f=[]; for i=1:num x(i,:)=sub+(up-sub).*rand(1,n); %初始化位置 f(i,1)=Target(p,x); end %以最小化为例 if opt==-1 [bestf,a]=min(f);%记录历史最优值 bestx=x(a,:);%记录历史最优解 elseif opt==1 [bestf,a]=max(f);%记录历史最优值 bestx=x(a,:);%记录历史最优解 end trace=[]; trace(1)=bestf; xx=[];ff=[]; for ii=1:det %猎物躲避,蒙特卡洛模拟周围1000次，并选择最佳的点作为下一逃跑点 d=[]; d=pdist2(bestx,x); d=sort(d); z=exp(-d(2)/mean(Mc));%猎物急躁系数 z=max(z,0.1); dx=[]; yx=[]; bestt=[]; bestc=[]; for i=1:1000 m=[]; for k=1:n m=[m,(-1)^randi([1,2])]; end xd=[]; xd=bestx+Mc.*z.*((det-ii)/det).*rand(1,n).*m; xd=max(sub,xd); xd=min(up,xd); dx=[dx;xd];%(det-ii)/det表示随着追捕，猎物可逃窜的范围越来越小 yx=[yx;Target(p,xd)]; end if opt==-1 [bestt,a]=min(yx);%选择最佳逃跑点 bestc=dx(a,:); if bestt<bestf bestf=bestt; bestx=bestc; end elseif opt==1 [bestt,a]=max(yx);%选择最佳逃跑点 bestc=dx(a,:); if bestt>bestf bestf=bestt; bestx=bestc; end end %鲸鱼追捕 for i=1:num %更新公式 if sqrt(sum((x(i,:)-bestx).^2))<=R xx(i,:)=x(i,:)+rand.*(x(i,:)-bestx); xx(i,:)=max(sub,xx(i,:)); xx(i,:)=min(up,xx(i,:)); ff(i,1)=Target(p,xx(i,:)); else xx(i,:)=x(i,:)+real(shayu(x(i,:)-bestx,k)); xx(i,:)=max(sub,xx(i,:)); xx(i,:)=min(up,xx(i,:)); ff(i,1)=Target(p,xx(i,:)); end end %引入上一代进行排序,并重新分配角色 F=[bestf;f;ff]; X=[bestx;x;xx]; if opt==-1 [F,b]=sort(F,'ascend');%按小到大排序 elseif opt==1 [F,b]=sort(F,'descend');%按大到大排序 end X=X(b,:); f=F(1:num,:); x=X(1:num,:); if opt==-1 [bestf,a]=min(f);%记录历史最优值 elseif opt==1 [bestf,a]=max(f);%记录历史最优值 end bestx=x(a,:);%记录历史最优解 trace=[trace;bestf]; end bestx=bestx./sum(bestx); disp(["投影寻踪计算权重结果为：",string(bestx)]) f1=[p*bestx']; |
| function f=shayu(X,K) %X为鲸鱼坐标减去猎物坐标 %k为最大围绕圈数 n=length(X); Y=[]; y=[]; costheta=[]; theta=[]; k=[]; for i=1:n-1 Y=[X(1:i),0]; y(i)=sqrt(sum(X(1:i+1).^2)); %计算角度（圈数） costheta(i)=(X(1:i+1)*Y')/(sqrt(sum(X(1:i+1).^2))*sqrt(sum(Y.^2))); if isnan(costheta(i))==1 costheta(i)=1; end if X(i+1)>=0 theta(i)=acos(costheta(i))/(2*pi); else theta(i)=2-acos(costheta(i))/(2*pi); end theta(i)=theta(i)*2*pi; %自适应调节k if y(i)>=10 k(i)=K*exp(-2); else k(i)=K*exp(-y(i)/5); end end %位置更新公式，左包围或右包围 f=[]; l=[]; yy=y; ttheta=theta; l=k(1).*rand; 以上是关于2021数学建模国赛思路 已修订的主要内容，如果未能解决你的问题，请参考以下文章 2021数学建模国赛A题思路 “FAST”主动反射面的形状调节 第一版思路 思路开源 已修订 2021数学建模国赛C题思路 生产企业原材料的订购与运输 第一版思路 思路开源 已修订 程序 2021数学建模国赛 || 高教社杯 || 赛题 || 思路 |