C++ implementation of MOEA\D algorithm, the code is simple and logically clear

Article directory

• Preface
• 1. Weight matrix and neighborhood?
• 2. Aggregation function
• 3. Complete source code
• Summarize

Foreword

I have been studying multi-objective optimization recently. I was confused when I first read the MOEA\D paper, so I decided to take a look at the author’s source code. Because the author’s C source code has too many packages, the code is not intuitive enough. There are also many C++ implementations of MOEA\D on the Internet. Few, so I decided to write a C++ version of MOEA\D. This version is simple and direct, and clearly breaks down each step of the algorithm. I hope it can help people like me who are new to the MOEA\D algorithm and are a little confused.

1. Weight matrix and neighborhood?

When I first read the paper, I didn’t understand what the weight matrix and neighborhood were at all. The writing in the paper was also very unclear. After reading the source code, some other blogs and papers, I realized that this thing is actually not complicated at all.

The weight matrix actually divides the multi-objective into multiple single-objectives and the weight of each single-objective. This weight is only used for whether to retain when generating new individual solutions during population evolution, and the fitness of each individual for each single-objective problem. The value doesn’t matter, that’s the point! ! ! This is also the reason why I was most confused at the beginning, because this weight is only used in the last step to update the individual’s use of the aggregation function, and is not used in the evaluation of function fitness. What does it mean? Look at the code

//Initialization of target weight vector matrix
void init_weights()
{<!-- -->
for (int i = 0; i <= H; i + + )
{<!-- -->
weights[i][0] = i*1.0 / H;
weights[i][1] = 1 - weights[i][0];
}
}

The above is a uniform setting of weights. It is as simple as this. This H is customized. For a two-objective problem, each individual only needs two weights. The population size can be determined through H and the number of targets. There is a formula for this. You can read the paper

The test function is the one below. Look, this test function has nothing to do with the weight. It only matters when it comes to the aggregation function.

//Test function ZDT1
void evaluate_ind(individual* ind)
{<!-- -->
ind->fvalue[0] = ind->xreal[0];
double g = 1, sum = 0;
for (int i = 1; i < Dimension; i + + )
{<!-- -->
sum + = ind->xreal[i];
}
sum = 9 * (sum / (Dimension - 1));
g + = sum;
ind->fvalue[1] = g * (1 - sqrt(ind->xreal[0] / g));
evaluations + + ;
}

Neighborhood is easy to understand. In fact, it is the weight matrix you set, and then set a neighborhood size T. The first T elements with the closest Euclidean distance between each weight matrix are the neighborhoods.

//Neighborhood matrix
void init_neighbors()
{<!-- -->
double distance[pop_size][pop_size];
memset(distance, 0, sizeof(double));
double sum;
for (int i = 0; i < pop_size; i + + )
{<!-- -->
for (int j = i + 1; j < pop_size; j + + )
{<!-- -->
sum = 0;
for (int t = 0; t < Target; t + + )
{<!-- -->
sum + = (weights[i][t] - weights[j][t]) * (weights[i][t] - weights[j][t]);
}
distance[i][j] = sum;
distance[j][i] = sum;
}
}


for (int i = 0; i < pop_size; i + + )
{<!-- -->
for (int j = 0; j < T_size; j + + )
{<!-- -->
double min_temp = 0x3f3f3f3f;
int sp;
for (int k = 0; k < pop_size; k + + )
{<!-- -->
if (distance[i][k] < min_temp)
{<!-- -->
min_temp = distance[i][k];
sp = k;
}
}
neighbors[i][j] = sp;
distance[i][sp] = 0x3f3f3f3f;
}
}
}

2. Aggregation function

Use Chebyshev aggregation. This is the last step to update the solution. Only at this time are the weight matrix and Chebyshev aggregation function used. The other steps are the most basic genetic or differential evolution. There is no difference. MOEA\D There are more things like weight matrices, neighborhoods, and aggregation functions.

//Update neighborhood solution
void update_neighbors(int i)
{<!-- -->
// Set the weight vector of 0 to a smaller value
for (int j = 0; j < T_size; j + + )
{<!-- -->
int k = neighbors[i][j];
for (int t = 0; t < Target; t + + )
{<!-- -->
if (weights[k][t] < 1e-6)
{<!-- -->
weights[k][t] = 0.00001;
}
}
}
// Chebyshev aggregation
double oldobj, newobj;
double max_temp = -0x3f3f3f3f, part[Target];
for (int j = 0; j < T_size; j + + )
{<!-- -->
int k = neighbors[i][j];
//New individual aggregate function value
for (int t = 0; t < Target; t + + )
{<!-- -->
part[t] = abs(new_ind.fvalue[t] - Z[t]);
if (max_temp < weights[k][t] * part[t])
{<!-- -->
max_temp = weights[k][t] * part[t];
}
}
newobj = max_temp;
max_temp = -0x3f3f3f3f;
// Neighborhood individual aggregation function value
for (int t = 0; t < Target; t + + )
{<!-- -->
part[t] = abs(parent_pop[k].fvalue[t] - Z[t]);
if (max_temp < weights[k][t] * part[t])
{<!-- -->
max_temp = weights[k][t] * part[t];
}
}
oldobj = max_temp;
if (newobj<oldobj)
{<!-- -->
parent_pop[k] = new_ind;
}
}
}

