前言

今天更新较晚主要还是学业繁忙，学习素材也不是很好找，可能很多同学们都在做数学建模以及应用统计时都会涉及到函数参数拟合的问题，一般最常用的方法是最小二乘法，但是当函数参数很多时，往往去普通最小二乘法迭代次数过多，容易陷入局部最优解的情况，因此，在参数拟合部分会用到现代优化算法，现代优化算法常用的有遗传算法，模拟退火，蚁群算法，粒子群算法等等由于之前地统计实验遗传写的拟合变异函数比较容易陷入局部最优，也找了不少资料去优化，但是仍然找不到比较好的方法，因此还是换一个比较好的算法来进行拟合效果会比较好，通常情况下，粒子群的收敛要比遗传收敛性要好，这句话不记得是在哪儿说的了。

1基本原理介绍

1.1最小二乘拟合

给定输入输出数列x,yx,yx,y求参量ccc使得
min∑i(F(c,xi)−yi)2min \sum_{i}(F(c,x_i)-y_i)^2 mini∑(F(c,xi)−yi)2
其实就是要求得在参数ccc的情况下，函数拟合值与真实值的残差平方和最小，即为最小二乘拟合参数的原理

1.2粒子群算法

粒子群优化算法(PSO：Particle swarm optimization) 是一种进化计算技术（evolutionary computation）。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解．
PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域
算法流程如下：
1、初始化
首先，我们设置最大迭代次数，目标函数的自变量个数，粒子的最大速度，位置信息为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为M，每个粒子随机初始化一个飞翔速度。
2、个体极值与全局最优解
定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。
3、更新速度和位置的公式
更新公式为：
Vid=wVid+C1random(0,1)(pid−Xid)+C2random(0,1)(Pgd−Xid)Xid=Xid+VidV_{id}=wV_{id}+C_1random(0,1)(p_{id}-X_{id})+C_2random(0,1)(P_{gd}-X_{id})\\ X_{id}=X_{id}+V_{id} Vid=wVid+C1random(0,1)(pid−Xid)+C2random(0,1)(Pgd−Xid)Xid=Xid+Vid
其中，www为惯性因子，C1C_1C1和C2C_2C2称为加速常数，一般取C1=C2∈[0,4]C_1=C_2\in[0,4]C1=C2∈[0,4],random(0,1)random(0,1)random(0,1)表示在区间[0,1][0,1][0,1]上的随机数，PidP_{id}Pid表示第iii个变量的个体极值的第ddd维，PgdP_{gd}Pgd表示全局最优解的第ddd维
4、终止条件
有两种终止条件可以选择，一时最大代数：TmaxT_{max}Tmax;二是相邻两代之间的偏差在一个设定的范围内停止，这里选择第一种
，具体流程如下图，该图是引用一篇博客中见参考链接，

2案例实践

2.1案例介绍

此次案例来源于数学建模算法与应用(第二版)中的第五章插值与拟合部分，观测数据如下：

要求如下：
根据上述观测数据，拟合函数
y=e−k1x1sin(k2x2)+x32y=e^{-k_1x_1}sin(k_2x_2)+x^2_3 y=e−k1x1sin(k2x2)+x32
可以看到要求拟合的参数有两个,即k1,k2k_1,k_2k1,k2

2.2最小二乘优化

1.编写拟合参数的函数m文件

function f=fun(Pos,x)
%其中Pos(1)=k1,Pos(2)=k2;x1,x2,x3分别为x(:,1)，x(:,2)，x(:,3)
f=exp(-Pos(1)*x(:,1)).*sin(Pos(2).*x(:,2))+x(:,3).^2;
end

2.调用lsqcurvefit函数命令

clear;clc
a=xlsread('C:\Users\asus\Desktop\data1.xlsx');
y0=a(:,2);x0=a(:,3:5);%准备原始数据，x0第1,2,3列代表x1,x2,x3
canshu0=rand(2,1);%拟合参数的初始值任意选取
lb=[-5,-5];%参数下界
ub=[5;5];%参数上界
canshu=lsqcurvefit(@fun,canshu0,x0,y0,lb,ub);%调用命令
yp=fun(canshu,x0);%计算拟合值
plot(yp);
hold on;
plot(y0,'*');%绘图查看拟合效果
RSS=sum((yp-y0).^2);%计算残差平方和

