在博客上看到了一些函数，为了练习python画图能力，所以就都绘制了部分函数

文章目录
函数
代码和图像
画图部分
获取x和y的值
Ackely's function
Sphere function
Rosebrock function
Beale's function
GoldsteinPrice function
Booth's function
Bukin function
Matyas function
Levi function
Three-hump camel function
Eason function
Cross-in-tray function
Holder table function
McCormick function
StyblinskiTang function

函数





代码和图像
画图部分
def draw_pic_3D(x, y, z, title, z_min, z_max, offset):
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap = plt.get_cmap('rainbow'),color='orangered')
    # 绘制等高线
    ax.contour(x,y,z,offset=offset,colors='green')
    ax.set_zlim(z_min, z_max)
    ax.set_title(title)
    plt.savefig("image")
    plt.show()

获取x和y的值
def get_x_and_y(x_min, x_max, y_min, y_max):
    x = np.arange(x_min, x_max, 0.1)
    y = np.arange(y_min, y_max, 0.1)
    x, y = np.meshgrid(x, y)  # 生成网格点坐标矩阵
    return x, y

Ackely’s function
def Ackelys_function(z_min = 0,z_max = 15, offset = 0):
    x, y = get_x_and_y(-5,5,-5,5)
    z = -20 * np.exp(-0.2 * np.sqrt(0.5 * (x ** 2 + y ** 2))) - np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y))) + 20 + np.e
    return x, y, z, 'Ackely function', z_min, z_max, offset


Sphere function
def Sphere(z_min = 0,z_max = 20,offset = 0):
    x,y = get_x_and_y( -3,3,-3,3)
    z = x ** 2 + y ** 2
    return x,y,z, "Sphere function", z_min, z_max, offset


Rosebrock function
def Rosebrock_function(z_min = 0,z_max = 1500,offset = 0):
    x, y = get_x_and_y(-1, 1, -1, 1)
    z = 100 * ((y - x ** 2) ** 2 + (x -1) ** 2)
    return x,y,z, "Rosebrock function", z_min, z_max, offset


Beale’s function
def Beale_function(z_min = 0,z_max = 2000,offset=0):
    x, y = get_x_and_y(-4.5, 4.5, -4.5, 4.5)
    z = (1.5 - x + x * y) ** 2 +  (2.25 - x + x * y) ** 2 + (2.625 -x + x * y) ** 2
    return x, y, z, "Beale function", z_min, z_max, offset


GoldsteinPrice function
def GoldsteinPrice_function(z_min =0 ,z_max =1000000 ,offset=0):
    x, y = get_x_and_y(-2,2,-2,2)
    z = (1 + ((x + y + 1) ** 2) * (19 - 14 * x + 3 * (x ** 2) - 14 * y + 6 * x * y + 3 * (y **2))) * (30 + ((2 * x - 3 * y) ** 2) * (18 - 32 * x + 12 * (x ** 2) + 48 * y - 36 * x * y + 27 * (y ** 2)))
    return x, y, z, "GoldsteinPrice function", z_min, z_max, offset


Booth’s function
def Booth_function(z_min=0, z_max=2500, offset=0):
    x, y = get_x_and_y(-10,10,-10,10)
    z = (x + 2 * y -7) ** 2 + (2 * x + y - 5) ** 2
    return x, y, z, "Booth function", z_min, z_max, offset


Bukin function
def  Bukin_function(z_min =0 ,z_max =250 ,offset=0):
    x, y = get_x_and_y(-15,-5,-3,3)
    z = 100 * np.sqrt(np.fabs(y - 0.01 * (x ** 2))) + 0.01 * np.fabs(x + 10)
    return x, y, z, "Bukin function", z_min, z_max, offset


Matyas function
def MAtyas_function(z_min =0 ,z_max =100 ,offset=0):
    x, y = get_x_and_y(-10,10,-10,10)
    z = 0.26 * (x ** 2 + y ** 2) - 0.48 * x * y
    return x, y, z, "Matyas function", z_min, z_max, offset


