数学建模竞赛算法全集附matlab实现

一.蒙特卡洛算法

01

1，定义：蒙特卡洛算法是以概率和统计的理论、方法为基础的一种数值计算方法，将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解，故又称随机抽样法或统计实验法。

02

2.适用范围：可以较好的解决多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题。

03

3.特点：蒙特卡洛算法可以应用在很多场合，但求的是近似解，在模拟样本越大的情况下，越接近于真实值，单样本数增加会带来计算量的大幅上升。对于一些简单问题来说，蒙特卡洛是个笨办法，但对于许多问题来说，它往往是个有效，有时甚至是唯一可行的方法。

matlab代码

%作图的代码
x = 0:0.25:12;
y1 = x.^2;
y2 = 12 - x;
plot(x, y1, x, y2)
xlabel('x');ylabel('y');
%产生图例
legend('y1=x^2', 'y2=12-x');
title('马驰绘制');
%图中x轴和y轴的范围，中括号前面是y轴范围，中括号后面是x轴范围
axis([0 15 0 15]);
text(3, 9, '交点');
%画图时加上网格线
grid on

%蒙特卡洛算法的具体实现
%产生一个1行10000000列的矩阵，矩阵中每个数是从0到12之间随机取
x = unifrnd(0, 12, [1, 10000000]);
y = unifrnd(0, 9, [1, 10000000]);
frequency = sum(y=3);
area = 12*9*frequency/10^7;
disp(area);
上述代码是求一重积分的例子

二.数据拟合

01

1.定义：不要求近似函数通过所有的数据点，而是要求他能较好的反应数据的整体变化趋势。

02

2.常用方法：最小二乘拟合方法

matlab代码

%读取表格
A = xlsread('E:\表格\1.xls', 'Sheet1', 'A1:AN2');
B = A;
[I, J] = size(B);

%数据拟合
%x为矩阵的第一行，y为矩阵的第二行
x = A(1,:);
y = A(2,:);
%polyfit为matlab中的拟合函数，第一个参数是数据的横坐标
%第二个参数是数据的纵坐标，第三个参数是多项式的最高阶数
%返回值p中包含n+1个多项式系数
p = polyfit(x, y, 2);
disp(p);
%下面是作图的代码
x1 = 300:10:600;
%polyval是matlab中的求值函数，求x1对应的函数值y1
y1 = polyval(p,x1);
plot(x,y,'*r',x1,y1,'-b');
%plot(x,'DisplayName','x','YDataSource','x');
%figure(gcf);

使用matlab的图形化拟合包（推荐）

（1）将数据导入工作区并通过cftool命令打开matlab的图形化拟合包
[object Object]
（2）选择x、y变量
[object Object]
（3）选择拟合方式和最高项次数
[object Object]
（4）得到拟合结果
[object Object]
使用图形化拟合工具不仅简单快捷，还可以使用多种拟合方式，寻找到最好的拟合曲线。

三.参数估计

01

1.定义：参数估计是统计推断的一种。根据从总体中抽取的随机样本来估计总体分布中未知参数的过程。

02

2.标准特点：无偏性、有效性、一致性。

03

3.主要分类：点估计和区间估计。

04

4.要处理的问题：求出未知参数的估计量和在一定信度下指出所求的估计量的精度。

四.插值

01

1.定义：在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。

02

2.作用：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值情况，估算出函数在其他点处的近似值。

03

3.区别：从定义上看，插值和拟合有一定的相似度，但插值要求近似函数通过给定的所有离散数据，而拟合并不要求这样，只要近似函数能较好的反映数据变化的趋势即可（近似含义不同），当测量值是准确的，没有误差时，一般用插值；当测量值与真实值有误差时，一般用数据拟合。

matlab代码

%years、service和wage是原始数据
years = 1950:10:1990;
service = 10:10:30;
wage = [ 150.697 199.592 187.625 179.323 195.072; 250.287 203.212 179.092 322.767 226.505;153.706 426.730 249.633 120.281 598.243];
[X, Y] = meshgrid(years, service);
% % 三维曲线
% plot3(X, Y, wage)
% 三维曲面
figure
surf(X, Y, wage)

%interp2是matlab中的二维插值函数，前两个参数是已知未知，后两个是未知位置，w是未知位置的插值结果

w = interp2(service,years,wage,15,1975);

五.线性规划

01

定义：线性规划是研究线性约束条件下线性目标函数的极值问题的数学理论和方法。

02

步骤：
（1）根据影响所要达到目的的因素找到决策变量。
（2）由决策变量和所在达到目的之间的函数关系确定目标函数。
（3）由决策变量所受的限制条件确定决策变量所要满足的约束条件。

03

特点：目标函数是决策变量的线性函数。根据具体问题可以是最大化或最小化，二者统称为最优化。约束条件也是决策变量的线性函数。

4.代码

