在科学计算中的应用9非传统解法选讲.ppt

《在科学计算中的应用9非传统解法选讲.ppt》由会员分享，可在线阅读，更多相关《在科学计算中的应用9非传统解法选讲.ppt（53页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、在科学计算中的应用9非传统解法选讲 Still waters run deep.流静水深流静水深,人静心深人静心深 Where there is life,there is hope。有生命必有希望。有生命必有希望9.1遗传算法 9.1.1遗传算法及其在最优化问题中的应用遗传算法是基于进化论，在计算机上模拟生命进化机制而发展起来的一门新学科，它根据适者生存、优胜劣汰等自然进化规则搜索和计算问题的解。美国 Michigen 大学的 John Holland 于 1975 年提出的。遗传算法最优化工具箱MATLAB 7.0的遗传算法与直接搜索工具箱遗传算法的基本思想从一个代表最优化问题解的一组初值

2、开始进行搜索，这组解称为一个种群，这里种群由一定数量的、通过基因编码的个体组成，其中每一个个体称为染色体，不同个体通过染色体的复制、交叉或变异又生成新的个体，依照适者生存的规则，个体也在一代一代进化，通过若干代的进化最终得出条件最优的个体。简单遗传算法的一般步骤选择 n 个个体构成初始种群 ，并求出种群内各个个体的函数值。设置代数为 i=1，即设置其为第一代。计算选择函数的值，所谓选择即通过概率的形式从种群中选择若干个个体的方式。通过染色体个体基因的复制、交叉、变异等创造新的个体，构成新的种群 。i=i+1，若终止条件不满足，则继续进化。遗传算法和传统优化算法比较不同于从一个点开始搜索最优解的

3、传统的最优化算法，遗传算法从一个种群开始对问题的最优解进行并行搜索，所以更利于全局最优化解的搜索。遗传算法并不依赖于导数信息或其他辅助信息来进行最优解搜索。遗传算法采用的是概率型规则而不是确定性规则，所以每次得出的结果不一定完全相同，有时甚至会有较大的差异。9.1.2 遗传算法在求解最优化问题中的应用举例GAOT 工具箱（目标求最大）bound=xm,xM为求解上下界构成的矩阵。a由最优解与目标构成，b为搜索的最终种群，c中间过程参数表。MATLAB 7.0GA工具箱界面，gatool()例：绘制目标函数曲线：ezplot(x*sin(10*pi*x)+2,-1,2)测试不同的初值：f=inl

4、ine(-x.*sin(10*pi*x)-2,x);v=;for x0=-1:0.8:1.5,1.5:0.1:2 x1=fmincon(f,x0,-1,2);v=v;x0,x1,f(x1);end vv=-1.0000 -1.0000 -2.0000 -0.2000 -0.6516 -2.6508 0.6000 0.6516 -2.6508 1.4000 1.4507 -3.4503 1.5000 0.2540 -2.2520 1.6000 1.6506 -3.6503 1.7000 1.2508 -3.2504 1.8000 1.8505 -3.8503 1.9000 0.4522 -2.4

5、511 2.0000 2.0000 -2.0000编写函数：function sol,y=c10mga1(sol,options)x=sol(1);y=x.*sin(10*pi*x)+2;调用gaopt()函数 a,b,c,d=gaopt(-1,2,c10mga1);a,ca=1.85054746606888 3.85027376676810c=1.0e+002*0.01000000000000 0.01644961385548 0.03624395818177 0.02000000000000 0.01652497353988 0.03647414028140 0.1600000000000