Levi function
def Levi_function(z_min =0 ,z_max =450 ,offset=0):
    x, y = get_x_and_y(-10,10,-10,10)
    z = np.sin(3 * np.pi * x) ** 2 + ((x - 1) ** 2) * (1 + np.sin(3 * np.pi * y) ** 2) + ((y - 1) ** 2) * (1 + np.sin(2 * np.pi * y) ** 2)
    return x, y, z, "Levi function", z_min, z_max, offset


Three-hump camel function
def Three_hump_camel_function(z_min =0 ,z_max =2000 ,offset=0):
    x, y = get_x_and_y(-5,5,-5,5)
    z = 2 * (x ** 2) - 1.05 * (x ** 4) + (x ** 6) / 6 + x * y + (y ** 2)
    return x, y, z, "Three-hump camel function", z_min, z_max, offset


Eason function
def Eason_function(z_min =-1.5 ,z_max =0 ,offset=-1.5):
    x, y = get_x_and_y(-10,10,-10,10)
    z = -np.cos(x) * np.cos(y) * np.exp(-((x - np.pi) ** 2 + (y - np.pi) ** 2))
    return x, y, z, "Eason function", z_min, z_max, offset


Cross-in-tray function
def Cross_in_tray_function(z_min =-4 ,z_max =0 ,offset=-4):
    x, y = get_x_and_y(-10,10,-10,10)
    z = -0.0001 * (np.fabs(np.sin(x) * np.sin(y) * np.exp(np.fabs(100 - (np.sqrt((x ** 2) + (y ** 2)) / np.pi)))) + 1) ** 0.1
    return x, y, z, "Cross-in-tray function", z_min, z_max, offset


Holder table function
def Holder_table_function(z_min =-30 ,z_max =0 ,offset=-30):
    x, y = get_x_and_y(-10,10,-10,10)
    z = -np.fabs(np.sin(x) * np.cos(y) * np.exp(np.fabs(1 - (np.sqrt(x ** 2 + y ** 2))/np.pi)))
    return x, y, z, "Holder table function", z_min, z_max, offset


McCormick function
def McCormick_function(z_min =0 ,z_max =50 ,offset=0):
    x, y = get_x_and_y(-1.5,4,-3,4)
    z = np.sin(x + y) + (x - y) ** 2 - 1.5 * x + 2.5 * y + 1
    return x, y, z, "McCormick function", z_min, z_max, offset


StyblinskiTang function
def StyblinskiTang_function(z_min =-200 ,z_max =250 ,offset=-200):
    x, y = get_x_and_y(-5,5,-5,5)
    z1 = x ** 4 -16 * (x ** 2) + 5 * x
    z2 = y ** 4 -16 * (y ** 2) + 5 * y
    z = (z1 + z2) / 2
    return x, y, z, "StyblinskiTang function", z_min, z_max, offset


欢迎指正，共同进步😎

写了算法当然要测试, 传统的有CUTEr但是似乎没有函数的解析形式
http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml 
Dataset Name  Level of 
Difficulty  Model 
Classification  Number of 
Parameters  Number of 
Observations 
Source 

---------------------------------------------------------------------------­----- 

Misra1a  Lower  Exponential   2    14      Observed 
Chwirut2  Lower  Exponential   3    54      Observed 
Chwirut1  Lower  Exponential   3    214      Observed 
Lanczos3  Lower  Exponential   6    24      Generated 
Gauss1  Lower  Exponential   8    250      Generated 
Gauss2  Lower  Exponential   8    250      Generated 
DanWood  Lower  Miscellaneous   2    6      Observed 
Misra1b  Lower  Miscellaneous   2    14      Observed 

Kirby2  Average  Rational   5    151      Observed 
Hahn1  Average  Rational   7    236      Observed 
Nelson  Average  Exponential   3    128      Observed 
MGH17  Average  Exponential   5    33      Generated 
Lanczos1  Average  Exponential   6    24      Generated 
Lanczos2  Average  Exponential   6    24      Generated 
Gauss3  Average  Exponential   8    250      Generated 
Misra1c  Average  Miscellaneous   2    14      Observed 
Misra1d  Average  Miscellaneous   2    14      Observed 
Roszman1  Average  Miscellaneous   4    25      Observed 
ENSO  Average  Miscellaneous   9    168      Observed 