结果如下：

普通最小二乘的拟合结果，k1=−0.0915,k2=0.3169k_1=-0.0915,k_2=0.3169k1=−0.0915,k2=0.3169,其中残差平方和RSS=297.0634RSS=297.0634RSS=297.0634

2.3粒子群优化

1.数据准备与参数设置

clear;clc
a=xlsread('C:\Users\asus\Desktop\data1.xlsx');
y0=a(:,2);x0=a(:,3:5);%准备原始数据，x0第1,2,3列代表x1,x2,x3
c1=2;%学习因子
c2=2;%学习因子
Dimension=2;%因为拟合的是两个参数k1,k2
Size=50;%设置粒子规模为50个
Tmax=500;%最大迭代次数
Vmaximum=10;%粒子最大速度
k1max=5;k1min=-5;
k2max=5;k2min=-5;%设置拟合参数的区间范围

2.粒子群初始化

position=zeros(Dimension,Size);velocity=zeros(Dimension,Size);
%初始化例子的位置与速度
vmax(1:Dimension)=Vmaximum;vmin(1:Dimension)=-Vmaximum;
%设置速度上下界
xmax=[k1max,k2max];xmin=[k1min,k2min];%设置位置即拟合参数上下限
[position,velocity]=Initial_position_velocity(Dimension,Size,xmax,xmin,vmax,vmin);function [Position,Velocity] = Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin)for i=1:DimensionPosition(i,:)=Xmin(i)+(Xmax(i)-Xmin(i))*rand(1,Size); % 产生合理范围内的随机位置，rand(1,Size)用于产生一行Size列个随机数Velocity(i,:)=Vmin(i)+(Vmax(i)-Vmin(i))*rand(1,Size);end
end

3.设置初始迭代前当前最优以及历史最优

pbest_position=position;%粒子的历史最优位置，初始值为粒子的起始位置，存储每个粒子的历史最优位置
gbest_position=zeros(Dimension,1);%全局最优的那个粒子所在位置，初始值认为是第1个粒子
for j=1:SizePos=position(:,j);%取第j列，即第j个粒子的位置fz(j)=Fitness_Function(Pos,x0,y0);%计算第j个粒子的适应值
end
[gbest_fitness,I]=min(fz);%求出所有适应值中最小的那个适应值，并获得该粒子的位置
gbest_position=position(:,I);%取最小适应值的那个粒子的位置，即I列

在这里适应度函数的适应度值为各个粒子群下每个个体对应的参数拟合的残差平方和，见如下代码：

function fit=Fitness_Function(Pos,x,y)yp=fun(Pos,x);%计算参数拟合下的值Pos为两个参数k1,k2的值,x为x1,x2,x3得值%yp为在该粒子下计算得y值，记为yp与y同维度，因此计算残差平方和视为适应度，使其最小fit=sum((yp-y).^2);
end

这里的fun函数就是上面编写的计算拟合值的函数m文件
4.开始迭代找到历史最优粒子群
在这里本例采用粒子群算法中惯性因子的参数设置尤为重要，为了更好地平衡算法的全局搜索能力与局部搜索能力，Shi.Y提出了线性递减惯性权重(LDIW)
即：
w(k)=wend+(wstart−wend)∗(Tmax−k)/Tmaxw(k) = w_{end} + (w_{start}- w_{end})*(T_{max}-k)/T_{max} w(k)=wend+(wstart−wend)∗(Tmax−k)/Tmax。*其中wstartw_{start}wstart 为初始惯性权重，wendw_{end}wend 为迭代至最大次数时的惯性权重;kkk为当前迭代次数， TmaxT_{max}Tmax为最大迭代次数。一般来说，wstart=0.9w_{start}=0.9wstart=0.9,wend=0.4w_{end}=0.4wend=0.4时，算法的性能最好。这样随着迭代的进行，惯性权重从0.9递减到0.4,迭代初期较大的惯性权重使算法保持了较强的全局搜索能力。而迭代后期较小的惯性权重有利于算法进行更精确的局部搜索

