Re: Cant reproduce results from Article. Close but not close enough

[Fabian Nick]

Hi,

did they use DEAP in the paper? If so, which version?

If they didn't use DEAP, I'd say your chances of reproducing the exact same result are fairly slim, cause the result will depend on the generation of random numbers (e.g. for the mutation), i.e. reproducibility across implementations is very hard if not impossible.

Cheers

Fabian

1768...@gmail.com schrieb am Mittwoch, 10. März 2021 um 13:07:12 UTC+1:

Hey everyone,

I have been trying to reproduce the results from a paper for quite some time and have gone through every equation and even re-derived their objective function. I have done everything I can think of and this has led me back to think that maybe I'm using the DEAP framework incorrectly.

I've gone through the guides,etc and still cant really see any issue.

I want to know whether someone can possibly have a look at my use. 

I've included the paper in which they provide "pseudo-code" and thorough explanations. I'm pretty sure I'm implementing everything correctly. 

I've also copied and pasted the relevant equations into my jupyter notebook and I've added a log file called LOG.txt which runs through each individual, mutation, etc and logs it.

https://github.com/Ian-Re/DE

Here is the link to the repo. Any questions, help or insight will be extremely appreciated.

Thanks in advance

[Fabian Nick]

Hi Ian,

maybe I should clarify one thing first: In my previous post I was talking about the EXACT reproducibility of the paper's results, which I think will be very hard to achieve (I will elaborate this point further down below)

You should of course be able to reproduce their results in terms of their general propositions (i.e. you should be able to at least get "similar" results).

Here's why I wouldn't expect you to get the exact same results / numbers with the DEAP framework:

As the paper states in Section 3.2, mutation introduces a random element into the DE algorithm. This happens mainly when creating the mutated vectors.

To achieve reproducibility within their own experiements, they would have used a random number generator with a fixed "seed", i.e. when they ran their experiments twice in a row, they would have got the exact same results. That is of course because the randomly generated vectors aren't truly random, but "pseudo-random".

However, the random number generator in DEAP will almost certainly give you a different - but still pseudo-random - behaviour.

And since they didn't use the DEAP framework, I don't think there would be a way to "align" those two sources of randomness.

With that being said, you should still be able to reproduce their results to a certain degree of accuracy. After all, at least if the problem is well posed, the randomness should not have too much of an impact (be aware of local optima, though! The mutation idea is trying to avoid those, but I don't think it can guarantee you don't run into them...)

May I ask how much your results differ from theirs?

Fabian

1768...@gmail.com schrieb am Samstag, 13. März 2021 um 10:43:47 UTC+1:

Hi Fabian,

Thank you for the reply.

They didn't use DEAP, however I'm not quite sure why it wouldn't be possible to reproduce the results. If you could maybe explain a little more I would highly appreciate it.

I attached the article just for reference. They use differential evolution and provide great detail of the mutation step as well as every other step. If you scroll through it, you'll see that they go in great depth on each step, therefore to me it looks like I should be able to reproduce their results. I've also implemented every step in the article as they provide it. It just seems that my issue might be using the DEAP framework incorrectly.

Any help would be highly appreciated.

Thanks again,

Ian

[Derek Tishler]

Are you monitoring the performance of the evolution?
https://deap.readthedocs.io/en/master/tutorials/basic/part3.html?highlight=plotting#some-plotting-sugar

This can help you monitor the fitness and std and other things visually to ensure you are not stopping too early, or having some odd convergence to a poor solution via an error in the operators or fitness functions. 

"I setup an extremely large amount of individuals just for the sake of debugging at this point: NP = 100000 and run the simulation for 100 generations."
This makes me worry, trying to brute force things in so few generations depends on many things(visualize to see if fitness is done convergine, may need Log scale to see). You also may want to print and explore your individuals over the first couple gens. If not set up correctly or if there is a mistake in an operator or something then the evo will just carry on with the problem so focusing on the final results may not get you anywhere as the problem can occur very early on.

Also in the paper they mention for selection they use equation 11, but i see in the notebook the use of:
toolbox.register("select", tools.selRandom, k=3)
This low selection pressure with such a large pop(in the paper do they mention using only 30?!) may be an issue depending on if i skimmed the paper correctly for this info.

On Friday, April 2, 2021 at 6:03:36 AM UTC-4 1768...@gmail.com wrote:

