数学建模常用算法
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;
for i=1:N
if (x(i)^2+y(i)^2
count=count+1;
end
end
PI=4*count/N
试给出下面赌博中的蒙特卡洛模拟
在一次旅游途中，小王看到有人用20枚签(其中10枚标有5分分值，10枚标有10分分值) 设赌。让游客从中抽出10枚，以10枚签的分值总和为奖罚金额，见表1
表1
分值 50，100 55，95 60，65，85，90 70，75，80
奖罚金额 奖100元 奖10元 不奖不罚 罚1元
你看，有奖有罚，在11个分值中有4个分值可以获奖，且最高奖额为100元；只有3个分值要受罚，而罚额仅为1元，很有吸引力吧？怪不得有些游客摩拳擦掌，跃跃欲试。那么这些奖是不是这么好拿呢？
试分析此游戏中，谁是真正的赢家？
%%假设前10个分值为5，后10个分值为10
income=0; %% 收入
n=10000; %% 模拟次数, 即有n 个人参加游戏
for i=1:n
a=randperm(20);
a=a(1:10);
b=find(a>10); %%10分分值的
sumb=length(b)*10+(10-length(b))*5;
if sumb==50||sumb==100
income=income-100;
elseif sumb==55||sumb==95
income=income-10;
elseif sumb==70||sumb==75||sumb==80
income=income+1;
end
end
Income
2. 数据拟合、参数估计、插值等算法
数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是98 年美国赛A 题，生物组织切片的
三维插值处理，94 年A 题逢山开路，山体海拔高度的插值计算，还有吵的沸沸扬扬可能会考的“非典”问题也要用到数据拟合算法，观察
数据的走向进行处理。此类问题在MA TLAB 中有很多现成的函数可
以调用，熟悉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 个节点：(x i , y j , z ij ) (i=1,2.....m;j=1,2.....n)并且
(x 1
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=[636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
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 个节点：(x i , y i , z i ) (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=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162
162 117.5];
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5
-66.5 84 -33.5];
z=-[4 8 6 8 6 8 8 9 9 8 8 9
4 9];
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 解方程组方法
在上面的记号下， J (a 1,..., a m ) =RA -Y
Matlab 中的线性最小二乘的标准型为 min RA -Y A 222,
命令为:A=R\Y。
例 4 用最小二乘法求一个形如:
y =a +bx 2
的经验公式，使它与如下所示的数据拟合。
x 19 25 31 38 44
y 19.0 32.3 49.0 73.3 97.8
解 编写程序如下
x=[19 25 31 38 44]';
y=[19.0 32.3 49.0 73.3 97.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 多项式拟合方法
如果取:{r 1(x ), r 2(x )..... r m +1(x )}={1, x 1..... x m }
即用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 中线性规划的标准型为:
min z =C T X x
⎧Ax ≤b ⎪s . t ⎨Aeq . x =beq
⎪lb ≤x ≤ub ⎩
基本函数形式为 linprog(c,A,b)， 它的返回值是向量 x 的值。 还有其它的一
些函数调用形式（在 Matlab 指令窗运行 help linprog 可以看到所有的函数调
用形式), 如：
[x,fval]=linprog(c,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
这里 fval 返回目标函数的值， LB 和 UB 分别是变量 x 的下界和上界，x0是
x 的初始值，OPTIONS 是控制参数。
例 2 求解下列线性规划问题
max z =2x 1+3x 2-5x 3
⎧x 1+x 2+x 3=7⎪2x -5x +x ≥10 ⎪223s . t ⎨⎪x 1+3x 2+x 3≤12
⎪⎩x 1, x 2, x 3≥0
解 （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 已知非线性整数规划为：
max z =x 1+x 2+3x 3+4x 4+2x 5-8x 1-2x 2-3x 3-x 4-2x 5
⎧0≤x i ≤99(i =1... 5)⎪x +x +x +x +x ≤40012345 ⎪⎪s . t ⎨x 1+2x 2+2x 3+x 4+6x 5≤800
⎪2x +x +6x ≤2003⎪12
⎪⎩x 3+x 4+5x 5≤20022222
如果用显枚举法试探，共需计算1010个点，其计算量非常之大。然而应用蒙特卡
洛去随机计算106个点，便可找到满意解，那么这种方法的可信度究竟怎样呢？
下面就分析随机取样采集106个点计算时，应用概率理论来估计一下可信度。不
失一般性，假定一个整数规划的最优点不是孤立的奇点。假设目标函数落在高值
区的概率分别为 0.01，0.00001，则当计算 个点后，有任一个点能落在高值区
的概率分别为:
1-0. 991000000≈0. 99..... 99(100多位)
1-0. [1**********]0 ≈0. 999954602。
解 （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
for i=1:10^6
x=99*rand(5,1);
x1=floor(x);x2=ceil(x);
[f,g]=mengte(x1);
if sum(g
if p0
x0=x1;p0=f;
end
end
[f,g]=mengte(x2);
if sum(g
if p0
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=1 1 1 1 1
1 2 2 1 6
2 1 6 0 0
0 0 1 1 5;
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 定义目标函数
function f=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,...)
function f=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 指令中运行 help quadprog 后的帮助） 。 例 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 定义目标函数如下： function f=fun6(x,s); f=sum((x-0.5).^2);
（2）编写 M 文件 fun7.m 定义约束条件如下： function [c,ceq,k1,k2,s]=fun7(x,s); c=[];ceq=[];
if isnan(s(1,1)) s=[0.2,0;0.2 0]; 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 定义向量函数如下： function f=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.4 Matlab 优化工具箱的用户图形界面解法:
Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。 optimtool 可应用到所有优化问题的求解，计算结果可以输出到 Matlab 工作空间中。
2.4 图论问题
98 年B 题、00 年B 题、95 年锁具装箱等问题体现了图
论问题的重要性，这类问题算法有很多，包括：Dijkstra 、Floyd 、Prim 、Bellman-Ford ，最大流，二分匹配等问题。每一个算法都应该实现一遍，否则到比赛时再写就晚了。 2.4.1 两个指定点间的最短路径
Dijkstra 算法（单源最短路径边权非负）
w=[0 8 inf inf inf inf 7 inf inf inf inf;inf 0 3 inf inf inf 6 inf inf inf inf;...
inf inf 0 5 6 inf inf inf inf inf inf;inf inf inf 0 1 inf inf inf inf inf 12;...
inf inf inf inf 0 2 inf inf inf 9 inf;inf inf inf inf inf 0 9 inf 3 inf inf;...
inf inf 5 inf inf inf 0 10 inf inf inf;8 inf inf inf inf inf inf 0 inf inf inf;...
inf inf inf inf 7 inf inf 9 0 inf inf;inf inf inf inf inf inf inf inf 2 0 2;...
inf inf inf inf 10 inf inf inf inf inf 0] n=size(w,1); w1=w(1,:); for i=1:n
dist(i)=w1(i); path(i)=1; end s=[]; s(1)=1; u=s(1); k=1 dist path
while k
for j=1:k
if i~=s(j)
if dist(i)>dist(u)+w(u,i) dist(i)=dist(u)+w(u,i); path(i)=u; end end end end dist path
distdist=dist; for i=1:n
for j=1:k
if i~=s(j)
distdist(i)=distdist(i); else
distdist(i)=inf; end end end
distv=inf; for i=1:n
if distdist(i)
s(k+1)=v k=k+1 u=s(k) end dist path
Linggo 代码：
编写 LINGO 程序如下： model: sets:
cities/A,B1,B2,C1,C2,C3,D/;
roads(cities,cities)/A B1,A B2,B1 C1,B1 C2,B1 C3,B2 C1, B2 C2,B2 C3,C1 D,C2 D,C3 D/:w,x; endsets data:
w=2 4 3 3 1 2 3 1 1 3 4; enddata
n=@size(cities); !城市的个数; min=@sum(roads:w*x);
@for(cities(i)|i #ne#1 #and# i #ne#n:
@sum(roads(i,j):x(i,j))=@sum(roads(j,i):x(j,i))); @sum(roads(i,j)|i #eq#1:x(i,j))=1; @sum(roads(i,j)|j #eq#n:x(i,j))=1; end
编写 LINGO 程序如下： model: sets:
cities/1..11/;
roads(cities,cities):w,x; endsets data: w=0; enddata calc:
w(1,2)=2;w(1,3)=8;w(1,4)=1; w(2,3)=6;w(2,5)=1;
w(3,4)=7;w(3,5)=5;w(3,6)=1;w(3,7)=2; w(4,7)=9;
w(5,6)=3;w(5,8)=2;w(5,9)=9; w(6,7)=4;w(6,9)=6; w(7,9)=3;w(7,10)=1; w(8,9)=7;w(8,11)=9;
w(9,10)=1;w(9,11)=2;w(10,11)=4;
@for(roads(i,j):w(i,j)=w(i,j)+w(j,i));
@for(roads(i,j):w(i,j)=@if(w(i,j) #eq# 0, 1000,w(i,j))); endcalc
n=@size(cities); !城市的个数; min=@sum(roads:w*x);
@for(cities(i)|i #ne#1 #and# i #ne#
n:@sum(cities(j):x(i,j))=@sum(cities(j):x(j,i))); @sum(cities(j):x(1,j))=1;
@sum(cities(j):x(j,1))=0; !不能回到顶点1; @sum(cities(j):x(j,n))=1; @for(roads:@bin(x)); end
与有向图相比较，在程序中只增加了一个语句@sum(cities(j):x(j,1))=0，即 从顶点 1 离开后，再不能回到该顶点。
求得的最短路径为 1→2→5→6→3→7→10→9→11，最短路径长度为 13。
2.4.2 Flord算法
%a=[0 50 inf inf inf;inf 0 inf inf 80;inf 30 0 20 inf;inf inf inf 0 70;65 inf 100 inf 0];
function [D,path]=floyd(a) % Floyd’s Algorithm n=size(a,1);
%设置D 和Path 的初值 D=a;path=zeros(n,n); for i=1:n for j=1:n if D(i,j)~=inf
path(i,j)=j; %j是i 的后继点 end end end
%做n 次迭代, 每次迭代均更新D(i,j)和path(i,j) for k=1:n for i=1:n for j=1:n
if D(i,k)+D(k,j)
D(i,j)=D(i,k)+D(k,j) %修改长度 path(i,j)=path(i,k) %修改路径 end end end end
%a=[0 50 inf inf inf;inf 0 inf inf 80;inf 30 0 20 inf;inf inf inf 0 70;65 inf 100 inf 0]; %[D,path]=floyd(a)
2.4.3 prim算法（最优选线问题，最小生成树）
clear;
a=[0 8 inf 1 5; 8 0 6 inf 7; inf 6 0 9 10; 1 inf 9 0 3; 5 7 10 3 0];
T=[];c=0;v=1;n=5;sb=2:n;
for j=2:n b(1,j-1)=1; b(2,j-1)=j; b(3,j-1)=a(1,j); end
while size(T,2)
[tmin,i]=min(b(3,:)); T(:,size(T,2)+1)=b(:,i); c=c+b(3,i); v=b(2,i);
temp=find(sb == b(2,i)); sb(temp)=[];b(:,i)=[]; for j=1:length(sb) d=a(v,b(2,j)); if d
b(1,j)=v;b(3,j)=d; end end end T,c
2.4.4 Bellman-Ford算法(单源最短路径边权可正可负)
function [D,R]= floydwarshall(A)
% %采用floyd 算法计算图中任意两点之间最短路程，可以有负权。 %参数D 为连通图的权矩阵 %R是路由矩阵
D=A; n=length(D); %赋初值
for (i=1:n)for (j=1:n)R(i,j)=j;end ; end %赋路径初值 for (k=1:n) for (i=1:n) for (j=1:n)
if (D(i,k)+D(k,j)
R(i,j)=k;end ; end ; end %更新rij ，需要通过k k; %显示迭代步数 D; %显示每步迭代后的路长 R; %显示每步迭代后的路径 pd=0;for i=1:n %含有负权时
if (D(i,i)
if (pd==1) fprintf(' 有负回路' ); break ; end %存在一条负回路, 跳出最外层循环 终止程序
end %程序结束
2.4.5 最大流算法
function [flow,val]=mincostmaxflow(rongliang,cost,flowvalue); %第一个参数：容量矩阵；第二个参数：费用矩阵； %前两个参数必须在不通路处置零
%第三个参数：指定容量值（可以不写，表示求最小费用最大流） %返回值flow 为可行流矩阵，val 为最小费用值 global M
flow=zeros(size(rongliang));allflow=sum(flow(1,:)); if nargin
while allflow
w=(flow0).*cost)'; path=floydpath(w);%调用floydpath 函数 if isempty(path)
val=sum(sum(flow.*cost)); return ; end
theta=min(min(path.*(rongliang-flow)+(path.*(rongliang-flow)==0).*M));
theta=min([min(path'.*flow+(path'.*flow==0).*M),theta]); flow=flow+(rongliang>0).*(path-path').*theta; allflow=sum(flow(1,:)); end
val=sum(sum(flow.*cost));
2.4.6 二分匹配
二分图最大匹配 匈牙利算法
#include #include const int N=502; int map[N][N];
int link[N],useif[N]; int n,m; int can(int t)
{
for (int i=1;i
if (useif[i]==0 && map[t][i]) {
useif[i]=1;
if (link[i]==-1 || can(link[i])) {
link[i]=t; return 1; } } }
return 0; }
int maxmatch() {
int i,num; num=0;
memset(link,-1, sizeof (link)); for (i=1;i
memset(useif,0, sizeof (useif)); if (can(i)) num++; }
return num; }
int main() {
int k,a,b;
while (scanf("%d",&k),k) {
scanf("%d%d",&n,&m); memset(map,0, sizeof (map)); while (k--) {
scanf("%d%d",&a,&b); map[a][b]=1; }
printf("%d\n",maxmatch()); }
return 0; }
2.5 计算机算法设计中的问题
计算机算法设计包括很多内容：动态规划、回溯搜索、分治算法、分支定界。比如92 年B 题用分枝定界法，97 年B 题是典型的动态规划问题，此外98 年B 题体现了分治算法。这方面问题和ACM 程序设计竞赛中的问题类似，推荐看一下《计算机算法设计与分析》（电子工业出版社）等与计算机算法有关的书。 2.5.1 动态规划
如果一个问题能用动态规划方法求解，那么，我们可以按下列步骤首先建立起动态规划的数学模型：
（i ）将过程划分成恰当的阶段。 （ii ）正确选择状态变量时确定允许状态集合
X k
x k
，使它既能描述过程的状态，又满足无后效性，同
()（iii ）选择决策变量u k ，确定允许决策集合U k x k
（iv ）写出状态转移方程。
（v ）确定阶段指标v k (x k , u k )及指标函数V kn 的形式（阶段指标之和，阶段指标之积，阶段指标之极大或极小等）。
（vi ）写出基本方程即最优值函数满足的递归方程，以及端点条件。
model:
Title Dynamic Programming; sets:
vertex/A,B1,B2,C1,C2,C3,C4,D1,D2,D3,E1,E2,E3,F1,F2,G/:L;
road(vertex,vertex)/A B1,A B2,B1 C1,B1 C2,B1 c3,B2 C2,B2 C3,B2 C4, C1 D1,C1 D2,C2 D1,C2 D2,C3 D2,C3 D3,C4 D2,C4 D3, D1 E1,D1 E2,D2 E2,D2 E3,D3 E2,D3 E3,
E1 F1,E1 F2,E2 F1,E2 F2,E3 F1,E3 F2,F1 G,F2 G/:D; endsets data:
D=5 3 1 3 6 8 7 6 6 8 3 5 3 3 8 4 2 2 1 2 3 3
3 5 5 2 6 6 4 3; L=0,,,,,,,,,,,,,,,; enddata
@for(vertex(i)|i#GT#1:L(i)=@min(road(j,i):L(j)+D(j,i))); End
2.5.2分支定界
2.6.2 神经网络
下面利用上文所叙述的网络结构及方法，对蠓虫分类问题求解。编写 Matlab 程序 如下： clear
p1=[1.24,1.27;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90; 1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08]; p2=[1.14,1.82;1.18,1.96;1.20,1.86;1.26,2.00 1.28,2.00;1.30,1.96]; p=[p1;p2]'; pr=minmax(p);
goal=[ones(1,9),zeros(1,6);zeros(1,9),ones(1,6)]; plot(p1(:,1),p1(:,2),'h',p2(:,1),p2(:,2),'o') net=newff(pr,[3,2],{'logsig','logsig'});
net.trainParam.show = 10; net.trainParam.lr = 0.05; net.trainParam.goal = 1e-10; net.trainParam.epochs = 50000; net = train(net,p,goal);
x=[1.24 1.80;1.28 1.84;1.40 2.04]'; y0=sim(net,p) y=sim(net,x)
层次分析法：
根据层次总排序权值，该生最满意的工作为工作 1。 计算的 Matlab 程序如下：
clc,clear
fid=fopen('txt3.txt','r'); n1=6;n2=3; a=[];
for i=1:n1
tmp=str2num(fgetl(fid));
a=[a;tmp]; %读准则层判断矩阵 end
for i=1:n1
str1=char(['b',int2str(i),'=[];']);
str2=char(['b',int2str(i),'=[b',int2str(i),';tmp];']); eval(str1); for j=1:n2
tmp=str2num(fgetl(fid));
eval(str2); %读方案层的判断矩阵 end end
ri=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45]; %一致性指标 [x,y]=eig(a);
lamda=max(diag(y));
num=find(diag(y)==lamda); w0=x(:,num)/sum(x(:,num)); cr0=(lamda-n1)/(n1-1)/ri(n1) for i=1:n1
[x,y]=eig(eval(char(['b',int2str(i)]))); lamda=max(diag(y));
num=find(diag(y)==lamda); w1(:,i)=x(:,num)/sum(x(:,num)); cr1(i)=(lamda-n2)/(n2-1)/ri(n2); end
cr1, ts=w1*w0, cr=cr1*w0
纯文本文件txt3.txt 中的数据格式如下： 1 1 1 4 1 1/2 1 1 2 4 1 1/2 1 1/2 1 5 3 1/2 1/4 1/4 1/5 1 1/3 1/3 1 1 1/3 3 1 1 2 2 2 3 3 1 1 1/4 1/2 4 1 3 2 1/3 1 1 1/4 1/5 4 1 1/2 5 2 1
1 3 1/3
1/3 1 1/7
3 7 1
1 1/3 5
3 1 7
1/5 1/7 1
1 1 7
1 1 7
1/7 1/7 1
1 7 9
1/7 1 1
1/9 1 1
模糊算法(折衷型多目标决策优化)
mohu.txt
290 85 90 100 85 90 100 75 80 85 75 80 85 288 85 90 100 75 80 85 85 90 100 60 70 75
288 75 80 85 85 90 100 50 55 60 60 70 75
285 85 90 100 75 80 85 75 80 85 75 80 85
283 75 80 85 85 90 100 75 80 85 75 80 85
283 75 80 85 50 55 60 85 90 100 75 80 85
280 85 90 100 75 80 85 60 70 75 75 80 85
280 75 80 85 85 90 100 85 90 100 60 70 75
280 75 80 85 75 80 85 85 90 100 75 80 85
280 50 55 60 75 80 85 85 90 100 60 70 75
278 50 55 60 60 70 75 75 80 85 85 90 100
277 85 90 100 75 80 85 60 70 75 85 90 100
275 75 80 85 60 70 75 50 55 60 85 90 100
275 50 55 60 75 80 85 85 90 100 75 80 85
274 85 90 100 75 80 85 60 70 75 75 80 85
273 75 80 85 85 90 100 75 80 85 60 70 75
%数据存在纯文本文件 mohu.txt 中
clc,clear
load mohu.txt
sj=[repmat(mohu(:,1),1,3),mohu(:,2:end)];
%首先进行归一化处理
m=size(sj,1); n=size(sj,2)/3;
w=[0.5*ones(1,3),0.125*ones(1,12)];
w=repmat(w,m,1);
y=[];
%收益型指标对应的模糊指标值归一化处理
for i=1:n
tm=sj(:,3*i-2:3*i);
max_t=max(tm);
max_t=repmat(max_t,m,1);
max_t=max_t(:,3:-1:1);
yt=tm./max_t; yt(:,3)=min([yt(:,3)',ones(1,m)]); y=[y,yt];
end
%成本型指标对应的模糊指标值归一化处理
% for i=1:n
% tm=sj(:,3*i-2:3*i);
% min_t=min(tm);
% min_t=repmat(min_t,m,1);
% tm=tm(:,3:-1:1);
% yt=min_t./tm; yt(:,3)=min([yt(:,3)',ones(1,m)]); % y=[y,yt];
% end
%下面求模糊决策矩阵
r=[];
for i=1:n
tm1=y(:,3*i-2:3*i);tm2=w(:,3*i-2:3*i);
r=[r,tm1.*tm2];
end
%求 M ＋、M －和距离
mplus=max(r);mminus=min(r);
dplus=dist(mplus,r');dminus=dist(mminus,r');
%求隶属度
mu=dminus./(dplus+dminus);
[mu_sort,ind]=sort(mu,'descend')