MGH09  Higher  Rational   4    11      Generated 
Thurber  Higher  Rational   7    37      Observed 
BoxBOD  Higher  Exponential   2    6      Observed 
Rat42  Higher  Exponential   3    9      Observed 
MGH10  Higher  Exponential   3    16      Generated 
Eckerle4  Higher  Exponential   3    35      Observed 
Rat43  Higher  Exponential   4    15      Observed 
Bennett5  Higher  Miscellaneous   3    154      Observed 
 

I leave no trace of wings in the air but I am glad I have had my flight. 
                                                                            ——Rabindranath Tagore，Stray Birds
粒子群优化算法（PSO）
Point 1 粒子群优化算法的自然界原型
鸟群在食物存在的空间内飞行觅食，食物在空间内的分布是不均匀的。每只鸟可能会凭借着自己的经验或是直觉，飞往它所觉得比较好的地方寻找食物。当一些鸟发现了更好的觅食地点时，鸟群间会有某种类似广播的沟通行为，渐渐地将整个鸟群引领至较佳的地点。
鸟类的这种通过群体协作而进行的觅食行为，是利用了社会中存在的相互影响的概念，来引领所有个体朝向最佳位置。还有其他生物群体（鱼群）都有这种社会性的行为，这种行为被抽象成了粒子群优化算法。
Point 2 粒子群优化算法的基本原理
粒子群中的每一个粒子（可以直接当成鸟）都是n维解空间中的一个点
每个粒子有自己的位置、速度
每个粒子在不同的位置都具有对应于目标函数的个体适应度（个体适应度相当于一个评价标准，对于鸟来说就是这个地方食物的多少，对于函数极值问题就是函数值的大小）。
对于每个粒子来说，在其所经历过的所有位置中，会有一个位置为最佳位置（个体最优解），记为 
对于粒子群体来说，所有粒子所处的位置中，会有一个位置为最佳位置（全局最优解），记为 
一开始先给定各粒子的初始位置和速度，并计算出个体最优解
与全局最优解 
。
不断更新按照以下公式更新粒子的速度和位置：
每次更新完位置之后，再次确定每个粒子的适应度，更新各粒子的个体最优解
，更新粒子群体的全局最优解 
每个粒子根据不断更新的个体最优解和全局最优解来更新自己的位置。 
 可理解为粒子的惯性，反映了粒子有维持自己先前运动的趋势。 
 可理解为粒子基于自我认知做出的判断， 
 可理解为粒子在协同合作中共享信息做出的判断。
Point 3 粒子群优化算法的实现
Point 4 粒子群优化解决函数优化问题代表性栗子——Rastrigin测试函数
Rastrigin测试函数
其函数图像如下：
使用PSO求解其极值，过程如下：
在1~80次迭代中，由上图可以看出收敛速度较快，但是在80次之后收敛速度就比较缓慢。这是因为PSO算法的趋同性，导致后期粒子之间距离接近，越往后速度越慢。
Rastrigin函数定义如下：
function f = Rastrigin(x, D)
    f = 0;
    for i = 1 : D
        f = f+(x(i)^2-10*cos(2*pi*x(i)));
    end 
    f = 10 * D + f;
作二维Rastrigin图像程序如下：
n=100;
AA = 5.12;
x=linspace(-AA,AA,n);
y=linspace(-AA,AA,n);
[x1,y1]=meshgrid(x,y);
for kki = 1:n
    for j = 1:n
        zz = [x1(kki,j); y1(kki,j)];
        z(kki,j) = Rastrigin(zz,2);
    end
end
meshc(x1,y1,z);
%contour(x1,y1,z); 
xlabel('x');
ylabel('y');
zlabel('z');
Point 5 粒子群优化解决函数优化问题栗子
利用粒子群优化算法求： 
 的极值点， 
目标函数定义如下：
function f = func(x, r)
    f = (x(1)^2 + x(2)^2)/1000;
