六月粒子群算法作业

问题：1.利用MATLAB编写课件中求解函数极值问题的代码。 2.利用MATLAB实现PSO求解多参数拟合问题。 3.编程实现基于PSO优化的灰度预测模型的建立。

 0.1：第一问图

 

0.2： 第二问图

注：第二问核心函数编程如下：

function MSE=fitness(para,x,y)

y_est=para(1) *sin( para(2) *x^2)+exp( para(3) *x)+ln( para(4) *sin(para(5)*x^3));

MSE=0;

for i=1:length(y)

MSE=MSE+(y_est(i)-y(i))^2;

end

MSE=MSE/length(y);

---

答案：

1：以下是利用MATLAB编写求解函数极值问题的代码：

% 定义参数
N = 50; % 粒子数量
c1 = 2; % 自身经验在计算粒子飞翔速度上的权重
c2 = 2; % 群体经验在计算粒子飞翔速度过程中的权重
w = 0.5; % 惯性因子
M = 1000; % 迭代次数
D = 2; % 维数，即自变量的个数
x = zeros(N, D); % 粒子群中粒子的位置
v = zeros(N, D); % 粒子群中粒子的飞行速度
p = zeros(N, 1); % 适应度函数值
y = zeros(N, D); % 粒子群中粒子的所有路程上的最优解的位置
Pbest = zeros(M, 1); % 每次迭代的最优解

% 初始化粒子位置和速度
for i = 1:N
    for j = 1:D
        x(i, j) = rand(); % 随机初始化粒子位置
        v(i, j) = rand(); % 随机初始化粒子速度
    end
end

% 迭代
for t = 1:M
    % 计算适应度函数值
    for i = 1:N
        p(i) = fitness(x(i, :)); % 计算粒子位置的适应度函数值
        if p(i) < fitness(y(i, :)) % 更新粒子的最优位置
            y(i, :) = x(i, :);
        end
    end
    
    % 更新全局最优位置
    [~, idx] = min(p);
    Gbest = y(idx, :);
    Pbest(t) = fitness(Gbest); % 记录每次迭代的最优解
    
    % 更新粒子速度和位置
    for i = 1:N
        v(i, :) = w * v(i, :) + c1 * rand() * (y(i, :) - x(i, :)) + c2 * rand() * (Gbest - x(i, :)); % 更新粒子速度
        x(i, :) = x(i, :) + v(i, :); % 更新粒子位置
    end
end

% 绘制迭代过程中最优解的变化曲线
plot(Pbest);
xlabel('迭代次数');
ylabel('最优解');
title('粒子群算法求解函数极值问题');

function f = fitness(x)
    f = (1 - x(1))^2 + 100 * (x(2) - x(1)^2)^2;
end

运行结果：

2 matlab运行代码：

function MSE = fitness(para, x, y)
    y_est = para(1)*sin(para(2)*x.^2) + exp(para(3)*x) + log(para(4)*sin(para(5)*x.^3));
    MSE = sum((y_est - y).^2)/length(y);
end

% 参数的取值范围
para_min = [0, 0, -10, 0, 0];
para_max = [10, 10, 10, 10, 10];

% 准备数据
x = linspace(0, 2*pi, 100);
y = cos(x) + randn(size(x))*0.1;

% 初始化粒子群
n_particles = 50;
n_iterations = 100;
w = 0.8;
c1 = 2;
c2 = 2;
particles = struct();
particles.position = rand(n_particles, 5) .* (para_max - para_min) + para_min;
particles.velocity = rand(n_particles, 5) .* (para_max - para_min) + para_min;
particles.best_position = particles.position;
particles.best_fitness = inf(n_particles, 1);

