粒子群算法求解非线性方程组matlab简单实现_WealdPoem

http://blog.sina.com.cn/u/3608600963

首页博文目录关于我

个人资料

微博

加好友发纸条

写留言加关注

博客等级：
博客积分：

博客访问：
关注人气：
获赠金笔：0支
赠出金笔：0支
荣誉徽章：

正文字体大小：大中小

粒子群算法求解非线性方程组matlab简单实现

(2019-03-12 15:57:37)

标签：

智能算法

it

matlab

粒子群优化算法

非线性方程组求解

非线性方程组的一般形式为

f 1 (x 1 , x 2 , …, x n )=0

f 2 (x 1 , x 2 , …, x n )=0

…

f m (x 1 , x 2 , …, x n )=0

式中, x 1 ...x n为所要求解的 n 个未知变量, f 1...f n 为定义在 m 维空间中的实值函数。用 PSO 求解非线性方程组时,设方程组中的未知数 x 为n 维搜索空间中的粒子Xi =(xi1 , xi2 , …, xin ), 定义粒子的适应值函数为

设 Fibest (x)表示第i 个粒子经历过的最好适应值 , Fgbest (x)表示全局经历过的最好适应值, 求解非线性方程组的过程等价于求解minF(x)的最优化问题 。PSO 求解非线性方程组的流程如下 :

Step1:输入方程组, 变量的数值范围和要求的精度

Step2:确定粒子的维数,根据 step1中的变量的范围 ,随机初始化粒子的位置和速度, 每个粒子的初始位置作为每个个体的个体极值 pbesti 。

Step3:对于第 i 个粒子,利用式计算粒子的适应值Fi (x),若其当前适应值小于它的历史最好适应值,则用当前粒子位置替换该粒子的个体极值 ,该粒子的适应值替换该粒子的最好适应值;若其当前适应值小于全局经历过的最好适应值, 则用当前粒子的位置替换全局极值 ,该粒子的适应值替换全局最好适应值。

Step4:根据(2)、(3)式更新粒子的速度和位置。

Step5:如果粒子的适应值达到预先设定的精度,则停止迭代 ,否则返回步骤 step3

V(n+1) = V(n) + c1*rand2*(pi – x(n)) + c2* rand1*(qi – x(n)) ------（2）

X(n+1) = x(n) + V(n+1) ------（3）

 MATLAB Code 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

%% fitness.m 文件
function F=fitness(x)
F = 0;
    f1 = x(1).^2+2*x(2).^2-0.3*cos(3*pi*x(1)) - 0.4*cos(4*pi*x(2)) + 0.7;
    f2 = 4*x(1).*exp(-x(1).^2-x(2).^2);
F = abs(f1)+ abs(f2);
end

------------------------------------------------------
%% PSO.m 文件
function[xm,fv] = PSO(fitness,N,c1,c2,w,M,D,a,b)
% c1,c2：学习因子
% w：惯性权重
% M：最大迭代次数
% D：搜索空间维数
% N：初始化群体个体数目

% 初始化种群的个体（可以在这里限定位置和速度的范围）
format long;
for i = 1:N
    for j=1:D
        x(i,j)  = unifrnd (a,b) ;% 随机初始化位置
        %unifrnd (a,b) ; 产生一个[a,b]的均匀随机数
        v(i,j) = randn; % 随即初始化速度       rand是0-1的均匀分布，randn是均值为0方差为1的正态分布；
    end
end

% 先计算各个粒子的适应度，并初始化pi和pg
for i=1:N
    p(i) = fitness(x(i,:));
    y(i,:) = x(i,:);
end
pg = x(N,:);  % pg为全局最优
for i=1:(N-1)
    if(fitness(x(i,:))
        pg = x(i,:);
    end
end

% 进入主要循环，按照公式依次迭代，直到满足精度要求
for t=1:M
    for i=1:N  % 更新速度、位移
        v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:));
        x(i,:)=x(i,:)+v(i,:);
        if fitness(x(i,:)) < p(i)
            p(i)=fitness(x(i,:));
            y(i,:)=x(i,:);
        end
        if p(i) < fitness(pg)
            pg=y(i,:);
        end

    end
    Pbest(t)=fitness(pg);

end

% 输出结果
disp(['PSO粒子群算法目标函数取最小值时的自变量：' num2str(pg) ]);
xm=pg';
disp(['PSO粒子群算法目标函数的最小值为：' num2str( Pbest(t))] );
fv=fitness(pg);

plot(1:t, Pbest);

---------------------------------------------
%% main.m 文件
clear
clc

[X ,Y] = meshgrid(-3:0.1:3, -3:0.1:3);
L = length(X);

for i = 1:L
    for j = 1:L
        Z(i,j) = fitness([X(i,j),Y(i,j)]);
    end
end

%作出适应度函数的图像  可见有两个零点  即方有两个解
surf(X,Y,Z)
xlabel('X(1)');
ylabel('X(2)');
zlabel('fitess(x(1),x(2))');

%利用粒子群算法计算方程组   可知方程组有一个解
%指定搜索范围 -1 - 1 之间  使用的是绝对值和作为fitness
figure
[xm1,fv1] = PSO(@fitness,100,2,2,0.8,40,2,-1,1);

legend('粒子群算法');
xlabel('迭代次数');ylabel('适应度');
title('粒子群算法求解方程组');

阅读┊ 收藏 ┊ 喜欢 ▼ ┊打印┊举报/Report

前一篇：A*寻路算法的简单实现

新浪BLOG意见反馈留言板　欢迎批评指正