主程序：
clc;
clear all;
close all;

%-------------------设置初始参数-------------------%
w = 0.6;    %设置权重
c1 = 2;   %设置加速度因子
c2 = 2;   %设置加速度因子

M = 60;        %设置粒子个数
D = 2;         %设置维数
N_max = 100;   %最大迭代次数

p_ub = 300;             %任意指定一个生成初始坐标的上界
p_lb = -300;            %任意指定一个生成初始坐标的下界
v_max_abs = 5;          %绝对值速度上界


%一些存储矩阵(元胞数组)
Position = [];          %粒子当前位置矩阵，D行M列矩阵，每一列向量对应一个粒子坐标
Fitness = [];           %粒子当前适应度矩阵，1行M列矩阵，每一个元素为Position矩阵对应坐标的函数值（适应度）
Velocity = [];          %粒子当前速度矩阵，D行M列矩阵，每一列向量对应一个粒子速度

Pbest_Fitness = [];             %当前个体最优解适应度的矩阵，1行M列矩阵，每一个元素为每个粒子个体最优解
Gbest_Fitness = 0;              %当前全局最优解适应度
Pbest_Position = [];            %当前个体最优解矩阵，D行M列矩阵，每一列向量对应一个坐标
Gbest_Position = [];            %当前全局最优解               
%-------------------生成粒子的初始位置、速度-------------------%
for i = 1:M
    for ii = 1:D
        rand_position(ii, 1) = p_lb + rand * (p_ub - p_lb);
        rand_velocity(ii, 1) = 2 * rand;
    end
    Position(:,i) = rand_position;         %初始位置
    Fitness(i) = func(rand_position);      %初始位置对应的适应度（函数值）
    Velocity(:,i) = rand_velocity;         %初始速度
end

%-------------------由初始坐标找到初始Pbest、Gbest-------------------%
Pbest_Position = Position;                           %个体最优解初始化为初始位置
Pbest_Fitness = Fitness;                             %个体最优适应度初始化为初始位置对应的适应度
[Gbest_Fitness, Gbest_n] = min(Pbest_Fitness);       %在个体最优适应度中找出全局最优适应度
Gbest_Position = Pbest_Position(:, Gbest_n);         %找出全局最优适应度对应的全局最优解的坐标

for N = 1:N_max
    %---------------将上一次的个体最优解、全局最优解录入历史数据---------------%
    History_Pbest_Position_x(N, :) = Pbest_Position(1, :);
    History_Pbest_Position_y(N, :) = Pbest_Position(2, :);
    History_Pbest_Fitness(N, :) = Pbest_Fitness;
    History_Pbest_Gbest_Position(N, 1) = Gbest_Position(1, :);
    History_Pbest_Gbest_Position(N, 2) = Gbest_Position(2, :);
    History_Gbest_Fitness(N, :) = Gbest_Fitness;
  
    %---------------更新位置和速度---------------%
    for i = 1:M                                         %遍历每一个粒子
        Velocity(:, i) = w * Velocity(:, i) + c1 * rand * (Pbest_Position(:, i) - Position(:, i)) + c2 * rand * (Gbest_Position - Position(:, i));
        for ii = 1:D                            %速度纠正
            if Velocity(D, ii) > v_max_abs
                Velocity(D, ii) = v_max_abs;
            elseif Velocity(D, ii) < -v_max_abs
                Velocity(D, ii) = -v_max_abs;
            end
        end
        Position(:, i) = Position(:, i) + Velocity(:, i);
        for ii = 1:D                            %位置纠正
            if Position(D, ii) > p_ub
                Position(D, ii) = p_ub;
            elseif Position(D, ii) < -p_lb
                Position(D, ii) = -p_lb;
            end
        end
    end
    %---------------更新个体适应度和全局适应度---------------%
    for i = 1:M                                         %遍历每一个粒子
        Fitness(i) = func(Position(:, i));              %计算更新位置后粒子的适应度
        if Fitness(i) < Pbest_Fitness(i)                %如果i粒子的新适应度小于个体最优解的适应度
            Pbest_Fitness(i) = Fitness(i);              %更新个体最优适应度 
            Pbest_Position(:, i) = Position(:, i);      %更新个体最优解
        end
    [Gbest_Fitness, Gbest_n] = min(Pbest_Fitness);
    Gbest_Position = Pbest_Position(:,Gbest_n);         %全局最优解的坐标
    end