model:
min = 2*x1 + 3*x2;
x1 + x2 >= 350;
x1 >= 100;
2*x1 + x2 <= 600;
end
Lingo实现0-1规划
model:
max = 2*x1 + 5*x2 + 3*x3 + 4*x4;
-4*x1 + x2 + x3 + x4 >= 0;
-2*x1 + 4*x2 + 2*x3 + 4*x4 >= 1;
x1 + x2 - x3 + x4 >= 1;
！@bin(变量)限定变量只能取0或1
@bin(x1);
@bin(x2);
@bin(x3);
@bin(x4);
end

七.多目标线性规划

01

1.定义：在相同的条件下，要求多个目标函数都得到最好的满足，这便是多目标规划。若目标函数和约束条件都是线性的，则为多目标线性规划。

02

2.求解方法：
（1)化多为少的方法，即把多目标华为比较容易求解的单目标或双目标，如主要目标法、线性加权法、理想点法等。
（2）分层序列法，即把目标按其重要性给出一个序列，每次都在前一个目标最优解集内求下一个目标最优解，知道求出共同的最优解。
（3）层次分析法，这是一种定性与定量相结合的多目标决策与分析方法，对于目标结构复杂且缺乏必要的数据的情况更为实用。

01

3.代码：由于解多目标规划问题比较复杂，等之后抽时间会更新这部分的代码。

八.二次规划

01

1.定义：二次规划是非线性规划中的一类特殊数学规划问题，在很多方面都有应用，如投资组合、约束最小二乘问题的求解、序列二次规划在非线性优化问题中应用等。

02

2.特点：目标函数是二次的，并且约束是线性的问题。在非线性约束最优化问题中非常重要，通常作为其他问题的子步骤存在。

九.最优化方法

1.定义：最优化方法，是指解决最优化问题的方法。所谓最优化问题，指在某些约束条件下，决定某些可选择的变量应该取何值，使所选定的目标函数达到最优的问题。即运用最新科技手段和处理方法，使系统达到总体最优。
2.数学意义：最优化方法是一种求极值的方法，即在一组约束为等式或不等式的条件下，使系统的目标函数达到极值，即最大值或最小值。
3.基本要素：最优化模型一般包括变量、约束条件和目标函数三要素。
4.无约束优化：
（1）在数值优化中，一般采用迭代法求解无约束优化问题。
（2）无约束优化的算法框架

3）关键问题：搜索方向、终止条件、步长

（4）线搜索技术:通过公式求出每次迭代的步长，分为精确线搜索和非精确线搜索。

（5）Armijo准则（非精确线搜索的准则）

5.无约束优化算法

原始的牛顿迭代法的缺点：

1）初值点选择（敏感）：基本牛顿法初始点需要足够靠近极小点，否则有可能导致算法不收敛。

2）Hessen矩阵正定

代码：

function [x, val, k] = niudun1(x0)
%基本牛顿法求函数最小值
%x是最小值点，val是最小值，k是迭代次数
epsilon = 1e-5; %终止误差值
k = 0; %迭代次数
%如果梯度的范数小于终止误差值则停止迭代
while norm(g(x0)) > epsilon
d = -G(x0)\g(x0);
x0 = x0 + d';
k = k+1;
end
x = x0;
val = f(x0);
end

%%原函数
function y = f(x)
y = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end

%%原函数的梯度,要弄成列向量的形式
function y = g(x)
y(1,1) = 400*x(1)*(x(1)^2 - x(2)) + 2*(x(1) - 1);
y(2,1) = -200*(x(1)^2 - x(2));
end

%%原函数的Hessen矩阵，注意矩阵每个元素的值
function y = G(x)
y(1, 1) = 1200*x(1)^2 - 400*x(2) + 2;
y(1, 2) = -400*x(1);
y(2, 1) = -400*x(1);
y(2, 2) = 200;
end

（2）阻尼牛顿法（

优点：阻尼牛顿法在基本牛顿法的基础上增加对步长的线搜索，克服了基本牛顿法对于初始点敏感的问题。

算法：

代码：（对于x0要输入列向量）