3. Complete source code

#include
#include
#include

using namespace std;

#define URAND ((double)rand()/((double)RAND_MAX + 1.0)) // [0.0,1.0) uniformly distributed random numbers

double rndreal(double low, double high)
{
double val;
val = URAND * (high - low) + low;
return val;
}

int rndint(int low, int high)
{
int res;
if (low >= high)
{
res = low;
}
else
{
res = low + (int)(URAND * (high - low + 1));
if (res > high)
{
res = high;
}
}
return res;
}

#define Dimension 30 // Decision variable dimension

#define H 100 // Custom positive integer, the population number N can be determined based on the positive integer and the problem dimension
#define Target 2 // Target quantity
#define pop_size 101 //The population size is determined by H and the number of targets
#define T_size 20 //The number of weight vectors in the field

#define max_evaluations 150000 // Maximum number of evaluations



typedef struct
{
double xreal[Dimension]; // value of each dimension
double fvalue[2]; // objective function value
}individual;

// Weight matrix
double weights[pop_size][Target];
// Area
int neighbors[pop_size][T_size];
//Reference point (minimum value of each sub-objective function)
double Z[Target];


// Parent population and new individual
individual parent_pop[pop_size];
individual new_ind;




double lowerBound[Dimension]; // lower bound
double upperBound[Dimension]; // upper bound


long int evaluations = 0; // Current number of evaluations
int gen_no = 0; // Current iteration number


// Target weight vector matrix initialization
void init_weights()
{
for (int i = 0; i <= H; i + + )
{
weights[i][0] = i*1.0 / H;
weights[i][1] = 1 - weights[i][0];
}
}

//neighborhood matrix
void init_neighbors()
{
double distance[pop_size][pop_size];
memset(distance, 0, sizeof(double));
double sum;
for (int i = 0; i < pop_size; i + + )
{
for (int j = i + 1; j < pop_size; j + + )
{
sum = 0;
for (int t = 0; t < Target; t + + )
{
sum + = (weights[i][t] - weights[j][t]) * (weights[i][t] - weights[j][t]);
}
distance[i][j] = sum;
distance[j][i] = sum;
}
}


for (int i = 0; i < pop_size; i + + )
{
for (int j = 0; j < T_size; j + + )
{
double min_temp = 0x3f3f3f3f;
int sp;
for (int k = 0; k < pop_size; k + + )
{
if (distance[i][k] < min_temp)
{
min_temp = distance[i][k];
sp = k;
}
}
neighbors[i][j] = sp;
distance[i][sp] = 0x3f3f3f3f;
}
}

}

//Initialize decision boundary
void init_variables()
{
for (int i = 0; i < Dimension; i + + )
{
lowerBound[i] = 0.0;
upperBound[i] = 1.0;
}
}

// Random initialization of population
void init_pop_uniform()
{
double low, up;

for (int i = 0; i < pop_size; i + + )
for (int j = 0; j < Dimension; j + + )
{
low = lowerBound[j];
up = upperBound[j];

parent_pop[i].xreal[j] = rndreal(low, up);
}
}

// Get reference point
void get_Z(individual* pop)
{
for (int i = 0; i < Target; i + + )
{
double min_temp = 0x3f3f3f3f;
for (int j = 0; j < pop_size; j + + )
{
if (pop[j].fvalue[i] < min_temp)
{
min_temp = pop[j].fvalue[i];
}
}
Z[i] = min_temp;
}
}


//Test function ZDT1
void evaluate_ind(individual* ind)
{<!-- -->
ind->fvalue[0] = ind->xreal[0];
double g = 1, sum = 0;
for (int i = 1; i < Dimension; i + + )
{<!-- -->
sum + = ind->xreal[i];
}
sum = 9 * (sum / (Dimension - 1));
g + = sum;
ind->fvalue[1] = g * (1 - sqrt(ind->xreal[0] / g));
evaluations + + ;
}

// Group fitness evaluation
void evaluate_pop(individual* pop, int size)
{
for (int i = 0; i < size; i + + )
{
evaluate_ind( & amp;pop[i]);
}
}

