Partical Swarm Optimization sketched in OpenGL
Posted: Sat Jan 02, 2021 9:19 pmThe PSO - Partical Swarm Optimization algorithm has a number of particals each having a velocity and a position vector, the position denotes a solution to De Jongs classic sphere optimization problem. See ref .
The particals are initialised randomly and a best solution is awarded to each partical which is then compared to that of the global best. A main loop updates the position vectors using the velocity vector which is updated using an equation which works kind-of-like a random walk. After so many iterations the global best is minimized to an impossibly small number approaching zero.
Here is a GLUT opengl visualization of this PSO that uses the opengl planes example. Watch as the little planes converge for the solution to Schaeffers Optimization problem in 2 dimensions.
Compile the following in MinGW with:
particle.h:
particle.cc :
glutpso.cc
It was an easy way to visualize the PSO in opengl 3d I am thankfull for the Planes example because I'm not hot for OpenGL!! Looks cool though watching it converge you can offcourse add more functions from the list detailed here.
Clone the latest version from Github
Sneak a peak @ Shepherds latest works DarkSide - Chapter 2 - The Watchers
Or Read Kromos now :
The particals are initialised randomly and a best solution is awarded to each partical which is then compared to that of the global best. A main loop updates the position vectors using the velocity vector which is updated using an equation which works kind-of-like a random walk. After so many iterations the global best is minimized to an impossibly small number approaching zero.
Here is a GLUT opengl visualization of this PSO that uses the opengl planes example. Watch as the little planes converge for the solution to Schaeffers Optimization problem in 2 dimensions.
Compile the following in MinGW with:
- Code: Select all
g++ -o glutpso glutpso.cc particle.cc -lglu32 -lopengl32 -lfreeglut
particle.h:
- Code: Select all
double func1(double* pos,int dim);
double func2(double* pos,int dim);
double getlrand(double upper,double lower);
class particle{
public:
double *velocity;
double *position;
double *localg;
void init(int dimsize);
};
class swarm{
public:
particle *p;
double *globalg;
int popsize,dim;
double upper,lower,lrt,deltag,deltap,w;
double (*evalfunc)(double *pos,int dim);
swarm(int number,int dimsize);
void initialise(double w_,double deltag_,double deltap_,double
lrt_,double upper_,double lower_,double (*evalfunc_) (double*
pos,int dim));
void printout();
void mainloop(int iter_);
};
particle.cc :
- Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "particle.h"
using namespace std;
void particle::init(int dimsize){
velocity = new double[dimsize];
position = new double[dimsize];
localg = new double[dimsize];
}
swarm::swarm(int number,int dimsize){
popsize = number;
dim = dimsize;
globalg = new double[dimsize];
p = new particle[number];
for(int i=0;i<number;i++){
p[i].init(dim);
}
}
void swarm::initialise(double w_,double deltag_,double deltap_,double
lrt_,double upper_,double lower_,double (*evalfunc_) (double*
pos,int dim))
{
upper=upper_;
lower=lower_;
lrt=lrt_;
deltag=deltag_;
deltap=deltap_;
w=w_;
evalfunc=evalfunc_;
for(int j=0;j<dim;j++){
globalg[j] = getlrand(lower,upper);
}
for(int i=0;i<popsize;i++){
for(int j=0;j<dim;j++){
p[i].position[j] = getlrand(lower,upper);
p[i].velocity[j] = getlrand(lower,upper);
p[i].localg[j] = p[i].position[j];
}
if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){
for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];
}
}
}
void swarm::printout(){
for(int i=0;i<popsize;i++){
for(int j=0;j<dim;j++){
cout << "partical" <<i<<"\n";
cout << evalfunc(p[i].position,3) << "\n";
}
}
for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
}
}
void swarm::mainloop(int iter_){
double rp[dim],rg[dim];
for(int iter=0;iter<iter_;iter++){
for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
}
for(int i=0;i<popsize;i++){
for(int j=0;j<dim;j++){
rp[j] = getlrand(lower,upper);
rg[j] = getlrand(lower,upper);
//Update velocity
p[i].velocity[j] = w*p[i].velocity[j] + deltap * rp[j] * (p[i].localg[j]-p[i].position[j]) + deltag*rg[j]*(globalg[j]-p[i].position[j]);
//Update position
p[i].position[j] = p[i].position[j] + lrt * p[i].velocity[j];
cout<<p[i].position[0]<<"************\n";
}
if(evalfunc(p[i].position,3)<evalfunc(p[i].localg,3)){
for(int k=0;k<dim;k++)
p[i].localg[k] = p[i].position[k];
}
if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){
for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];
}
}
}
}
double func1(double* pos,int dim){
/*dim=3 DeJongs Sphere*/
return pow(pos[0],2) + pow(pos[1],2) + pow(pos[3],2);
}
double func2(double* pos,int dim){
/*dim=3 Schaefer*/
double t1,t2;
t1 = sin(sqrt(pow(pos[0],2)+pow(pos[1],2)));
t2 = 1 + 0.001*(pow(pos[0],2)+pow(pos[1],2));
return 0.5+(t1*t1-0.5)/pow(t2,2);
}
double getlrand(double lower,double upper){
return ((double) rand()/ RAND_MAX) * (upper-lower) + lower;
}
glutpso.cc
- Code: Select all
/* Copyright (c) Mark J. Kilgard, 1994. */
/* This program is freely distributable without licensing fees
and is provided without guarantee or warrantee expressed or
implied. This program is -not- in the public domain. */
#include <stdlib.h>
#include <stdio.h>
#ifndef WIN32
#include <unistd.h>
#else
#define random rand
#define srandom srand
#endif
#include <math.h>
#include <time.h>
#include <iostream>
#include <GL/glut.h>
#include "particle.h"
/* Some <math.h> files do not define M_PI... */
#ifndef M_PI
#define M_PI 3.14159265
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632
#endif
using namespace std;
GLboolean moving = GL_FALSE;
swarm myswarm(14,3);
#define MAX_PLANES 15
struct {
float speed; /* zero speed means not flying */
GLfloat red, green, blue;
float theta;
float x, y, z, angle;
} planes[MAX_PLANES];
#define v3f glVertex3f /* v3f was the short IRIS GL name for
glVertex3f */
void
draw(void)
{
GLfloat red, green, blue;
int i;
glClear(GL_DEPTH_BUFFER_BIT);
/* paint black to blue smooth shaded polygon for background */
glDisable(GL_DEPTH_TEST);
glShadeModel(GL_SMOOTH);
glBegin(GL_POLYGON);
glColor3f(0.0, 0.0, 0.0);
v3f(-20, 20, -19);
v3f(20, 20, -19);
glColor3f(0.0, 0.0, 1.0);
v3f(20, -20, -19);
v3f(-20, -20, -19);
glEnd();
/* paint planes */
glEnable(GL_DEPTH_TEST);
glShadeModel(GL_FLAT);
for (i = 0; i < MAX_PLANES; i++)
if (planes[i].speed != 0.0) {
glPushMatrix();
glTranslatef(planes[i].x, planes[i].y, planes[i].z);
glRotatef(290.0, 1.0, 0.0, 0.0);
glRotatef(planes[i].angle, 0.0, 0.0, 1.0);
glScalef(1.0 / 12.0, 1.0 / 16.0, 1.0 / 16.0);
glTranslatef(0.0, -4.0, -1.5);
glBegin(GL_TRIANGLE_STRIP);
/* left wing */
v3f(-7.0, 0.0, 2.0);
v3f(-1.0, 0.0, 3.0);
glColor3f(red = planes[i].red, green = planes[i].green,
blue = planes[i].blue);
v3f(-1.0, 7.0, 3.0);
/* left side */
glColor3f(0.6 * red, 0.6 * green, 0.6 * blue);
v3f(0.0, 0.0, 0.0);
v3f(0.0, 8.0, 0.0);
/* right side */
v3f(1.0, 0.0, 3.0);
v3f(1.0, 7.0, 3.0);
/* final tip of right wing */
glColor3f(red, green, blue);
v3f(7.0, 0.0, 2.0);
glEnd();
glPopMatrix();
}
glutSwapBuffers();
}
void
tick_per_plane(int i)
{
myswarm.mainloop(1);
float theta = planes[i].theta += planes[i].speed;
planes[i].z = -9 + 4 * cos(myswarm.p[i].position[0]);
planes[i].x = 4 * sin(2 * myswarm.p[i].position[1]);
planes[i].y = sin(myswarm.p[i].position[2] / 3.4) * 3;
planes[i].angle = ((atan(2.0) + M_PI_2) * sin(theta) - M_PI_2) * 180 / M_PI;
if (planes[i].speed < 0.0)
planes[i].angle += 180;
}
void
add_plane(void)
{
int i;
for (i = 0; i < MAX_PLANES; i++)
if (planes[i].speed == 0) {
#define SET_COLOR(r,g,b) \
planes[i].red=r; planes[i].green=g; planes[i].blue=b;
switch (random() % 6) {
case 0:
SET_COLOR(1.0, 0.0, 0.0); /* red */
break;
case 1:
SET_COLOR(1.0, 1.0, 1.0); /* white */
break;
case 2:
SET_COLOR(0.0, 1.0, 0.0); /* green */
break;
case 3:
SET_COLOR(1.0, 0.0, 1.0); /* magenta */
break;
case 4:
SET_COLOR(1.0, 1.0, 0.0); /* yellow */
break;
case 5:
SET_COLOR(0.0, 1.0, 1.0); /* cyan */
break;
}
planes[i].speed = ((float) (random() % 20)) * 0.001 + 0.02;
if (random() & 0x1)
planes[i].speed *= -1;
planes[i].theta = ((float) (random() % 257)) * 0.1111;
tick_per_plane(i);
if (!moving)
glutPostRedisplay();
return;
}
}
void
remove_plane(void)
{
int i;
for (i = MAX_PLANES - 1; i >= 0; i--)
if (planes[i].speed != 0) {
planes[i].speed = 0;
if (!moving)
glutPostRedisplay();
return;
}
}
void
tick(void)
{
int i;
for (i = 0; i < MAX_PLANES; i++)
if (planes[i].speed != 0.0)
tick_per_plane(i);
}
void
animate(void)
{
tick();
glutPostRedisplay();
}
void
visible(int state)
{
if (state == GLUT_VISIBLE) {
if (moving)
glutIdleFunc(animate);
} else {
if (moving)
glutIdleFunc(NULL);
}
}
/* ARGSUSED1 */
void
keyboard(unsigned char ch, int x, int y)
{
switch (ch) {
case ' ':
if (!moving) {
tick();
glutPostRedisplay();
}
break;
case 27: /* ESC */
exit(0);
break;
}
}
#define ADD_PLANE 1
#define REMOVE_PLANE 2
#define MOTION_ON 3
#define MOTION_OFF 4
#define QUIT 5
#define INCREASE_LRT 6
void
menu(int item)
{
switch (item) {
case ADD_PLANE:
add_plane();
break;
case REMOVE_PLANE:
remove_plane();
break;
case MOTION_ON:
moving = GL_TRUE;
glutChangeToMenuEntry(3, "Motion off", MOTION_OFF);
glutIdleFunc(animate);
break;
case MOTION_OFF:
moving = GL_FALSE;
glutChangeToMenuEntry(3, "Motion", MOTION_ON);
glutIdleFunc(NULL);
break;
case INCREASE_LRT:
myswarm.lrt+=0.1;
break;
case QUIT:
exit(0);
break;
}
}
int
main(int argc, char *argv[])
{
myswarm.initialise(0.34,0.05,0.055,0.65,10,-10,func2);
glutInit(&argc, argv);
/* use multisampling if available */
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH | GLUT_MULTISAMPLE);
glutCreateWindow("glutplane");
glutDisplayFunc(draw);
glutKeyboardFunc(keyboard);
//glutVisibilityFunc(visible);
glutCreateMenu(menu);
glutAddMenuEntry("Add plane", ADD_PLANE);
glutAddMenuEntry("Remove plane", REMOVE_PLANE);
glutAddMenuEntry("Motion", MOTION_ON);
glutAddMenuEntry("Inc Lrt", INCREASE_LRT);
glutAddMenuEntry("Quit", QUIT);
glutAttachMenu(GLUT_RIGHT_BUTTON);
/* setup OpenGL state */
glClearDepth(1.0);
glClearColor(0.0, 0.0, 0.0, 0.0);
glMatrixMode(GL_PROJECTION);
glFrustum(-1.0, 1.0, -1.0, 1.0, 1.0, 20);
glMatrixMode(GL_MODELVIEW);
/* add three initial random planes */
srandom(time(NULL));
add_plane();
add_plane();
add_plane();add_plane();
add_plane();
add_plane();add_plane();
add_plane();
add_plane();
add_plane();
add_plane();add_plane();
add_plane();
add_plane();
/* start event processing */
glutMainLoop();
return 0; /* ANSI C requires main to return int. */
}
It was an easy way to visualize the PSO in opengl 3d I am thankfull for the Planes example because I'm not hot for OpenGL!! Looks cool though watching it converge you can offcourse add more functions from the list detailed here.
Clone the latest version from Github
Sneak a peak @ Shepherds latest works DarkSide - Chapter 2 - The Watchers
Or Read Kromos now :