w_start=0.9;w_end=0.4;
for itrtn=1:Tmax%Weight=1;%惯性因子Weight=w_end+(w_start-w_end)*(Tmax-itrtn)/Tmax;%Weight=((0.95-0.4)/(1-Tmax))*(itrtn-1)+0.95;r1=rand(1);r2=rand(1);%进行速度更新for i=1:Sizevelocity(:,i)=Weight*velocity(:,i)+c1*r1*(pbest_position(:,i)-position(:,i))+c2*r2*(gbest_position-position(:,i));end%限制速度边界for i=1:Sizefor row=1:Dimensionif velocity(row,i)>vmax(row)velocity(row,i)=vmax(row);elseif velocity(row,i)<vmin(row)velocity(row,i)=vmin(row);elseendendend%大于最大速度的用最大速度代替，小于最小速度的用最小速度代替position=position+velocity;%位置更新%限制位置边界for i=1:Sizefor row=1:Dimensionif position(row,i)>xmax(row)position(row,i)=xmax(row);elseif position(row,i)<xmin(row)position(row,i)=xmin(row);elseendendendfor j=1:Sizep_position=position(:,j);%取一个粒子的位置fitness_p(j)=Fitness_Function(p_position,x0,y0);%计算第j个粒子适应度的值if fitness_p(j)< fz(j) %粒子的适应值比运动之前的适应值要好，更新原来的适应值pbest_position(:,j)=position(:,j);fz(j)=fitness_p(j);endif fitness_p(j)<gbest_fitnessgbest_fitness=fitness_p(j);%如果该粒子比当前全局适应度的值还好，则代替endend[gbest_fitness_new,I]=min(fz);%更新后的所有粒子的适应值，取最小的那个，并求出其编号best_fitness(itrtn)=gbest_fitness_new; %记录每一代的最好适应值gbest_position=pbest_position(:,I);%最好适应值对应的个体所在位置
end

5.结果检验
RSS=297.0634,k1=−0.0915,k2=0.3169RSS=297.0634,k_1=-0.0915,k_2=0.3169RSS=297.0634,k1=−0.0915,k2=0.3169,
历代迭代次数的最优适应度如下

在200次迭代之后就会收敛最后变得相对比较平稳，拟合效果如下：

3总结

1.本次粒子群算法拟合函数参数的博客编写给我的提升主要是在于对于matlab的代码编写，在矩阵的基础上来去实习那些伪代码，希望同学们可以借鉴方法来对其他的优化问题进行套用，比较复杂的优化问题主要是复杂在适应度函数的编写这一部分，虽然体验过数模比赛二维装箱问题那么一大长串的适应度函数代码，写代码写一天结果还找不出错误出来，在matlab实现算法的过程中要弄清楚每个矩阵变量的维数是怎么变化的，这一层面是需要我们写代码实现算法过程中最重要的核心问题，如果觉得算法太复杂，可以现在matlab交互的命令行下一行一行运行，对每个变量进行记录，其矩阵的大小在for循环，乘积开方等数学运算中是怎么变化的。
2.粒子群算法本身是一个全局优化算法，特别是在参数设置方面调整时比遗传算法的交叉概率，变异概率调整起来要流畅一些，也体现了各种优化算法的特点。
3.从算法的维度上来看，粒子群算法特别适用于多维参数以及自变量的模拟，大家可以想想如果参数过多，遗传算法就必须要设置多个种群，多个估值精度等等各项参数，有时还需要用到元胞数组等来批量调用，但是粒子群算法就不同，其将各个不同的维度全部都集中于一个速度与位置矩阵当中，无论维数有多少，我只需要更高矩阵的规模即可，这也是为什么粒子群算法在优化算法中比较受欢迎

4参考链接

5附录代码