// Differential evolution
void evolution(int i)
{
double CR = 0.9, F = 0.4;
int r1, r2, r3;
do {
r1 = rndint(0, T_size-1);
r1 = neighbors[i][r1];
} while (r1 == i);
do {
r2 = rndint(0, T_size-1);
r2 = neighbors[i][r2];
} while (r2 == i || r2 == r1);
do {
r3 = rndint(0, T_size-1);
r3 = neighbors[i][r3];
} while (r3 == i || r3 == r2 || r3 == r1);

// Differential mutation
for (int j = 0; j < Dimension; j + + )
{
new_ind.xreal[j] = parent_pop[r1].xreal[j] +
F * (parent_pop[r2].xreal[j] - parent_pop[r3].xreal[j]);

// Boundary judgment
if (new_ind.xreal[j] < lowerBound[j])
{
new_ind.xreal[j] = lowerBound[j];
}
if (new_ind.xreal[j] > upperBound[j])
{
new_ind.xreal[j] = upperBound[j];
}
}

// crossover operator
int j_rand = rndint(0, Dimension - 1);
for (int j = 0; j < Dimension; j + + )
{
if (rndreal(0, 1) < CR || j == j_rand)
{
new_ind.xreal[j] = new_ind.xreal[j];
}
else
{
new_ind.xreal[j] = parent_pop[i].xreal[j];
}
}

}

// Gaussian mutation increases population diversity
void gaussian_mutation()
{
double pr = 1 / Dimension;
for (int i = 0; i < Dimension; i + + )
{
double rnd;
double U1 = URAND, U2 = URAND;
double Z = sqrt(-2 * log(U1)) * cos(2 * 3.1415926 * U2); // Normal distribution with mean 0 and variance 1
rnd = new_ind.xreal[i] + 1 / 20 * Z; // Normal distribution with mean 1 and variance 1, 20

double newparm = min(max(rnd, 0.0), 1.0);
if (URAND
new_ind.xreal[i] = newparm;
}
}
}

//Update neighborhood solution
void update_neighbors(int i)
{<!-- -->
// Set the weight vector of 0 to a smaller value
for (int j = 0; j < T_size; j + + )
{<!-- -->
int k = neighbors[i][j];
for (int t = 0; t < Target; t + + )
{<!-- -->
if (weights[k][t] < 1e-6)
{<!-- -->
weights[k][t] = 0.00001;
}
}
}
// Chebyshev aggregation
double oldobj, newobj;
double max_temp = -0x3f3f3f3f, part[Target];
for (int j = 0; j < T_size; j + + )
{<!-- -->
int k = neighbors[i][j];
//New individual aggregate function value
for (int t = 0; t < Target; t + + )
{<!-- -->
part[t] = abs(new_ind.fvalue[t] - Z[t]);
if (max_temp < weights[k][t] * part[t])
{<!-- -->
max_temp = weights[k][t] * part[t];
}
}
newobj = max_temp;
max_temp = -0x3f3f3f3f;
// Neighborhood individual aggregation function value
for (int t = 0; t < Target; t + + )
{<!-- -->
part[t] = abs(parent_pop[k].fvalue[t] - Z[t]);
if (max_temp < weights[k][t] * part[t])
{<!-- -->
max_temp = weights[k][t] * part[t];
}
}
oldobj = max_temp;
if (newobj<oldobj)
{<!-- -->
parent_pop[k] = new_ind;
}
}
}

void moead_method()
{
for (int i = 0; i < pop_size; i + + )
{
evolution(i); // Differential evolution
gaussian_mutation(); // Gaussian mutation
evaluate_ind( & amp;new_ind); // Update fitness
// Update Z
for (int j = 0; j < Target; j + + )
{
if (new_ind.fvalue[j] < Z[j])
{
Z[j] = new_ind.fvalue[j];
}
}
update_neighbors(i);
}
}

int main()
{
srand((unsigned)time(NULL));
\t
init_weights();
init_neighbors();

init_variables();
init_pop_uniform();
evaluate_pop(parent_pop, pop_size);

get_Z(parent_pop);

gen_no = 1;
while (evaluations <= max_evaluations & amp; & amp; gen_no <= 500)
{
moead_method();
cout << "Current iteration number:" << gen_no << " " << "Current solution:" << endl;
for (int i = 0; i < pop_size; i + + )
{
cout << "F1: " << parent_pop[i].fvalue[0] << " F2: " << parent_pop[i].fvalue[1] << endl;
}

gen_no + + ;
}

return 0;
}

Summary

The code can run perfectly. I didn't explain the principle of the MOEA\D algorithm very clearly, but the code logic I wrote is still very clear. Regarding the principle, you can read the original paper and other blogs, but you can use mine for the code. .