Fabian,

Sorry for the delayed reply. 

I understand that exact results are highly unlikely and not what I'm going for. I plan on obtaining a solution, there can be slight variation in the results obtained however the current output I'm obtaining is not usable.

I think this would be the first start to showing the difference. I calculate the curve from the obtained data:

I'm happy with the shape as it follows the same profile of an IV curve.

At x = 0, I should be getting the y value corresponding to the short circuit current of 3.45, this is a given though and I'll always obtain this.

at x = 17.4 i should be a the maximum power point, which will correspond to Vmp - 17.4 and Imp = 3.15. This is also relatively close as I enter the x value = Vmp = 17.4 and I obtain the y-value = 3.2065.

The issue comes close to the end. Here we can see that at an x value = 21.7 V, I get a y value of 1.1065. It should in fact be 0, so it is unacceptably far from the desired point.

Here is a printed list of each point I calculate:

Imp = 3.20651033729589  

0.0 [3.44651759191854] 

0.5 [3.44590833228358] 

1.0 [3.44529906826181] 

1.5 [3.44468979708421] 

2.0 [3.44408051423389] 

2.5 [3.44347121234283] 

3.0 [3.44286187939219] 

3.5 [3.44225249577662] 

4.0 [3.44164302951559] 

4.5 [3.44103342844195] 

5.0 [3.44042360745987] 

5.5 [3.43981342776003] 

6.0 [3.43920266291541] 

6.5 [3.43859094357738] 

7.0 [3.43797766726633] 

7.5 [3.43736185122895] 

8.0 [3.43674189243940] 

8.5 [3.43611517616955] 

9.0 [3.43547743764310] 

9.5 [3.43482172118358] 

10.0 [3.43413668349887] 

10.5 [3.43340382900209] 

11.0 [3.43259300807938] 

11.5 [3.43165509519443] 

12.0 [3.43051010207966] 

12.5 [3.42902793815189] 

13.0 [3.42699742354360] 

13.5 [3.42407678527865] 

14.0 [3.41971562203327] 

14.5 [3.41303459460073] 

15.0 [3.40264674431231] 

15.5 [3.38640872009538] 

16.0 [3.36111308275005] 

16.5 [3.32219241968673] 

17.0 [3.26361061163081] 

17.5 [3.17821464309342] 

18.0 [3.05875118377565] 

18.5 [2.89936567452677] 

19.0 [2.69690972299289] 

19.5 [2.45136226549314] 

20.0 [2.16527103919394] 

20.5 [1.84270794378545] 

21.0 [1.48828879489440] 

21.7 [1.10651681079745]  

The result is also relatively consistent. I think the lowest y-value I have obtained at this point is 0.7

The fitness of the selected individual is the following = fitness = (mpf('9.2722737821127055e-10'),)

I setup an extremely large amount of individuals just for the sake of debugging at this point: NP = 100000 and run the simulation for 100 generations.

Here is the output of the last generation in the logbook :

gen evals std min avg max  

99 100000 5.43212132428282e-5 9.27227378211271e-10 5.52332140756996e-5 0.00067915266363619 

Finally comparing the parameters they solve for I obtain the following at 25 degrees Celsius (converted to kelvin in the simulation of 298.15) :

io = 6.722687751087442e-10 

Rs = 0.8283842440483365

Rp = 819.8491640723421 

a = 1.1036259660321337

Vt = 0.9249338197129674 

r = 0.9796420772941139  

Where they obtain :

Rp  = 2340

Rs = 0.3

a = 1.37

The designated limits in the simulation for these parameters are :

Rs : [0.1, 1], Rp : [100, 3000], a : [1, 2].

I think from this it can be seen that I do not really obtain a solution and thus my implementation of their algorithm is failing.

I have however debugged every step of my implementation and ensured that I am correctly implementing each single step of their algorithm. 

Within the github link I posted there is a jupyter notebook file which you can see how I went through the process of listing their algorithm and showing my implementation:

Therefore I'm not sure whether I'm implementing the deap framework incorrectly? 

Any help will really be appreciated.

Thanks again Fabian

[Derek Tishler]

