-
Notifications
You must be signed in to change notification settings - Fork 1
/
HS.py
132 lines (119 loc) · 4.01 KB
/
HS.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2022/11/9 9:50
# @Author : Xavier Ma
# @Email : xavier_mayiming@163.com
# @File : HS.py
# @Statement : Harmony Search
# @Reference : A new heuristic optimization algorithm: Harmony search[J]. Simulation, 2001, 2(2):60-68
import random
import math
import matplotlib.pyplot as plt
def obj(x):
"""
The objective function of pressure vessel design
:param x:
:return:
"""
x1 = x[0]
x2 = x[1]
x3 = x[2]
x4 = x[3]
g1 = -x1 + 0.0193 * x3
g2 = -x2 + 0.00954 * x3
g3 = -math.pi * x3 ** 2 - 4 * math.pi * x3 ** 3 / 3 + 1296000
g4 = x4 - 240
if g1 <= 0 and g2 <= 0 and g3 <= 0 and g4 <= 0:
return 0.6224 * x1 * x3 * x4 + 1.7781 * x2 * x3 ** 2 + 3.1661 * x1 ** 2 * x4 + 19.84 * x1 ** 2 * x3
else:
return 1e10
def boundary_check(value, lb, ub):
"""
The boundary check
:param value:
:param lb: the lower bound (list)
:param ub: the upper bound (list)
:return:
"""
for i in range(len(value)):
value[i] = max(value[i], lb[i])
value[i] = min(value[i], ub[i])
return value
def main(hms, iter, hmcr, par, bw, nnew, lb, ub):
"""
The main function of the HS
:param hms: harmony memory size
:param iter: the number of iterations
:param hmcr: harmony memory consideration rate
:param par: pitch adjustment rate
:param bw: bandwidth
:param nnew: the number of new harmonies created in each iteration
:param lb: the lower bound (list)
:param ub: the upper bound (list)
:return:
"""
# Step 1. Initialization
pos = [] # the set of harmonies
score = [] # the score of harmonies
dim = len(lb) # dimension
for i in range(hms):
temp_pos = [random.uniform(lb[j], ub[j]) for j in range(dim)]
pos.append(temp_pos)
score.append(obj(temp_pos))
iter_best = []
gbest = min(score) # the score of the best-so-far harmony
gbest_pos = pos[score.index(gbest)].copy() # the best-so-far harmony
con_iter = 0
# Step 2. The main loop
for t in range(iter):
new_pos = []
new_score = []
# Step 2.1. Create new harmonies
for _ in range(nnew):
temp_pos = []
for j in range(dim):
if random.random() < hmcr: # utilize harmony memory
ind = random.randint(0, hms - 1)
temp_pos.append(pos[ind][j])
if random.random() < par: # pitch adjustment
temp_pos[j] += random.normalvariate(0, 1) * bw * (ub[j] - lb[j])
else:
temp_pos.append(random.uniform(lb[j], ub[j]))
temp_pos = boundary_check(temp_pos, lb, ub)
new_pos.append(temp_pos)
new_score.append(obj(temp_pos))
# Step 2.2. Update harmony memory
new_pos.extend(pos)
new_score.extend(score)
sorted_score = sorted(new_score)
pos = []
score = []
for i in range(hms):
score.append(sorted_score[i])
pos.append(new_pos[new_score.index(sorted_score[i])])
# Step 2.3. Update the global best
if score[0] < gbest:
gbest = score[0]
gbest_pos = pos[0].copy()
con_iter = t + 1
iter_best.append(gbest)
# Step 3. Sort the results
x = [i for i in range(iter)]
plt.figure()
plt.plot(x, iter_best, linewidth=2, color='blue')
plt.xlabel('Iteration number')
plt.ylabel('Global optimal value')
plt.title('Convergence curve')
plt.show()
return {'best score': gbest, 'best solution': gbest_pos, 'convergence iteration': con_iter}
if __name__ == '__main__':
# Parameter settings
hms = 30
iter = 1000
hmcr = 0.9
par = 0.1
bw = 0.02
nnew = 20
lb = [0, 0, 10, 10]
ub = [99, 99, 200, 200]
print(main(hms, iter, hmcr, par, bw, nnew, lb, ub))