前言
最近想学一些玄学算法，比如说退火和粒子群。退火网上教程有很多，我学得比较舒服。但是粒子群...我在网上真的没有找到几篇适合我的，要么使用的语言不是c++，要么长得像工程文件，要么就是过于简练...但我在经过一番自己的努力摸索后，还是写出来了。所以在此我写一篇博客，希望对后来有兴趣学习的人有所帮助。
本次的代码是带可视化的，使用了hdc画图和一些控制台的操作。看不懂没关系，删了一样能跑。
算法针对的问题
粒子群算法可以解决在一个很大的范围内短时间内找到全局最优解的问题，但要求该问题是没有约束的（但其实可以将有约束的问题转为无约束问题）。就比如说我这里代码解决的样例问题：给定一个点，要求求出离这个点最近的点。我知道你想说什么：这不就是这个点本身吗！但是要知道，正式的问题可能还有不同的权值计算方式。对于某些问题，粒子群算法将十分有优势。
算法原理
算法的原型是鸟群觅食活动：有一群鸟在广阔的林场中觅食，它们都想要去到食物最多的地方。这些鸟有某种奇特的方式，可以共享当前找到的食物最多的地方（假设每只鸟都不想恰独食）。对于每只鸟，它下一步都会考虑两个东西：当前所有鸟的记忆中最好的地方和它的记忆中最好的地方。又因为它自己飞行有一个“惯性”，它下一步的飞行也会受到前一步飞行的影响。
该原理经过数学家们的研究、提炼，变成了一个数学模型，并得到了如下的公式：
（图来自百度百科）
该公式中,Vx是一个向量，表示了微粒x的速度(实际上我认为Vx前面还应该加上一个惯性系数ω)；‘2’是参数，实际运用中可以根据情况改动；rand与Rand都是一个随机数，范围是0~1；pbest是这个微粒之前所搜索到的最优解，pbest-x是一个由x指向pbest的向量；gbest是全局最优解
来解释一下上面那个公式。该公式大体可以分为三个部分：ω*Vx 、p*rand()*(pbest-x) 、q*Rand()*(gbest-x)。第一个部分就是该粒子前一次运动的惯性的影响；第二部分是该粒子记忆最优解对该粒子的“吸引力”；第三部分是全局最优解对该粒子的“吸引力”。
公式中的可调参数有三个：ω、p、q。ω减小会使得粒子更快地改变方向，但可能会因此失去寻找更优解的能力，陷入局部最优；p>q会使粒子更趋向于去到全局最优，但可能导致陷入局部最优；p<q会使粒子更趋向于去到个体最优，但可能导致解的最优性下降或最优解精度丢失；p、q的大小都会影响粒子收敛的速度，p、q值越大速度越快，但是精度低；值越小速度越慢，但精度高。
有了这些原理，我们就可以写出代码啦
代码
下面两个代码都加上了可视化
二维
#include<bits/stdc++.h>
#include<windows.h>
#include<conio.h>
using namespace std;
//选择窗口，新建画笔
HWND hWnd=GetConsoleWindow();
HDC hdc=GetDC(GetConsoleWindow());
HPEN hpenwhite = CreatePen( PS_SOLID, 1, RGB(255,255,255) );
HPEN hpenblack = CreatePen( PS_SOLID, 1, RGB(0,0,0) );
//
//下面是PSO算法主要内容
struct Vector{
double x,y;
double px,py;//px,py是原点为零的方向向量坐标
};
struct Point{
double x,y,val;
}glb;
double Bx,By;
double val(double x,double y){//val函数是对一个点权值的计算，是实际运用中最重要的地方
int len=sqrt((Bx-x)*(Bx-x)+(By-y)*(By-y));
return -len;
}
mt19937 mrand(time(0));
int rnd(int x){
return mrand()%x;
}
struct P{
Vector dir;
Point invb;
};
Vector operator +(Vector a,Vector b){//这边默认a是底子
Vector c;
c.x=a.x,c.y=a.y;
c.px=a.px+b.px,c.py=a.py+b.py;
return c;
}
Vector operator *(double a,Vector b){
Vector c;
c.x=b.x,c.y=b.y;
c.px=a*b.px,c.py=a*b.py;
return c;
}
vector<P> pts;
void Prework(){
for(int i=1;i<=300;i++){//调整i的范围可以改变离子的数量
int x=rnd(20000)-10000,y=rnd(20000)-10000;
pts.push_back({{x,y,rnd(2000)-1000,rnd(2000)-1000},{x,y,val(x,y)}});
}
}
//
//Print是画图函数
void Print(){
while(1){
SelectObject(hdc,hpenwhite);
Rectangle(hdc,0,0,1980,1080);
SelectObject(hdc,hpenblack);
cout<<glb.x<<' '<<glb.y<<endl;
for(int i=0;i<pts.size();i++){
MoveToEx(hdc,pts[i].dir.x+990,1080-(pts[i].dir.y+540),NULL);
LineTo(hdc,pts[i].dir.x+pts[i].dir.px+990,1080-(pts[i].dir.y+pts[i].dir.py+540));
}
Sleep(25);
}
}
int main(){
cin>>Bx>>By;
int cnt=0;
Prework();
thread t1(Print);//专门开一个线程搞画图
t1.detach();
glb.val=pts[0].invb.val;
for(int i=0;i<pts.size();i++)
if(glb.val<pts[i].invb.val)
glb=pts[i].invb;
while(cnt<=100){
cnt++;
for(int i=0;i<pts.size();i++){
pts[i].dir=0.9*pts[i].dir
+
(0.05*(double)((rnd(100)/100.0))*
(Vector){pts[i].dir.x,pts[i].dir.y,pts[i].invb.x-pts[i].dir.x,pts[i].invb.y-pts[i].dir.y})
+
(0.1*(double)((rnd(100)/100.0))*
(Vector){pts[i].dir.x,pts[i].dir.y,glb.x-pts[i].dir.x,glb.y-pts[i].dir.y});//调整这个式子里的三个参数可以改变粒子运动方式
int tempx=pts[i].dir.x,tempy=pts[i].dir.y;
pts[i].dir.x+=pts[i].dir.px,pts[i].dir.y+=pts[i].dir.py;
double v=val(pts[i].dir.x,pts[i].dir.y);
if(pts[i].invb.val<v)
pts[i].invb={pts[i].dir.x,pts[i].dir.y,v};
if(glb.val<pts[i].invb.val)
glb=pts[i].invb;
}
Sleep(25);
}
cout<<glb.x<<' '<<glb.y<<' '<<glb.val<<endl;
return 0;
}三维
三维表现的模板也是我自己瞎琢磨的，纯净的3D模板可以看看我的洛谷博客
#include<bits/stdc++.h>
#include<windows.h>
#include<conio.h>
#define KeyDown(VK_NONAME) ((GetAsyncKeyState(VK_NONAME)&0x8000) ? 1:0)
using namespace std;
double D=900,sen=0.2;
const double PI=acos(-1),g=9.8;
double cx,cy,cz;
double dirxz=0,diryz=0,mvdir;
double vy=0,ay=0,vx=0,ax=0;
bool ismove=false,isfly=false;
HANDLE Out=GetStdHandle(STD_OUTPUT_HANDLE),In=GetStdHandle(STD_INPUT_HANDLE);
HWND hWnd=GetConsoleWindow();
HDC hdc=GetDC(GetConsoleWindow());
HPEN hpenwhite = CreatePen( PS_SOLID, 1, RGB(255,255,255) );
HPEN hpenblack = CreatePen( PS_SOLID, 1, RGB(0,0,0) );
struct Vector{
double x,y,z;
double px,py,pz;
};
struct Point{
double x,y,z,val;
}glb;
double Bx,By,Bz;
double val(double x,double y,double z){
int len=sqrt((Bx-x)*(Bx-x)+(By-y)*(By-y)+(Bz-z)*(Bz-z));
return -len;
}
mt19937 mrand(time(0));
int rnd(int x){
return mrand()%x;
}
struct P{
Vector dir;
Point invb;
};
Vector operator +(Vector a,Vector b){//这边默认a是底子
Vector c;
c.x=a.x,c.y=a.y,c.z=a.z;
c.px=a.px+b.px,c.py=a.py+b.py,c.pz=a.pz+b.pz;
return c;
}
Vector operator *(double a,Vector b){
Vector c;
c.x=b.x,c.y=b.y,c.z=b.z;
c.px=a*b.px,c.py=a*b.py,c.pz=a*b.pz;
return c;
}
struct line{
double x1,y1,z1,x2,y2,z2;
};
vector<line> Linelist;
vector<P> pts;
void Prework(){
for(int i=1;i<=500;i++){
int x=rnd(200000)-100000,y=rnd(200000)-100000,z=rnd(200000)-100000;
pts.push_back({{x,y,z,rnd(200)-100,rnd(200)-100,rnd(200)-100},{x,y,z,val(x,y,z)}});
}
}
double Hu(double dir){
return dir*PI/180;
}
void Turn(double &x1,double &y1,double ox,double oy,double dir){
double tmpx=x1-ox,tmpy=y1-oy;
double cdir=cos(Hu(dir)),sdir=sin(Hu(dir));
x1=tmpx*cdir-tmpy*sdir;
y1=tmpx*sdir+tmpy*cdir;
}
void CursorControl()
{
while(1)
{
POINT p;
GetCursorPos(&p);
double dx=p.x-990,dy=p.y-540;
dirxz+=dx*sen;
diryz-=dy*sen;
SetCursorPos(990,540);
Sleep(10);
}
}
int cnt=0;
void TurnAndPrint(double x1,double y1,double z1,double x2,double y2,double z2){
Turn(x1,z1,cx,cz,dirxz); Turn(x2,z2,cx,cz,dirxz);
z1+=cz,z2+=cz;
Turn(z1,y1,cz,cy,-diryz); Turn(z2,y2,cz,cy,-diryz);
if(z1<0&&z2<0)
return;
if(z2<0){
x2=x1-z1*((x1-x2)/(z1-z2));
if(y1!=y2) y2=y2-z2*((y1-y2)/(z1-z2));
z2=0;
}
else if(z1<0){
x1=x1-z1*((x1-x2)/(z1-z2));
if(y1!=y2) y1=y1-z1*((y1-y2)/(z1-z2));
z1=0;
}
double X1,Y1,X2,Y2;
if(z1!=0) X1=x1*(D/z1)+989,Y1=-y1*(D/z1)+540;
else X1=x1*D+989,Y1=-y1*D+540;
if(z2!=0) X2=x2*(D/z2)+989,Y2=-y2*(D/z2)+540;
else X2=x2*D+989,Y2=-y2*D+540;
MoveToEx(hdc,X1,Y1,NULL);
LineTo(hdc,X2,Y2);
}
void print()
{
while(1)
{
SelectObject(hdc,hpenwhite);
Rectangle( hdc,0,0,1980,1080);
SelectObject(hdc,hpenblack);
for(int o=0;o<int(Linelist.size());o++){
SelectObject(hdc,hpenblack);
TurnAndPrint(Linelist[o].x1,Linelist[o].y1,Linelist[o].z1,Linelist[o].x2,Linelist[o].y2,Linelist[o].z2);
}
for(int o=0;o<int(pts.size());o++){
SelectObject(hdc,hpenblack);
TurnAndPrint(pts[o].dir.x,pts[o].dir.y,pts[o].dir.z,pts[o].dir.x+pts[o].dir.px,pts[o].dir.y+pts[o].dir.py,pts[o].dir.z+pts[o].dir.pz);
}
Sleep(25);
}
}
void PSO(){
glb.val=pts[0].invb.val;
for(int i=0;i<pts.size();i++)
if(glb.val<pts[i].invb.val)
glb=pts[i].invb;
while(cnt<=100){
cnt++;
for(int i=0;i<pts.size();i++){
pts[i].dir=0.9*pts[i].dir
+
(0.2*(double)((rnd(100)/100.0))*
(Vector){pts[i].dir.x,pts[i].dir.y,pts[i].dir.z,pts[i].invb.x-pts[i].dir.x,pts[i].invb.y-pts[i].dir.y,pts[i].invb.z-pts[i].dir.z})
+
(0.15*(double)((rnd(100)/100.0))*
(Vector){pts[i].dir.x,pts[i].dir.y,pts[i].dir.z,glb.x-pts[i].dir.x,glb.y-pts[i].dir.y,glb.z-pts[i].dir.z});
pts[i].dir.x+=pts[i].dir.px,pts[i].dir.y+=pts[i].dir.py,pts[i].dir.z+=pts[i].dir.pz;
double v=val(pts[i].dir.x,pts[i].dir.y,pts[i].dir.z);
if(pts[i].invb.val<v)
pts[i].invb={pts[i].dir.x,pts[i].dir.y,pts[i].dir.z,v};
if(glb.val<pts[i].invb.val)
glb=pts[i].invb;
}
Sleep(75);
}
}
void Move()
{
while(1){
if(ismove) ax=5;
else ax=-5;
if(vx<=30&&vx>=0) vx+=0.1*ax;
if(vx<0) vx=0;
if(vx>30) vx=30;
cx+=0.1*vx*sin(Hu(mvdir));
cz+=0.1*vx*cos(Hu(mvdir));
Sleep(10);
}
}
void fv(){
while(1){
if(!isfly){
cy+=0.1*vy;
Sleep(10);
}
}
}
void jump(){
while(1){
if(KeyDown(0x20)&&cy==0&&!isfly){
vy=50,ay=0;
Sleep(500);
}
else if(KeyDown(0x20)&&isfly){
cy+=5;
Sleep(10);
}
else if(KeyDown(0x10)&&isfly){
cy-=5;
Sleep(10);
}
}
}
void Gravity(){
while(1){
if(!isfly){
if(cy<0) vy=0,ay=0,cy=0;
vy+=0.1*ay;
ay-=g*0.1;
Sleep(10);
}
}
}
int main()
{
ShowWindow(hWnd,SW_SHOWMAXIMIZED);
DWORD mode;
GetConsoleMode(In,&mode);
mode &= ~ENABLE_QUICK_EDIT_MODE;
SetConsoleMode(In,mode);
cin>>Bx>>By>>Bz;
CONSOLE_CURSOR_INFO A;
GetConsoleCursorInfo(Out,&A);
A.bVisible=false;
SetConsoleCursorInfo(Out,&A);
thread t1(print);
thread t2(CursorControl);
thread t3(fv),t4(Gravity),t5(jump);
thread t6(Move);
for(int i=-1000;i<=1000;i+=50) Linelist.push_back({i,0,-1000,i,0,1000});
for(int i=-1000;i<=1000;i+=50) Linelist.push_back({-1000,0,i,1000,0,i});
t1.detach();
t2.detach();
t3.detach(),t4.detach(),t5.detach();
t6.detach();
Prework();
thread t7(PSO);t7.detach();
while(1)
{
if(KeyDown(0x41)||KeyDown(0x44)||KeyDown(0x57)||KeyDown(0x53)) ismove=true;
else ismove=false;
if(KeyDown(0x41)) mvdir=dirxz-90.0;
if(KeyDown(0x44)) mvdir=dirxz+90.0;
if(KeyDown(0x57)) mvdir=dirxz;
if(KeyDown(0x53)) mvdir=dirxz-180.0;
if(KeyDown(0x46)) isfly=!isfly,Sleep(100);
if(KeyDown(0x4F)) sen+=0.01;
if(KeyDown(0x50)) sen-=0.01;
if(KeyDown(0xDD)) D+=10;
if(KeyDown(0xDB)) D-=10;
if(KeyDown(0x1B)) exit(0);
Sleep(40);
}
return 0;
}
本文含有隐藏内容，请 开通VIP 后查看