/**
* PSO algorithm
*/
// f1 - Best: x0 = -3 e f(x0) = 0
// f2 - Best: x0 = 2,5 e f(x0) = -2,25
// f3 - Best: x0 = -1,414213562, x1 = 0 e f(x0, x1) = -1,414213562
// f4 - Best: x0 = 1, x1 = 1 e f(x0, x1) = 0
// f5 - Best: x0 = 1, x1 = 1 e f(x0, x1) = 1
// G4 - Best: -30665.53867178332
// G6 - Best: -6961.81387558015
// G12 - Best: -1
// G15 - Best: 961.715022289961
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N_VAR 5 // f1 e f2: 1, f3: 2, f4: 2, f5: 2, G4: 5, G6: 2, G12: 3, G15: 3
#define N_CON 6 // f3: 1, f5: 2, G4: 6, G6: 2, G12: 729, G15: 2
#define MAX_PAR 40
#define MAX_TIME 100
typedef struct best_particle_s {
double * position;
double fitness;
} best_particle_t;
typedef struct particle_s {
double * position;
double * velocity;
double fitness;
best_particle_t pBest;
} particle_t;
typedef struct swarm_s {
particle_t particles[MAX_PAR];
best_particle_t gBest;
} swarm_t;
/**
* Allocate function - allocates the swarm that will be used.
*/
swarm_t * allocate () {
int i;
swarm_t * swarm = malloc (sizeof(swarm_t));
swarm->gBest.position = malloc (N_VAR * sizeof(double));
for (i = 0; i < MAX_PAR; i++) {
swarm->particles[i].position = malloc (N_VAR * sizeof(double));
swarm->particles[i].velocity = malloc (N_VAR * sizeof(double));
swarm->particles[i].pBest.position = malloc (N_VAR * sizeof(double));
}
return swarm;
}
/**
* Boundary function - boundary condition.
*/
void boundary (double lowerBound[], double upperBound[]) {
int m;
// f1
//lowerBound[0] = -4;
//upperBound[0] = 0;
// f2
//lowerBound[0] = 1;
//upperBound[0] = 5;
// f3
/*lowerBound[0] = -2;
upperBound[0] = 0;
lowerBound[1] = 0;
upperBound[1] = 2;*/
// f4
/*lowerBound[0] = -1;
upperBound[0] = 2;
lowerBound[1] = -1;
upperBound[1] = 2;*/
// f5
/*lowerBound[0] = -2;
upperBound[0] = 5;
lowerBound[1] = -2;
upperBound[1] = 5;*/
// G4
lowerBound[0] = 78;
upperBound[0] = 102;
lowerBound[1] = 33;
upperBound[1] = 45;
for (m = 2; m < N_VAR; m++) {
lowerBound[m] = 27;
upperBound[m] = 45;
}
// G6
/*lowerBound[0] = 13;
upperBound[0] = 100;
lowerBound[1] = 0;
upperBound[1] = 100;*/
// G12
/*for (m = 0; m < N_VAR; m++) {
lowerBound[m] = 0;
upperBound[m] = 10;
}*/
// G15
/*for (m = 0; m < N_VAR; m++) {
lowerBound[m] = 0;
upperBound[m] = 10;
}*/
}
/**
* Initialize function - initialize the position and velocity randomly.
*/
void initialize (swarm_t * swarm, double lowerBound[], double upperBound[]) {
int i, j;
float rand1, rand2;
swarm->gBest.fitness = 99999999999;
for (i = 0; i < MAX_PAR; i++) {
swarm->particles[i].fitness = 0;
swarm->particles[i].pBest.fitness = 9999999999;
for (j = 0; j < N_VAR; j++) {
rand1 = (rand() % 10000);
rand1 /= 10000;
swarm->particles[i].position[j] = lowerBound[j] + rand1 * (upperBound[j] - lowerBound[j]);
//printf(" %d = Position: %.15lf Velocity: ", j, swarm->particles[i].position[j]);
rand2 = (rand() % 10000);
rand2 /= 10000;
swarm->particles[i].velocity[j] = lowerBound[j] + rand2 * (upperBound[j] - lowerBound[j]);
//printf("%.15lf", swarm->particles[i].velocity[j]);
}
//printf("\n");
}
}
/**
* Fitness function - executes the objective function.
*/
void fitness (swarm_t * swarm) {
int i, j, m, cont1, cont2, cont3;
double x[N_VAR], g[N_CON], function = 0;
double penaltyCoef;
for (i = 0; i < MAX_PAR; i++) {
for (j = 0; j < N_VAR; j++) {
x[j] = swarm->particles[i].position[j];
}
// Function
// f1
//function = x[0] * x[0] + 6 * x[0] + 9;
// f2
//function = x[0] * x[0] - 5 * x[0] + 4;
// f3
/*function = x[0] + x[1];
// Constraints
penaltyCoef = 0.35;
g[0] = -2 + x[0] * x[0] + x[1] * x[1];
if (g[0] > 0) {
function = function + penaltyCoef * g[0] * g[0];
}*/
// f4
//function = 100 * pow((x[1] - x[0]), 2) + pow((1 - x[0]), 2);
// f5
/*function = pow((x[0] - 2), 2) + pow((x[1] - 1), 2);
// Constraints
penaltyCoef = 100000;
g[0] = x[0] * x[0] - x[1];
if (g[0] > 0) {
function = function + penaltyCoef * g[0] * g[0];
}
g[1] = x[0] + x[1] - 2;
if (g[1] > 0) {
function = function + penaltyCoef * g[1] * g[1];
}*/
// G4
function = 5.3578547 * x[2] * x[2] + 0.8356891 * x[0]* x[4] + 37.293239 * x[0] - 40792.141;
// Constraints
penaltyCoef = 1000000;
g[0] = 85.334407 + 0.0056858 * x[1] * x[4] + 0.0006262 * x[0] * x[3] - 0.0022053 * x[2] * x[4] - 92;
g[1] = -85.334407 - 0.0056858 * x[1] * x[4] - 0.0006262 * x[0] * x[3] + 0.0022053 * x[2] * x[4];
g[2] = 80.51249 + 0.0071317 * x[1] * x[4] + 0.0029955 * x[0] * x[1] + 0.0021813 * x[2] * x[2] - 110;
g[3] = -80.51249 - 0.0071317 * x[1] * x[4] - 0.0029955 * x[0] * x[1] - 0.0021813 * x[2] * x[2] + 90;
g[4] = 9.300961 + 0.0047026 * x[2] * x[4] + 0.0012547 * x[0] * x[2] + 0.0019085 * x[2] * x[3] - 25;
g[5] = -9.300961 - 0.0047026 * x[2] * x[4] - 0.0012547 * x[0] * x[2] - 0.0019085 * x[2] * x[3] + 20;
for (m = 0; m < N_CON; m++) {
if (g[m] > 0)
function = function + penaltyCoef * g[m] * g[m];
}
// G6
/*x[0] =14.09500000000000064;
x[1] =0.8429607892154795668;*/
/*function = pow((x[0] - 10), 3) + pow((x[1] - 20), 3);
// Constraints
penaltyCoef = 500;
g[0] = -pow((x[0] - 5), 2) - pow((x[1] - 5), 2) + 100;
//printf("g0 = %e\t", g[0]);
if (g[0] > 0)
function = function + penaltyCoef * g[0] * g[0];
g[1] = pow((x[0] - 6), 2) + pow((x[1] - 5), 2) - 82.81;
//printf("g1 = %e\n", g[1]);
if (g[1] > 0)
function = function + penaltyCoef * g[1] * g[1];*/
// G12
/*function = -(100 - pow((x[0] - 5), 2) - pow((x[1] - 5), 2) - pow((x[2] - 5), 2)) / 100;
// Constraints
penaltyCoef = 0.0004;
for (m = 0; m < N_CON; m++) {
for (cont1 = 0; cont1 < 9; cont1++)
for (cont2 = 0; cont2 < 9; cont2++)
for (cont3 = 0; cont3 < 9; cont3++)
g[m] = -pow((x[0] - cont1), 2) - pow((x[1] - cont2), 2) - pow((x[2] - cont3), 2) + 0.0625;
}
for (m = 0; m < N_CON; m++) {
if (g[m] > 0)
function = function + penaltyCoef * g[m] * g[m];
}*/
// G15
/*function = 1000 - (pow(x[0], 2)) - (2 * pow(x[1], 2)) - (pow(x[2], 2)) - (x[0] * x[1]) - (x[0] * x[2]);
// Constraints
penaltyCoef = 0.7;
g[0] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] - 25.0;
if (g[0] > 0)
function = function + penaltyCoef * g[0] * g[0];
g[1] = 8.0 * x[0] + 14.0 * x[1] + 7.0 * x[2] - 56.0;
if (g[1] > 0)
function = function + penaltyCoef * g[1] * g[1];*/
/*for (j = 0; j < N_VAR; j++)
printf("x[%d] = %.15lf\t", j, x[j]);
printf("function = %.15lf\n", function);*/
swarm->particles[i].fitness = function;
}
}
/**
* Calc_best function - calculate pBest and gBest.
*/
void calc_best (swarm_t * swarm) {
int j, m;
for (j = 0; j < MAX_PAR; j++) {
if (swarm->particles[j].fitness < swarm->particles[j].pBest.fitness) {
swarm->particles[j].pBest.fitness = swarm->particles[j].fitness;
for (m = 0; m < N_VAR; m++) {
swarm->particles[j].pBest.position[m] = swarm->particles[j].position[m];
}
}
//printf("pBest: %.15lf\n", swarm->particles[j].pBest.fitness);
if (swarm->particles[j].pBest.fitness < swarm->gBest.fitness) {
swarm->gBest.fitness = swarm->particles[j].pBest.fitness;
for (m = 0; m < N_VAR; m++) {
swarm->gBest.position[m] = swarm->particles[j].pBest.position[m];
}
}
}
//printf(" \ngBest: %.15lf\n\n", swarm->gBest.fitness);
}
/**
* Update function - update velocity and position.
*/
void update (swarm_t * swarm, double lowerBound[], double upperBound[], double alfa) {
int j, m, c1, c2;
double rand1, rand2, newVelocity, newPosition;
c1 = c2 = 2;
for (j = 0; j < MAX_PAR; j++) {
for (m = 0; m < N_VAR; m++) {
rand1 = rand () % 10000;
rand1 /= 10000;
rand2 = rand () % 10000;
rand2 /= 10000;
newVelocity = alfa * swarm->particles[j].velocity[m] + c1 * rand1 * (swarm->particles[j].pBest.position[m] - swarm->particles[j].position[m]) + c2 * rand2 * (swarm->gBest.
