4. Optimization (Python vs C++)
4.1. Differential Evolution Algorithm
In this section we implement a powerful optimization algorithm which is called Differential Evolution (Storn & Price, 1997). This algorithm is one of the most popular population-based algorithms that is very efficient in optimization of complex and multi-dimensional problems. We explain how it works and then implement it in Python:
Lets say we have a 2-dimensional problem meaning that there are 2 parameters that we want to find in a way that the cost function is minimized. First, we need to set up a list of population. The population list include vectors \(\vec{x_1} , ... , \vec{x_n}\). Each vector has two components corresponding to the unknowns of interest. The number of the population should be set by the user:
Then we need to define a trial vector \(\vec{u}\) for each member in the population list. The components of the trial vector are determined based on the below equation:
Note
In the above equation the \(F\) is called Mutation Factor that should be in the range of [0,2]
It should be noted that we need to set a limit for each parameter for the optimization. In the trial vector, if one of the components corresponding to a parameter is higher than the maximum limit for that parameter, it will be pushed to the maximum value of the limit. Similarly, if it is less than the minimum of the limit it should be pushed back to the minimum value of the limit that has been set for that parameter.
Now, we check in the trial vector to know if we should keep that particular component or replace it with the corresponding component from the member in the population list that the trial vector has been produced based on. In this regard, we should define a cross-over value (\(C.R\) )which is a number we should set between 0 and 1. A random number should be produced 2 times to determine each component of the trial vector. If the \(C.R\) is larger than the produced random number, then the component is kept inside the trial vector and if it is less than the produced random number it should be replaced with the component of the corresponding member in the population list:
In the above, the \(\vec{X}\) is called Child and corresponding member in the population list is called Parent. Now we should decide if we should keep the child or the parent. This should be done by replacing the component of the Child Vector and Parent Vector in the function that we want to minimize and evaluate the function. We should choose the one resulting in the smaller value of the function. Based on that, only one of the Child or Parent will survive in the population list.
After the above step, the population list has been updated. Finally, we should compare all members in the population list to see which one leads to a smaller value of the function that we want to minimize. Then we should select that particular member as the best vector containing the best optimum parameters.
Now we are done with the first generation. The same procedure should be repeated based on the updated population list.
4.2. Python Implementation
The python code is implemented for a simple function \(f = (parameter.1)^3+ (parameter.2)^3\) with the limit of \([-3,3]\) for both parameters to find the minimum of the function \(f\)
import random
import numpy as np
def urand():
return random.random()
def func(X):
#print ("aaa")
return pow(X[0],3)+pow(X[1],3)
up_bound = [3,3]
low_bound = [-3,-3]
dim= 2
NP = 10
cr = 0.90 # crossover probability
F5 = 0.50 # Scaling factor
Iteration_Number=20
parent_start = []
parent_val_start = []
def start():
for j in range(NP):
x = []
for k in range (dim):
x.append(random.uniform(low_bound[k],up_bound[k]))
parent_start.append(x)
parent_val_start.append(func(x))
start()
all_ind = []
for j in range (len(parent_start)):
all_ind.append(j)
def DE():
generation = 0
iteration = []
param1 = []
param2 = []
cost_func = []
for z in range(0, Iteration_Number):
generation=generation+1
for i in range (len(parent_start)):
trials = []
p1 = random.choice([e for k,e in zip(all_ind,all_ind) if k != i])
p2 = random.choice([e for k,e in zip(all_ind,all_ind) if k != i and k != p1])
p_target = random.choice([e for k,e in zip(all_ind,all_ind) if k != i and k != p1 and k != p2])
for k in range(0, dim):
if urand() <= cr:
trials.append((parent_start[p_target])[k] + F5 * ((parent_start[p1])[k] - (parent_start[p2])[k]))
else:
trials.append((parent_start[i])[k])
for s in range(dim):
if trials[s] > up_bound[s]:
trials[s]=up_bound[s]
elif trials[s] < low_bound[s]:
trials[s]=low_bound[s]
elif low_bound[s] <= trials[s] <= up_bound[s]:
trials[s]=trials[s]
child_val = func(trials)
if child_val <= parent_val_start[i]:
parent_start[i] = trials
parent_val_start[i] = child_val
soretd_val = sorted(parent_val_start,reverse=False)
ind = parent_val_start.index(soretd_val[0])
iteration.append(generation)
param1.append(parent_start[ind][0])
param2.append(parent_start[ind][1])
cost_func.append(soretd_val[0])
col_format = "{:<30}" * 4 + "\n"
with open("RESULTS-Python.txt", 'w') as of:
for x in zip(iteration, param1, param2,cost_func):
of.write(col_format.format(*x))
DE()
The results are saved in a text file (RESULTS-Python.txt). In this file the columns 1-4 corresponds to the number of generation, best value found for the first parameter, best value found for the second parameter and the minimum value of the objective function respectively.
1 -3 0.13575177113689574 -26.997498292598483
2 -3 -0.32372692704626427 -27.033926298141374
3 -2.8804034576214943 -3 -50.897912723155216
4 -2.8914790672098825 -3 -51.17464792180381
5 -2.91210613749475 -3 -51.69571468276654
6 -2.91210613749475 -3 -51.69571468276654
7 -2.923664426142372 -3 -51.99093876783297
8 -2.9752511297830933 -3.0 -53.33727790449086
9 -2.9752511297830933 -3.0 -53.33727790449086
10 -2.9778505435169755 -3.0 -53.40636919427237
11 -2.9792931692582787 -3.0 -53.444765647001205
12 -2.9832801955402197 -3.0 -53.55107657228403
13 -2.9832801955402197 -3.0 -53.55107657228403
14 -2.9909008960491836 -3.0 -53.75506858321391
15 -2.9979827647658714 -3.0 -53.945571263611825
16 -2.9979827647658714 -3.0 -53.945571263611825
17 -2.9979827647658714 -3.0 -53.945571263611825
18 -2.9979827647658714 -3.0 -53.945571263611825
19 -3 -3.0 -54.0
20 -3 -3.0 -54.0
4.3. C++ Implementation
#include <random>
#include <vector>
#include <fstream>
#include <iostream>
using namespace std;
float func(float x, float y){
return pow(x,3)+pow(y,3);
}
vector<float> up_bound {3,3};
vector<float> low_bound {-3,-3};
int dim= 2;
int NP = 10;
float cr = 0.90; // crossover probability
float F5 = 0.50; // Scaling factor
int Iteration_Number=20 ;
vector<vector<float>> parent_start;
vector<float> parent_val_start;
random_device rd;
default_random_engine eng(rd());
void start() {
for (int j = 0; j < NP; j++ ){
vector< float > x;
for (int k = 0; k < dim; k++ ){
uniform_real_distribution<> distr{low_bound[k], up_bound[k]};
x.push_back(distr(eng));
}
parent_start.push_back(x);
parent_val_start.push_back(func(x[0],x[1]));
}
};
int find_index (int size, float target, vector<float> val_array ){
int iterator;
for (int q = 0 ; q < size; q++){
if (val_array[q] == target){
iterator = q;
}
}
return iterator;
};
int generation = 0;
vector<int> iteration;
vector<float> param1;
vector<float> param2;
vector<float> cost_func;
uniform_real_distribution<> distr_cr {0, 1};
int main(){
srand(unsigned(time(NULL)));
start();
for (int z = 0; z < Iteration_Number; z++){
generation = generation+1;
for (int i = 0; i < parent_start.size(); i++){
vector<float> trials;
int p1;
int p2;
int p_target;
do {
p1 = rand() % parent_start.size();
p2 = rand() % parent_start.size();
p_target = rand() % parent_start.size();
//cout << "p1 :" << p1 << endl;
} while (p1 == i || p1 == p2 || p1 == p_target || p2 == p_target);
for (int k = 0; k < dim; k++){
float rand_0_1 = distr_cr (eng);
if (rand_0_1 <= cr){
trials.push_back((parent_start[p_target])[k] + F5 * ((parent_start[p1])[k] - (parent_start[p2])[k]));
}
else {
trials.push_back((parent_start[i])[k]);
}
};
for (int s = 0; s < dim; s++){
if (trials[s] > up_bound[s]){
trials[s]=up_bound[s];
}
else if (trials[s] < low_bound[s]){
trials[s]=low_bound[s];
}
else if (low_bound[s] <= trials[s] <= up_bound[s]){
trials[s]=trials[s];
}
};
float child_val = func(trials[0],trials[1]);
if (child_val <= parent_val_start[i]){
parent_start[i] = trials;
parent_val_start[i] = child_val;
}
};
vector<float> soretd_val = parent_val_start;
sort(soretd_val.begin(), soretd_val.end());
int ind = find_index (parent_val_start.size(),soretd_val[0],parent_val_start);
iteration.push_back(generation);
param1.push_back(parent_start[ind][0]);
param2.push_back(parent_start[ind][1]);
cost_func.push_back(soretd_val[0]);
}
string filename("Results-CPP.txt");
fstream file_out;
file_out.open(filename, std::ios::out);
for (int r=0 ; r < iteration.size(); r++){
file_out << iteration[r] << " " << param1[r]<< " " << param2[r]<< " " << cost_func[r]<< endl;
};
file_out.close();
return 0;
};
The results are saved in a text file (RESULTS-CPP.txt). In this file the columns 1-4 corresponds to the number of generation, best value found for the first parameter, best value found for the second parameter and the minimum value of the objective function respectively.
1 1.18459 -3 -25.3377
2 -1.54879 -2.90724 -28.2871
3 -1.83923 -2.95362 -31.9886
4 -1.83923 -2.95362 -31.9886
5 -2.30088 -2.9091 -36.8003
6 -2.56735 -2.93374 -42.1723
7 -2.56735 -2.93374 -42.1723
8 -2.56735 -2.93374 -42.1723
9 -2.6831 -2.90672 -43.8746
10 -2.73581 -2.88452 -44.4773
11 -2.81982 -2.93374 -47.6718
12 -2.90422 -2.90321 -48.9656
13 -2.91029 -2.93327 -49.8878
14 -2.91029 -2.96987 -50.8442
15 -3 -2.96256 -53.0016
16 -3 -2.96256 -53.0016
17 -3 -3 -54
18 -3 -3 -54
19 -3 -3 -54
20 -3 -3 -54
Like the previous Implementation in Python, the C++ code minimizes the objective function to -54 where the best value for each parameter was found equal to -3.
Note
The limits of the parameter that we want to find as well as the number of populations, \(C.R\) , \(F\) and number of iterations should be set in way a that it leads to best possible results. It usually needs some try and error and good understanding of the physics of the problem.