C++ 実行結果が想定通りにならない

実現したいこと

次に示す内容を実装したい。

[適応的実数値交叉AREXの提案と評価](https://www.jstage.jst.go.jp/article/tjsai/24/6/24_6_446/_pdf

発生している問題・分からないこと

実行結果として、0.00000に収束することを想定しているが、現在、0.2程度に収束してしまう。収束までは正常に数値が改善されていくため、原因がわからない。

エラーメッセージ

error
1エラーメッセージはありません

該当のソースコード

C++
1#include <cstdio>
2#include <cstdlib>
3#include <climits>
4#include <cfloat>
5#include <iostream>
6#include <fstream>
7#include <cmath>
8#include <unistd.h>
9#include <sys/time.h>
10#include <cassert>
11#include <algorithm>
12#include <random>
13
14#define dimension 20
15#define population_size 100
16#define ncalcs 2000
17#define sigmanum sqrt(1.0 / dimension)
18#define number_child_times 4
19#define maxvalue 5.0
20#define PI 3.141592
21#define c_alpha (1.0 / (5.0 * dimension))
22using namespace std;
23
24random_device rd;
25mt19937 mt(rd());
26uniform_real_distribution<double> dist(-maxvalue, maxvalue);
27uniform_real_distribution<double> distone(0.0, 1.0);
28normal_distribution<double> distn(0.0,sigmanum);
29
30typedef struct{
31    vector<double> edge = vector<double>(dimension);
32    double fitness;
33   
34} Indiv;
35
36vector<int> randomparent(const int pop){    //pop個、ランダムに親選択
37    vector<int> number(population_size);
38    vector<int> parent(pop);
39    //第1, 2の親の選択
40    for(int i = 0; i < population_size; ++i){
41        number[i] = i;
42    }
43    shuffle(number.begin(), number.end(), mt);
44    for(int i = 0; i < pop; ++i){
45        parent[i] = number[i];
46    }
47    return parent;
48}
49
50double evaluate(const Indiv& indiv){
51    double sum = 0;
52    //sphere
53    for(int i = 0 ; i < dimension; ++i){
54        sum += pow(indiv.edge[i], 2);
55    }
56    return sum;
57}
58
59vector<double> middle(const vector<Indiv>& group){
60    vector<double> m(dimension, 0.0);
61    for(int i = 0; i < dimension; ++i){
62        for(int j = 0; j < dimension+1; ++j){
63            m[i] += group[j].edge[i];
64        }
65        m[i] /= dimension+1.0;
66    }
67    return m;
68}
69
70double update_alpha(double abefore, const vector<Indiv>& upchild){
71    //更新前のα=abefore, 交叉中心ベクトルmiddle, 上位子個体upchild
72
73    //上位子個体の平均
74    vector<double> mid(dimension, 0.0);
75    mid = middle(upchild);
76
77    //ε　イプシロン
78    vector<vector<double>> epsilon(dimension+1, vector<double>(dimension));
79    for(int i = 0; i < dimension+1; ++i){
80        for(int j = 0; j < dimension; ++j){
81            epsilon[i][j] = upchild[i].edge[j] - mid[j];
82        }
83    }
84
85    //Lcdp, Lavg
86    double sigma1 = 0.0;
87    double sigma2 = 0.0;
88    for (int i = 0; i < dimension+1; ++i){
89        double term = 0.0;
90        for(int j = 0; j < dimension; ++j){
91            term += upchild[i].edge[j] / (dimension + 1.0);
92        }
93        sigma1 += pow(term, 2);
94        sigma2 += term;
95    }
96    sigma2 = pow(sigma2, 2) / (dimension + 1.0);
97    double L_cdp = (pow(abefore, 2) * dimension) * (sigma1 - sigma2);
98    double L_avg = pow(abefore, 2) * dimension * (dimension + 1.0);
99    double updatedalpha = abefore * sqrt((1.0 - c_alpha) + (c_alpha * (L_cdp / L_avg)));
100
101    // //αを更新
102    // double updatedalpha = abefore * sqrt((1 - c_alpha) + c_alpha * (L_cdp / L_avg));
103    //printf("L_cdp = %lf, L_avg = %lf, alpha = %lf\n", L_cdp, L_avg, updatedalpha);
104    if (updatedalpha < 1.0){
105        updatedalpha = 1.0;
106    }
107    if (updatedalpha > 2.0){
108        updatedalpha = 1.0;
109    }
110    return updatedalpha;
111}
112
113vector<Indiv> arex(const vector<Indiv>& parent, double alpha, const vector<double>& middle){
114    //生成子個体
115    vector<Indiv> children(number_child_times*(dimension+1));
116    vector<Indiv> sortpa = parent;
117    //親個体を昇順ソート
118    sort(sortpa.begin(), sortpa.end(), [](const Indiv& a, const Indiv& b){
119        return a.fitness < b.fitness;
120    });
121    //加重平均
122    vector<double> center(dimension, 0.0);
123    double sumone;
124    for(int i = 0; i < dimension; ++i){
125        sumone = 0;
126        for(int j = 0; j < dimension+1; ++j){
127            center[i] += sortpa[j].edge[i] * (2.0 * (dimension+2.0-(j+1)) / ((dimension+1.0) * (dimension+2.0)));
128            sumone += (2.0 * (dimension+2.0-(j+1)) / ((dimension+1.0) * (dimension+2.0)));
129        }
130        //printf("1 = %lf\n", sumone);
131    }
132    //子個体生成
133    for(int k = 0; k < number_child_times * (dimension+1); ++k){
134        for(int i = 0; i < dimension+1; ++i){
135            Indiv c;
136            double epsilon = distn(mt) * (-2.0);
137            //double epsilon = distn(mt) *(-2.0) * cos(2 * PI * distone(mt));
138            for(int j = 0; j < dimension; ++j){
139                double diff = parent[i].edge[j] - middle[j];
140                c.edge[j] = center[j] + alpha * epsilon * diff;
141                c.edge[j] = max(-maxvalue, min(maxvalue, c.edge[j]));
142            }
143            c.fitness = evaluate(c);
144            children[k] = c;
145        }
146    }
147    return children;
148}
149
150//JGGモデル
151vector<Indiv> JGG(const vector<Indiv>& children){
152    vector<Indiv> selectchild(number_child_times*(dimension+1));
153    for(int i = 0; i < number_child_times*(dimension+1); ++i){
154        selectchild[i] = children[i];
155    }
156
157    //並び替えておく
158    sort(selectchild.begin(), selectchild.end(), [](const Indiv& a, const Indiv& b){
159        return a.fitness < b.fitness;
160    });
161    vector<Indiv> upperchild(dimension+1);
162    for(int i = 0; i < dimension+1; ++i){
163        upperchild[i] = selectchild[i];
164    }
165    return upperchild;
166}
167
168int main(){
169    clock_t start = clock();
170    // vector<Indiv> population(population_size);
171    //初期化
172    vector<Indiv> population(population_size);
173    for (int i = 0; i < population_size; ++i){
174        for(int j = 0; j < dimension; ++j){
175            population[i].edge[j] = dist(mt);
176            population[i].fitness = evaluate(population[i]);
177        }
178    }
179    printf("世代数 : 最良値 : 平均値\n");
180    double alpha = 1;
181
182    // FILE *fp;
183    // fp = fopen("allchild.csv", "w");
184    for(long i=0;i<ncalcs;i++){
185        //double best_fitness = population[0].fitness;
186        double sum_fitness = 0.0;
187        double best_fitness = population[0].fitness;
188
189        //最良値と平均値
190        for(int j = 0; j < population_size; ++j){
191            if(population[j].fitness < best_fitness){
192                best_fitness = population[j].fitness;
193            }
194            sum_fitness += population[j].fitness;
195        }
196
197        long idx = static_cast<long>(ncalcs / 2000.0);
198        if (i % idx == 0){
199            printf("%d\t%lf\t%lf\n", i, best_fitness, sum_fitness / population_size);  
200        }
201        //推移表示
202        // for(int j = 0; j <= 50; ++j){
203        //     long idx = static_cast<long>(j * ncalcs / 5000.0);
204        //     if (i == idx){
205        //         printf("%d\t%lf\t%lf\n", i, best_fitness, sum_fitness / population_size);
206        //         //fprintf(fp, "%lf, %lf\n", best_fitness, sum_fitness / population_size);
207        
208        //     }
209        // }
210
211        //親選択
212        vector<int> parentc = randomparent(dimension+1);        //選択された親個体のインデックス
213        vector<Indiv> parents(dimension+1);     //選択された親個体
214        vector<Indiv> upper(dimension+1);   //生成子個体の上位, 上位順
215
216        //親群parent 作成 dimension+1個の個体
217        for(int j = 0; j < dimension+1; ++j){
218            parents[j] = population[parentc[j]];
219        }
220
221        //選択された親個体の平均
222        vector<double> midpa(dimension);
223        midpa = middle(parents);
224
225        //交叉
226        vector<Indiv> child(number_child_times*dimension);
227        child = arex(parents, alpha, midpa);
228
229        //printf("child.edge[0] = %lf\n", child[0].edge[0]);
230
231        //選択
232        upper = JGG(child);
233
234        //α更新
235        alpha = update_alpha(alpha, upper);
236        // if(i == 50){
237        //     printf("選択した親個体たち = ");
238        //     for(int k = 0; k < dimension+1; ++k){
239        //         printf("%lf, ", population[parentc[k]].fitness);
240        //     }
241        //     printf("\n代入する数値たち = ");
242        //     for(int k = 0; k < dimension+1; ++k){
243        //         printf("%lf, ", parents[k].fitness);
244        //     }
245        //     printf("\n生成個体の一部 = ");
246        //     for(int k = 0; k < dimension*number_child_times; ++k){
247        //         printf("%lf ", child[k].fitness);
248        //     }
249        //     printf("\n");
250        // }
251
252        //代入
253        for(int j = 0; j < dimension+1; ++j){
254            population[parentc[j]] = upper[j];
255        }
256    }
257    // fclose(fp);
258    clock_t end = clock();
259    double duration = static_cast<double>(end - start) / CLOCKS_PER_SEC;
260    printf("time : %lf\n", duration);
261}