matlab在科學(xué)計(jì)算中的應(yīng)用9數(shù)學(xué)問(wèn)題的非傳統(tǒng)解法選講課件

1、9.1遺傳算法 9.1.1遺傳算法及其在最優(yōu)化問(wèn)題中的應(yīng)用遺傳算法是基于進(jìn)化論，在計(jì)算機(jī)上模擬生命進(jìn)化機(jī)制而發(fā)展起來(lái)的一門新學(xué)科，它根據(jù)適者生存、優(yōu)勝劣汰等自然進(jìn)化規(guī)則搜索和計(jì)算問(wèn)題的解。美國(guó) Michigen 大學(xué)的 John Holland 于 1975 年提出的。遺傳算法最優(yōu)化工具箱MATLAB 7.0的遺傳算法與直接搜索工具箱第1頁(yè)，共52頁(yè)。遺傳算法的基本思想從一個(gè)代表最優(yōu)化問(wèn)題解的一組初值開(kāi)始進(jìn)行搜索，這組解稱為一個(gè)種群，這里種群由一定數(shù)量的、通過(guò)基因編碼的個(gè)體組成，其中每一個(gè)個(gè)體稱為染色體，不同個(gè)體通過(guò)染色體的復(fù)制、交叉或變異又生成新的個(gè)體，依照適者生存的規(guī)則，個(gè)體也在一代一代進(jìn)

2、化，通過(guò)若干代的進(jìn)化最終得出條件最優(yōu)的個(gè)體。第2頁(yè)，共52頁(yè)。簡(jiǎn)單遺傳算法的一般步驟選擇 n 個(gè)個(gè)體構(gòu)成初始種群 ，并求出種群內(nèi)各個(gè)個(gè)體的函數(shù)值。設(shè)置代數(shù)為 i=1，即設(shè)置其為第一代。計(jì)算選擇函數(shù)的值，所謂選擇即通過(guò)概率的形式從種群中選擇若干個(gè)個(gè)體的方式。通過(guò)染色體個(gè)體基因的復(fù)制、交叉、變異等創(chuàng)造新的個(gè)體，構(gòu)成新的種群 。i=i+1，若終止條件不滿足，則繼續(xù)進(jìn)化。第3頁(yè)，共52頁(yè)。遺傳算法和傳統(tǒng)優(yōu)化算法比較不同于從一個(gè)點(diǎn)開(kāi)始搜索最優(yōu)解的傳統(tǒng)的最優(yōu)化算法，遺傳算法從一個(gè)種群開(kāi)始對(duì)問(wèn)題的最優(yōu)解進(jìn)行并行搜索，所以更利于全局最優(yōu)化解的搜索。遺傳算法并不依賴于導(dǎo)數(shù)信息或其他輔助信息來(lái)進(jìn)行最優(yōu)解搜索。遺傳

3、算法采用的是概率型規(guī)則而不是確定性規(guī)則，所以每次得出的結(jié)果不一定完全相同，有時(shí)甚至?xí)休^大的差異。第4頁(yè)，共52頁(yè)。9.1.2 遺傳算法在求解最優(yōu)化問(wèn)題中的應(yīng)用舉例GAOT 工具箱（目標(biāo)求最大） bound=xm,xM為求解上下界構(gòu)成的矩陣。a由最優(yōu)解與目標(biāo)構(gòu)成，b為搜索的最終種群，c中間過(guò)程參數(shù)表。MATLAB 7.0GA工具箱界面， gatool()第5頁(yè)，共52頁(yè)。例：繪制目標(biāo)函數(shù)曲線： ezplot(x*sin(10*pi*x)+2,-1,2)第6頁(yè)，共52頁(yè)。測(cè)試不同的初值： f=inline(-x.*sin(10*pi*x)-2,x); v=; for x0=-1:0.8:1.5,

4、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.4511 2.0000 2.0000 -2.0000第7頁(yè)，共52頁(yè)。編寫函數(shù)：func

5、tion sol,y=c10mga1(sol,options) x=sol(1); y=x.*sin(10*pi*x)+2;調(diào)用gaopt( )函數(shù) 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.16000000000000 0.01850468596975 0.038