6、0 0.01850468596975 0.03850268083951 0.23000000000000 0.01850553961009 0.03850273728228 1.00000000000000 0.01850547466069 0.03850273766768比较：ff=optimset;ff.Display=iter;x0=1.8;x1=fmincon(f,x0,-1,2,ff);f(x1)ans=-3.85027376676808 f(a(1)%遗传算法结果ans=-3.85027376676810 ezplot(x*sin(10*pi*x)+2,-1,20)改变求解区间 a

7、,b,c,d=gaopt(-1,20,c10mga1);a,ca=19.45005206632863 21.45002604650601c=1.0e+002*0.01000000000000 0.17243264358456 0.18858649532480 0.02000000000000 0.19253552639304 0.21133759487918 0.25000000000000 0.19450021530572 0.21450017081177 0.27000000000000 0.19450024961756 0.21450018981219 0.29000000000000

8、0.19450055493368 0.21450025935531 1.00000000000000 0.19450052066329 0.21450026046506 ezplot(x*sin(10*pi*x)+2,12,20)放大区间 a,b,c,d=gaopt(12,20,c10mga1);a,ca=19.85005104334383 21.85002552164857c=1.0e+002*0.01000000000000 0.17647930304626 0.19610637643594 0.03000000000000 0.17648091337382 0.1961637407469

9、7 0.05000000000000 0.18841858256128 0.202288599115410.21000000000000 0.19850064250944 0.21850023812862 0.23000000000000 0.19850055906254 0.21850025289993 1.00000000000000 0.19850051043344 0.21850025521649例：求最小值编写函数：function sol,f=c10mga3(sol,options)x=sol(1:4);f=-(x(1)+x(2)2-5*(x(3)-x(4)2-(x(2)-2*x(

10、3)4-10*(x(1)-x(4)4;a,b,c,d=gaopt(-1,1;-1 1;-1 1;-1 1,c10mga3);a,ca=-0.0666 0.0681 -0.0148 -0.0154 -0.0002c=1.0000 -0.3061 0.2075 -0.2235 -0.1206 -0.25805.0000 -0.2294 0.2076 0.0352 -0.1217 -0.125393.0000 -0.0666 0.0682 -0.0148 -0.0154 -0.0002100.0000 -0.0666 0.0681 -0.0148 -0.0154 -0.0002 求解区域太小，有误差

11、GAOT 的最优化函数 其中：p可给目标函数增加附加参数,v为精度及显示控制向量,P0为初始种群,fun1为终止函数的名称,默认值maxGenTerm,n为最大的允许代数。例：求最小值 tic,xmM=-ones(4,1),ones(4,1)*1000;a,b,c,d=gaopt(xmM,c10mga3,maxGenTerm,2000);a(1:4),dd=c(1:100:end,:);c(end,:),tocans=-0.0049 0.0049 -0.0081 -0.0081dd=1.0e+009*0.0000 0.0000 -0.0000 -0.0000 0.0000 -5.9663 0.

12、0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000elapsed_time=76.5200描述函数：matlab7.0function f=c10mga3a(x)f=(x(1)+x(2)2+5*(x(3)-x(4)2+(x(2)-2*x(3)4+10*(x(1)-x(4)4;x,f=ga(c10mga3a,4)%四个自变量Optimization terminated:maximum number of generations exceeded.x=0.06976151754582 -0.05491931584170 0.04952579333589 0.0

13、6130810339402f=0.00147647985822 ff=gaoptimset;ff.Generations=2000;ff.PopulationSize=80;ff.CrossoverFcn=crossoverheuristic;x=ga(c10mga3a,4,ff)Optimization terminated:maximum number of generations exceeded.x=-0.00216363106525 0.00216366042770 -0.00039985387788 -0.00039996677375f=1.739330597649231e-010

14、 f=inline.%目标函数描述(x(1)+x(2)2+5*(x(3)-x(4)2+(x(2)-2*x(3)4+10*(x(1)-x(4)4,x);时间少，精度高 ff=optimset;ff.MaxIter=10000;ff.TolX=1e-7;tic,x,f1=fminsearch(f,10*ones(4,1),ff);toc;x,f1Elapsed time is 0.595406 seconds.ans=1.0e-006*0.03039572499758 -0.03039585246164 -0.75343487601326 -0.75343518285272f1=9.0140528

15、14563438e-024例：求下面的最优化问题 x,y=meshgrid(-1:0.1:3,-3:0.1:3);z=sin(3*x.*y+2)+x.*y+x+y;surf(x,y,z);shading interp%用光滑曲面表示目标函数函数描述：传统方法 function y=c10mga5(x)y=sin(3*x(1)*x(2)+2)+x(1)*x(2)+x(1)+x(2);x0=1,3;x=fmincon(c10mga5,x0,-1;-3,3;3)x=-1.00000000000000 1.19031291227215函数描述：function sol,y=c10mga6(sol,op

16、tions)x=sol(1:2);y=-sin(3*x(1)*x(2)+2)-x(1)*x(2)-x(1)-x(2);xmM=-1 3;-3 3;a,b,c,d=gaopt(xmM,c10mga6,maxGenTerm,500);aa=2.51604948433614 -3.00000000000000 9.00709500762913遗传算法优化中间结果（40代即可，无需500代，可用默认100）9.1.3 遗传算法在有约束最优化问题中的应用不能直接用于有约束最优化问题求解需通过变换处理划为无约束最优化问题对等式约束可通过等式求解将若干个自变量用其它自变量表示。不等式约束可用惩罚函数方法转移

17、到目标函数中。仍采用 gaopt()或 ga()函数求解例：描述函数：function sol,y=c10mga4(sol,options)x=sol(1:2);x=x(:);x(3)=(6+4*x(1)-2*x(2)/3;y1=-2 1 1*x;y2=-1 1 0*x;if(y19|y2-4|x(3)a,b,c=gaopt(-1000 0;-1000 0,c10mga4,maxGenTerm,1000);c=c(1:15:end,:);c(end,:);a,ca=-6.99981015633155-10.99962347934527 28.99905078165773c=1.0e+003*0

18、.00100000000000 -0.32769544124065 -0.20423049398177 -0.10000000000000 0.05900000000000 -0.00146223175991 0 0.00131115879955 0.10200000000000 -0.00416116639726 -0.00666729713459 0.01480583198631 0.84900000000000 -0.00689401645967 -0.01080365682806 0.02847008229837 0.89200000000000 -0.00694511749224 -

19、0.01089232545085 0.02872558746118 0.93200000000000 -0.00698531391213 -0.01097813084259 0.02892656956064 0.96800000000000 -0.00699692906988 -0.01099399300138 0.02898464534940 1.00000000000000 -0.00699981015633 -0.01099962347935 0.02899905078166可用线性规划得出更精确的结果 f=1 2 3;A=-2 1 1;1-1 0;B=9;4;Aeq=4-2-3;Beq

20、=-6;x=linprog(f,A,B,Aeq,Beq,-inf;-inf;0,0;0;inf);xOptimization terminated successfully.ans=-6.99999999999967-10.99999999999935 0.00000000000000 f*xans=-28.99999999999836 建议求解方法：用 GA 找出全局最优解的大致位置,以其为初值调用最优化函数求精确解。9.2神经网络及其在数据拟合中的应用9.2.1神经网络基础知识单个人工神经元的数学表示形式例:常用传输函数曲线 x=-2:0.01:2;y=tansig(x);plot(x,y

21、)x=-2:0.01:2;y=logsig(x);plot(x,y)BP 神经网络结构示意图 其中:xm,xM分别为列向量,为各样本数据的最大最小值。hi为一行向量，各隐层节点数。fi每层传输函数，同一层应使用相同的传输函数。例：考虑一个前馈网络，2个隐层，第一个有8个节点，采用Sigmoid传输函数，第二层节点个数应该等于输出信号的路数，故节点数为1，传输函数为对数Sigmoid函数。net=newff(0,1;-1,5,8,1,tansig,logsig);3个隐层，1层4个点，线性函数；2层6个点，Sigmoid函数；3层1个点，logsig函数。net=newff(0,1;-1,5,4

22、 6 1,purelin,tansig,logsig);可用下面的语句格式设定其它参数。net.trainParam.epochs=300;net.trainFcn=trainlm;9.2.2 神经网络的训练与泛化神经网络训练函数 X为n*M，n为输入变量的路数，M为样本的组数，Y为m*M，m为输出变量的路数。tr为结构体数据，返回训练的相关跟踪信息。Y1和E为输出和误差矩阵。可多次训练，原加权矩阵为初值。目标值曲线函数神经网络泛化例:由前面最小拟合的例子中的数据进行曲线拟合，2个隐层，隐层节点选择为5。x=0:.5:10;y=0.12*exp(-0.213*x)+0.54*exp(-0.17

23、*x).*sin(1.23*x);x0=0:0.1:10;y0=0.12*exp(-0.213*x0)+0.54*exp(-0.17*x0).*sin(1.23*x0);net=newff(0,10,5,1,tansig,tansig);net.trainParam.epochs=1000;%设置最大步数 net=train(net,x,y);%训练神经网络 net.IW1 net.LW2,1%隐层权值和输出层权值ans=0.4765 -1.9076 0.5784 0.9450 -0.2888 -2.7916 0.3052 -2.9388 0.9780 1.1814可改变求解算法 net=ne

24、wff(0,10,5,1,tansig,tansig);net.trainParam.epochs=100;net.trainFcn=trainlm;net,b1=train(net,x,y);net=newff(0,10,5,1,tansig,tansig);net.trainParam.epochs=100;net.trainFcn=traincgf;net,b2=train(net,x,y);net=newff(0,10,5,1,tansig,tansig);net.trainParam.epochs=100;net.trainFcn=traingdx;net,b3=train(net,

25、x,y);可改变各层传输函数 net=newff(0,10,5,1,tansig,logsig);net.trainParam.epochs=100;net.trainFcn=trainlm;net,b2=train(net,x,y);net=newff(0,10,5,1,logsig,tansig);net,b3=train(net,x,y);net=newff(0,10,5,1,logsig,logsig);net,b4=train(net,x,y);可改变结构，选择隐层15个节点 net=newff(0,10,15,1,tansig,tansig);net.trainParam.epoc

26、hs=100;net.trainFcn=trainlm;net,b2=train(net,x,y);figure;y1=sim(net,x0);plot(x0,y0,x0,y1,x,y,o)例：二元函数的拟合 x,y=meshgrid(-3:.6:3,-2:.4:2);x=x(:);y=y(:);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);%这三个变量均应为行向量 net=newff(-3 3;-2 2,10,10,1,tansig,tansig,tansig);net.trainParam.epochs=1000;net.trainFcn=trainlm;net,b=tr

27、ain(net,x;y,z);%训练神经网络 x2,y2=meshgrid(-3:.1:3,-2:.1:2);x1=x2(:);y1=y2(:);figure;z1=sim(net,x1;y1);z2=reshape(z1,size(x2);surf(x2,y2,z2)改变第二层节点数 net=newff(-3 3;-2 2,10,20,1,tansig,tansig,tansig);net,b=train(net,x;y,z);%训练神经网络 z1=sim(net,x1;y1);z2=reshape(z1,size(x2);surf(x2,y2,z2)效果不好给出密集一点的的样本点 x,y=

28、meshgrid(-3:.2:3,-2:.2:2);x=x(:);y=y(:);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);net=newff(-3 3;-2 2,10,10,1,tansig,tansig,tansig);net.trainParam.epochs=100;net.trainFcn=trainlm;net=train(net,x;y,z);x1,y1=meshgrid(-3:.1:3,-2:.1:2);a=x1;x1=x2(:);y1=y2(:);z1=sim(net,x1;y1);z2=reshape(z1,size(a);surf(x2,y2,z2)n

29、et=newff(-3 3;-2 2,10,20,1,tansig,tansig,tansig);net=train(net,x;y,z);%修改节点个数后的泛化效果 figure;z1=sim(net,x1;y1);z2=reshape(z1,size(a);surf(x2,y2,z2)9.2.3 神经网络界面启动神经网络界面 nntool 遗传算法的MATLAB编程包括如下几个程序文件：genetic.m (主程序文件）gen_encode.m(二进制数组编码P）gen_decode.m(将二进制数组P解码为状态矩阵)crossover.m(两个染色体间的交叉)mutation.m (变异

30、)shuffle.m (打乱染色体次序)function xo,fo=genetic(f,x0,l,u,Np,Nb,Pc,Pm,eta,kmax)%基因算法求f(x）最小值 s.t.l=x=u%f为待求函数,x0初值，l，u上下限，Np群体大小，Nb每一个变量的基因值（二进制数）%Pc交叉概率，Pm变异概率，eta学习率，kmax最大迭代次数N=length(x0);%确定各变量缺省值if nargin 10 kmax=100;%最大迭代次数缺省为100endif nargin 1|eta=0 eta=1;%学习率eta,(0 eta 1)endif nargin 8 Pm=0.01;%变异概

31、率缺省0.01endif nargin 7 Pc=0.5;%交叉概率缺省0.5endif nargin 6 Nb=8*ones(1,N);%每一变量的基因值（二进制数）endif nargin 5 Np=10;%群体大小（染色体数）end%生成初始群体NNb=sum(Nb);xo=x0(:);l=l(:);u=u(:);fo=feval(f,xo);X(1,:)=xo;for n=2:Np X(n,:)=l+rand(size(x0).*(u-l);%初始群体随机数组endP=gen_encode(X,Nb,l,u);%编码为2进制字串for k=1:kmax X=gen_decode(P,N

32、b,l,u);%解码为10进制数 for n=1:Np fX(n)=feval(f,X(n,:);end fxb,nb=min(fX);%选择最适合的，函数值最小的 if fxb fo fo=fxb;xo=X(nb,:);end fX1=max(fxb)-fX;%将函数值转化为非负的适合度值 fXm=fX1(nb);if fXm eps%如果所有的染色体值相同，终止程序 return;end%复制下一代 for n=1:Np X(n,:)=X(n,:)+eta*(fXm-fX1(n)/fXm*(X(nb,:)-X(n,:);%复制准则 end P=gen_encode(X,Nb,l,u);%对

33、下一代染色体编码%随机配对/交叉得新的染色体数组 is=shuffle(1:Np);for n=1:2:Np-1 if rand Pc P(is(n:n+1),:)=crossover(P(is(n:n+1),:),Nb);end end%变异 P=mutation(P,Nb,Pm);endfunction P=gen_encode(X,Nb,l,u)%将群体X的状态编码为二进制数组PNp=size(X,1);%群体大小N=length(Nb);%变量（状态）维数for n=1:Np b2=0;for m=1:N b1=b2+1;b2=b2+Nb(m);Xnm=(2Nb(m)-1)*(X(n,

34、m)-l(m)/(u(m)-l(m);%编码方程 P(n,b1:b2)=dec2bin(Xnm,Nb(m);%10进制转换为2进制 endendfunction X=gen_decode(P,Nb,l,u)%将二进制数组P解码为群体X的状态矩阵Np=size(P,1);%群体大小N=length(Nb);%变量维数for n=1:Np b2=0;for m=1:N b1=b2+1;b2=b1+Nb(m)-1;X(n,m)=bin2dec(P(n,b1:b2)*(u(m)-l(m)/(2Nb(m)-1)+l(m);%解码方程 endendfunction chrms2=crossover(chr

35、ms2,Nb)%两个染色体间的交叉Nbb=length(Nb);b2=0;for m=1:Nbb b1=b2+1;bi=b1+mod(floor(rand*Nb(m),Nb(m);b2=b2+Nb(m);tmp=chrms2(1,bi:b2);chrms2(1,bi:b2)=chrms2(2,bi:b2);chrms2(2,bi:b2)=tmp;endfunction P=mutation(P,Nb,Pm)%变异Nbb=length(Nb);for n=1:size(P,1)b2=0;for m=1:Nbb if rand Pm b1=b2+1;bi=b1+mod(floor(rand*Nb(

36、m),Nb(m);b2=b2+Nb(m);P(n,bi)=P(n,bi);end endendfunction is=shuffle(is)%打乱染色体次序N=length(is);for n=N:-1:2 in=ceil(rand*(n-1);tmp=is(in);is(in)=is(n);is(n)=tmp;%将第n个元素与第in个元素交换end%例:ex.m 遗传算法计算最优解f=inline(x(1)4-16*x(1)2-5*x(1)*x(2)+x(2)4-16*x(2)2-5*x(2),x);l=-5-5;%下限u=5 5;%上限x0=0 0;Np=30;%群体大小Nb=12 12;%代表每个变量的二进制位数Pc=0.5;%交叉概率Pm=0.01;%变异概率eta=0.8;%学习率kmax=200;%最大迭代次数xos,fos=fminsearch(f,x0)xo_gen,fo_gen=genetic(f,x0,l,u,Np,Nb,Pc,Pm,eta,kmax)