Download Comparative Analysis of Numerical Optimization Algorithms and more Cheat Sheet Computer science in PDF only on Docsity!
HW
Yun-Chen, Tsai (R13222026)
November 2024
1 Problem 1
1.1 (a)
Take linear congruential generator(LCG) as an example: Let X be a sequence of pseudo-random numbers. with X 0 as a seed. Then the generator is defined as a recurrence relation:
Xn+1 = (a · Xn + b) mod m
Set a = 1664525, b = 1013904223, m = 2^32 , seed = 20 and generate 10000 pseudo-random numbers and plot the results using histogram:
Figure 1: 10000 pseudo-random numbers generated by linear congruential generator and numpy
As seen from above, it has great uniformity. The code is implemented in random number.py.
1.2 (b) Random sampling (hw2 1b.py)
Here I use np.random.uniform(low=0, high=1, size=3) to generate 3 uniform random numbers between 0 and 1. Initially the position is at (0.5, 0.5, 0.5) and the objective function value is -0.628022. These 2 values are set as currently best (xb, fb). At each iterations, 3 uniform random numbers are generated to be the position within the unit hypercube. Then the correponding objective function value is calculated. If the new objective function value is less than fb, then (xb, fb) will be updated. The random sampling algorithm is implemented in hw2 1b.py.
Iteration Point (xb) Function Value (fb) 1
Scipy Optimization Result: [0.114589, 0.555649, 0.852547], with function value: -3.862780.
1.3 (c) Mutation function
Here I define same mutation function as taught in the lecture. First generate 3 random numbers between -1 and 1 to be the random vectoru⃗ , then normalized it ˆu =u/⃗ u|⃗ |. The new point is calculated as follows: −x−→ new =^ −x−→ old +^ step^ ·^ uˆ.^ If one of the coordinate of the new point exceed either lower or upper bound, the correponding coordinate of the new point is set to the bound value. For example, if the new point is (1.2, 0.8, 0.7), x coordinate exceeds 1.0, so the new point will be (1.0, 0.8, 0.7). This critirion is the same for hill climbing and random walk. And the mutation function is called ”update” in hw2 1d.py and hw2 1e.py.
1.4 (d) Hill climbing (hw2 1d.py)
Step Iteration Point (xb) Function Value (fb)
- 00000 e − 01 1
- 00000 e − 01 4
- 00000 e − 01 4
- 00000 e − 01 5
- 00000 e − 01 5
- 00000 e − 01 5
- 00000 e − 01 5
- 00000 e − 01 5
- 00000 e − 01 5
- 00000 e − 01 22
- 00000 e − 01 31
- 00000 e − 01 34
- 00000 e − 01 1000
- 00000 e − 02 1000
- 50000 e − 02 67
- 50000 e − 02 104
- 50000 e − 02 180
- 50000 e − 02 437
- 50000 e − 02 1000
- 25000 e − 02 32
- 25000 e − 02 34
- 25000 e − 02 1000
- 25000 e − 03 367
- 25000 e − 03 411
- 25000 e − 03 655
- 25000 e − 03 1000
- 12500 e − 03 45
- 12500 e − 03 223
- 12500 e − 03 283
- 12500 e − 03 1000
- 56250 e − 03 379
- 56250 e − 03 1000
- 81250 e − 04 1000
- 90625 e − 04 1000
- 95313 e − 04 1000
- 76563 e − 05 1000
- 88281 e − 05 1000
- 44141 e − 05 1000
- 22070 e − 05 1000
- 10352 e − 06 1000
- 05176 e − 06 1000
- 52588 e − 06 1000
2.2 (b) Cauchy method (hw2 2b.py)
New direction in Cauchy method is simply: S 1 = −∇(f (X))|X 0 = [− 1. 046557 , 2. 742792 , 3 .643449]. Next use direction exploitation to find optimal step along this direction. Given a maximal step λ 1 = 0.1, slice it into 100 equally spaced points and perform a complete search. The result is λ∗ 1 = 0.080808, so the new point is X 1 = X 0 + λ∗ 1 S 1 = [0. 4154300. 7216400 .794420], with correponding function value: f (X 1 ) = − 2 .750658.
2.3 (c) Fletcher-Reeves (hw2 2d.py)
Start with X 1 calculated from 2(b), first calculate −∇f (X)|X 1 = [− 0. 730195 , − 7. 537728 , 5 .682158]. So
the ratio: r = |∇f^ (X)|X^1 |
2 |∇f (X)|X 0 |^2 = 4.094360. New direction: S 2 = −∇f (X)|X 1 +r·(−∇f (X)|X 0 ) = [− 5. 015177 , 3. 692249 , 20 .599747]. Set the maximal step to be λ 1 = 0.1 and use the same direction exploitation technique to find λ∗ 2 = 0.002020. So the new point is: X 2 = X 1 + λ∗ 2 S 2 = [0. 405298 , 0. 729099 , 0 .836036], with function value: f (X 2 ) : − 2 .866772.
2.4 (d) Fletcher-Reeves continued (hw2 2d.py)
Step (λ∗) Iteration Point (xb) Function Value (fb)
- 02020 e − 03 2
- 01010 e − 02 3
- 08081 e − 03 4
- 05051 e − 03 5
- 71717 e − 02 6
- 87879 e − 02 7
- 74747 e − 02 8
- 05051 e − 03 9
- 41414 e − 02 10
- 11111 e − 02 11
- 28283 e − 02 12
- 31313 e − 02 13
- 09091 e − 03 14
- 21212 e − 02 15
- 53535 e − 02 16
- 83838 e − 02 17
- 07071 e − 03 18
- 11111 e − 02 19
- 00000 e − 01 20
- 21212 e − 02 21
- 07071 e − 03 22
- 91919 e − 02 23
- 33333 e − 02 24
- 22222 e − 02 25
- 07071 e − 03 26
3 Problem 3 (hw2 3.py)
3.1 (a) Initial temperature
10 random points are the columns of the following matrix:
The initial temperature is the mean function values of these 10 points: -1.319623.
3.2 (b) First iteration of SA
New point: pnew.x = [0. 536213 , 0. 479619 , 0 .457870] with new function value: pnew.y = − 0 .429740. dE = pcur.y − pnew.y = − 0. 628022 − (− 0 .429740) = − 0. 198282 < 0 (Note that I changed the definition of dE since the initial temperature is negative) Although dE < 0, there’s a chance to accept it using the following procedure: Generate a uniform random number between 0 and 1: 0.805223. It is less then exp−dE/T^ = 0.860487, so this point is accepted as pcur.
3.3 (c) Second iteration of SA
Next point: pnew.x = [0. 540542 , 0. 561349 , 0 .421717] with new function value: pnew.y = − 0 .471208. dE = pcur.y − pnew.y = − 0. 429740 − (− 0 .471208) = 0. 041468 > 0. Since pnew is better then pcur, it is accepted as pcur.
3.4 (d) Remaining iterations of SA
Iteration Point (xb) Function Value (fb) 3
It takes 196843 iterations for the function value to converge. For the sake of curiosity, I tried another value of search range = ± 0 .01. It shows better convergence behavior:
Using 0-based indexing, 0th, 1st, 3rd, 4th, 5th, 6th, 7th and 8th particles have lower function values. So their local best positions should be updated acoordingly.
−−−−−−→ Plocalbest =
- 415046 0. 056009 0. 633648 0. 404275 0. 284544 0. 153614 0. 140826 0. 409884 0. 126977 0. 088340
- 700600 0. 828815 0. 003948 0. 617765 0. 761476 0. 665760 0. 718010 0. 512376 0. 811776 0. 714576
- 573058 0. 381538 0. 373341 0. 636004 0. 524882 0. 521142 0. 609041 0. 558991 0. 624650 0. 601039
Also, 8th particle has lowest function value now, so global best position is updated as:
−−−−−−→ Pglobalbest =
- 126977
- 811776
- 624650
With global best function value: -2.843445.
4.3 (c) Next iteration of PSO (hw2 4abc.py)
Using similar procedure, new random numbers, velocities and positions are:
r 1 , r 2 = 0. 908649 , 0. 319236
V 2 =
h− 0. 400000 0. 057913 − 0. 096924 − 0. 400000 − 0. 264264 − 0. 079686 − 0. 061658 − 0. 400000 − 0. 042134 0. 012334
- 050732 − 0. 130018 0. 173412 0. 167508 − 0. 035087 0. 099848 0. 026188 0. 316079 − 0. 105998 0. 031030
- 046983 0. 316978 0. 053159 − 0. 041754 0. 114899 0. 120171 − 0. 003743 0. 066814 − 0. 025748 0. 007537
i
X 2 =
h 0. 015046 0. 113921 0. 252266 0. 004275 0. 020281 0. 073928 0. 079168 0. 009884 0. 084843 0. 100674
- 751332 0. 698797 0. 548057 0. 785273 0. 726389 0. 765608 0. 744198 0. 828456 0. 705778 0. 745606
- 620041 0. 698515 0. 545278 0. 594250 0. 639781 0. 641314 0. 605298 0. 625805 0. 598902 0. 608576
i
Their function values are:
[− 2. 676097 − 2. 579238 − 1. 306032 − 2. 831252 − 2. 566089 − 2. 688510 − 2. 738507 − 2. 773287 − 2. 563926 − 2. 744254 ]
−−−−−−→ Plocalbest =
- 015046 0. 113921 0. 252266 0. 004275 0. 020281 0. 073928 0. 079168 0. 009884 0. 126977 0. 100674
- 751332 0. 698797 0. 548057 0. 785273 0. 726389 0. 765608 0. 744198 0. 828456 0. 811776 0. 745606
- 620041 0. 698515 0. 545278 0. 594250 0. 639781 0. 641314 0. 605298 0. 625805 0. 624650 0. 608576
−−−−−−→ Pglobalbest =
- 126977
- 811776
- 624650
Global best function value is: − 2 .843445.
4.4 (d) Remaining iterations of PSO (hw2 4d.py)
Next set the maximum iterations = 200000. Convergence critirion: If the variance of particle positions is less than 1e − 6 for 3 coordinates, then break the while loop. Result: iterations break at 3609.
−−−−−−→ Pglobalbest =
- 103705
- 553259
- 851089
With function value: − 3 .862313.