% 迭代更新粒子群
for i = 1:n_iterations
    % 计算每个粒子的适应度值
    fitness_values = zeros(n_particles, 1);
    for j = 1:n_particles
        fitness_values(j) = fitness(particles.position(j,:), x, y);
        if fitness_values(j) < particles.best_fitness(j)
            particles.best_fitness(j) = fitness_values(j);
            particles.best_position(j,:) = particles.position(j,:);
        end
    end
    
    % 更新全局最优解
    [global_best_fitness, global_best_index] = min(particles.best_fitness);
    global_best_position = particles.best_position(global_best_index,:);
    
    % 更新粒子的速度和位置
    for j = 1:n_particles
        r1 = rand(1,5);
        r2 = rand(1,5);
        particles.velocity(j,:) = w*particles.velocity(j,:) ...
            + c1*r1.*(particles.best_position(j,:) - particles.position(j,:)) ...
            + c2*r2.*(global_best_position - particles.position(j,:));
        particles.position(j,:) = particles.position(j,:) + particles.velocity(j,:);
        
        % 限制粒子的位置在参数的取值范围内
        particles.position(j,:) = max(particles.position(j,:), para_min);
        particles.position(j,:) = min(particles.position(j,:), para_max);
    end
end

% 计算最优参数
[best_fitness, best_index] = min(particles.best_fitness);
best_position = particles.best_position(best_index,:);

% 绘制拟合曲线
x_fit = linspace(0, 2*pi, 1000);
y_fit = best_position(1)*sin(best_position(2)*x_fit.^2) ...
    + exp(best_position(3)*x_fit) + log(best_position(4)*sin(best_position(5)*x_fit.^3));

figure;
plot(x, y, 'o');
hold on;
plot(x_fit, y_fit);
legend('原始数据', '拟合曲线');

运行结果：

 

3matlab运行代码：

clc;
clear;
close all;

% 数据准备
data = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100];
n = length(data);

% 灰度预测模型
x0 = data(1); % 初始值
alpha = 0.5; % 灰度生成参数
for i = 2:n
    x(i) = (data(i) - alpha * x0) / (1 - alpha);
    x0 = x(i);
end

% PSO优化
max_iter = 100; % 最大迭代次数
pop_size = 50; % 粒子群大小
c1 = 2; % 学习因子1
c2 = 2; % 学习因子2
w = 0.8; % 惯性权重
v_max = 5; % 最大速度限制
dim = 2; % 变量维度
lb = [0, 0]; % 变量下界
ub = [1, 1]; % 变量上界

% 初始化粒子群
for i = 1:pop_size
    pop(i, :) = lb + rand(1, dim) .* (ub - lb);
    vel(i, :) = -v_max + rand(1, dim) .* (2 * v_max);
    fitness(i) = inf;
    pbest(i, :) = pop(i, :);
    pbest_fitness(i) = inf;
end

% 迭代优化
for iter = 1:max_iter
    % 计算适应度
    for i = 1:pop_size
        a = pbest(i, 1);
        b = pbest(i, 2);
        y = x(1);
        for j = 2:n
            y = (1 - a) * (y + b * x(j));
        end
        fitness(i) = abs(y - data(n));
        if fitness(i) < pbest_fitness(i)
            pbest(i, :) = pop(i, :);
            pbest_fitness(i) = fitness(i);
        end
    end
    [gbest_fitness, gbest_index] = min(pbest_fitness);
    gbest = pbest(gbest_index, :);
    % 更新速度和位置
    for i = 1:pop_size
        vel(i, :) = w * vel(i, :) + c1 * rand(1, dim) .* (pbest(i, :) - pop(i, :)) + c2 * rand(1, dim) .* (gbest - pop(i, :));
        vel(i, :) = min(max(vel(i, :), -v_max), v_max);
        pop(i, :) = pop(i, :) + vel(i, :);
        pop(i, :) = min(max(pop(i, :), lb), ub);
    end
end

% 绘制预测结果
a = gbest(1);
b = gbest(2);
y = x(1);
for i = 2:n
    y = (1 - a) * (y + b * x(i));
end
figure;
plot(data, 'b-o');
hold on;
plot(y, 'r-*');
legend('原始数据', '预测结果');
xlabel('时间');
ylabel('数值');
title('灰度预测模型预测结果');

运行结果：