end
    

fprintf('蚁群算法求解函数极值：n')
fprintf('最小值点为：[%f; %f]n', Gbest_Position)
fprintf('最小值为：%fn',Gbest_Fitness)

plot(1:1:N_max,History_Gbest_Fitness);


% for i = 1:N_max
%     figure(i)
%     plot(History_Pbest_Position_x(i,:),History_Pbest_Position_y(i,:),'g.','markersize',30);
%     axis([-p_ub p_ub -p_ub p_ub])
%     saveas(figure(i),['C:UsersAirHeDesktop1',int2str(i),'.jpg']);
% end
计算结果如下：

下面是下载的demo源码（运行于Linux环境下），作者加上自己的注释方便理解。
#include "Header.h"
#include <sys/time.h>
#include <cstdio>
#include <unistd.h>

int main(){
  /*  Test the basic benchmark function */
  double* X;                                 # 定义粒子指针
  Benchmarks* fp=NULL;                       # 定义函数对象指针（不确定）
  unsigned dim = 1000;                       # 粒子维度，表示每个粒子所包含的参数个数
  unsigned funToRun[] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};     # 要运行的函数编号
  // unsigned funToRun[] = {1};
  // unsigned funToRun[] = {15};
  unsigned funNum = 15;                      # 要运行的函数个数
  unsigned run = 1;                          # 运行次数（不确定）

  vector<double> runTimeVec;                 # 时间相关，windows下可删除
  struct timeval start, end;
  long seconds, useconds;    
  double mtime;
        
  X = new double[dim];                       # 定义粒子并初始化
  for (unsigned i=0; i<dim; i++){ 
    X[i]=0;
  }

  for (unsigned i=0; i<funNum; i++){         # 遍历函数并生成函数对象，打印编号和计算值；嵌套for循环
    fp = generateFuncObj(funToRun[i]); 
    printf("F %d value = %1.20E\n", fp->getID(), fp->compute(X));
    gettimeofday(&start, NULL);
    for (unsigned j=0; j < run; j++){
      fp->compute(X);
    }
    gettimeofday(&end, NULL);
    
    seconds  = end.tv_sec  - start.tv_sec;
    useconds = end.tv_usec - start.tv_usec;

    mtime = (((seconds) * 1000 + useconds/1000.0) + 0.5)/1000;

    runTimeVec.push_back(mtime);
    printf ( "F %d, Running Time = %f s\n\n", fp->getID(), mtime);
    
    delete fp;
  }

  delete []X;

  // for (unsigned i=0; i<runTimeVec.size(); i++){
  // 	printf ( "%f\n", runTimeVec[i] );
  // }

  return 0;
}

// create new object of class with default setting
Benchmarks* generateFuncObj(int funcID){                   # 根据输入函数ID创建类对象
  Benchmarks *fp;
  // run each of specified function in "configure.ini"
  if (funcID==1){
    fp = new F1();
  }else if (funcID==2){
    fp = new F2();
  }else if (funcID==3){
    fp = new F3();
  }else if (funcID==4){
    fp = new F4();
  }else if (funcID==5){
    fp = new F5();
  }else if (funcID==6){
    fp = new F6();
  }else if (funcID==7){
    fp = new F7();
  }else if (funcID==8){
    fp = new F8();
  }else if (funcID==9){
    fp = new F9();
  }else if (funcID==10){
    fp = new F10();
  }else if (funcID==11){
    fp = new F11();
  }else if (funcID==12){
    fp = new F12();
  }else if (funcID==13){
    fp = new F13();
  }else if (funcID==14){
    fp = new F14();
  }else if (funcID==15){
    fp = new F15();
  }else{
    cerr<<"Fail to locate Specified Function Index"<<endl;
    exit(-1);
  }
  return fp;
}