数值最优化方法

算法来源：《数值最优化方法~高立》

算法目的：实现函数的局部最优化寻找，以二元函数为例，展示了最速下降法和牛顿寻优的算法过程

主要Python模块：numpy,sympy

（1）python实现

（2）MATLAB实现

（3）比较

（1）Python实现

A 最速下降法

 

[python] view plain copy
 
  1. # -*- coding: utf-8 -*-  
  2. """ 
  3. Created on Sat Oct 01 15:01:54 2016 
  4. @author: zhangweiguo 
  5. """  
  6. import sympy,numpy  
  7. import math  
  8. import matplotlib.pyplot as pl  
  9. from mpl_toolkits.mplot3d import Axes3D as ax3  
  10. #最速下降法，二维实验  
  11. def SD(x0,G,b,c,N,E):  
  12.     f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c  
  13.     f_d = lambda x: numpy.dot(G, x)+b  
  14.     X=x0;Y=[];Y_d=[];  
  15.     xx=sympy.symarray('xx',(2,1))  
  16.     n = 1  
  17.     ee = f_d(x0)  
  18.     e=(ee[0]**2+ee[1]**2)**0.5  
  19.     Y.append(f(x0)[0,0]);Y_d.append(e)  
  20.     a=sympy.Symbol('a',real=True)  
  21.     print '第%d次迭代：e=%d' % (n, e)  
  22.     while n<N and e>E:  
  23.         n=n+1  
  24.         yy=f(x0-a*f_d(x0))  
  25.         yy_d=sympy.diff(yy[0,0],a,1)  
  26.         a0=sympy.solve(yy_d)  
  27.         x0=x0-a0*f_d(x0)  
  28.         X=numpy.c_[X,x0]  
  29.         Y.append(f(x0)[0,0])  
  30.         ee = f_d(x0)  
  31.         e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5)  
  32.         Y_d.append(e)  
  33.         print '第%d次迭代：e=%s'%(n,e)  
  34.     return X,Y,Y_d  
  35. if __name__=='__main__':  
  36.     G=numpy.array([[21.0,4.0],[4.0,15.0]])  
  37.     b=numpy.array([[2.0],[3.0]])  
  38.     c=10.0  
  39.     f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + c  
  40.     f_d = lambda x: numpy.dot(G, x) + b  
  41.     x0=numpy.array([[-30.0],[100.0]])  
  42.     N=40  
  43.     E=10**(-6)  
  44.     X, Y, Y_d=SD(x0,G,b,c,N,E)  
  45.     X=numpy.array(X)  
  46.     Y=numpy.array(Y)  
  47.     Y_d=numpy.array(Y_d)  
  48.     figure1=pl.figure('trend')  
  49.     n=len(Y)  
  50.     x=numpy.arange(1,n+1)  
  51.     pl.subplot(2,1,1)  
  52.     pl.semilogy(x,Y,'r*')  
  53.     pl.xlabel('n')  
  54.     pl.ylabel('f(x)')  
  55.     pl.subplot(2,1,2)  
  56.     pl.semilogy(x,Y_d,'g*')  
  57.     pl.xlabel('n')  
  58.     pl.ylabel("|f'(x)|")  
  59.     figure2=pl.figure('behave')  
  60.     x=numpy.arange(-110,110,1)  
  61.     y=x  
  62.     [xx,yy]=numpy.meshgrid(x,y)  
  63.     zz=numpy.zeros(xx.shape)  
  64.     n=xx.shape[0]  
  65.     for i in xrange(n):  
  66.         for j in xrange(n):  
  67.             xxx=numpy.array([xx[i,j],yy[i,j]])  
  68.             zz[i,j]=f(xxx.T)  
  69.     ax=ax3(figure2)  
  70.     ax.contour3D(xx,yy,zz)  
  71.     ax.plot3D(X[0,:],X[1,:],Y,'ro--')  
  72.     pl.show()  

 

结果：

第1次迭代：e=1401
第2次迭代：e=285.423850209
第3次迭代：e=36.4480186834
第4次迭代：e=7.42196747556
第5次迭代：e=0.947769462918
第6次迭代：e=0.192995789132
第7次迭代：e=0.0246451518433
第8次迭代：e=0.00501853110317
第9次迭代：e=0.000640855749362
第10次迭代：e=0.000130498466038
第11次迭代：e=1.66643765921e-05
第12次迭代：e=3.39339326369e-06
第13次迭代：e=4.33329103766e-07

 

B 牛顿方法

 