The convergence is pretty fast in the plot(not a bad thing in itself but is that expected? The paper didn't share their fitness curves but you could be rapidly falling into a local minima).

I am not quite understanding how the paper is intending to performing the 'perturbation' evolution based on the last bit of code you posted: is this intending to perform everything(mutation AND crossover AND 'selection') on an individual case while looping over the pop? The paper suggests this process but the code above appears to be changing the population in-place as it loops through the pop(while still taking from the pop via toolbox.select). Will this further encourage rapid convergence as the 'better' individuals in the early gens are continually selected-for-mutation and reinserted? 'Thus, if the new trial vector acquires a lower value of the objective function, it swaps the corresponding target vector in the next generation' (after eq 11).  I can be way off on this as its a quick thought inline with more standard evolutions vs what the paper is doing.

Can see what I mean with comparing operations(mutation & crossover in varAnd, then population wide selection in a separate step) in the example code which use aSimple: 
https://github.com/DEAP/deap/blob/d328fe6b68e7528b2d2d990bb2ab1ad1786e6f58/deap/algorithms.py#L85

One solution to a possible issue:
"such that the indices i, r1, r2 and r3 are distinct" and the line:
a, b, c = [toolbox.clone(ind) for ind in toolbox.select(pop)]
Wont this have a chance to grab duplicate objects when forming each trio? And since looping over the pop while modifying it; later gens may have 2 or 3 matching individuals in the trio easily. You could do something like this to ensure unique but random index at this step:
unique_idx = np.random.choice(np.random.permutation(len(pop)), 3)
a,b,c = [toolbox.clone(pop[i]) for i in unique_idx]

On Tuesday, April 6, 2021 at 4:35:41 PM UTC-4 1768...@gmail.com wrote:

Hi!

Thank you for the reply.

"Are you monitoring the performance of the evolution?

https://deap.readthedocs.io/en/master/tutorials/basic/part3.html?highlight=plotting#some-plotting-sugar

This can help you monitor the fitness and std and other things visually to ensure you are not stopping too early, or having some odd convergence to a poor solution via an error in the operators or fitness functions. "

I've changed my monitoring to display the graph as well. Here's the text output and the graph:

                        fitness              

          -----------------------------------

gen evals min               avg  max

0   30    0.107498520522443 +inf inf

1   30    0.103424940918559 +inf inf

2   30    0.0478141140837091 0.161194489331119 0.184258873366956

3   30    0.0138565027954866 0.151395257994144 0.181583428945098

4   30    0.0138565027954866 0.131337471673565 0.181456262531419

5   30    0.0138565027954866 0.114684379539387 0.18136454794196 

6   30    0.0138565027954866 0.112935809229708 0.18113715866916 

7   30    0.00160144612612353 0.0941891082772629 0.179872481615692

8   30    0.00160144612612353 0.0739546648145661 0.17287135502154 

9   30    0.00160144612612353 0.0729480083821712 0.17287135502154 

10  30    0.00160144612612353 0.064845833262041  0.17287135502154 

11  30    0.00160144612612353 0.0555559135551622 0.170680085524535

12  30    0.000898255922016561 0.0443383809110422 0.153808297234372

13  30    0.000898255922016561 0.0362329009588852 0.153808297234372

14  30    0.000898255922016561 0.0322831476948328 0.153808297234372

15  30    0.000898255922016561 0.0287909348642001 0.153808297234372

16  30    0.000898255922016561 0.0277770924732048 0.153808297234372

17  30    0.000898255922016561 0.0262199428179946 0.153808297234372

18  30    0.000898255922016561 0.0244565216698966 0.153808297234372

19  30    0.000898255922016561 0.0218977764875409 0.0791575270800119

20  30    0.000898255922016561 0.0215672432777107 0.0784961050489076

21  30    0.000898255922016561 0.0202560448455632 0.0784961050489076

22  30    0.000898255922016561 0.0176044888081345 0.0784575728440888

23  30    0.000898255922016561 0.0173417079952114 0.0784575728440888

24  30    0.000898255922016561 0.0173417079952114 0.0784575728440888

25  30    0.000898255922016561 0.0163371177808751 0.0784575728440888

26  30    0.000898255922016561 0.0130015769654964 0.0770156686403184

27  30    0.000898255922016561 0.0115369573789649 0.0770156686403184

28  30    0.000898255922016561 0.0113053759800573 0.0770156686403184

29  30    0.000898255922016561 0.0103192204883294 0.0557207433062525

30  30    0.000898255922016561 0.00919338086195996 0.0557207433062525

31  30    0.000898255922016561 0.0081287228141644  0.0557207433062525

32  30    0.000898255922016561 0.00742330888202007 0.0557207433062525

33  30    0.000898255922016561 0.00553372233251289 0.0180774526731676

34  30    6.70016465277312e-5  0.00548642427047371 0.0180774526731676

35  30    2.87010107429542e-5  0.00508708501345932 0.0180774526731676

36  30    2.87010107429542e-5  0.00489478943711212 0.0180774526731676

37  30    2.87010107429542e-5  0.00465394059122298 0.0180774526731676

38  30    2.87010107429542e-5  0.00463273471435232 0.0180774526731676

39  30    2.87010107429542e-5  0.00374453015303631 0.0180774526731676

40  30    2.87010107429542e-5  0.00370761172510751 0.0180774526731676

41  30    2.87010107429542e-5  0.00370134733190036 0.0180774526731676

42  30    2.87010107429542e-5  0.00362299394569295 0.0180774526731676

43  30    2.87010107429542e-5  0.0035756894675224  0.0180683969057817

44  30    2.87010107429542e-5  0.00346744904971329 0.0180683969057817

45  30    2.87010107429542e-5  0.00304807855927434 0.0180683969057817

46  30    2.87010107429542e-5  0.00269002847829674 0.0180683969057817

47  30    2.87010107429542e-5  0.00267245317927223 0.0180683969057817

48  30    2.87010107429542e-5  0.00266993252211907 0.0180683969057817

49  30    2.87010107429542e-5  0.00208644781620565 0.00775852843955674

50  30    2.87010107429542e-5  0.00192933578330644 0.00775852843955674

51  30    2.87010107429542e-5  0.0019238165298625  0.00775852843955674

52  30    2.87010107429542e-5  0.00172779190213034 0.00692456633944916

53  30    2.87010107429542e-5  0.00170348884804259 0.00619547471681664

54  30    2.87010107429542e-5  0.00170214472464467 0.00619547471681664

55  30    2.87010107429542e-5  0.00169818588580015 0.00619547471681664

56  30    2.87010107429542e-5  0.00165652520147283 0.00619547471681664

57  30    2.87010107429542e-5  0.00165652520147283 0.00619547471681664

58  30    2.87010107429542e-5  0.00165610385379056 0.00619547471681664

59  30    2.87010107429542e-5  0.00160902393767292 0.00619547471681664

60  30    2.87010107429542e-5  0.00160863543995544 0.00619547471681664

61  30    2.87010107429542e-5  0.00160863543995544 0.00619547471681664

62  30    2.87010107429542e-5  0.00160863543995544 0.00619547471681664

63  30    2.87010107429542e-5  0.00160863543995544 0.00619547471681664

64  30    2.87010107429542e-5  0.00159070969086693 0.00619547471681664

65  30    2.87010107429542e-5  0.00159070969086693 0.00619547471681664

66  30    2.87010107429542e-5  0.00152404259873593 0.00619547471681664

67  30    2.87010107429542e-5  0.00152404259873593 0.00619547471681664

68  30    2.87010107429542e-5  0.0014496884932955  0.00619547471681664

69  30    2.87010107429542e-5  0.00139899625170319 0.00619547471681664

70  30    2.87010107429542e-5  0.00139899625170319 0.00619547471681664

71  30    2.87010107429542e-5  0.00130855558909627 0.00619547471681664

72  30    2.87010107429542e-5  0.00128833673387474 0.00619547471681664

73  30    2.87010107429542e-5  0.00126357065607264 0.00619547471681664

74  30    2.87010107429542e-5  0.00122539942318871 0.00619547471681664

75  30    2.87010107429542e-5  0.00122539942318871 0.00619547471681664

76  30    1.90983323070104e-6  0.00116049611523825 0.00619547471681664

77  30    1.90983323070104e-6  0.000979161705418197 0.00340759734774973

78  30    1.90983323070104e-6  0.000913133787835428 0.00340759734774973

79  30    1.90983323070104e-6  0.000835082794640667 0.00241359143105091

80  30    1.90983323070104e-6  0.000834269606226184 0.00241359143105091

81  30    1.90983323070104e-6  0.000832929164699387 0.00241359143105091

82  30    1.90983323070104e-6  0.000832929164699387 0.00241359143105091

83  30    1.90983323070104e-6  0.000829783440046698 0.00241359143105091

84  30    1.90983323070104e-6  0.000829783440046698 0.00241359143105091

85  30    1.90983323070104e-6  0.000829299606437929 0.00241359143105091

86  30    1.90983323070104e-6  0.000817588644822336 0.00241359143105091

87  30    1.90983323070104e-6  0.000817276023109249 0.00241359143105091

88  30    1.90983323070104e-6  0.000783697321914271 0.00241359143105091

89  30    1.90983323070104e-6  0.000783697321914271 0.00241359143105091

90  30    1.90983323070104e-6  0.000783484831885969 0.00241359143105091

91  30    1.90983323070104e-6  0.000778126899100021 0.00241359143105091

92  30    1.90983323070104e-6  0.000739906989456947 0.00241359143105091

93  30    1.90983323070104e-6  0.000692672872693339 0.00241359143105091

94  30    1.90983323070104e-6  0.000664634971484744 0.00241359143105091

95  30    1.90983323070104e-6  0.000598849920978774 0.00208491900316071

96  30    1.90983323070104e-6  0.000597463439761913 0.00208491900316071

97  30    1.90983323070104e-6  0.000557235199953356 0.00208491900316071

98  30    1.90983323070104e-6  0.000557235199953356 0.00208491900316071

99  30    1.90983323070104e-6  0.00054807016227624  0.00208491900316071

To me it looks like it converges to a solution? 

""I setup an extremely large amount of individuals just for the sake of debugging at this point: NP = 100000 and run the simulation for 100 generations."

This makes me worry, trying to brute force things in so few generations depends on many things(visualize to see if fitness is done convergine, may need Log scale to see). You also may want to print and explore your individuals over the first couple gens. If not set up correctly or if there is a mistake in an operator or something then the evo will just carry on with the problem so focusing on the final results may not get you anywhere as the problem can occur very early on."

I've decreased the amount of individuals to 30 and the generations are still 100 as per the paper. My assumption was that maybe I was getting stuck in a local minima(unlikely although I'm struggling to figure this out) so I drastically increased the number of individuals initially.

 I did print out the individuals over a few generations and monitored them to ensure the code was working as expected (none go out of bounds, mutation and evaluation work, etc).

