2024高教社杯数学建模竞赛A题B题C题D题E题完整成品文章和全部问题的解题代码完整版本更新如下:https://www.yuque.com/u42168770/qv6z0d/rytbc1nelty1mu4o
要知道国赛的国奖不是按照比例来给的,每年国赛国一ABC题每题固定100个队左右,所以选的人越少越容易获奖!
推荐大家选A题冲奖,获奖机会大,而且推荐国奖不容易被刷,去年我们赛区A题推荐国奖没有一个被刷,B题刷了将近25%国奖,C题刷了15%以上国奖。本人打了20余场比赛,国赛和美赛做的都是A题,拿了21国赛国一,22美赛的F,今年我又带队参加24年的数学建模国赛了
很多人觉得A题难,其实A题只是看起来难,其实也有很多套路。下面我就来分享一下我的A题建模经验。
2024数学建模国赛常见模型汇总
2024年高教社杯数学建模国赛各种评价方法的适用条件
国赛做题套路
不管是A题B题C题,优化问题是国赛和美赛中最重要也是最常见的问题,基本上大部分模型都可以建成优化模型,所以,优化问题需要重点关注
不管是国赛还是美赛,优化模型都可以分成外部的优化模型和内部的具体模型
外部模型目标函数:(,,,)控制方程(内部模型):(,,)约束条件:;;;初始条件:;;;外部模型{目标函数:minS=G(a,b,x,y)控制方程(内部模型):y=F(a,b,x)约束条件:a;b;x;y初始条件:a;b;x;y
第一问建模方法
分析过去5年的国赛A题,第一问几乎都是参数优化问题,一般会给附件数据,我们把它抽象出来,也就是需要先建立一个模型,通过附件数据拟合模型中的未定参数。
首先,什么是参数优化?举个例子,对于模型y=f(x)=ax+b,附件中会给定x和y的值。当然,模型y=ax+b不会告诉你,这是自己建立的模型。建立了模型之后,通过数据去拟合,求解参数a和b
怎么求a和b,这就是关于算法的问题。对于y=ax+b这样简单的模型,拟合本质上是优化,我们可以使用各种成熟的优化算法,如:单步长(变步长)搜索法、二分搜索法、黄金分割搜索法、最小二乘法、遗传粒子群模拟退火算法
以下是2021国赛训练做的2020国赛A题炉温曲线的第一问的参数拟合求解
第二、三问建模方法
第二、三问是正问题和反问题的优化求解。通过建立的模型和标定好的参数,将模型中的某些变量作为自变量,再将这些量作为若干个优化目标函数的自变量,求解目标函数对应的解。
这里的目标函数可能是以下几种情况,函数关系式记为F:
-
模型中的参数或者参数的函数关系式,比如a,b,或者F(a;b)
-
模型的自变量或者自变量的函数关系式,比如x,或者F(x)
-
模型的因变量或者因变量的函数关系式,比如y,或者F(y)
-
关于模型中的若干个变量的函数关系式,比如F(a;b;x;y)
注意上面的a,b,x,y可以推广为向量或者矩阵形式,a=[a1,a2,a3...an];b=[b1,b2,b3...bn];x=[x11,x12,x13,...x1n;x21,x22,x23,...x2n;...;xn1,xn2,xn3...,xnn]
同时,参数a,b也可以分为恒定参数和可变参数。例如2020A题中,各温区的温度([T1-5, T6, T7, T8-9])和过炉速度(v)是可变参数,热传导因子( )是恒定参数,()已经由第一问拟合得到。
上面第1、2种情况是反问题,第3种情况是正问题,第4种情况是正反问题的结合
以y=ax+b为例,第一问已经建立好了y=ax+b这个模型,并且根据附件数据求解了最优的参数a和b的值。现在,整个模型已经确定。第二问的目标函数的自变量为y;x;a;b,目标函数的因变量(记为y2),y2=f(y,x,a,b),求解使得y2最优的(y、x、a,b)的其中一个或者某几个。
以2020A题为例,问题二给出了各温区温度的设定值,其中小温区1-5 为182ºC,小温区6 为203ºC,小温区7 为237ºC、小温区89为254ºC,需要求解在满足制程界限的前提下,允许的最大传送带过炉速度。对应到上面的y=ax+b的模型,目标函数量是传送带过炉速度v,其中v是模型中的参数,即符合第一种情况
问题3需要求解[T1-5, T6, T7, T8-9, v],使得超过217ºC 到峰值温度所覆盖的面积最小。这里的优化自变量是[T1-5, T6, T7, T8-9, v],这些都是模型中的参数,目标函数面积S是温度曲线的函数,温度曲线可以由参数[T1-5, T6, T7, T8-9, v]和自变量(时间t)等通过模型推导得到。因此这一问的目标函数是第四种情况,即F(a;b;x;y)
添加图片注释,不超过 140 字(可选)
问题4的是多目标函数,面积S最小化和对称性指标D最小化。同第三问,目标函数是第四种情况
如何判别是否是优化问题
-
凡事关于最小最大相关或者尽量相关的词汇,一般都是优化问题。收益最大,风险最小,尽量对称,尽量小,尽量大等等
-
另外还有一种最后趋于稳定的问题,一般也是优化问题,例如2019A高压油管压力控制,油压最后稳定在100,优化目标是最后一秒的油压曲线与100的偏差最小
优化问题的组成部分
注意优化问题特别是微分方程问题,一般由5个主要部分构成,它们分别是:
-
优化变量(自变量)
-
优化目标(单目标、多目标)(因变量)
-
控制方程(微分方程模型主体及其他附属方程)
-
边界条件
-
初始条件
在论文写作时,我们需要做的是把优化目标、控制方程、约束条件等用数学语言表达出来。需要注意:应该明确写出上面5个部分,如果某些部分没有,应该至少写出优化目标和控制方程
以我们2021国赛暑期训练为例,我们做的是2020国赛A题炉温曲线,各位可以自行阅读题目,此处不赘述题目
模型的建立部分
在写作时,明确写出以xx为优化变量、以xx为控制方程(约束条件)、以xx为优化目标的单目标(多目标)优化模型
上面写完了模型的引入,下面开始写模型的公式:
其中优化变量和优化目标可以写一起,举例如下。其中排第一位的是优化变量的估计,表示求解出来的值只是全局最优的估计;argmin(argmax)下标优化变量,表示使得后面的优化目标方程最小时的优化变量的精确值;最后是优化目标方程
上面的是单目标的写法,对于多目标优化,写法也一样,如下:
下面开始写控制方程以及约束,注意有下面两点。
第一点、使用大括号把控制方程和约束条件、边界条件写一起。
第二点、后面需要有变量符号的解释(其中,a为:,b为:...)
为什么要这么写?原因是这样评委一看就知道所有的东西,一目了然!
当然,你也可以将上面5个部分分开写,我2022美赛A题就是这样写的,如下,也就是每个部分写个小标题并加粗:
第一、决策变量(优化变量))
第二、目标方程(优化目标)
第三、控制方程
第四、约束条件
第五、初始条件
不管是国赛还是美赛,优化模型都可以分成外部的优化模型和内部的具体模型
下面将外部模型的求解和内部模型的求解分开讲述,首先是外部优化模型的求解:
外部优化模型的求解
上面建立完了模型,我们就开始写模型的求解,这部分就是算法了
对于一个优化问题,如何选择算法?当然算法不是乱选的,你得体现你的思考与创新
比如对于一个参数搜索问题,变步长搜索就比单步长搜索要好,为什么?因为变步长搜索在趋于最优解的时候缩小步长,能够减少运算次数,这就是它的优点。
再进一步,二分搜索法比变步长搜索法更好。你用二分,别人用穷举搜索,你的算法就比别人的更优秀。
注意,不是越复杂的算法就越好,而是使用简单的算法在最短的时间内解决问题才是最优的算法。
下面是一些常见优化问题的求解算法
参数拟合或者参数标定问题
利用最小二乘法,优化目标是误差最小,拟合函数一般是复杂的常微分方程(2019A)或者偏微分方程(2020A,2018A),例如2020A题需要拟合带10个未知参数的偏微分方程,这时候就需要用遗传算法等智能算法了
求复杂函数的最小\最大值\零点
-
求解最大最小值一般,先用变范围搜索的方法大致搜索一遍,然后:
-
对于单调函数,用二分法
-
对太复杂的,一般用智能算法
内部具体模型的求解
-
求解常微分方程(组)的方法有欧拉法,改进欧拉法,龙格库塔法,算法最好自己写,尽量不调库。例子:2019A高温油管压力控制
-
求解偏微分方程(组)一般是有限差分法,一般需要自己写算法。例子:2018A高温服装设计、2020A炉温曲线
-
多变量、多约束优化问题求解,使用罚函数法,例子:2021A题FAST反射面节点调节
常微分方程的求解
为了编程时更好的与我们外部优化模型的迭代算法匹配,常微分方程我们一般使用手写的改进欧拉法代码迭代求解。
利用欧拉法求解微分方程 ,将时间等分,区间长度是h, ,则 ,精度为O(h)。
利用改进欧拉法求解微分方程 ,将时间等分,区间长度是h, ,则 ,精度为O(h2)。
利用改进欧拉法求解混合微分方程 ,将时间等分,区间长度是h, 。则y1的迭代公式为 ,y2的迭代公式为
微分方程组求解要点:
-
因变量个数与方程个数相同
-
尽量使用最基础的理论(牛顿第二定理等)
-
只采用一阶微分方程组的形式
-
尽量采用向量方式处理
-
尽量做到可进行数组运算(并行计算, 存贮空间换时间)
-
如无法理论求解,不做化简(只做计算机数值计算)
-
保留混合微分方程的形式
1.建立混合微分方程的函数代码,(对复杂的微分方程组往往比较困难)
2.利用库函数求解(精度满足条件),或用欧拉、改进欧拉自编程求解(时间换精度)
在无法理论求解的情况下,高阶理论往往只有理论上的优势,对数值求解帮助不大,基础理论形式更简洁明了
以2019年A题高压油管控制为例,第一问建立的模型如下:
取定参数a, T,将时间等分,区间长度是h,初始化: ,根据改进欧拉法,由 计算 。依次循环计算下去,可得所有时间节点的压力及密度值
编程的代码如下:
function [p,t]=f(time,T,p1) %time 总时长, T 开阀时间
global p0 A C V
h=0.001; n=floor(time*1000/h)+1; a=0;
t=(0:n-1)*h; %时间列
p=zeros(n,1); %对应的压强
midu=zeros(n,1); %对应的密度
rho0=dens(p0); p(1)=p1; midu(1)=dens(p1);
tmp1=h*C*A*eta(t,a,T)*sqrt(2*rho0)/V;
tmp2=h*q(t)/V;
%改进欧拉法迭代求解
for i=1:n-1
tmp=tmp1(i)*sqrt(p0-p(i))-midu(i)*tmp2(i);
midu(i+1)=midu(i)+tmp;
p(i+1)=p(i)+tmp*sc(p(i))/midu(i);
end
end
第二问建立的模型如下,求解方法同理:
偏微分方程的求解
偏微分方程的求解使用有限差分法,有前向差分、后项差分、CN隐格式差分。使用差分法的时候,需要考虑稳定性和收敛性的问题。Crank-Nicolson 隐格式收敛性与稳定性都优于古典显格式和古典隐格式。
同时,需要明确是第几类边界条件!!选错边界条件类型,直接全盘出错!!!
以热方程的求解为例(2018A题高温服装设计):
添加图片注释,不超过 140 字(可选)
添加图片注释,不超过 140 字(可选)
建立好模型之后,编程求解也不是很困难,核心就是3个矩阵:
多变量、多约束优化问题求解 (罚函数法)
对于多变量、多约束优化问题: ,转化为新无约束优化问题, 。
求解步骤如下:
-
令x1=x0,计算 ,如果 ,跳到3,否则2
-
找 ,设最小点为 ,令k=k+1,
-
计算等式、不等式平方和及|x0-x1|, 若最大值<esp, 结束,否则 1
其中esp为误差精度; 为初始点,初始取k=1
要点:
-
多变量、无约束优化问题有很多算法,选合适的
-
变量太多,粒子群、遗传、蚁群等算法应该不适应
-
要注意用向量并行计算,避免循环语句
以2021A题FAST节点调节为例,建立以下模型:
部分求解代码如下
while tol>eps
x1=x0;
while 1
y0=dfun(x0)+c*(dineq(x0)+deqf(x0));
tmp=norm(y0);
if tmp<eps
break;
end
h=-y0/tmp;
fhandle=@(t) fun(x0+t*h)+c*(ineq(x0+t*h)+eqf(x0+t*h));
range=max(tmp,range);
[t,g0]=funm(fhandle,range);
range=t*50;
g1=fun(x0);
g2=ineq(x0);
g3=eqf(x0);
x0=x0+t*h;
end
if c>1e8
break;
else
c=r*c;
end
tol=max([err1,err2,err3]);
end
function y=dineq(x) %x(n,4) 主索节点坐标及促动器伸缩量
%y 不等式约束在x的梯度值
global side_ind_square index
r=0.0007; % 边长最大容许伸缩比
E=x(:,1:3); % 主索节点坐标
s=x(:,4); %促动器伸长量
y=zeros(size(x));
tmp=E(side_ind_square(:,1),:)-E(side_ind_square(:,2),:);
tmp1=sum(tmp.*tmp,2);
tmp2=tmp1-side_ind_square(:,3)*(1+r)^2;
tmp3=side_ind_square(:,3)*(1-r)^2-tmp1;
tmp2=max(0,tmp2);
tmp3=max(0,tmp3);
tmp4=4*(tmp2-tmp3).*tmp;
tmp5=max(0,s.^2-0.36);
y(:,1:3)=index*tmp4;
y(:,4)=4*s.*tmp5;
end
function g=ineq(x) %x(n,4) 主索节点坐标及促动器伸缩量
%g 不等式约束在x值
global side_ind_square err1 err3
r=0.0007; % 边长最大容许伸缩比
E=x(:,1:3); % 主索节点坐标
s=x(:,4); %促动器伸长量
tmp=E(side_ind_square(:,1),:)-E(side_ind_square(:,2),:);
tmp1=sum(tmp.*tmp,2);
tmp2=tmp1-side_ind_square(:,3)*(1+r)^2;
tmp3=side_ind_square(:,3)*(1-r)^2-tmp1;
tmp2=max(0,tmp2);
tmp3=max(0,tmp3);
tmp4=max(0,s.^2-0.36);
g=sum(tmp2.^2)+sum(tmp3.^2)+sum(tmp4.^2);
err3=max(tmp4/0.72);
err1=max(max(tmp2,tmp3)./side_ind_square(:,3));
end
function g=eqf(x) %x(n,4) 主索节点坐标及促动器伸缩量
%g 等式约束在x值
global L l P e err2
E=x(:,1:3); % 主索节点坐标
s=x(:,4); %促动器伸长量
tmp=P-E+e(:,1:3).*(L+s);
tmp1=sum(tmp.*tmp,2)-l;
t1=tmp1.^2;
g=sum(t1);
err2=max(abs(tmp1)./l)/2;
end
function y=deqf(x) %x(n,4) 主索节点坐标及促动器伸缩量
%y 等式约束在x的梯度值
global L l P e
E=x(:,1:3); % 主索节点坐标
s=x(:,4); %促动器伸长量
y=zeros(size(x));
tmp=P-E+e(:,1:3).*(L+s);
tmp1=sum(tmp.*tmp,2)-l;
tmp2=sum((P-E).*e,2);
y(:,1:3)=4*tmp1.*(E-P-(L+s).*e);
y(:,4)=4*tmp1.*(L+s+tmp2);
end
function g=fun(x) %x(n,4) 主索节点坐标及促动器伸缩量
%g 优化函数在x值
global k R p
f=0.466*R;
s=-0.3361;
F=f-s; %焦距
E=x(k,1:3);
tmp=sum(E.^2,2);
tmp1=E(:,1:3)*p';
tmp2=tmp1.^2;
tmp3=tmp-tmp2-4*F*(tmp1-s+R);
g=tmp3'*tmp3;
end
function y=dfun(x) %x(n,4) 主索节点坐标及促动器伸缩量
%y 优化函数在x的梯度值
global k R p
f=0.466*R;
s=-0.3361;
F=f-s; %焦距
y=zeros(size(x));
E=x(k,1:3);
tmp=sum(E.^2,2);
tmp1=E(:,1:3)*p';
tmp2=tmp1.^2;
tmp3=tmp-tmp2-4*F*(tmp1-s+R);
y(k,1:3)=4*tmp3.*(E-tmp1*p-2*F*ones(length(k),1)*p);
end