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
| import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
G = 100
fun = '5 * sin(x * y) + x ** 2 + y ** 2'
choice = 'max'
dim = 2
limit_x = [-4,4]
limit_y = [-4,4]
pop_size = 10
mutate_rate = 0.7
delta = 0.2
beta = 1
colone_num = 10
if choice == 'max':
alpha = 2
else:
alpha = -2
def f(fun,x,y):
return eval(fun)
def init_pop(dim,pop_size,*limit):
pop = np.random.rand(dim,pop_size)
for i in range(dim):
pop[i,:] *= (limit[i][1] - limit[i][0])
pop[i,:] += limit[i][0]
return pop
def calc_density(pop,delta):
density = np.zeros([pop.shape[1],])
for i in range(pop.shape[1]):
density[i] = np.sum(len(np.ones([pop.shape[1],])[np.sqrt(np.sum((pop - pop[:,i].reshape([2,1])) ** 2,axis = 0)) < delta]))
return density / pop.shape[1]
def calc_simulation(simulation,density):
return (alpha * simulation - beta * density)
def mutate(x,gen,mutata_rate,dim,*limit):
for i in range(dim):
if np.random.rand() <= mutata_rate:
x[i] += (np.random.rand() - 0.5) * (limit[i][1] - limit[i][0]) / (gen + 1)
if (x[i] > limit[i][1]) or (x[i] < limit[i][0]):
x[i] = np.random.rand() * (limit[i][1] - limit[i][0]) + limit[i][0]
pop = init_pop(dim,pop_size,limit_x,limit_y)
init_simulation = f(fun,pop[0,:],pop[1,:])
density = calc_density(pop,delta)
simulation = calc_simulation(init_simulation,density)
index = np.argsort(-simulation)
pop = pop[:,index]
new_pop = pop.copy()
simulation = simulation[index]
for gen in range(G):
best_a = np.zeros([dim,int(pop_size / 2)])
best_a_simulation = np.zeros([int(pop_size / 2),])
for i in range(int(pop_size / 2)):
a = new_pop[:,i].reshape([2,1])
b = np.tile(a,(1,colone_num))
for j in range(colone_num):
mutate(a,gen,mutate_rate,dim,limit_x,limit_y)
b[:,0] = pop[:,i]
b_simulation = f(fun,b[0,:],b[1,:])
index = np.argsort(-b_simulation)
best_a_simulation = b_simulation[index][0]
best_a[:,i] = b[:,index][:,0]
a_density = calc_density(best_a,delta)
best_a_simulation = calc_simulation(best_a_simulation,a_density)
new_a = init_pop(2,int(pop_size / 2),limit_x,limit_y)
new_a_simulation = f(fun,new_a[0,:],new_a[1,:])
new_a_density = calc_density(new_a,delta)
new_a_simulation = calc_simulation(new_a_simulation,new_a_density)
pop = np.hstack([best_a,new_a])
simulation = np.hstack([best_a_simulation,new_a_simulation])
index = np.argsort(-simulation)
pop = pop[:,index]
simulation = simulation[index]
new_pop = pop.copy()
figure = plt.figure(figsize = (10,8),dpi = 80)
ax = Axes3D(figure)
plt.xlabel("x")
plt.ylabel("y")
for i in range(int(pop_size / 2)):
ax.scatter(pop[0,i],pop[1,i],f(fun,pop[0,i],pop[1,i]),color = 'red')
print('最优解:','x = ',pop[0,i],'y = ',pop[1,i],end = '\n')
print('结果:','z = ',f(fun,pop[0,i],pop[1,i]))
x = np.arange(limit_x[0],limit_x[1],(limit_x[1] - limit_x[0]) / 50)
y = np.arange(limit_y[0],limit_y[1],(limit_y[1] - limit_y[0]) / 50)
x,y = np.meshgrid(x,y)
z = f(fun,x,y)
ax.plot_surface(x,y,z, rstride=1, cstride=1, color = 'green',alpha = 0.5)
plt.show()
|
