从零开始的数模（十七）粒子群算法

目录

一、概念

1.1背景

​编辑

1.2基本形式 

1.3基本流程

1.4流程图

1.5比较

​编辑 1.6代码

 二、基于Python的粒子群算法

 2.1工具包的使用

​编辑 2.2例子说明

例子1有限制的粒子群

例子3python 实现粒子群 

三、基于MATLAB 的粒子群算法

3.1例题

---

一、概念

1.1背景

1.2基本形式 

1.3基本流程

1.4流程图

 看不懂了吧

上例题

1.5比较

 1.6代码

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2

from sko.PSO import PSO

pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)

import matplotlib.pyplot as plt

plt.plot(pso.gbest_y_hist)
plt.show()

 二、基于Python的粒子群算法

 2.1工具包的使用

sko.PSO 工具包讲解

​​​​​​​
pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)

参数说明

 2.2例子说明

例子1有限制的粒子群

第一步，定义问题

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2

第二步，做粒子群算法

from sko.PSO import PSO

pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)

 第三步，画出结果

import matplotlib.pyplot as plt

plt.plot(pso.gbest_y_hist)
plt.show()

 例子2无限制的粒子群

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2


# %% Do PSO
from sko.PSO import PSO

pso = PSO(func=demo_func, dim=3)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)

# %% Plot the result
import matplotlib.pyplot as plt

plt.plot(pso.gbest_y_hist)
plt.show()

 数据批量做粒子群优化

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2

from sko.PSO import PSO
import numpy as np
import pandas as pd



data={'lb':[[0,-1,0.5],[1,1,1],[2,3,4]],'ub':[[1,1,1],[2,2,2],[4,5,6]]}
data=pd.DataFrame(data)

print(data.shape[0])

def pso(lb,ub):
    pso = PSO(func=demo_func, dim=3, pop=40, max_iter=150, lb=lb, ub=ub, w=0.8, c1=0.5, c2=0.5)
    pso.run()
    print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)


for i in range(data.shape[0]):
    pso(data['lb'][i], data['ub'][i])

例子3python 实现粒子群 

这是一个简单算例，通过该算例简单感受下粒子群
求取 函数x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x) 的最大值

import numpy as np
import matplotlib.pyplot as plt
# 粒子（鸟）
class particle:
    def __init__(self):
        self.pos = 0  # 粒子当前位置
        self.speed = 0
        self.pbest = 0  # 粒子历史最好位置


class PSO:
    def __init__(self):
        self.w = 0.5  # 惯性因子
        self.c1 = 1  # 自我认知学习因子
        self.c2 = 1  # 社会认知学习因子
        self.gbest = 0  # 种群当前最好位置
        self.N = 20  # 种群中粒子数量
        self.POP = []  # 种群
        self.iter_N = 100  # 迭代次数

    # 适应度值计算函数
    def fitness(self, x):
        return x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x)

    # 找到全局最优解
    def g_best(self, pop):
        for bird in pop:
            if bird.fitness > self.fitness(self.gbest):
                self.gbest = bird.pos

    # 初始化种群
    def initPopulation(self, pop, N):
        for i in range(N):
            bird = particle()#初始化鸟
            bird.pos = np.random.uniform(-10, 10)#均匀分布
            bird.fitness = self.fitness(bird.pos)
            bird.pbest = bird.fitness
            pop.append(bird)

        # 找到种群中的最优位置
        self.g_best(pop)

    # 更新速度和位置
    def update(self, pop):
        for bird in pop:
            # 速度更新
            speed = self.w * bird.speed + self.c1 * np.random.random() * (
                bird.pbest - bird.pos) + self.c2 * np.random.random() * (
                self.gbest - bird.pos)

            # 位置更新
            pos = bird.pos + speed

            if -10 < pos < 10: # 必须在搜索空间内
                bird.pos = pos
                bird.speed = speed
                # 更新适应度
                bird.fitness = self.fitness(bird.pos)

                # 是否需要更新本粒子历史最好位置
                if bird.fitness > self.fitness(bird.pbest):
                    bird.pbest = bird.pos

    # 最终执行
    def implement(self):
        # 初始化种群
        self.initPopulation(self.POP, self.N)

        # 迭代
        for i in range(self.iter_N):
            # 更新速度和位置
            self.update(self.POP)
            # 更新种群中最好位置
            self.g_best(self.POP)


pso = PSO()
pso.implement()

best_x=0
best_y=0
for ind in pso.POP:
    #print("x=", ind.pos, "f(x)=", ind.fitness)
    if ind.fitness>best_y:
        best_y=ind.fitness
        best_x=ind.pos
print(best_y)
print(best_x)




x = np.linspace(-10, 10, 100000)


def fun(x):
    return x + 16 * np.sin(5 * x) + 10 * np.cos(4 * x)
y=fun(x)
plt.plot(x, y)

plt.scatter(best_x,best_y,c='r',label='best point')
plt.legend()
plt.show()

三、基于MATLAB 的粒子群算法

3.1例题

 步骤

​​​​​​​

%% 清空环境
clc;
clear all;
close all;
 
%% 参数初始化
%粒子群算法中的两个参数
c1 = 1.5; %学习因子1
c2 = 1.5; %学习因子2
 
T=100;   % 进化次数  
N=50;   %种群规模
D=2; %粒子堆数
 
Vmax=1; %速度最大数
Vmin=-1; %速度最小数
Xmax=3; %位置最大数
Xmin=-3; %位置最小数
[zbest,fitnesszbest]=pso_train(T,N,D,c1,c2,Vmax,Vmin,Xmax,Xmin);
%% 产生初始粒子和速度
X = rand(N,D)*(Xmax-Xmin)-Xmax;
V = rand(N,D)*(Vmax-Vmin)-Vmax;
for i=1:N
    %随机产生一个种群
%     pop(i,:)=rands(1,2)*(Xmax-Xmin)-Xmax;    %初始种群
%     V(i,:)=rands(1,2)*(Vmax-Vmin)-Vmax;  %初始化速度
%     计算适应度
    fit(i)=fun(X(i,:));   %染色体的适应度
end
 
%% 个体极值和群体极值
[bestfitness bestindex]=min(fit);
zbest=X(bestindex,:);   %全局最佳
gbest=X;    %个体最佳
fitnessgbest=fit;   %个体最佳适应度值
fitnesszbest=bestfitness;   %全局最佳适应度值
 
%% 迭代寻优
for i=1:T
    for j=1:N
        %速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - X(j,:)) + c2*rand*(zbest - X(j,:));
        V(j,find(V(j,:)>Vmax))=Vmax;
        V(j,find(V(j,:)<Vmin))=Vmin;
        %种群更新
        X(j,:)=X(j,:)+0.5*V(j,:);
        X(j,find(X(j,:)>Xmax))=Xmax;
        X(j,find(X(j,:)<Xmin))=Xmin;
        %适应度值
        fit(j)=fun(X(j,:)); 
    end
    for j=1:N
        
        %个体最优更新
        if fit(j) < fitnessgbest(j)
            gbest(j,:) = X(j,:);
            fitnessgbest(j) = fit(j);
        end
        
        %群体最优更新
        if fit(j) < fitnesszbest
            zbest = X(j,:);
            fitnesszbest = fit(j);
        end
    end 
    yy(i)=fitnesszbest;    
end
%% 结果分析
% 绘图
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
fprintf('最优个体是：%8.7f%8.7f\n',zbest(1),zbest(2))
fprintf('最小值是：%8.7f\n',fun(zbest))