clear;clc
a=xlsread('C:\Users\asus\Desktop\data1.xlsx');
y0=a(:,2);x0=a(:,3:5);%准备原始数据，x0第1,2,3列代表x1,x2,x3
c1=2;%学习因子
c2=2;%学习因子
Dimension=2;%因为拟合的是两个参数k1,k2
Size=50;%设置粒子规模为50个
Tmax=500;%最大迭代次数
Vmaximum=10;%粒子最大速度
k1max=5;k1min=-5;
k2max=5;k2min=-5;%设置拟合参数的区间范围
position=zeros(Dimension,Size);velocity=zeros(Dimension,Size);
%初始化例子的位置与速度
vmax(1:Dimension)=Vmaximum;vmin(1:Dimension)=-Vmaximum;
%设置速度上下界
xmax=[k1max,k2max];xmin=[k1min,k2min];%设置位置即拟合参数上下限
[position,velocity]=Initial_position_velocity(Dimension,Size,xmax,xmin,vmax,vmin);
pbest_position=position;%粒子的历史最优位置，初始值为粒子的起始位置，存储每个粒子的历史最优位置
gbest_position=zeros(Dimension,1);%全局最优的那个粒子所在位置，初始值认为是第1个粒子
for j=1:SizePos=position(:,j);%取第j列，即第j个粒子的位置fz(j)=Fitness_Function(Pos,x0,y0);%计算第j个粒子的适应值
end
[gbest_fitness,I]=min(fz);%求出所有适应值中最小的那个适应值，并获得该粒子的位置
gbest_position=position(:,I);%取最小适应值的那个粒子的位置，即I列
%开始循环找到最优解
w_start=0.9;w_end=0.4;
for itrtn=1:Tmax%Weight=1;%惯性因子Weight=w_end+(w_start-w_end)*(Tmax-itrtn)/Tmax;%Weight=((0.95-0.4)/(1-Tmax))*(itrtn-1)+0.95;r1=rand(1);r2=rand(1);%进行速度更新for i=1:Sizevelocity(:,i)=Weight*velocity(:,i)+c1*r1*(pbest_position(:,i)-position(:,i))+c2*r2*(gbest_position-position(:,i));end%限制速度边界for i=1:Sizefor row=1:Dimensionif velocity(row,i)>vmax(row)velocity(row,i)=vmax(row);elseif velocity(row,i)<vmin(row)velocity(row,i)=vmin(row);elseendendend%大于最大速度的用最大速度代替，小于最小速度的用最小速度代替position=position+velocity;%位置更新%限制位置边界for i=1:Sizefor row=1:Dimensionif position(row,i)>xmax(row)position(row,i)=xmax(row);elseif position(row,i)<xmin(row)position(row,i)=xmin(row);elseendendendfor j=1:Sizep_position=position(:,j);%取一个粒子的位置fitness_p(j)=Fitness_Function(p_position,x0,y0);%计算第j个粒子适应度的值if fitness_p(j)< fz(j) %粒子的适应值比运动之前的适应值要好，更新原来的适应值pbest_position(:,j)=position(:,j);fz(j)=fitness_p(j);endif fitness_p(j)<gbest_fitnessgbest_fitness=fitness_p(j);%如果该粒子比当前全局适应度的值还好，则代替endend[gbest_fitness_new,I]=min(fz);%更新后的所有粒子的适应值，取最小的那个，并求出其编号best_fitness(itrtn)=gbest_fitness_new; %记录每一代的最好适应值gbest_position=pbest_position(:,I);%最好适应值对应的个体所在位置
end
plot(best_fitness,'b','Linewidth',2);
xlabel('Number of iterations','FontSize',14);
ylabel('best_fitness','FontSize',14);
title('Optimal fitness change','FontSize',20);
set(gca,'Fontname','times new Roman');
figure;
yp=fun(gbest_position,x0);
plot(yp,'b','Linewidth',2);
hold on;
plot(y0,'r*');
xlabel('Serial number','FontSize',14);
ylabel('value','FontSize',14);
legend('fitted value','actual value');%添加图例
legend('boxoff','FontSize',14);%删除图例的轮廓宇背景
title('Fitting effect','FontSize',20);
set(gca,'Fontname','times new Roman');
RSS=sum((yp-y0).^2);%计算残差平方和

function [Position,Velocity] = Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin)for i=1:DimensionPosition(i,:)=Xmin(i)+(Xmax(i)-Xmin(i))*rand(1,Size); % 产生合理范围内的随机位置，rand(1,Size)用于产生一行Size列个随机数Velocity(i,:)=Vmin(i)+(Vmax(i)-Vmin(i))*rand(1,Size);end
end

function fit=Fitness_Function(Pos,x,y)yp=fun(Pos,x);%计算参数拟合下的值Pos为两个参数k1,k2的值,x为x1,x2,x3得值%yp为在该粒子下计算得y值，记为yp与y同维度，因此计算残差平方和视为适应度，使其最小fit=sum((yp-y).^2);
end