[python] view plain copy
 
  1. # -*- coding: utf-8 -*-  
  2. """ 
  3. Created on Sat Oct 01 15:01:54 2016 
  4. @author: zhangweiguo 
  5. """  
  6. import numpy  
  7. import math  
  8. import matplotlib.pyplot as pl  
  9. from mpl_toolkits.mplot3d import Axes3D as ax3  
  10. #Newton寻优，二维实验  
  11. def NE(x0,y0,N,E):  
  12.     X1=[];X2=[];Y=[];Y_d=[]  
  13.     n = 1  
  14.     ee = g(x0,y0)  
  15.     e=(ee[0,0]**2+ee[1,0]**2)**0.5  
  16.     X1.append(x0)  
  17.     X2.append(y0)  
  18.     Y.append(f(x0,y0))  
  19.     Y_d.append(e)  
  20.     print '第%d次迭代：e=%s' % (n, e)  
  21.     while n<N and e>E:  
  22.         n=n+1  
  23.         d=-numpy.dot(numpy.linalg.inv(G(x0,y0)),g(x0,y0))  
  24.         #d=-numpy.dot(numpy.linalg.pinv(G(x0,y0)),g(x0,y0))  
  25.         x0=x0+d[0,0]  
  26.         y0=y0+d[1,0]  
  27.         ee = g(x0, y0)  
  28.         e = (ee[0,0] ** 2 + ee[1,0] ** 2) ** 0.5  
  29.         X1.append(x0)  
  30.         X2.append(y0)  
  31.         Y.append(f(x0, y0))  
  32.         Y_d.append(e)  
  33.         print '第%d次迭代：e=%s'%(n,e)  
  34.     return X1,X2,Y,Y_d  
  35. if __name__=='__main__':  
  36.     f = lambda x,y: 3*x**2+3*y**2-x**2*y                    #原函数  
  37.     g = lambda x,y: numpy.array([[6*x-2*x*y],[6*y-x**2]])   #一阶导函数向量  
  38.     G = lambda x,y: numpy.array([[6-2*y,-2*x],[-2*x,6]])    #二阶导函数矩阵  
  39.     x0=-2;y0=4  
  40.     N=10;E=10**(-6)  
  41.     X1,X2,Y,Y_d=NE(x0,y0,N,E)  
  42.     X1=numpy.array(X1)  
  43.     X2=numpy.array(X2)  
  44.     Y=numpy.array(Y)  
  45.     Y_d=numpy.array(Y_d)  
  46.     figure1=pl.figure('trend')  
  47.     n=len(Y)  
  48.     x=numpy.arange(1,n+1)  
  49.     pl.subplot(2,1,1)  
  50.     pl.semilogy(x,Y,'r*')  
  51.     pl.xlabel('n')  
  52.     pl.ylabel('f(x)')  
  53.     pl.subplot(2,1,2)  
  54.     pl.semilogy(x,Y_d,'g*')  
  55.     pl.xlabel('n')  
  56.     pl.ylabel("|f'(x)|")  
  57.     figure2=pl.figure('behave')  
  58.     x=numpy.arange(-6,6,0.1)  
  59.     y=x  
  60.     [xx,yy]=numpy.meshgrid(x,y)  
  61.     zz=numpy.zeros(xx.shape)  
  62.     n=xx.shape[0]  
  63.     for i in xrange(n):  
  64.         for j in xrange(n):  
  65.             zz[i,j]=f(xx[i,j],yy[i,j])  
  66.     ax=ax3(figure2)  
  67.     ax.contour3D(xx,yy,zz,15)  
  68.     ax.plot3D(X1,X2,Y,'ro--')  
  69.     pl.show()  


 

 

（2）MATLAB实现

A 最速下降法

先建立调用函数，输入不同的系数矩阵G和不同的初始迭代点会有不同的表现

再调用一下即可

 