6、50268083951 0.23000000000000 0.01850553961009 0.03850273728228 1.00000000000000 0.01850547466069 0.03850273766768第8頁(yè)，共52頁(yè)。比較： ff=optimset; ff.Display=iter; x0=1.8; x1=fmincon(f,x0,-1,2,ff); f(x1)ans = -3.85027376676808 f(a(1) % 遺傳算法結(jié)果ans = -3.85027376676810 ezplot(x*sin(10*pi*x)+2,-1,20) 改變求解區(qū)間 a,b,

7、c,d=gaopt(-1,20,c10mga1); a,ca = 19.45005206632863 21.45002604650601第9頁(yè)，共52頁(yè)。c = 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.29

8、000000000000 0.19450055493368 0.21450025935531 1.00000000000000 0.19450052066329 0.21450026046506第10頁(yè)，共52頁(yè)。 ezplot(x*sin(10*pi*x)+2,12,20) 放大區(qū)間 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

9、.17648091337382 0.19616374074697 0.05000000000000 0.18841858256128 0.202288599115410.21000000000000 0.19850064250944 0.21850023812862 0.23000000000000 0.19850055906254 0.21850025289993 1.00000000000000 0.19850051043344 0.21850025521649第11頁(yè)，共52頁(yè)。例：求最小值編寫函數(shù)：function sol,f=c10mga3(sol,options) x=sol(1:

10、4); f=-(x(1)+x(2)2-5*(x(3)-x(4)2-(x(2)-2*x(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.000

11、0 -0.0666 0.0681 -0.0148 -0.0154 -0.0002 求解區(qū)域太小，有誤差第12頁(yè)，共52頁(yè)。第13頁(yè)，共52頁(yè)。GAOT 的最優(yōu)化函數(shù) 其中：p可給目標(biāo)函數(shù)增加附加參數(shù), v為精度及顯示控制向量, P0為初始種群, fun1為終止函數(shù)的名稱,默認(rèn)值maxGenTerm, n為最大的允許代數(shù)。第14頁(yè)，共52頁(yè)。例：求最小值 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,:), tocan

12、s = -0.0049 0.0049 -0.0081 -0.0081dd = 1.0e+009 * 0.0000 0.0000 -0.0000 -0.0000 0.0000 -5.9663 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000elapsed_time = 76.5200第15頁(yè)，共52頁(yè)。第16頁(yè)，共52頁(yè)。描述函數(shù)：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) %四個(gè)自

13、變量Optimization terminated: maximum number of generations exceeded.x = 0.06976151754582 -0.05491931584170 0.04952579333589 0.06130810339402f = 0.00147647985822 ff=gaoptimset; ff.Generations=2000; ff.PopulationSize=80;ff.CrossoverFcn=crossoverheuristic; x=ga(c10mga3a,4,ff)Optimization terminated: maxi

14、mum number of generations exceeded.x = -0.00216363106525 0.00216366042770 -0.00039985387788 -0.00039996677375f = 1.739330597649231e-010第17頁(yè)，共52頁(yè)。 f=inline. % 目標(biāo)函數(shù)描述(x(1)+x(2)2+5*(x(3)-x(4)2+(x(2)-2*x(3)4+10*(x(1)-x(4)4,x); 時(shí)間少，精度高 ff=optimset; ff.MaxIter=10000; ff.TolX=1e-7; tic, x,f1=fminsearch(f,1

15、0*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.014052814563438e-024第18頁(yè)，共52頁(yè)。例：求下面的最優(yōu)化問(wèn)題 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 % 用光滑曲面表示目標(biāo)函數(shù)第19頁(yè)，

16、共52頁(yè)。函數(shù)描述：傳統(tǒng)方法 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函數(shù)描述：function sol,y=c10mga6(sol,options) 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,maxGen

17、Term,500); aa = 2.51604948433614 -3.00000000000000 9.00709500762913第20頁(yè)，共52頁(yè)。遺傳算法優(yōu)化中間結(jié)果（40代即可，無(wú)需500代，可用默認(rèn)100）第21頁(yè)，共52頁(yè)。9.1.3 遺傳算法在有約束最優(yōu)化問(wèn)題中的應(yīng)用不能直接用于有約束最優(yōu)化問(wèn)題求解需通過(guò)變換處理劃為無(wú)約束最優(yōu)化問(wèn)題對(duì)等式約束可通過(guò)等式求解將若干個(gè)自變量用其它自變量表示。不等式約束可用懲罰函數(shù)方法轉(zhuǎn)移到目標(biāo)函數(shù)中。仍采用 gaopt() 或 ga() 函數(shù)求解第22頁(yè)，共52頁(yè)。例：描述函數(shù)：function sol,y=c10mga4(sol,options)

18、 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.00100000000000 -0.32769544124065 -0.2042304939

19、8177 -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 -0.01089232545085 0.02872558746118 0.93200000000

20、000 -0.00698531391213 -0.01097813084259 0.02892656956064 0.96800000000000 -0.00699692906988 -0.01099399300138 0.02898464534940 1.00000000000000 -0.00699981015633 -0.01099962347935 0.02899905078166第24頁(yè)，共52頁(yè)。可用線性規(guī)劃得出更精確的結(jié)果 f=1 2 3; A=-2 1 1; 1 -1 0; B=9; 4; Aeq=4 -2 -3; Beq=-6; x=linprog(f,A,B,Aeq,Beq

21、,-inf;-inf;0,0;0;inf); xOptimization terminated successfully.ans = -6.99999999999967 -10.99999999999935 0.00000000000000 f*xans = -28.99999999999836 建議求解方法：用 GA 找出全局最優(yōu)解的大致位置,以其為初值調(diào)用最優(yōu)化函數(shù)求精確解。第25頁(yè)，共52頁(yè)。9.2神經(jīng)網(wǎng)絡(luò)及其在數(shù)據(jù)擬合中的應(yīng)用9.2.1神經(jīng)網(wǎng)絡(luò)基礎(chǔ)知識(shí)單個(gè)人工神經(jīng)元的數(shù)學(xué)表示形式第26頁(yè)，共52頁(yè)。例:常用傳輸函數(shù)曲線 x=-2:0.01:2; y=tansig(x); plot(x,

22、y) x=-2:0.01:2; y=logsig(x); plot(x,y)第27頁(yè)，共52頁(yè)。BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖第28頁(yè)，共52頁(yè)。 其中:xm,xM分別為列向量,為各樣本數(shù)據(jù)的最大最小值。 hi為一行向量，各隱層節(jié)點(diǎn)數(shù)。fi每層傳輸函數(shù)，同一層應(yīng)使用相同的傳輸函數(shù)。第29頁(yè)，共52頁(yè)。例： 考慮一個(gè)前饋網(wǎng)絡(luò)，2個(gè)隱層，第一個(gè)有8個(gè)節(jié)點(diǎn)，采用Sigmoid傳輸函數(shù)，第二層節(jié)點(diǎn)個(gè)數(shù)應(yīng)該等于輸出信號(hào)的路數(shù)，故節(jié)點(diǎn)數(shù)為1，傳輸函數(shù)為對(duì)數(shù)Sigmoid函數(shù)。 net=newff(0,1; -1,5,8,1,tansig,logsig);3個(gè)隱層，1層4個(gè)點(diǎn)，線性函數(shù)；2層6個(gè)點(diǎn)， Sigmoi

23、d函數(shù)；3層1個(gè)點(diǎn)， logsig函數(shù)。 net=newff(0,1; -1,5,4 6 1,purelin,tansig,logsig);可用下面的語(yǔ)句格式設(shè)定其它參數(shù)。 net.trainParam.epochs=300; net.trainFcn=trainlm;第30頁(yè)，共52頁(yè)。9.2.2 神經(jīng)網(wǎng)絡(luò)的訓(xùn)練與泛化神經(jīng)網(wǎng)絡(luò)訓(xùn)練函數(shù) X為n*M，n為輸入變量的路數(shù)，M為樣本的組數(shù)，Y為m*M，m為輸出變量的路數(shù)。tr為結(jié)構(gòu)體數(shù)據(jù)，返回訓(xùn)練的相關(guān)跟蹤信息。Y1和E為輸出和誤差矩陣。 可多次訓(xùn)練，原加權(quán)矩陣為初值。目標(biāo)值曲線函數(shù)神經(jīng)網(wǎng)絡(luò)泛化第31頁(yè)，共52頁(yè)。例:由前面最小擬合的例子中的數(shù)據(jù)進(jìn)

24、行曲線擬合，2個(gè)隱層，隱層節(jié)點(diǎn)選擇為5。 x=0:.5:10; y=0.12*exp(-0.213*x)+0.54*exp(-0.17*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; % 設(shè)置最大步數(shù) net=train(net,x,y); % 訓(xùn)練神經(jīng)網(wǎng)絡(luò)第32頁(yè)，共52頁(yè)。 net.IW1 net.LW2,1 % 隱層權(quán)值和輸出層權(quán)值ans

25、= 0.4765 -1.9076 0.5784 0.9450 -0.2888 -2.7916 0.3052 -2.9388 0.9780 1.1814可改變求解算法 net=newff(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); ne

26、t=newff(0,10,5,1,tansig,tansig); net.trainParam.epochs=100; net.trainFcn=traingdx; net,b3=train(net,x,y);第33頁(yè)，共52頁(yè)。可改變各層傳輸函數(shù) 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(

27、0,10,5,1,logsig,logsig); net,b4=train(net,x,y);可改變結(jié)構(gòu)，選擇隱層15個(gè)節(jié)點(diǎn) net=newff(0,10,15,1,tansig,tansig); net.trainParam.epochs=100; net.trainFcn=trainlm; net,b2=train(net,x,y); figure; y1=sim(net,x0); plot(x0,y0,x0,y1,x,y,o)第34頁(yè)，共52頁(yè)。例：二元函數(shù)的擬合 x,y=meshgrid(-3:.6:3, -2:.4:2); x=x(:); y=y(:); z=(x.2-2*x).*e

28、xp(-x.2-y.2-x.*y); % 這三個(gè)變量均應(yīng)為行向量 net=newff(-3 3; -2 2,10,10,1,tansig,tansig,tansig); net.trainParam.epochs=1000; net.trainFcn=trainlm; net,b=train(net,x; y,z); % 訓(xùn)練神經(jīng)網(wǎng)絡(luò) 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)第35頁(yè)，共52頁(yè)。改

29、變第二層節(jié)點(diǎn)數(shù) net=newff(-3 3; -2 2,10,20,1,tansig,tansig,tansig); net,b=train(net,x; y,z); % 訓(xùn)練神經(jīng)網(wǎng)絡(luò) z1=sim(net,x1; y1); z2=reshape(z1,size(x2); surf(x2,y2,z2) 效果不好第36頁(yè)，共52頁(yè)。給出密集一點(diǎn)的的樣本點(diǎn) x,y=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,tans

30、ig,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)第37頁(yè)，共52頁(yè)。 net=newff(-3 3; -2 2,10,20,1, tansig,tansig,tansig); net=train(net,x; y,z); % 修改節(jié)點(diǎn)個(gè)數(shù)后的泛化

31、效果 figure; z1=sim(net,x1; y1); z2=reshape(z1,size(a); surf(x2,y2,z2)第38頁(yè)，共52頁(yè)。9.2.3 神經(jīng)網(wǎng)絡(luò)界面啟動(dòng)神經(jīng)網(wǎng)絡(luò)界面 nntool第39頁(yè)，共52頁(yè)。 遺傳算法的MATLAB編程包括如下幾個(gè)程序文件： genetic.m (主程序文件） gen_encode.m (二進(jìn)制數(shù)組編碼P） gen_decode.m ( 將二進(jìn)制數(shù)組P解碼為狀態(tài)矩陣) crossover.m (兩個(gè)染色體間的交叉) mutation.m (變異) shuffle.m (打亂染色體次序)第40頁(yè)，共52頁(yè)。function xo,fo =

32、genetic(f,x0,l,u,Np,Nb,Pc,Pm,eta,kmax)% 基因算法求f(x）最小值 s.t. l = x = u%f為待求函數(shù),x0初值，l，u上下限，Np群體大小，Nb每一個(gè)變量的基因值（二進(jìn)制數(shù)）%Pc交叉概率，Pm變異概率，eta學(xué)習(xí)率，kmax最大迭代次數(shù)N = length(x0);%確定各變量缺省值if nargin 10 kmax = 100; %最大迭代次數(shù)缺省為100end第41頁(yè)，共52頁(yè)。if nargin 1|eta = 0 eta = 1; %學(xué)習(xí)率eta,(0 eta 1)endif nargin 8 Pm = 0.01; %變異概率缺省0.0

33、1endif nargin 7 Pc = 0.5; %交叉概率缺省0.5endif nargin 6 Nb = 8*ones(1,N);%每一變量的基因值（二進(jìn)制數(shù)）endif nargin 5 Np = 10; %群體大?。ㄈ旧w數(shù)）end第42頁(yè)，共52頁(yè)。%生成初始群體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); %初始群體隨機(jī)數(shù)組endP = gen_encode(X,Nb,l,u);

34、%編碼為2進(jìn)制字串第43頁(yè)，共52頁(yè)。for k = 1:kmax X = gen_decode(P,Nb,l,u); %解碼為10進(jìn)制數(shù) for n = 1:Np fX(n) = feval(f,X(n,:); end fxb,nb = min(fX); %選擇最適合的，函數(shù)值最小的 if fxb fo fo = fxb; xo = X(nb,:); end fX1 = max(fxb) - fX; %將函數(shù)值轉(zhuǎn)化為非負(fù)的適合度值 fXm = fX1(nb); if fXm eps %如果所有的染色體值相同，終止程序 return; end 第44頁(yè)，共52頁(yè)。%復(fù)制下一代 for n =

35、1:Np X(n,:) = X(n,:) + eta*(fXm - fX1(n)/fXm*(X(nb,:) - X(n,:); %復(fù)制準(zhǔn)則 end P = gen_encode(X,Nb,l,u); %對(duì)下一代染色體編碼 %隨機(jī)配對(duì)/交叉得新的染色體數(shù)組 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);end第45頁(yè)，共52頁(yè)。function P = gen_encode(X,Nb,l,u)%將群體X的狀態(tài)編碼為二進(jìn)制數(shù)組PNp=size(X,1); %群體大小N = length(Nb); %變量（狀態(tài)）維數(shù)for n = 1:Np b2 = 0; for m = 1:N b1 = b2+1; b2 = b2 + Nb(m); Xnm =(2Nb(m)- 1)*(X(n,m) - l(m)/(u(m) - l(m); %編碼方程 P(n,b1:b2) = dec2bin(Xnm,Nb(m); %10進(jìn)制轉(zhuǎn)換為2進(jìn)制 en