问题: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('灰度预测模型预测结果');
运行结果: