数学建模常用算法

1、数学建模常用算法数学建模常用算法数学数学建模校内培训建模校内培训2011年年8月月31日日目录：一.数据拟合数据拟合二二. 规划问题规划问题三三. 灰色预测灰色预测四四. 遗传算法遗传算法五五. 模拟退火算法模拟退火算法六六. 人工神经网络算法人工神经网络算法七七. 小波分析小波分析八八. 蒙特卡罗方法蒙特卡罗方法一一. 数据拟合、参数估计、插值等算法数据拟合、参数估计、插值等算法 数据拟合在很多赛题中有应用，与图形处理有关的问题很多与拟合有关系，一个例子就是98 年美国赛A 题，生物组织切片的三维插值处理，94 年A 题逢山开路，山体海拔高度的插值计算，观察数据的走向进行处理。此类问题在MA

2、TLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。曲线拟合的线性最小二乘法曲线拟合的线性最小二乘法 曲线拟合问题的提法是，已知一组（二维）数据，即平面上的n 个点寻求一个函数（曲线） y = f (x) ，使f (x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。线性最小二乘法是解决曲线拟合最常用的方法，基本思路是，令 其中是事先选定好的一组线性无关的函数，)()()()(21xraxraxraxfmm 是待定系数，拟合的准则是使 与 的距离的平方和最小，即为最小二乘准则。记为使 与 的距离的平方和最小，只需利用极值的必要条件 得到关于 的线性方程组，

3、该方程组的解即为所要求的待定系数而函数 的选取也是首要、关键的一步。一般需要经验分析与判断。常用的曲线有直线 多项式 （一般m=2,3, 不宜太高）双曲线（一支） 指数函数 已知一组数据，用什么样的曲线拟合最好，可以在直观判断的基础上，选几种曲线分别拟合，然后比较，看哪条曲线的最小二乘指标J 最小。ka, 2 , 1,niyi), 1(0mkaJkmaa,1niiiniimyxfaaJ12121)(),()(ixfmaa,1, 2 , 1,niyi)(ixf21axay11mmmaxaxay21axayxaeay21)(xrk多项式拟合方法多项式拟合方法如果取 即用m 次多项式拟合给定数据，M

4、atlab中有现成的函数 a=polyfit(x0,y0,m)其中输入参数x0,y0 为要拟合的数据，m 为拟合多项式的次数，输出参数a 为拟合多项式 y=a_mx_m+a_1x+a_0 系数a= a_m, , a_1, a_0。多项式在x 处的值y 可用下面的函数计算 y=polyval(a,x)x=1 2 3 4 5 6 7 8 9;y=9 7 6 3 -1 2 5 7 20;P=polyfit(x,y,3);xi=0:.2:10;yi=polyval(P,xi);plot(xi,yi,x,y,r*);, 1,),(11mmxxrxrsyms tx=0;0.4;1.2;2;2.8;3.6;

5、4.4;5.2;6;7.2;8;9.2;10.4;11.6;12.4;13.6;14.4;15;y=1;0.85;0.29;-0.27;-0.53;-0.4;-0.12;0.17;0.28;0.15;-0.03;-0.15;-0.071;0.059;0.09;0.032;-0.015;-0.02;%指定函数形式为f(t)=acos(kt)e(wt),进行拟合f=fittype(a*cos(k*t)*exp(w*t),independent,t,coefficients,a,k,w);cfun=fit(x,y,f) %显示拟合函数xi=0:.1:20;yi=cfun(xi);plot(x,y,r

6、*,xi,yi,b-);clearclc%读人口数据Y=33815 33981 34004 34165 34212 34327 34344 34458 34498 34476, 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513, 34515 34517 34519 34519 34521 34521 34523 34525 34525 34527T=1:30%因为logistic 曲线模型的基本形式为y=1/(a+be(-t),%所以只要令y=1/y,x=e(-t),就可以转化为线性模型 y=a+bx%下面利用Matlab进

7、行回归分析拟合计算for t=1:30 x(t)=exp(-t); y(t)=1/Y(t);end% 计算，并输出回归系数 Bc=zeros(30,1)+1;X=c,x;B=inv(X*X)*X*yfor i=1:30 %计算回归拟合值 z(i)=B(1,1)+B(2,1)*x(i); %计算离差 s(i)=y(i)-sum(y)/30; %计算误差 w(i)=z(i)-y(i);end %计算离差平方和SS=s*s;%回归误差平方和QQ=w*w;%计算回归平方和UU=S-Q;%计算 输出F检验值F=28*U/Q%计算 非线性回归模型的拟合值for j=1:30 Y(j)=1/(B(1,1)+

8、B(2,1)*exp(-j);endplot(T,Y)min E(k,a,b)=sum_j=110a+be-0.02Kt_j-C_j2tdata=linspace(100,1000,10);cdata=1e-05.*454 499 535 565 590 610 626 639 650 659;x0=0.2,0.05,0.05;x=lsqcurvefit(curvefun,x0,tdata,cdata) %对曲线使用自定义函数进行最小二乘拟合f=curvefun(x,tdata)plot(tdata,cdata,o,tdata,f,r-)function f=curvefun(x,tdata)

9、f=x(1)+x(2)*exp(-0.02*x(3)*tdata); %x(1)=a,x(2)=b,x(3)=k数据的可视化数据的可视化解决地形地貌问题的关键是要将未测量地点的高度用数据插值的方法求出来，然后就可以用Matlab 绘图方法画出来。x,y=meshgrid(1:10) ； %构造测量网络h=rand(10,10); %测量数据点 (此时为随机产生的数据)xi,yi=meshgrid(1:0.1:10) ； %构造插值网络hi=interp2(x,y,h,xi,yi,spline); %二维插值命令surf(hi); %画出地貌图xlabel(x),ylabel(y),zlabel

10、(h);2002 年年A 题是要绘制投影区域。需要建立车灯投影的数学模型，然后再题是要绘制投影区域。需要建立车灯投影的数学模型，然后再根据模型汇出投影效果图根据模型汇出投影效果图根据得到的线光源长度，用投点法可画出测试屏上反射光的亮区根据得到的线光源长度，用投点法可画出测试屏上反射光的亮区。p=0.03;x=25.0216;for y1=-0.002:0.0004:0.002 y0=(-0.036:0.001:0.036)*ones(1,73); z0=ones(73,1)*(-0.036:0.001:0.036); x0=(y0.2+z0.2)/(2*p); xn=(p3+4*x0*2*p.

11、*x0+p*(-4*y1*y0+3*2*p*x0)./(2*(p2+2*p*x0); yn=(2*p*x0.*y0+p2*(-y1+y0)+y1*(y0.2-z0.2)./(p2+2*p*x0); zn=(p2+2*p*x0+2*y1*y0).*z0./(p2+2*p*x0); y=y0+(yn-y0).*(x-x0)./(xn-x0); z=z0+(zn-z0).*(x-x0)./(xn-x0); plot(y,z,r.) xlabel(y) ylabel(z) hold onend 二二. 规划问题规划问题线性规划线性规划线性规划是在一组线性约束条件的限制下， 求一线性目标函数最大或最小的

12、问题。在解决实际问题时，把问题归结成一个线性规划数学模型是很重要的一步，但往往也是困难的一步，模型建立得是否恰当，直接影响到求解。而选适当的决策变量，是我们建立有效模型的关键之一。线性规划的Malab 标准形式 其中c 和 为n 为列向量，A、Aeq 为适当维数矩阵， b,beq 为适当维数的列向量。基本函数形式为linprog(c,A,b)，它的返回值是向量x 的值。还有其它的一些函数调用形式（在 Matlab 指令窗运行 help linprog 可以看到所有的函数调用形式），如：x,fval=linprog(c,A,b,Aeq,beq,VLB,VUB,X0,OPTIONS)这里fval

13、返回目标函数的值，LB 和UB 分别是变量x 的下界和上界，x0 是x 的初始值，OPTIONS 是控制参数。VUBXVLBbeqAeqXbAXtsXczTx.min非线性规划非线性规划如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。（i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。（ii）提出追求目标：经过资料分析，

14、根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。（iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。（iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。Matlab 中非线性规划的数学模型写成以下形式其

15、中f (x) 是标量函数，A, B, Aeq, Beq, LB,UB 是相应维数的矩阵和向量，C(x),Ceq(x) 是非线性向量函数。VUBXVLBXCeqXCBeqAeqXBAXtsxf0)(0)(.)(minMatlab 中的命令是X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)它的返回值是向量x ，其中FUN 是用M 文件定义的函数f (x)； 是x 的初始值；A,B,Aeq,Beq 定义了线性约束A* X B, Aeq * X = Beq，如果没有线性约束，则A=,B=,Aeq=,Beq=；LB 和UB 是变量x 的下界和上界，如

16、果上界和下界没有约束，则LB=，UB=，如果x 无下界，则LB 的各分量都为-inf，如果x 无上界，则UB的各分量都为inf；NONLCON 是用M 文件定义的非线性向量函数C(x),Ceq(x) ；OPTIONS定义了优化参数，可以使用Matlab 缺省的参数设置。0 x线性规划线性规划f=-5;-4;-6;A=1 -1 1; 3 2 4; 3 2 0;b=20; 42; 30;lb=zeros(3,1);x,fval,exifflag, output, lambda=linprog(f,A,b,lb)%非线性规划非线性规划options=optimset; x,y=fmincon(fun

17、1,rand(2,1),zeros(2,1),fun2,options)function f=fun1(x);f=x(1)2+x(2)2+8;function g,h=fun2(x);g=-x(1)2+x(2);h=-x(1)-x(2)2+2;x0=1;1;1;1;lb=0;0;0;0;ub=;A=;b=;Aeq=;beq=;x,fval=fmincon(fun44,x0,A,b,Aeq,beq,lb,ub,mycon1)function g,ceq=mycon1(x)g(1)=x(1)-400;g(2)=1.1*x(1)+x(2)-440;g(3)=1.21*x(1)+1.1*x(2)+x(

18、3)-484;g(4)=1.331*x(1)+1.218*x(2)+1.1*x(3)+x(4)-532.4;ceq=0;function f=fun44(x)f=-(sqrt(x(1)+sqrt(x(2)+sqrt(x(3)+sqrt(x(4);整数规划整数规划规划中的变量（部分或全部）限制为整数时，称为整数规划。若在线性规划模型中，变量限制为整数，则称为整数线性规划。目前所流行的求解整数规划的方法，往往只适用于整数线性规划。目前还没有一种方法能有效地求解一切整数规划。随机取样计算方法：随机取样计算方法：尽管整数规划由于限制变量为整数而增加了难度，但由于整数解是有限个，于是为枚举法提供了方便。

19、应用概率理论可以证明在有一定计算量的情况下，随机取样计算方法完全可以得到一个满意解。例：function f,g=mengte(x);f=x(1)2+x(2)2+3*x(3)+4*x(4)2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-x(4)-2*x(5);g(1)=sum(x)-400;g(2)=x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-800;g(3)=2*x(1)+x(2)+6*x(3)-200;g(4)=x(3)+x(4)+5*x(5)-200;rand(state,sum(clock); %初始化随机数产生序列的办法,每次产生随机数的时候，随机数生成

20、器触发器的状态都会翻转一次。p0=0;ticfor i=1:105 x=99*rand(5,1); x1=floor(x);x2=ceil(x); f,g=mengte(x1); if sum(g=0)=4 if p0=f x0=x1;p0=f; endendf,g=mengte(x2);if sum(g=0)=4 if p0=f x0=x2;p0=f; endendendx0,p0toc灰色预测灰色预测灰色预测则是应用灰色模型对灰色系统进行分析、建模、求解、预测的过程。由于灰色建模理论应用数据生成手段，弱化了系统的随机性，使紊乱的原始序列呈现某种规律，规律不明显的变得较为明显，建模后还能进行

21、残差辨识，即使较少的历史数据，任意随机分布，也能得到较高的预测精度。因此，灰色预测在社会经济、管理决策、农业规划、气象生态等各个部门和行业都得到了广泛的应用。灰色预测通过鉴别系统因素之间发展趋势的相异程度，并对原始数据进行生成处理来寻找系统变动的趋势，生成有较强趋势的数据序列，然后建立相应的微分方程模型，从而预测事物的未来发展趋势。灰色预测的数据是通过生成数据的模型所得到的预测值的逆处理结果。在诸多灰色预测模型中，以灰色系统中单一序列一阶线性微分方程GM(1,1) 模型最为常用。下面简要介绍GM(1,1)模型。 已知某公司已知某公司1998-2008年利润为年利润为89677, 99215,

22、109655,120333,135823,159878,182321,209407,246619,300670，现要预测该，现要预测该公司未来几年的利润情况。公司未来几年的利润情况。clear syms a b;c=a b;A=89677, 99215, 109655, 120333,135823,159878,182321,209407,246619,300670;B=cumsum(A); %原始数据累加n=length(A);for i=1:(n-1) C(i)=(B(i)+B(i+1)/2; %生成累加矩阵end%计算待定参数的值D=A;D(1)=;D=D;E=-C;ones(1,n-1

23、);c=inv(E*E)*E*D;c=c;a=c(1);b=c(2);%预测后续数据F=;F(1)=A(1);for i=2:(n+10) F(i)=(A(1)-b/a)/exp(a*(i-1)+b/a;endG=;G(1)=A(1);for i=2:(n+10) G(i)=F(i)-F(i-1); %得出预测数据endt1=1999:2008;t2=1999:2018;G；plot(t1,A,o,t2,G)%2005A题题 长江水质的预测长江水质的预测clear syms a b;c=a b;A=174 179 183 189 207 234 220.5 256 270 285;B=cums

24、um(A);n=length(A);for i=1:(n-1) C(i)=(B(i)+B(i+1)/2;endD=A;D(1)=;D=D;E=-C;ones(1,n-1);c=inv(E*E)*E*D;c=c;a=c(1);b=c(2);F=;F(1)=A(1);for i=2:(n+10) F(i)=(A(1)-b/a)/exp(a*(i-1)+b/a;endG=;G(1)=A(1);for i=2:(n+10) G(i)=F(i)-F(i-1);endt1=1995:2004;t2=1995:2014;G;a,bplot(t1,A,o,t2,G)遗传算法 BEGIN t=0 %遗传代数 计

25、算P(t) 的适应值 while (不满足停止准则) do begin t=t+1 从P(t-1)中选择P (t) %选择(或复制) 重组 P(t) %交叉和变异 计算P(t) 的适应值 end END遗传算法的起源遗传算法的起源遗传算法是由美国的J. Holland教授于1975年在他的专著自然界和人工系统的适应性中首先提出的，它是一类借鉴生物界自然选择和自然遗传机制的随机化搜索算法 。遗传算法的搜索机制遗传算法的搜索机制 遗传算法模拟自然选择和自然遗传过程中发生的繁殖、交叉和基因突变现象，在每次迭代中都保留一组候选解，并按某种指标从解群中选取较优的个体，利用遗传算子(选择、交叉和变异)对这

26、些个体进行组合，产生新一代的候选解群，重复此过程，直到满足某种收敛指标为止。 基本遗传算法基本遗传算法（Simple Genetic Algorithms，简称SGA，又称简单遗传算法或标准遗传算法），是由Goldberg总结出的一种最基本的遗传算法，其遗传进化操作过程简单，容易理解，是其它一些遗传算法的雏形和基础。基本遗传算法的组成（1）编解码（产生初始种群）（2）适应度函数（3）遗传算子（选择、交叉、变异）（4）运行参数编码 GA是通过某种编码机制把对象抽象为由特定符号按一定顺序排成的串。正如研究生物遗传是从染色体着手，而染色体则是由基因排成的串。SGA使用二进制串进行编码。解码解码的目的

27、是为了将不直观的二进制数据还原为十进制。遗传算法的编码和解码在宏观上可以对应生物的基因型和表现型，在微观上可以对应DNA的转录和翻译两个过程。函数优化示例函数优化示例 求下列一元函数的最大值： x-1,2，求解结果精确到6位小数。0 .2)10sin()(xxxfSGA对于本例的编码 由于区间长度为3，求解结果精确到6位小数，因此可将自变量定义区间划分为3106等份。又因为221 3106 222 ，所以本例的二进制编码长度至少需要22位，本例的编码过程实质上是将区间-1，2内对应的实数值转化为一个二进制串（b21b20b0）。初始种群初始种群 SGA采用随机方法生成若干个个体的集合，该集合称

28、为初始种群。初始种群中个体的数量称为种群规模。适应度函数 遗传算法对一个个体（解）的好坏用适应度函数值来评价，适应度函数值越大，解的质量越好。适应度函数是遗传算法进化过程的驱动力，也是进行自然选择的唯一标准，它的设计应结合求解问题本身的要求而定。选择算子选择算子遗传算法使用选择运算来实现对群体中的个体进行优胜劣汰操作：适应度高的个体被遗传到下一代群体中的概率大；适应度低的个体，被遗传到下一代群体中的概率小。选择操作的任务就是按某种方法从父代群体中选取一些个体，遗传到下一代群体。SGA中选择算子采用轮盘赌选择方法。 轮盘赌选择方法轮盘赌选择方法 轮盘赌选择又称比例选择算子，它的基本思想是：各个个

29、体被选中的概率与其适应度函数值大小成正比。设群体大小为n ，个体i 的适应度为 Fi，则个体i 被选中遗传到下一代群体的概率为：niiiiFFP1/轮盘赌选择方法的实现步骤（1） 计算群体中所有个体的适应度函数值（需要解码）；（2） 利用比例选择算子的公式，计算每个个体被选中遗传到下一代群体的概率；（3） 采用模拟赌盘操作（即生成0到1之间的随机数与每个个体遗传到下一代群体的概率进行匹配）来确定各个个体是否遗传到下一代群体中。交叉算子所谓交叉运算，是指对两个相互配对的染色体依据交叉概率 Pc 按某种方式相互交换其部分基因，从而形成两个新的个体。交叉运算是遗传算法区别于其他进化算法的重要特征，它

30、在遗传算法中起关键作用，是产生新个体的主要方法。 SGA中交叉算子采用单点交叉算子。 单点交叉运算 交叉前：00000|0111000000001000011100|00000111111000101交叉后：00000|0000011111100010111100|01110000000010000交叉点交叉点变异算子所谓变异运算，是指依据变异概率 Pm 将个体编码串中的某些基因值用其它基因值来替换，从而形成一个新的个体。遗传算法中的变异运算是产生新个体的辅助方法，它决定了遗传算法的局部搜索能力，同时保持种群的多样性。交叉运算和变异运算的相互配合，共同完成对搜索空间的全局搜索和局部搜索。 SG

31、A中变异算子采用基本位变异算子。基本位变异算子 基本位变异算子是指对个体编码串随机指定的某一位或某几位基因作变异运算。对于基本遗传算法中用二进制编码符号串所表示的个体，若需要进行变异操作的某一基因座上的原有基因值为0，则变异操作将其变为1；反之，若原有基因值为1，则变异操作将其变为0 。基本位变异算子的执行过程 变异前：000001110000000010000变异后：000001110001000010000变异点变异点运行参数（1）M ： 种群规模 （2）T ： 遗传运算的终止进化代数 （3）Pc ： 交叉概率 （4）Pm ： 变异概率 SGA的框图 产生初始群体产生初始群体是否满足停止准

32、则是否满足停止准则是是输出结果并结束输出结果并结束计算个体适应度值计算个体适应度值比例选择运算比例选择运算单点交叉运算单点交叉运算基本位变异运算基本位变异运算否否产生新一代群体产生新一代群体执行执行M/2M/2次次运用运用GA 工具箱求解工具箱求解 cumcm 中多约束非线性规划问题中多约束非线性规划问题问题为问题为min f(x)=e(x_1)(4x_12+2x_22+4x_1x_2+2x_2+1)s.t. 1.5+x_1x_2-x_1-x_2=0; -x_1x_20|g210) f=100;else f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2

33、)+1);endBestx = -9.45946346276848 1.05596649888319BestFval = 0.02520168130380 模拟退火算法模拟退火算法算法简介模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程

34、。假设材料在状态i 之下的能量为E(i) ，那么材料在温度T 时从状态i 进入状态j 就遵循如下规律：（1）如果E( j) E(i) ，接受该状态被转换。（2）如果E( j) E(i) ，则状态转换以如下概率被接受： 其中K 是物理学中的波尔兹曼常数， T是材料温度。KTiEjEe)()(0TT 0 x算法基本步骤：算法基本步骤：（1）令 即开始退火的初始温度，随机生成一个初始解 ，并计算相应的目标函数值。（2）令 T 等于冷却进度表中的一个值 （3）根据当前解 进行扰动，产生一个新解 ，计算相应的目标函数值 ，得到 （4）若 ，则新解 被接受，作为新的当前解；若 ，则新解 按照概率 接受，

35、为当前温度。（5）在温度 下，重复Markov 链长度 次的扰动和接受过程，即执行（3）与（4）。（6）判断T 是否已到达 ，是，则终止算法；否，则转到（2）继续执行。ixjx0TT 0 x)()(ijxExEE)(jxE0Ejx0Ejx)/exp(iTEiTiTkLfTiT算法实质分为两循环，在任一温度随机扰动产生新解，并计算目标函数的变化，决定是否被接受。由于算法初始温度比较高，这样使得E增大的新解在初始是也可以被接受，因而能跳出局部最小值然后通过缓慢地减低温度，算法最终可能得到全局最优解。新解产生的机制的基本要求是能够尽量遍及解空间的各个区域，这样在某一恒定温度，不断产生新解时，就可能跳

36、出当前区域的极小以搜索其他区域，这是模拟退火算法能够进行广域搜索的一个重要条件。控制函数T的衰减函数有多种形式，一个常用的衰减函数是： 其中 是一个常数，可以取为0.50.99，它的取值，决定了降温的过程。Markov 链长度Markov 链长度的选取原则是：在控制参数T的衰减函数已选定的前提下， 应能使在控制参数T的每一取值上达到平衡。从经验上看，对简单的情况可以令 ，n 为问题规模。kkTT1kLnLk100旅行商问题(TSP)代笔一类组合优化问题，在物流、计算机网络、电子地图、交通疏导等方面都有重要的工程和理论价值，引起许多学者的关注。TSP简单描述为：一个商人要到n 个不同的地方去推销

37、商品，每2个城市i 和j 之间的距离为 ， 如何选择一条路径使得商人每个城市走一遍后回到起点，所走的路径最短。该问题当城市数目在100以上，一般很难精确的求出其全局最优解。说明：工程中许多实际优化问题的目标函数都是非凸的，存在许多局部最优解，特别是随着优化规模的增大，局部最优解的数目将迅速增加。有效地求出一般非凸目标函数的全局解至今仍是一个难题。一般确定性算法往往容易陷入局部而非全局最优。模拟退火算法是一种通用概率算法，用来在一个大的搜寻空间内寻找问题的最优解。具有高效、鲁棒、通用、灵活的优点。将模拟退火算法引入TSP求解，可以避免在求解过程中陷入TSP的局部最优。ijd算法设计步骤：1. 问

38、题的解空间和初始解问题的解空间是遍访每个城市恰好一次的所有回路，是所有城市排列的集合，即其中每一个排列 表示遍访n 个城市的一个路径， 表示第i 次访问城市j.初始解为随机生成一个 的随机排列作为 。2. 目标函数TSP 问题的目标函数即为访问所有城市的路径的总长度，也可以称为代价函数;现在TSP 问题的求解就是通过模拟退火算法求出目标函数 的最小值。相应地， 即为TSP 问题的最优解。), 2 , 1 (),( | ),(2121的排列为nccccccSnnisjcj, 2 , 1n0s111121),(),(),(niniinccdccdcccC),(21ncccC),(*2*1ncccS

39、)(iissCC0),/ exp(0, 1CTCP3. 新解的产生新解的产生对问题的求解非常重要。新解可通过分别或者交替使用以下两种方法来产生: (1) 二变换法：任选序号u,v (设uvn), 交换 u 和v 之间的访问顺序（2）三变换法：任选序号u,v，w (设u=vw), 将 u 和v 之间的路径插到w之后访问。4. 新解产生计算变换前的解和变换后目标函数的差值：5. Metropolis 接受准则以新解与当前解的目标函数差定义接受概率，即0-1 背包问题：有一个贼在偷窃一家商店时发现有背包问题：有一个贼在偷窃一家商店时发现有N件物品：第件物品：第i件物品件物品值值v_i 元，重元，重

40、w_i 磅（磅（1=i=tf for r=1:100 %产生随机扰动 tmp=ceil(rand.*num); sol_new(1,tmp)=sol_new(1,tmp); while 1 q=(sol_new*d=restriction); ifq p=p; %实现交错着逆转头尾的第一个1 tmp=find(sol_new=1); if p sol_new(1,tmp)=0; else sol_new(1,tmp(end)=0; end else break end end %计算背包中的物品价值 E_new=sol_new*k; if E_newE_current E_current=E_

41、new; sol_current=sol_new; if E_newE_best %把冷却过程中的解保存下来 E_best=E_new; sol_best=sol_new; end else if randexp(-(E_new-E_current)./t) E_current=E_new; sol_current=sol_new; else sol_new=sol_current; end end end t=t.*a;enddisp(最优解为：) sol_bestdisp(物品总价值等于;) val=-E_best;disp(val)disp(背包中物品重量是：)disp(sol_best

42、*d) 人工神经网络人工神经网络 人工神经网络具有人脑学习的记忆功能，主要应用到数据建模、预测、模式识别和函数优化等方向，在数学建模中也是经常被用到，如人工神经网络预测非典的传播、人口的增长等。 人工神经网络 是由大量简单的基本原件- 神经元相互连接，通过模拟人的大脑神经处理信息的方式，进行信息并行处理和非线性转换的复杂网络系统。其优点是多输入多输出实现了数据的并行处理及自学能力。它有三个基本要素：(i) 一组连接（对应于生物神经元的突触），连接强度由各连接上的权值表示，权值为正表示激活，为负表示抑制。（ii）一个求和单元，用于求取各输入信号的加权和（线性组合）。（iii）一个非线性激活函数，

43、起非线性映射作用并将神经元输出幅度限制在一定范围内（一般限制在(0,1) 或(1,1) 之间）。此外还有一个阈值 （或偏置 ）。以上作用可分别以数学式表达出来：式中 为输入信号，为神经元 k 之权值， 为线性组合结果， 为阈值(由刺激引起应激组织反应的最低值)， 为激活函数， 为神经元k 的输出。pxx,1kpkkwww,21ku)(kky激活函数 可以有以下几种：（i）阈值函数 即阶梯函数。这时，相应的输出 为 其中 ，常称此种神经元为M-P 模型。（ii）分段线性函数它类似于一个放大系数为1 的非线性放大器，当工作于线性区时它是一个线性组合器，放大系数趋于无穷大时变成一个阈值单元。（iii

44、）sigmoid 函数最常用的函数形式为 ，参数 可控制其斜率。 另一种常用的是双曲正切函数。这类函数具有平滑和渐进性，并保持单调性。0001)(vvv0001kkkvvyky)exp(11)(vv0Matlab 中的激活（传递）函数如下表所示： 各个函数的定义及使用方法，可以参看Matlab 的帮助。常见的神经网络理论常见的神经网络理论BP网络基本数学原理BP网络是一种多层前馈神经网络（各神经元接受前一层的输入，并输出给下一层，没有反馈）。BP 网络是一种具有三层或三层以上神经元的神经网络，包括输入层、中间层（隐含层）和输出层。隐含层或输出层将前一层所有的神经元传来的信息进行整合，通常还会在

45、整合过的信息中添加一个阈值，这主要是模仿生物学中神经元必须达到一定的阈值才会触发的原理，然后将整合过的信息作为该层神经元输入。当一对学习样本提供给输入神经后，神经元的激活值从输入层经过个隐含层向输出层传播，在输出层的个神经元获得网络输出相应，然后按照减少网络输出与实际输出样本之间误差方向，从输出层反向经过个隐含层，从而逐步修正个连接权值，这种算法交误差反向传播法，即BP算法。随着这种误差逆向传播修正的反复进行，网络对输入模式响应的正确率也不断上升。BP 算法的核心是数学中是“负梯度下降”理论，即BP网络的误差调整方向总是沿着误差下降最快的方向进行，常规三层BP网络权值和阈值调整公式如下：其中，

46、 E为网络输出与实际输出样本之间的误差平方和； 为网络的学习速率即权值调整幅度； 为t 时刻输入层第i 个神经元与隐含层第j 个神经元的连接权值； 为t+1 时刻输入层第i 个神经元与隐含层第j 个神经元的连接权值， 为t+1 时刻输隐含层第j 个神经元与输出层第k个神经元的连接权值; 为神经元的阈值，下标的意义与权值的相同。根据神经元具体激励函数表达式，可以具体地求出 和 。 最终网络通过负梯度 下降学习规则自行修正权值和阈值是误差平方和E逐步变小并最终达到理想误差。),() 1(twwEtwijijij),() 1(tEtijijij),()1(tEtjkjkjk),() 1(twwEtw

47、jkjkjkijwEjkwE)(twij) 1( twij) 1( twjk%公路运输主要包括公路客运量和公路客运量。公路运输主要包括公路客运量和公路客运量。%据研究，某地区公路运量主要与该地区的人数、机动车数量和公路面据研究，某地区公路运量主要与该地区的人数、机动车数量和公路面积有关。积有关。%数据如下（见下面程序）数据如下（见下面程序）%根据有关部门，该地区根据有关部门，该地区2010年年2011 年的人数分别为年的人数分别为73.39 和和75.55万万人，机动车数量为人，机动车数量为3.9635 和和4.0975万量，公路面积分别为万量，公路面积分别为0.9880和和1.0268 万平

48、方千米。万平方千米。%请利用请利用BP 网络预测该地区网络预测该地区2010 年和年和2011 年的公路客运量和货运量。年的公路客运量和货运量。function main()clcclear all;close all;SamNum=20; %输入样本数量为20TestSamNum=20; %测试样本也是20ForcastSamNum=2; %预测样本数量为2HiddenUnitNum=8; %中间层隐节点数量为8， 比工具箱多了1个InDim=3; %网络输入维度为3OutDim=2; %网络输出维度为2%原始数据%人数（单位：万人）sqrs=20.55 22.44 25.37 27.13

49、29.45 30.1 30.96 34.06 36.42 38.09 39.13 39.99 41.93 44.59 47.30 52.89 55.73 56.76 59.17 60.63;%机动车数（万辆）sqjdcs=0.6 0.75 0.85 0.90 1.05 1.35 1.45 1.60 1.70 1.85 2.15 2.2 2.25 2.35 2.5 2.6 2.7 2.85 2.95 3.1;%公路面积（万平方千米）sqglmj=0.09 0.11 0.11 0.14 0.20 0.23 0.23 0.32 0.32 0.34 0.36 0.36 0.38 0.49 0.56 0

50、.59 0.59 0.67 0.69 0.79;%公路客运量（万人）glhyl=5126 6217 7730 9145 10460 11387 12353 15750 18304 19836 21024 19490 20433 22598 25107 33442 36836 40548 42927 43462;%公路货运量（万吨）glkyl=1237 1379 1385 1399 1663 1714 1834 4322 8132 8936 11099 11203 10524 11115 13320 16762 18673 20724 20803 21804;p=sqrs;sqjdcs;sqgl

51、mj; %输入数据矩阵t=glkyl;glhyl; %目标数据矩阵SamIn, minp, maxp, tn, mint, maxt=premnmx(p,t); %原始样本对（输入和输出）初始化rand(state,sum(100*clock); %依据系统时钟种子产生随机数Noisevar=0.01; %噪音强度为0.01（添加噪音的目的是为了防止网络过度拟合）Noise=Noisevar*randn(2,SamNum); %生成噪音SamOut=tn + Noise; %将噪音添加到输出样本上TestSamIn=SamIn; %这里去输入样本与测试样本相同，因为样本容量偏少TestSamO

52、ut=SamOut; %也取输出样本与测试样本相同MaxEpochs=50000; %最多训练次数为50000lr=0.035; %学习速率为0.035E0=0.65*10(-3); %目标误差W1=0.5*rand(HiddenUnitNum,InDim)-0.1; %初始化输入层与隐含层之间的权值B1=0.5*rand(HiddenUnitNum,1)-0.1; %初始化输入层与隐含层之间的阈值W2=0.5*rand(OutDim,HiddenUnitNum)-0.1; %初始化输出层与隐含层之间的权值B2=0.5*rand(OutDim,1)-0.1; %初始化输出层与隐含层之间的阈值E

53、rrHistory=;for i=1:MaxEpochs HiddenOut=logsig(W1*SamIn+repmat(B1,1,SamNum); %隐含层网络输出 NetworkOut=W2*HiddenOut+repmat(B2,1,SamNum); %输出层网络输出 Error=SamOut-NetworkOut; %实际输出与网络输出之差 SSE=sumsqr(Error); %能量函数（误差平方和） ErrHistory=ErrHistory SSE; if SSEE0, break, end %如果达到误差要求，则跳出学习循环%以下6行是Bp网络的核心程序%它们是权值（阈值）依

54、据能量函数负梯度下降原理所作的每一步动态调整量Delta2=Error;Delta1=W2*Delta2.*HiddenOut.*(1-HiddenOut);dW2=Delta2*HiddenOut;dB2=Delta2*ones(SamNum,1);dW1=Delta1*SamIn;dB1=Delta1*ones(SamNum,1);%对输出层与隐含层之间的权值和估值进行修正W2=W2+lr*dW2;B2=B2+lr*dB2;%对输入层与隐含层之间的权值和估值进行修正W1=W1+lr*dW1;B1=B1+lr*dB1;endHiddenOut=logsig(W1*SamIn+repmat(B

55、1,1,TestSamNum); %隐含层输出最终结果NetworkOut=W2*HiddenOut+repmat(B2,1,TestSamNum); %输出层输出最终结果a=postmnmx(NetworkOut,mint,maxt); %还原网络输出层的结果%用原始的数据仿真的结果与已知数据进行对比测试x=1990:2009;newk=a(1,:); %网络输出客运量newh=a(2,:); %网络输出货运量figure;subplot(2,1,1);plot(x,newk,r-o,x,glkyl,b-+);legend(网络输出客运量,实际客运量);xlabel(年份);ylabel(客

56、运量/万人);title(源程序神经网络客运量学习和测试对比图);title(源程序神经网络客运量学习和测试对比图);subplot(2,1,2);plot(x,newh,r-o,x,glhyl,b-+);legend(网络输出货运量,实际货运量);xlabel(年份); ylabel(货运量/万吨);pnew=73.39 75.55 3.9635 4.0975 0.9889 1.0268; %2010年和2011 年的相关数据pnewn=tramnmx(pnew,minp,maxp); %利用原始输入数据的归一化参数对新数据进行归一化HiddenOut=logsig(W1*pnewn+rep

57、mat(B1,1,ForcastSamNum); %隐含层输出预测结果anewn=W2*HiddenOut+repmat(B2,1,ForcastSamNum); %输出层输出预测结果%把网络预测得到的数据还原为原始的数量级anew=postmnmx(anewn,mint,maxt);小波分析小波分析 小波分析是纯数学、应用数学和工程技术的完美结合。从数学来说是大半个世纪“调和分析”的结晶（包括傅里叶分析、函数空间等）。 小波变换是20世纪最辉煌科学成就之一。在计算机应用、信号处理、图象分析、非线性科学、地球科学和应用技术等已有重大突破，预示着小波分析进一步热潮的到来。 “小波分析” 是分析原

58、始信号各种变化的特性，进一步用于数据压缩、噪声去除、特征选择等。例如歌唱信号：是高音还是低音，发声时间长短、起伏、旋律等。从平稳的波形发现突变的尖峰。小波分析是利用多种 “小波基函数” 对 “原始信号” 进行分解、再重构。傅立叶变换缺点傅立叶变换缺点函数f(t)的连续傅立叶变换定义为函数f(t)的逆傅立叶变换定义为（1）用傅立叶变换提取信号的频谱需要利用信号的全部时域信息。（2）傅立叶变换没有反映出随着时间的变化信号频率成分的变化情况。（3）傅立叶变换的积分作用平滑了非平稳信号的突变成分。t(t)(-tj -defFdeFf-tj)(21(t)反映傅立叶变换缺点的一个例子：解决傅立叶变换缺点的

59、方法：小波的定义小波的定义 小波（Wavelet）是由“wave”和“let”构成，表示的是一种长度有限、平均值为0的波形。 确切定义：设 为一平方可积函数，即 若其傅立叶变换 满足条件： 则称为一个基本小波、母小波或者容许小波。)(t)()(2RLt )(dCR2)(小波有两个基本特点：小波有两个基本特点：(1)小。在 空间内，常选取紧支集或近似紧支集（具有时域局限性）且具有正则性（具有频域的局限性）的实数或复数作为小波母函数。这样的小波母函数在时频域都会有较好的局部特性，即快速收敛性。(2)正负交替的“波动性”，也即直流分量为零，即震荡性。)(2RL母小波伸缩平移母小波伸缩平移将母小波函数

60、 进行伸缩和平移，就可以得到子函数其中 称为依赖于参数a,b 的小波家族子函数。a 为伸缩因子，用来控制家族子小波 的体形（“肥胖”和“高矮”）；b 为平移因子，用来控制小波子函数的中心位置。 称为连续小波函数基两点说明：（1）相互正交；（2）具有相等的能量。)(t)(,tba0,),(1)(,aRbaabtatba)(,tba)(,tbaHaar小波基母函数小波基母函数（a）Haar “近似”基函数 （b）Haar “细节”基函数 低频滤波系数 高频滤波系数 H0= 1 1 q H1= 1 -1 q = q q = q -q 其中：7071.02 qRbafdtabttfattfbaWT)(