function [x,val,k]=damp_Newton(fun,gfun, Hess,x0)
%功能: 用阻尼牛顿法求解无约束问题: min f(x)
%输入: x0是初始点, fun, gfun, Hess 分别是求
% 目标函数,梯度,Hesse 阵的函数
%输出: x, val分别是近似最优点和最优值, k是迭代次数.
maxk=10000; %给出最大迭代次数
rho=0.55;sigma=0.5;
k=0; epsilon=1e-5;
tic
while(k gk=feval(gfun,x0); %计算梯度
Gk=feval(Hess,x0); %计算Hesse阵
dk=-Gk\gk; %解方程组Gk*dk=-gk, 计算搜索方向
if(norm(gk) m=0; mk=0;
while(m<20) % 用Armijo搜索求步长
if(feval(fun,x0+rho^m*dk) mk=m; break;
end
m=m+1;
end
x0=x0+rho^mk*dk;
k=k+1
end
x=x0;
val=feval(fun,x);
toc
%gval=norm(gfun(x));
%%%% 目标函数 %%%%%%%%%
%function f=fun(x)
%f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
%%%% 梯度 %%%%%%%%%%%%%%%%%%%
%function g=gfun(x)
%g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1), -200*(x(1)^2-x(2))]';
%%%% Hesse 阵 %%%%%%%%%%%%%%%%%%%
%function He=Hess(x)
%n=length(x);
%He=zeros(n,n);
%He=[1200*x(1)^2-400*x(2)+2, -400*x(1);
% -400*x(1), 200 ];

3）最速下降法

优点：最速下降法只需要求函数的梯度，并不需要求函数的Hessen矩阵，而且最速梯度法对于初始值的选取并不敏感，全局收敛。

缺点：迭代速度较慢。

代码：

function [x, val, k] = zuisuxiajiang(x0)
%使用最速下降法求解最小值
%x是迭代后的最小值点，val是最小值，k是迭代次数
%x0是初始值，最速下降法对于初始值不敏感
k = 0;
epsilon = 1e-10; %终止误差值
%rho和sigma是线搜索法的两个参数
rho = 0.55;
sigma = 0.5;
%当函数梯度的范数小于等于终止误差值时停止迭代
while norm(g(x0)) > epsilon
d = -g(x0);
m = 0;
while (f(x0 + rho^m*d) > (f(x0)+sigma*rho^m*g(x0)'*d))
m = m + 1;
end
x0 = x0 + rho^m*d;
k = k + 1;
end
x = x0;
val = f(x0);
end

%%原函数
function y = f(x)
y = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end

%%原函数的梯度,要弄成列向量的形式
function y = g(x)
y(1,1) = 400*x(1)*(x(1)^2 - x(2)) + 2*(x(1) - 1);
y(2,1) = -200*(x(1)^2 - x(2));
end

%%原函数的Hessen矩阵，注意矩阵每个元素的值
function y = G(x)
y(1, 1) = 1200*x(1)^2 - 400*x(2) + 2;
y(1, 2) = -400*x(1);
y(2, 1) = -400*x(1);
y(2, 2) = 200;
end

04

（4）信赖域方法
优点：信赖域方法是一种保证全局收敛的优化算法，并且可以直接确定位移，产生新的迭代点。
思想：（求出位移，并根据位移的相关信息调整信赖域半径）

算法：

代码：

function [xk,val,k]=trustm(x0)
%功能: 牛顿型信赖域方法求解无约束优化问题 min f(x)
%输入: x0是初始迭代点
%输出: xk是近似极小点, val是近似极小值, k是迭代次数
n=length(x0); x=x0; dta=1;
eta1=0.15; eta2=0.75; dtabar=2.0;
tau1=0.5; tau2=2.0; epsilon=1e-6;
k=0; Bk=Hess(x); %Bk=eye(n);
while(k<150)
gk=gfun(x); %%%% 梯度值
if(norm(gk) break;
end
[d,val,lam,ik]=trustq(gk,Bk,dta);%%%%%信頼域子问题
deltaq=-qk(x,d); %%%%预测下降
deltaf=fun(x)-fun(x+d);%%%实际下降
rk=deltaf/deltaq;%%%%%二者比值
if(rk<=eta1)
dta=tau1*dta;%%%%缩小信頼域半径，只有这种情况需要重新计算子问题
else if (rk>=eta2&norm(d)==dta)
dta=min(tau2*dta,dtabar); %%%%放大
else
dta=dta; %%%维持不动
end
end
if(rk>eta1)
x0=x; x=x+d; %%%产生新的迭代点
% sk=x-x0; yk=gfun(x)-gfun(x0);
%vk=sqrt(yk'*Bk*yk)*(sk/(sk'*yk)-Bk*yk/(yk'*Bk*yk));
%Bk=Bk-Bk*yk*yk'*Bk/(yk'*Bk*yk)+sk*sk'/(sk'*yk)+vk*vk'
%pause
Bk=Hess(x);
end
k=k+1;
end
xk=x;
val=fun(xk);
%%% 目标函数 %%%%%%%%%%%%%%%
function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
%%% 子问题目标函数 %%%%%%%%%%%%%
function qd=qk(x,d)
gk=gfun(x); Bk=Hess(x);
qd=gk'*d+0.5*d'*Bk*d;
%%% 目标函数的梯度 %%%%%%%%%%%%%%
function gf=gfun(x)
gf=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1), -200*(x(1)^2-x(2))]';
%%% 目标函数的Hesse阵 %%%%%%%%%%%%%%
function He=Hess(x)
n=length(x);
He=zeros(n,n);
He=[1200*x(1)^2-400*x(2)+2, -400*x(1);
-400*x(1), 200 ];

距离国赛仅剩2天，如果你仍然迷茫，不知道如何临阵磨枪，不妨看看两位国一学长备战国赛的心酸历程，点击下方小程序，让学长面对面与你交流竞赛流程、题目的选择技巧、需要掌握的基本知识、一起揣摩出题人思路：

国赛在即，数模乐园微小店正式上线，现小店已上架了50余种数学建模的相关产品：包括国赛真题讲解、超全优秀论文、必备模型总结、常用数学建模软件教程、国赛备赛大礼包等备赛资料一应俱全，各种备战数模好物点击下方图片立即进店挑选：

2022年新学期第一场保研综测加分的权威竞赛开放报名！