function f=fun(Pos,x)
%其中Pos(1)=k1,Pos(2)=k2;x1,x2,x3分别为x(:,1)，x(:,2)，x(:,3)
f=exp(-Pos(1)*x(:,1)).*sin(Pos(2).*x(:,2))+x(:,3).^2;
end

clear;clc
a=xlsread('C:\Users\asus\Desktop\data1.xlsx');
y0=a(:,2);x0=a(:,3:5);%准备原始数据，x0第1,2,3列代表x1,x2,x3
canshu0=rand(2,1);%拟合参数的初始值任意选取
lb=[-5,-5];%参数下界
ub=[5;5];%参数上界
canshu=lsqcurvefit(@fun,canshu0,x0,y0,lb,ub);%调用命令
yp=fun(canshu,x0);%计算拟合值
plot(yp,'b','Linewidth',2);
hold on;
plot(y0,'r*');%绘图查看拟合效果
xlabel('Serial number','FontSize',14);
ylabel('value','FontSize',14);
legend('fitted value','actual value');%添加图例
legend('boxoff','FontSize',14);%删除图例的轮廓宇背景
title('Fitting effect','FontSize',20);
set(gca,'Fontname','times new Roman');
RSS=sum((yp-y0).^2);%计算残差平方和

基于粒子群算法与最小二乘拟合函数参数相关推荐

基于粒子群算法的神经网络非线性函数拟合
基于粒子群算法的神经网络非线性函数拟合文章初心最近在学机器学习,自己的方向是智能算法,课程报告需要,于是试着把机器学习和粒子群算法相结合,写出来供大家参考,交流. 文末有这部分内容相关的代码,已开 ...
【回归预测-lssvm】基于粒子群算法优化最小二乘支持向量机lssvm实现数据回归预测附matlab代码
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
【机器学习】基于粒子群算法的非线性函数寻优
本微信图文介绍了基于粒子群算法的非线性函数寻优过程,并利用Matlab实现.
基于粒子群算法的组卷系统的研究与实现
摘要组卷系统的主要任务是根据用户的需要用当前数据库中的试题组成一套符合用户需求的试卷.随着数据库与题量增大,传统采用随机选取和回朔试探法的组卷抽提算法因其抽题时间长,占用的空间复杂度太大,容易陷入 ...
【ELM预测】基于粒子群算法PSO优化极限学习机预测含Matlab源码
1 模型为了提高空气质量预测精度,提出一种基于粒子群算法优化极限学习机的空气质量预测模型.运用粒子群算法优化极限学习机的初始权值和偏置,在保证预测误差最小的情况下实现空气质量最优预测.选择平均绝对百 ...
【回归预测-ELM预测】基于粒子群算法PSO优化极限学习机预测附matlab代码
1 内容介绍风电功率预测为电网规划提供重要的依据,研究风电功率预测方法对确保电网在安全稳定运行下接纳更多的风电具有重要的意义.针对极限学习机(ELM)回归模型预测结果受输入参数影响的问题,现将粒子群 ...
matlab 重叠峰分解算法,一种基于粒子群算法的光谱重叠峰分解方法与流程
本发明涉及一种基于粒子群算法的光谱重叠峰分解方法. 背景技术: 由于探测器能量分辨率等原因,峰位接近且峰宽较大的不同谱峰之间常常出现严重重叠干扰的现象,要对光谱作进一步较为准确.全面的成分定量和定性分 ...
基于粒子群算法优化的Elman神经网络数据预测-附代码
基于粒子群算法优化的Elman神经网络数据预测 - 附代码文章目录基于粒子群算法优化的Elman神经网络数据预测 - 附代码 1.Elman 神经网络结构 2.Elman 神经用络学习过程 3.电 ...
回归预测 | MATLAB实现PSO-LSSVM粒子群算法优化最小二乘支持向量机多输入单输出
回归预测 | MATLAB实现PSO-LSSVM粒子群算法优化最小二乘支持向量机多输入单输出目录回归预测 | MATLAB实现PSO-LSSVM粒子群算法优化最小二乘支持向量机多输入单输出预测效 ...

基于粒子群算法与最小二乘拟合函数参数

前言