数学建模算法整理.docx
《数学建模算法整理.docx》由会员分享，可在线阅读，更多相关《数学建模算法整理.docx（36页珍藏版）》请在冰豆网上搜索。
数学建模算法整理
数学建模常用算法
1.大多数建模赛题中都离不开计算机仿真，随机性模拟是非常常见的算法之一。
举个例子就是97年的A题，每个零件都有自己的标定值，也都有自己的容差等级，而求解最优的组合方案将要面对着的是一个极其复杂的公式和108种容差选取方案，根本不可能去求解析解，那如何去找到最优的方案呢？
随机性模拟搜索最优方案就是其中的一种方法，在每个零件可行的区间中按照正态分布随机的选取一个标定值和选取一个容差值作为一种方案，然后通过蒙特卡罗算法仿真出大量的方案，从中选取一个最佳的。
另一个例子就是去年的彩票第二问，要求设计一种更好的方案，首先方案的优劣取决于很多复杂的因素，同样不可能刻画出一个模型进行求解，只能靠随机仿真模拟。
1.1蒙特卡罗算法
蒙特卡罗模拟
就是随机数相关的东西，你只要知道随机数是怎么得到。
其它的事就要好办了。
rand（m,n）产生m*n均匀随机数。
ex:
用概率方法求pi
N=100000;
x=rand（N,1）;
y=rand（N,1）;
count=0;
fori=1:
N
if（x（i）^2+y（i）^2<=1）
count=count+1;
end
end
PI=4*count/N
试给出下面赌博中的蒙特卡洛模拟
在一次旅游途中，小王看到有人用20枚签（其中10枚标有5分分值，10枚标有10分分值）设赌。
让游客从中抽出10枚，以10枚签的分值总和为奖罚金额，见表1
表1
分值50，10055，9560，65，85，9070，75，80
奖罚金额奖100元奖10元不奖不罚罚1元
你看，有奖有罚，在11个分值中有4个分值可以获奖，且最高奖额为100元；只有3个分值要受罚，而罚额仅为1元，很有吸引力吧？
怪不得有些游客摩拳擦掌，跃跃欲试。
那么这些奖是不是这么好拿呢？
试分析此游戏中，谁是真正的赢家？
%%假设前10个分值为5，后10个分值为10
income=0;%%收入
n=10000;%%模拟次数,即有n个人参加游戏
fori=1:
n
a=randperm（20）;
a=a（1:
10）;
b=find（a>10）;%%10分分值的
sumb=length（b）*10+（10-length（b））*5;
ifsumb==50||sumb==100
income=income-100;
elseifsumb==55||sumb==95
income=income-10;
elseifsumb==70||sumb==75||sumb==80
income=income+1;
end
end
Income
2. 数据拟合、参数估计、插值等算法
数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是98年美国赛A题，生物组织切片的三维插值处理，94年A题逢山开路，山体海拔高度的插值计算，还有吵的沸沸扬扬可能会考的“非典”问题也要用到数据拟合算法，观察数据的走向进行处理。
此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。
2.1三次样条插值在Matlab中的实现
在Matlab中数据点称之为断点。
如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。
这个条件强迫第1个和第2个三多项式的三阶导数相等。
对最后一个和倒数第2个三次多项式也做同样地处理。
Matlab中三次样条插值也有现成的函数：
y=interp1（x0,y0,x,'spline'）；
y=spline（x0,y0,x）；
pp=csape（x0,y0,conds），y=ppval（pp,x）。
其中x0,y0是已知数据点,x是插值点,y是插值点的函数值。
对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。
pp=csape（x0,y0）：
使用默认的边界条件，即Lagrange边界条件。
pp=csape（x0,y0,conds）中的conds指定插值的边界条件，其值可为：
'complete'边界为一阶导数，即默认的边界条件
'not-a-knot'非扭结条件
'periodic'周期条件
'second'边界为二阶导数，二阶导数的值[0,0]。
'variational'设置边界的二阶导数值为[0,0]。
对于一些特殊的边界条件，可以通过conds的一个
1×2矩阵来表示，conds元素的取值为1，2。
此时，使用命令
pp=csape（x0,y0_ext,conds）
其中y0_ext=[left,y0,right]，这里left表示左边界的取值，right表示右边界的取值。
conds（i）=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条
件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶
导数，对应的值由left和right给出。
2.2二维插值
前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。
若节点是二维的，插值函数就是二元函数，即曲面。
如在某区域测量了若干点（节点）的高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）。
2.1.1插值节点为网格节点
已知m×n个节点：
（i=1,2.....m;j=1,2.....n）并且
。
求点（x,y）处的插值z。
Matlab中有一些计算二维插值的程序。
如
z=interp2（x0,y0,z0,x,y,'method'）
其中x0,y0分别为m维和n维向量，表示节点，z0为n×m维矩阵，表示节点值，x,y为一维数组，表示插值点，x与y应是方向不同的向量，即一个是行向量，另一个是列向量，z为矩阵，它的行数为y的维数，列数为x的维数，表示得到的插值，'method'的用法同上面的一维插值。
如果是三次样条插值，可以使用命令
pp=csape（{x0,y0},z0,conds,valconds），z=fnval（pp,{x,y}）
其中x0,y0分别为m维和n维向量，z0为
m×n维矩阵，z为矩阵，它的行数为x的维数，列数为y的维数，表示得到的插值，具体使用方法同一维插值。
例2在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程如2表，试插值一曲面，确定合适的模型，并由此找出最高点和该点的高程。
解：
编写程序如下：
clear,clc
x=100:
100:
500;
y=100:
100:
400;
z=[636697624478450
698712630478420
680674598412400
662626552334310];
pp=csape（{x,y},z'）
xi=100:
10:
500;yi=100:
10:
400
cz1=fnval（pp,{xi,yi}）
cz2=interp2（x,y,z,xi,yi','spline'）
[i,j]=find（cz1==max（max（cz1）））
x=xi（i）,y=yi（j）,zmax=cz1（i,j）
2.1.2插值节点为散乱节点
已知n个节点：
（i=1,2.....n）
求点（x，y）处的插值z。
对上述问题，Matlab中提供了插值函数griddata，其格式为：
ZI=GRIDDATA（X,Y,Z,XI,YI）
其中X、Y、Z均为n维向量，指明所给数据点的横坐标、纵坐标和竖坐标。
向量XI、YI是给定的网格点的横坐标和纵坐标，返回值ZI为网格（XI，YI）处的函数值。
XI与YI应是方向不同的向量，即一个是行向量，另一个是列向量。
例3在某海域测得一些点（x,y）处的水深z由下表给出，在矩形区域（75,200）
×（-50,150）内画出海底曲面的图形。
解：
编写程序如下：
x=[129140103.588185.5195105157.5107.57781162162117.5];
y=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];
z=-[48686889988949];
xi=75:
1:
200;
yi=-50:
1:
150;
zi=griddata（x,y,z,xi,yi','cubic'）
subplot（1,2,1）,plot（x,y,'*'）
subplot（1,2,2）,mesh（xi,yi,zi）
2.2最小二乘法的Matlab实现
2.2.1解方程组方法
在上面的记号下，
Matlab中的线性最小二乘的标准型为
命令为:
A=R\Y。
例4用最小二乘法求一个形如:
的经验公式，使它与如下所示的数据拟合。
x1925313844
y19.032.349.073.397.8
解编写程序如下
x=[1925313844]';
y=[19.032.349.073.397.8]';
r=[ones（5,1）,x.^2];
ab=r\y
x0=19:
0.1:
44;
y0=ab
（1）+ab
（2）*x0.^2;
plot（x,y,'o',x0,y0,'r'）
2.2.2多项式拟合方法
如果取:
{
}={
}
即用m次多项式拟合给定数据，Matlab中有现成的函数
a=polyfit（x0,y0,m）
其中输入参数x0,y0为要拟合的数据，m为拟合多项式的次数，输出参数a为拟合多项式y=amxm+…+a1x+a0系数a=[am,…,a1,a0]。
多项式在x处的值y可用下面的函数计算
y=polyval（a,x）。
2.3 规划类问题算法
竞赛中很多问题都和数学规划有关，可以说不少的模型都可以归结为一组不等式作为约束条件、几个函数表达式作为目标函数的问题，遇到这类问题，求解就是关键了，比如98年B题，用很多不等式完全可以把问题刻画清楚，因此列举出规划后用Lindo、Lingo等软件来进行解决比较方便，所以还需要熟悉这两个软件。
2.3.1求解线性规划的Matlab解法
单纯形法是求解线性规划问题的最常用、最有效的算法之一。
下面我们介绍线性规划的Matlab
解法。
Matlab中线性规划的标准型为:
基本函数形式为linprog（c,A,b），它的返回值是向量x的值。
还有其它的一些函数调用形式（在Matlab指令窗运行helplinprog可以看到所有的函数调用形式）,如：
[x,fval]=linprog（c,A,b,Aeq,beq,LB,UB,X0,OPTIONS）
这里fval返回目标函数的值，LB和UB分别是变量x的下界和上界，x0是x的初始值，OPTIONS是控制参数。
例2求解下列线性规划问题
解（i）编写M文件
c=[2;3;-5];
a=[-2,5,-1;1,3,1];b=[-10;12];
aeq=[1,1,1];
beq=7;
x=linprog（-c,a,b,aeq,beq,zeros（3,1））
value=c'*x
§4蒙特卡洛法（随机取样法）
前面介绍的常用的整数规划求解方法，主要是针对线性整数规划而言，而对于非线性整数规划目前尚未有一种成熟而准确的求解方法，因为非线性规划本身的通用有效解法尚未找到，更何况是非线性整数规划。
然而，尽管整数规划由于限制变量为整数而增加了难度；然而又由于整数解是有限个，于是为枚举法提供了方便。
当然，当自变量维数很大和取值范围很宽情况下，企图用显枚举法（即穷举法）计算出最优值是不现实的，但是应用概率理论可以证明，在一定的计算量的情况下，完全可以得出一个满意解。
例7已知非线性整数规划为：
如果用显枚举法试探，共需计算
个点，其计算量非常之大。
然而应用蒙特卡洛去随机计算
个点，便可找到满意解，那么这种方法的可信度究竟怎样呢？
下面就分析随机取样采集
个点计算时，应用概率理论来估计一下可信度。
不失一般性，假定一个整数规划的最优点不是孤立的奇点。
假设目标函数落在高值区的概率分别为0.01，0.00001，则当计算个点后，有任一个点能落在高值区的概率分别为:
解（i）首先编写M文件mente.m定义目标函数f和约束向量函数g，程序如下：
function[f,g]=mengte（x）;
f=x
（1）^2+x
（2）^2+3*x（3）^2+4*x（4）^2+2*x（5）-8*x
（1）-2*x
（2）-3*x（3）-...
x（4）-2*x（5）;
g=[sum（x）-400
x
（1）+2*x
（2）+2*x（3）+x（4）+6*x（5）-800
2*x
（1）+x
（2）+6*x（3）-200
x（3）+x（4）+5*x（5）-200];
（ii）编写M文件mainint.m如下求问题的解：
rand（'state',sum（clock））;
p0=0;
tic
fori=1:
10^6
x=99*rand（5,1）;
x1=floor（x）;x2=ceil（x）;
[f,g]=mengte（x1）;
ifsum（g<=0）==4
ifp0<=f
x0=x1;p0=f;
end
end
[f,g]=mengte（x2）;
ifsum（g<=0）==4
ifp0<=f
x0=x2;p0=f;
end
end
end
x0,p0
toc
本题可以使用LINGO软件求得精确的全局最有解，程序如下：
model:
sets:
row/1..4/:
b;
col/1..5/:
c1,c2,x;
link（row,col）:
a;
endsets
data:
c1=1,1,3,4,2;
c2=-8,-2,-3,-1,-2;
a=11111
12216
21600
00115;
b=400,800,200,200;
enddata
max=@sum（col:
c1*x^2+c2*x）;
@for（row（i）:
@sum（col（j）:
a（i,j）*x（j））@for（col:
@gin（x））;
@for（col:
@bnd（0,x,99））;
end
解（i）编写M文件fun1.m定义目标函数
functionf=fun1（x）;
f=sum（x.^2）+8;
（ii）编写M文件fun2.m定义非线性约束条件
function[g,h]=fun2（x）;
g=[-x
（1）^2+x
（2）-x（3）^2
x
（1）+x
（2）^2+x（3）^3-20];%非线性不等式约束
h=[-x
（1）-x
（2）^2+2
x
（2）+2*x（3）^2-3];%非线性等式约束
（iii）编写主程序文件example2.m如下：
options=optimset（'largescale','off'）;
[x,y]=fmincon（'fun1',rand（3,1）,[],[],[],[],zeros（3,1）,[],...
'fun2',options）
解：
编写M文件fun2.m如下：
function[f,g]=fun2（x）;
f=100*（x
（2）-x
（1）^2）^2+（1-x
（1））^2;
g=[-400*x
（1）*（x
（2）-x
（1）^2）-2*（1-x
（1））;200*（x
（2）-x
（1）^2）];
编写主函数文件example6.m如下：
options=optimset（'GradObj','on'）;
[x,y]=fminunc（'fun2',rand（1,2）,options）即可求得函数的极小值。
在求极值时，也可以利用二阶导数，编写M文件fun3.m如下：
function[f,df,d2f]=fun3（x）;
f=100*（x
（2）-x
（1）^2）^2+（1-x
（1））^2;
df=[-400*x
（1）*（x
（2）-x
（1）^2）-2*（1-x
（1））;200*（x
（2）-x
（1）^2）];
d2f=[-400*x
（2）+1200*x
（1）^2+2,-400*x
（1）
-400*x
（1）,200];
编写主函数文件example62.m如下：
options=optimset（'GradObj','on','Hessian','on'）;
[x,y]=fminunc（'fun3',rand（1,2）,options）
即可求得函数的极小值。
求多元函数的极值也可以使用Matlab的fminsearch命令，其使用格式为：
[X,FVAL,EXITFLAG,OUTPUT]＝FMINSEARCH（FUN,X0,OPTIONS,P1,P2,...）
functionf=fun4（x）;
f=sin（x）+3;
编写主函数文件example7.m如下：
x0=2;
[x,y]=fminsearch（@fun4,x0）
即求得在初值2附近的极小点及极小值。
Matlab中求解二次规划的命令是
[X,FVAL]=QUADPROG（H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS）
返回值X是决策向量x的值，返回值FVAL是目标函数在x处的值。
（具体细节可以参
看在Matlab指令中运行helpquadprog后的帮助）。
例8求解二次规划
h=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9];
[x,value]=quadprog（h,f,a,b,[],[],zeros（2,1））
解
（1）编写M文件fun6.m定义目标函数如下：
functionf=fun6（x,s）;
f=sum（（x-0.5）.^2）;
（2）编写M文件fun7.m定义约束条件如下：
function[c,ceq,k1,k2,s]=fun7（x,s）;
c=[];ceq=[];
ifisnan（s（1,1））
s=[0.2,0;0.20];
end
%取样值
w1=1:
s（1,1）:
100;
w2=1:
s（2,1）:
100;
%半无穷约束
k1=sin（w1*x
（1））.*cos（w1*x
（2））-1/1000*（w1-50）.^2-sin（w1*x（3））-x（3）-1;
k2=sin（w2*x
（2））.*cos（w2*x
（1））-1/1000*（w2-50）.^2-sin（w2*x（3））-x（3）-1;
%画出半无穷约束的图形
plot（w1,k1,'-',w2,k2,'+'）;
（3）调用函数fseminf
在Matlab的命令窗口输入
[x,y]=fseminf（@fun6,rand（3,1）,2,@fun7）
即可。
解
（1）编写M文件fun8.m定义向量函数如下：
functionf=fun8（x）;
f=[2*x
（1）^2+x
（2）^2-48*x
（1）-40*x
（2）+304
-x
（1）^2-3*x
（2）^2
x
（1）+3*x
（2）-18
-x
（1）-x
（2）
x
（1）+x
（2）-8];
（2）调用函数fminimax
[x,y]=fminimax（@fun8,rand（2,1））
解
（1）编写M文件fun9.m定义目标函数及梯度函数：
function[f,df]=fun9（x）;
f=exp（x
（1））*（4*x
（1）^2+2*x
（2）^2+4*x
（1）*x
（2）+2*x
（2）+1）;
df=[exp（x
（1））*（4*x
（1）^2+2*x
（2）^2+4*x
（1）*x
（2）+8*x
（1）+6*x
（2）+1）;exp（x
（1））*（4*x
（2）
+4*x
（1）+2）];
（2）编写M文件fun10.m定义约束条件及约束条件的梯度函数：
function[c,ceq,dc,dceq]=fun10（x）;
c=[x
（1）*x
（2）-x
（1）-x
（2）+1.5;-x
（1）*x
（2）-10];
dc=[x
（2）-1,-x
（2）;x
（1）-1,-x
（1）];
ceq=[];dceq=[];
（3）调用函数fmincon，编写主函数文件example13.m如下：
%采用标准算法
options=optimset（'largescale','off'）;
%采用梯度
options=optimset（options,'GradObj','on','GradConstr','on'）;
[x,y]=fmincon（@fun9,rand（2,1）,[],[],[],[],[],[],@fun10,options）
3.4Matlab优化工具箱的用户图形界面解法:
Matlab优化工具箱中的optimtool命令提供了优化问题的用户图形界面解法。
optimtool可应用到所有优化问题的求解，计算结果可以输出到Matlab工作空间中。
2.4 图论问题
98年B题、00年B题、95年锁具装箱等问题体现了图论问题的重要性，这类问题算法有很多，包括：
Dijkstra、Floyd、Prim、Bellman-Ford，最大流，二分匹配等问题。
每一个算法都应该实现一遍，否则到比赛时再写就晚了。
2.4.1两个指定点间的最短路径
Dijkstra算法（单源最短路径边权非负）
w=[08infinfinfinf7infinfinfinf;inf03infinfinf6infinfinfinf;...
infinf056infinfinfinfinfinf;infinfinf01infinfinfinfinf12;...
infinfinfinf02infinfinf9inf;infinfinfinfinf09i
升级会员 