"Also in the paper they mention for selection they use equation 11, but i see in the notebook the use of: toolbox.register("select", tools.selRandom, k=3)
This low selection pressure with such a large pop(in the paper do they mention using only 30?!) may be an issue depending on if i skimmed the paper correctly for this info. " 

From the population, I use the registered select function to select three random individuals for the mutation process as denoted in section 3.2:

"Mutation is a perturbation or change with a random element.

In literature, a parent vector from the current generation

is known as target vector, a mutant vector achieved

through the differential mutation operation is called a

donor vector and finally an offspring formed by recombining

the donor with the target vector is called trial vector.

For a given parameter vector Xi,G, three vectors (Xr1,G,

Xr2,G, Xr3,G) are randomly selected in the range [1, NP],

such that the indices i, r1, r2 and r3 are distinct."

All generations are however evaluated and "selected" by the algorithm given in the paper in this section (here you can also see the selection for mutation highlighted in red):

while (generation < max):

        for i, agent in enumerate(pop):

            # Mutation

            a, b, c = [toolbox.clone(ind) for ind in toolbox.select(pop)]

            mutant1 = toolbox.clone(agent)

            mutant = toolbox.clone(agent)

            mutant = toolbox.mutate(mutant, a, b, c)

            trialvec = toolbox.mate(mutant1, mutant)

            trialvec1 = cxPenalty(trialvec)

            trialeval = (toolbox.evaluate(trialvec1[0], trialvec1[1], trialvec1[2], trialvec1[3]))

            curreval = (toolbox.evaluate(agent[0], agent[1], agent[2], agent[3]))

            pop[i].fitness.values = curreval, 

            if ((trialeval < curreval)):

                #selection -> trialvector is better than original vector

                pop[i] = trialvec

                pop[i].fitness.values = trialeval,

        hof.update(pop)

        record = mstats.compile(pop)

        logbook.record(gen=generation, evals=len(pop), **record)

        print(logbook.stream)

        x.append(logbook.chapters["fitness"].select("min"))

        lengthy.append(generation)

        generation = generation + 1 

Where the if statement applies the selection process denoted in equation 11.