[plain] view plain copy
 
  1. function [X Y Y_d]=SD(x0)  
  2. %%%最速下降法  
  3. %%%采用精确搜索  
  4. tic  
  5. G=[21 4;4 15];  
  6. %G=[21 4;4 1];  
  7. b=[2 3]';  
  8. c=10;  
  9. %x=[x1;x2]  
  10. f=@(x)0.5*(x'*G*x)+b'*x+c;%目标函数  
  11. f_d=@(x)G*x+b;%一阶导数  
  12. %设定终止条件  
  13. N=100;E=10^(-6);  
  14.   
  15. n=1;  
  16. e=norm(abs(f_d(x0)));  
  17. X=x0;Y=f(x0);Y_d=e;  
  18. syms a real  
  19. while n<N && e>E  
  20.     n=n+1;  
  21.     f1=f(x0-a*f_d(x0));  
  22.     a0=solve(diff(f1,'a',1));  
  23.     a0=double(a0);  
  24.     x0=x0-a0*f_d(x0);  
  25.     X(:,n)=x0;Y(n)=f(x0);  
  26.     e=norm(abs(f_d(x0)));  
  27.     Y_d(n)=e;  
  28. end  
  29. toc  
  30. figure(1)  
  31. subplot(2,1,1)  
  32. semilogy(Y,'r*')  
  33. title(['函数最优值：' num2str(Y(n))])  
  34. ylabel('f(x)')  
  35. xlabel('迭代次数')  
  36. subplot(2,1,2)  
  37. semilogy(Y_d,'g<')  
  38. ylabel('f''(x)')  
  39. xlabel('迭代次数')  
  40. title(['导函数最优值：' num2str(Y_d(n))])  
  41. figure(2)  
  42. x1=-110:1:110;y1=x1;  
  43. [X1 Y1]=meshgrid(x1,y1);  
  44. nn=length(x1);  
  45. Z1=zeros(nn,nn);  
  46. for i=1:nn  
  47.     for j=1:nn  
  48.         Z1(i,j)=f([X1(i,j);Y1(i,j)]);  
  49.     end  
  50. end  
  51. hold on  
  52. contour(X1,Y1,Z1)  
  53. plot(X(1,:),X(2,:),'o-','linewidth',1)  
  54. end  


调用：

 

x0为初始点，可替换

 

[html] view plain copy
 
  1. x0=[-30;100];  
  2. [X Y Y_d]=SD(x0);  

 

 

结果：




B 牛顿方法

 

 

[plain] view plain copy
 
  1. function [X Y Y_d]=NT(x0)  
  2. %%%固定步长为1，Newton方法  
  3. %%%设定终止条件  
  4. tic  
  5. N=100;E=10^(-6);  
  6. %%%f为原函数，g为导函数，G为二阶导数  
  7. f=@(x)3*x(1).^2+3*x(2).^2-x(1).^2*x(2);  
  8. g=@(x)[6*x(1)-2*x(1)*x(2);6*x(2)-x(1)^2];  
  9. G=@(x)[6-2*x(2),-2*x(1);-2*x(1),6];  
  10. e=norm(g(x0));  
  11. X=x0;Y=f(x0);Y_d=e;  
  12. n=1;  
  13. while n<100 && e>E  
  14.     n=n+1;  
  15.     d=-inv(G(x0))*g(x0);  
  16.     %d=-pinv(G(x0))*g(x0);  
  17.     x0=x0+d;  
  18.     X(:,n)=x0;Y(n)=f(x0);  
  19.     e=norm(g(x0));  
  20.     Y_d(n)=e;  
  21. end  
  22. toc  
  23. figure(1)  
  24. subplot(2,1,1)  
  25. semilogy(Y,'r*')  
  26. title(['函数最优值：' num2str(Y(n))])  
  27. ylabel('f(x)')  
  28. xlabel('迭代次数')  
  29. subplot(2,1,2)  
  30. semilogy(Y_d,'g<')  
  31. ylabel('f''(x)')  
  32. xlabel('迭代次数')  
  33. title(['导函数最优值：' num2str(Y_d(n))])  
  34. figure(2)  
  35. x1=-6:0.05:6;y1=x1;  
  36. [X1 Y1]=meshgrid(x1,y1);  
  37. nn=length(x1);  
  38. Z1=zeros(nn,nn);  
  39. for i=1:nn  
  40.     for j=1:nn  
  41.         Z1(i,j)=f([X1(i,j);Y1(i,j)]);  
  42.     end  
  43. end  
  44. hold on  
  45. contour(X1,Y1,Z1,20)  
  46. plot(X(1,:),X(2,:),'go-','linewidth',1)  
  47. end  

 

[html] view plain copy
 
  1. %x0=[1.5;1.5];  
  2. %x0=[-2;4];  
  3. x0=[0;3];  
  4. %x0=[0;3]时矩阵为奇异的  
  5. %矩阵G为病态或奇异的时候，pinv能适当改善结果，更差的时候需采用正则化的方法  
  6. [X Y Y_d]=NT(x0);  
  7. disp('[1.5;1.5]的最优值：')  
  8. disp(min(Y))  

 

 

结果：



（3）比较

相同的算法，在不同的环境下，由于截断参数不同，会有些许的差异，MATLAB更易上手，Python却更严谨，但Python在符号运算时却有诸多不便，两者之间的差异，只能说各有所长各有所短，能结合起来用也许效果会更好。

 

posted @ 2017-08-02 16:27  于繁华求淡然  阅读(874)  评论(0)  编辑  收藏  举报