基于粒子群最小二乘算法的多项式参数拟合方法

最小二乘法误差平方和作为粒子好坏判据，粒子群算法（PSO）用来寻找极值（也就是最佳拟合参数），二者结合，基于粒子群最小二乘算法的多项式参数拟合方法，有问题请指出

1.问题提出

给定一组数据点$(x_i,y_i)$,目标是找到多项式$P(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$的系数
$a_0,a_1,\dots ,a_n$,使得误差最小。

例程中给出的点如下，且多项式的最高次数n为4,需要拟合的点如下

x = [1, 2, 3, 4, 5];
y = [2.1, 4.9, 9.8, 16.1, 25.2];

2.最小二乘法

最小二乘法通过最小化误差平方和来拟合多项式：

$$E=\sum _{i=1}^m(y_i-P(x_i){)}^2$$

其中，$m$是数据点数量。

3.粒子群算法介绍(PSO)

PSO 是一种优化算法，通过粒子群搜索最优解。每个粒子代表一个可能的解 (即一组多项式系数),并根
据个体和群体经验更新位置和速度。

4.结合过程

1.初始化粒子群：随机生成粒子，每个粒子包含一组多项式系数$a_0,a_1,\dots ,a_n$。
2.计算适应度：使用最小二乘法计算每个粒子的适应度(即误差$E$)。3.更新个体和群体最优：记录每个粒子的历史最优位置和群体最优位置。4. 更新粒子速度和位置：

$$v_{id}=w\cdot v_{id}+c_1\cdot r_1\cdot (p_{id}-x_{id})+c_2\cdot r_2\cdot (p_{gd}-x_{id})\text{\hspace{1em}(4.1)}$$
$$x_{id}=x_{id}+v_{id}\text{\hspace{1em}(4.2)}$$

其中，$v_{id}$是速度，$x_{i}d$是位置，$w$是惯性权重，$c_1$和$c_2$是学习因子，$r_1$和$r_2$是随机数。$p_{id}$是当前粒子最优值，$p_{gd}$是粒子群中最优值
5.迭代：重复计算适应度、更新最优位置和粒子状态，直到达到停止条件 (如最大迭代次数或误差阈

值)。

5代码算法初步讲解

5.1首先要搞清楚，粒子是什么？

不同的问题中，粒子不一样，本问题目的是为了拟合最好的参数$a_0$到$a_4$，所以每个粒子是一个向量[a0 a1 a2 a3 a4],如果问题是给定函数$y=x^2+x+1$，求函数在范围[-1 1]的极值，那么每个粒子是一个点。

5.2粒子群算法的核心思想是什么

是反馈和搜索

在本问题中，反馈的判据是对于一个粒子（特定的[a0 a1 a2 a3 a4]）组成的四阶多项式，带入所拟合的点的x,理想y值和实际y值的误差平方和err。举个简单的例子，有5个粒子，带入要拟合的点后，产生了5个方差和err1-err5,这里面err5最小，那么粒子1到粒子4的参数将会朝着粒子5的参数移动，同时由粒子群算法的公式可得，粒子1到粒子4的历史最优值，也会影响粒子移动的方向。粒子5因为速度不为0，所以也会移动，直到迭代次数完毕，选出了迭代过程中，最优的粒子。所以，粒子的初值很重要，迭代次数也很重要，粒子移动的速度也很重要

5.3简单理解一下4中的两个公式

4中第二个公式，位置是速度的积分。没啥好说的。第一个公式，速度和之前的速度有个比例分量，类似于PID中的比例调节，和粒子最优位置，和粒子群的最优位置，相乘一个随机分量再相加，类比PID中的积分调节。随机分量不至于让粒子陷入局部死区。

6.代码例程

以四阶多项式参数拟合为例，matlab代码如下

% 清空环境
clc;
clear;

% 数据点
x = [1, 2, 3, 4, 5];
y = [2.1, 4.9, 9.8, 16.1, 25.2];

% 多项式阶数
n = 4; % 四次多项式

% PSO 参数
num_particles = 30; % 粒子数量
max_iter = 100;     % 最大迭代次数
w = 0.7;            % 惯性权重
c1 = 1.5;           % 个体学习因子
c2 = 1.5;           % 群体学习因子

% 初始化粒子位置和速度
particle_pos = rand(num_particles, n+1); % 每个粒子代表一组系数 [a0, a1, a2, a3, a4]
particle_vel = rand(num_particles, n+1); % 粒子速度

% 初始化个体最优位置和适应度
particle_best_pos = particle_pos;
particle_best_fitness = inf(num_particles, 1);

% 初始化全局最优位置和适应度
global_best_pos = particle_pos(1, :);
global_best_fitness = inf;

% PSO 迭代
for iter = 1:max_iter
    for i = 1:num_particles
        % 提取当前粒子的系数
        a = particle_pos(i, :);
        
        % 计算多项式拟合值
        y_fit = a(1) + a(2) * x + a(3) * x.^2 + a(4) * x.^3 + a(5) * x.^4;
        
        % 计算误差（适应度）
        fitness = sum((y - y_fit).^2); % 最小二乘误差
        
        % 更新个体最优
        if fitness < particle_best_fitness(i)
            particle_best_fitness(i) = fitness;
            particle_best_pos(i, :) = a;
        end
        
        % 更新全局最优
        if fitness < global_best_fitness
            global_best_fitness = fitness;
            global_best_pos = a;
        end
    end
    
    % 更新粒子速度和位置
    for i = 1:num_particles
        r1 = rand(1, n+1);
        r2 = rand(1, n+1);
        particle_vel(i, :) = w * particle_vel(i, :) ...
            + c1 * r1 .* (particle_best_pos(i, :) - particle_pos(i, :)) ...
            + c2 * r2 .* (global_best_pos - particle_pos(i, :));
        particle_pos(i, :) = particle_pos(i, :) + particle_vel(i, :);
    end
    
    % 显示当前迭代结果
    fprintf('Iteration %d: Best Fitness = %.4f\n', iter, global_best_fitness);
end

% 输出最终结果
fprintf('Optimal Coefficients:\n');
fprintf('a0 = %.4f, a1 = %.4f, a2 = %.4f, a3 = %.4f, a4 = %.4f\n', ...
    global_best_pos(1), global_best_pos(2), global_best_pos(3), ...
    global_best_pos(4), global_best_pos(5));

% 绘制拟合结果
x_fit = linspace(min(x), max(x), 100);
y_fit = global_best_pos(1) + global_best_pos(2) * x_fit + global_best_pos(3) * x_fit.^2 ...
    + global_best_pos(4) * x_fit.^3 + global_best_pos(5) * x_fit.^4;

figure;
plot(x, y, 'bo', 'DisplayName', 'Data Points');
hold on;
plot(x_fit, y_fit, 'r-', 'DisplayName', 'Fitted Polynomial');
xlabel('x');
ylabel('y');
title('PSO + Least Squares 4th Degree Polynomial Fitting');
legend;
grid on;

7.代码运行结果

a0 = 0.6049, a1 = 1.2357, a2 = 0.1376, a3 = 0.1926, a4 = -0.0147