README
scikit-opt
[](https://pypi.org/project/scikit-opt/)
[](https://github.com/guofei9987/scikit-opt)
[](https://travis-ci.com/guofei9987/scikit-opt)
[](https://codecov.io/gh/guofei9987/scikit-opt)
[](https://pypi.org/project/scikit-opt/)
[](https://github.com/guofei9987/scikit-opt/stargazers)
[](https://github.com/guofei9987/scikit-opt/network/members)
[](https://gitter.im/guofei9987/scikit-opt?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)
Heuristic Algorithms in Python
(Genetic Algorithm, Particle Swarm Optimization, Simulated Annealing, Ant Colony Algorithm, Immune Algorithm,Artificial Fish Swarm Algorithm in Python)
- Documentation: https://scikit-opt.github.io/scikit-opt/#/en/
- 文档： https://scikit-opt.github.io/scikit-opt/#/zh/
- Source code: https://github.com/guofei9987/scikit-opt
---
En/ Coverpage
scikit-opt
Powerful Python module for Heuristic Algorithms
* [](https://pypi.org/project/scikit-opt/)
[](https://travis-ci.com/guofei9987/scikit-opt)
[](https://codecov.io/gh/guofei9987/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/blob/master/LICENSE)
[](https://pepy.tech/project/scikit-opt)
* Differential Evolution, Genetic Algorithm
* Particle Swarm Optimization, Simulated Annealing
* Ant Colony Algorithm, Immune Algorithm
* Artificial Fish Swarm Algorithm
---
En/ Sidebar
* Document
* About parameters
* More Genetic Algorithm
* More Particle Swarm Optimization
* More Simulated Annealing
* Curve fiting
* Speed Up
---
En/Args
input parameters
Use help (like help(GA), GA?) to see the details.
import skohelp(sko.GA.GA)
help(sko.GA.GA_TSP)
help(sko.PSO.PSO)
help(sko.DE.DE)
help(sko.SA.SA)
help(sko.SA.SA_TSP)
help(sko.ACA.ACA_TSP)
help(sko.IA.IA_TSP)
help(sko.AFSA.AFSA)
GA
| input parameters | default value | meaning |
|-----------------|--------|------------------------|
| func | \- | objective function |
| n\_dim | \- | dim of objective function (number of input parameters) |
| size\_pop | 50 | size of population |
| max\_iter | 200 | max iteration |
| prob\_mut | 0\.001 | probability of mutation |
| lb | \-1 | lower bound of variables. Can be int/float/list |
| ub | 1 | upper bound of variables. Can be int/float/list |
| constraint\_eq | tuple() | equal constraint |
| constraint\_ueq | tuple() | unequal constraint |
| precision | 1e\-7 | precision，int/float or list |
| n_processes | 0 | number of processes, 0 means use all cpu |
GA_TSP
| input parameters | default value | meaning |
|-----------|--------|--------|
| func | \- | objective function |
| n\_dim | \- | num of cities |
| size\_pop | 50 | size of population |
| max\_iter | 200 | max iteration |
| prob\_mut | 0\.001 | probability of mutation |
PSO
| input parameters | default value | meaning |
|-----------|------|----------|
| func | \- | objective function |
| n\_dim | \- | dim of objective function |
| size\_pop | 50 | size of population |
| max\_iter | 200 | max iteration |
| lb | None | lower bound of variables |
| ub | None | upper bound of variables |
| w | 0\.8 | inertia weight |
| c1 | 0\.5 | cognitive parameter |
| c2 | 0\.5 | social parameter |
| constraint\_ueq | tuple() | unequal constraint |
| n_processes | 0 | number of processes, 0 means use all cpu |
DE
| input parameters | default value | meaning |
|-----------|--------|----------|
| func | \- | objective function |
| n\_dim | \- | dim of objective function |
| size\_pop | 50 | size of population |
| max\_iter | 200 | max iteration |
| prob\_mut | 0\.001 | probability of mutation |
| F | 0\.5 | coefficient of mutation |
| lb | \-1 | lower bound of variables |
| ub | 1 | upper bound of variables |
| constraint\_eq | tuple() | equal constraint |
| constraint\_ueq | tuple() | unequal constraint |
| n_processes | 0 | number of processes, 0 means use all cpu |
SA
| input parameters | default value | meaning |
|--------------------|-------|-------|
| func | \- | objective function |
| x0 | \- | initial point |
| T\_max | 100 | max temperature |
| T\_min | 1e\-7 | min temperature |
| L | 300 | long of chain |
| max\_stay\_counter | 150 | cooldown time |
| lb | \-1 | lower bound of variables. Can be int/float/list |
| ub | 1 | upper bound of variables. Can be int/float/list |
SA_TSP
| input parameters | default value | meaning |
|--------------------|-------|-------|
| func | \- | objective function |
| x0 | \- | initial point |
| T\_max | 100 | max temperature |
| T\_min | 1e\-7 | min temperature |
| L | 300 | long of chain |
| max\_stay\_counter | 150 | cooldown time |
ACA_TSP
| input parameters | default value | meaning |
|------------------|------|----------------------|
| func | \- | objective function |
| n\_dim | \- | number of cities ( also called dim of objective function) |
| size\_pop | 10 | number of ants |
| max\_iter | 20 | max iteration |
| distance\_matrix | \- | distance matrix between cities |
| alpha | 1 | importance of pheromone |
| beta | 2 | importance of fitness |
| rho | 0\.1 | evaporation speed of pheromone |
IA_TSP
| input parameters | default value | meaning |
|-----------|--------|----------------------------------|
| func | \- | objective function |
| n\_dim | \- | number of cities ( also called dim of objective function) |
| size\_pop | 50 | size of population |
| max\_iter | 200 | max iteration |
| prob\_mut | 0\.001 | probability of mutation |
| T | 0\.7 | concentration of antibody |
| alpha | 0\.95 | importance of diversity compared to concentration of antibody |
AFSA
| input parameters | default value | meaning |
|---------------|-------|------------------|
| func | \- | objective function |
| n\_dim | \- | dim of objective function |
| size\_pop | 50 | size of population |
| max\_iter | 300 | max iteration |
| max\_try\_num | 100 | max try of prey in one movement |
| step | 0\.5 | max scale of movement |
| visual | 0\.3 | max range of perception |
| q | 0\.98 | perception of fish will go down every movement |
| delta | 0\.5 | fishes toleration of crowd |
outputs
(in this part,
x is also called value of objective funtion. x is also called inputs of objective funtion)GA&GA_TSP
- ga.generation_best_Y best Y of every generation
- ga.generation_best_X X for best Y of every generation
- ga.all_history_FitV fitness value of every generation
- ga.all_history_Y function value of every generation and every individual
- ga.best_y best y
- ga.best_x best x
DE
- de.generation_best_Y best Y of every generation
- de.generation_best_X X for best Y of every generation
- de.all_history_Y Y of every generation and every individual
- de.best_y best y
- de.best_x best x
PSO
-
pso.record_value location, velocity, function value of every generation and every particles. only when pso.record_mode = True.-
pso.gbest_y_hist best y every generation-
pso.best_y best y （In PSO, use pso.gbest_x, pso.gbest_y）-
pso.best_x best xSA
- sa.generation_best_Y best Y of every generation
- sa.generation_best_X X for best Y of every generation
- sa.best_x best x
- sa.best_y best y
ACA
-
sa.generation_best_Y best Y of every generation-
sa.generation_best_X X for best Y of every generation-
aca.best_y best y-
aca.best_x best xAFSA
-
afsa.best_x best x-
afsa.best_y best yIA
- ia.generation_best_Y best Y of every generation
- ia.generation_best_X X for best Y of every generation
- ia.all_history_FitV 每一代的每个个体的适应度
- ia.all_history_Y 每一代每个个体的函数值
- ia.best_y best y
- ia.best_x best x
---
En/Contributors
contributors
Thanks to those contributors:
- Bowen Zhang
- hranYin
---
En/Curve Fitting
curve fitting using GA
Generate toy train datasets
import numpy as np
import matplotlib.pyplot as plt
from sko.GA import GAx_true = np.linspace(-1.2, 1.2, 30)
y_true = x_true 3 - x_true + 0.4 * np.random.rand(30)
plt.plot(x_true, y_true, 'o')
Make up residuals
def f_fun(x, a, b, c, d):
return a x 3 + b x 2 + c * x + d
def obj_fun(p):
a, b, c, d = p
residuals = np.square(f_fun(x_true, a, b, c, d) - y_true).sum()
return residuals
Do GA
ga = GA(func=obj_fun, n_dim=4, size_pop=100, max_iter=500,
lb=[-2] 4, ub=[2] 4)best_params, residuals = ga.run()
print('best_x:', best_params, '\n', 'best_y:', residuals)
Plot the fitting results
y_predict = f_fun(x_true, *best_params)fig, ax = plt.subplots()
ax.plot(x_true, y_true, 'o')
ax.plot(x_true, y_predict, '-')
plt.show()
---
En/More Ga
Do integer programming with genetic algorithm
If you want some variables to be integer, then set the corresponding precision to an integer
For example, our objective function is demo_func. We want the variables to be integer interval 2, integer interval 1, float. We set precision=[2, 1, 1e-7]:
from sko.GA import GAdemo_func = lambda x: (x[0] - 1) 2 + (x[1] - 0.05) 2 + x[2] 2
ga = GA(func=demo_func, n_dim=3, max_iter=500, lb=[-1, -1, -1], ub=[5, 1, 1], precision=[2, 1, 1e-7])
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
Notice:
- If precision is an integer, the number of all possible value would better be $2^n$, in which case the performance is the best. It also works if the number is not $2^n$
- If precision is not an integer, but you still want this mode, manually deal with it. For example, your original precision=0.5, just make a new variable, multiplied by 2
How to fix start point and end point with GA for TSP
If it is not a cycle graph, no need to do this.
if your start point and end point is (0, 0) and (1, 1). Build up the object function :
- Start point and end point is not the input of the object function. If totally n+2 points including start and end points, the input is the n points.
- And build up the object function, which is the total distance, as actually they are.
import numpy as np
from scipy import spatial
import matplotlib.pyplot as pltnum_points = 20
points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points
start_point=[[0,0]]
end_point=[[1,1]]
points_coordinate=np.concatenate([points_coordinate,start_point,end_point])
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
num_points, = routine.shape
routine = np.concatenate([[num_points], routine, [num_points+1]])
return sum([distance_matrix[routine[i], routine[i + 1]] for i in range(num_points+2-1)])
And the same with others:
from sko.GA import GA_TSPga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
best_points, best_distance = ga_tsp.run()
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([[num_points],best_points, [num_points+1]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
ax[1].plot(ga_tsp.generation_best_Y)
plt.show()
For more information, click here
How to set up starting point or initial population
- For GA, after ga=GA(params), use codes like ga.Chrom = np.random.randint(0,2,size=(80,20)) to manually set the initial population.
- For DE, set de.X to your initial X.
- For SA, there is a parameter x0, which is the init point.
- For PSO, set pso.X to your initial X, and run pso.cal_y(); pso.update_gbest(); pso.update_pbest()
---
En/More Pso
demonstrate PSO with animation
step1:do pso
-> Demo code: examples/demo_pso_ani.py#s1
import numpy as np
from sko.PSO import PSO
def demo_func(x):
x1, x2 = x
return -20 np.exp(-0.2 np.sqrt(0.5 * (x1 2 + x2 2))) - np.exp(
0.5 (np.cos(2 np.pi x1) + np.cos(2 np.pi * x2))) + 20 + np.e
constraint_ueq = (
lambda x: (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2
,
)
max_iter = 50
pso = PSO(func=demo_func, n_dim=2, pop=40, max_iter=max_iter, lb=[-2, -2], ub=[2, 2]
, constraint_ueq=constraint_ueq)
pso.record_mode = True
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
step2: plot animation
-> Demo code: examples/demo_pso_ani.py#s2
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimationrecord_value = pso.record_value
X_list, V_list = record_value['X'], record_value['V']
fig, ax = plt.subplots(1, 1)
ax.set_title('title', loc='center')
line = ax.plot([], [], 'b.')
X_grid, Y_grid = np.meshgrid(np.linspace(-2.0, 2.0, 40), np.linspace(-2.0, 2.0, 40))
Z_grid = demo_func((X_grid, Y_grid))
ax.contour(X_grid, Y_grid, Z_grid, 30)
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
t = np.linspace(0, 2 * np.pi, 40)
ax.plot(0.5 np.cos(t) + 1, 0.5 np.sin(t), color='r')
plt.ion()
p = plt.show()
def update_scatter(frame):
i, j = frame // 10, frame % 10
ax.set_title('iter = ' + str(i))
X_tmp = X_list[i] + V_list[i] * j / 10.0
plt.setp(line, 'xdata', X_tmp[:, 0], 'ydata', X_tmp[:, 1])
return line
ani = FuncAnimation(fig, update_scatter, blit=True, interval=25, frames=max_iter * 10)
plt.show()
ani.save('pso.gif', writer='pillow')
---
En/More Sa
3 types of Simulated Annealing
In the ‘fast’ schedule the updates are:
u ~ Uniform(0, 1, size = d)
y = sgn(u - 0.5) T ((1 + 1/T)abs(2*u - 1) - 1.0)xc = y * (upper - lower)
x_new = x_old + xc
c = n exp(-n quench)
T_new = T0 exp(-c kquench)
In the ‘cauchy’ schedule the updates are:
u ~ Uniform(-pi/2, pi/2, size=d)
xc = learn_rate T tan(u)
x_new = x_old + xcT_new = T0 / (1 + k)
In the ‘boltzmann’ schedule the updates are:
std = minimum(sqrt(T) ones(d), (upper - lower) / (3learn_rate))
y ~ Normal(0, std, size = d)
x_new = x_old + learn_rate * yT_new = T0 / log(1 + k)
Do Simulated Annealing
#### 1.1 Fast Simulated Annealing
-> Demo code: examples/demo_sa.py#s4
from sko.SA import SAFastsa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_fast.run()
print('Fast Simulated Annealing: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)
#### 1.2 Fast Simulated Annealing with bounds
-> Demo code: examples/demo_sa.py#s5
from sko.SA import SAFastsa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
lb=[-1, 1, -1], ub=[2, 3, 4])
sa_fast.run()
print('Fast Simulated Annealing with bounds: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)
#### 2.1 Boltzmann Simulated Annealing
-> Demo code: examples/demo_sa.py#s6
from sko.SA import SABoltzmannsa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_boltzmann.run()
print('Boltzmann Simulated Annealing: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_boltzmann.best_y)
#### 2.2 Boltzmann Simulated Annealing with bounds
-> Demo code: examples/demo_sa.py#s7
from sko.SA import SABoltzmannsa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
lb=-1, ub=[2, 3, 4])
sa_boltzmann.run()
print('Boltzmann Simulated Annealing with bounds: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_boltzmann.best_y)
#### 3. Cauchy Simulated Annealing
-> Demo code: examples/demo_sa.py#s8
from sko.SA import SACauchysa_cauchy = SACauchy(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_cauchy.run()
print('Cauchy Simulated Annealing: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)
#### 3. Cauchy Simulated Annealing with bounds
-> Demo code: examples/demo_sa.py#s9
from sko.SA import SACauchysa_cauchy = SACauchy(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
lb=[-1, 1, -1], ub=[2, 3, 4])
sa_cauchy.run()
print('Cauchy Simulated Annealing with bounds: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)
---
En/README
scikit-opt
[](https://pypi.org/project/scikit-opt/)
[](https://travis-ci.com/guofei9987/scikit-opt)
[](https://codecov.io/gh/guofei9987/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/blob/master/LICENSE)
[](https://github.com/guofei9987/scikit-opt/fork)
[](https://pepy.tech/project/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/discussions)
Swarm Intelligence in Python
(Genetic Algorithm, Particle Swarm Optimization, Simulated Annealing, Ant Colony Algorithm, Immune Algorithm,Artificial Fish Swarm Algorithm in Python)
- Documentation: https://scikit-opt.github.io/scikit-opt/#/en/
- 文档： https://scikit-opt.github.io/scikit-opt/#/zh/
- Source code: https://github.com/guofei9987/scikit-opt
- Help us improve scikit-opt https://www.wjx.cn/jq/50964691.aspx
install
pip install scikit-optFor the current developer version:
git clone [email protected]:guofei9987/scikit-opt.git
cd scikit-opt
pip install .Features
Feature1: UDF
UDF (user defined function) is available now!
For example, you just worked out a new type of selection function.
Now, your selection function is like this:
-> Demo code: examples/demo_ga_udf.py#s1
step1: define your own operator:
def selection_tournament(algorithm, tourn_size):
FitV = algorithm.FitV
sel_index = []
for i in range(algorithm.size_pop):
aspirants_index = np.random.choice(range(algorithm.size_pop), size=tourn_size)
sel_index.append(max(aspirants_index, key=lambda i: FitV[i]))
algorithm.Chrom = algorithm.Chrom[sel_index, :] # next generation
return algorithm.ChromImport and build ga
-> Demo code: examples/demo_ga_udf.py#s2
import numpy as np
from sko.GA import GA, GA_TSPdemo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + (x[2] - 0.5) 2
ga = GA(func=demo_func, n_dim=3, size_pop=100, max_iter=500, prob_mut=0.001,
lb=[-1, -10, -5], ub=[2, 10, 2], precision=[1e-7, 1e-7, 1])
Regist your udf to GA
-> Demo code: examples/demo_ga_udf.py#s3
ga.register(operator_name='selection', operator=selection_tournament, tourn_size=3)scikit-opt also provide some operators
-> Demo code: examples/demo_ga_udf.py#s4
from sko.operators import ranking, selection, crossover, mutationga.register(operator_name='ranking', operator=ranking.ranking). \
register(operator_name='crossover', operator=crossover.crossover_2point). \
register(operator_name='mutation', operator=mutation.mutation)
Now do GA as usual
-> Demo code: examples/demo_ga_udf.py#s5
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)Until Now, the udf surportcrossover,mutation,selection,rankingof GA
scikit-opt provide a dozen of operators, see here
For advanced users:
-> Demo code: examples/demo_ga_udf.py#s6
class MyGA(GA):
def selection(self, tourn_size=3):
FitV = self.FitV
sel_index = []
for i in range(self.size_pop):
aspirants_index = np.random.choice(range(self.size_pop), size=tourn_size)
sel_index.append(max(aspirants_index, key=lambda i: FitV[i]))
self.Chrom = self.Chrom[sel_index, :] # next generation
return self.Chrom ranking = ranking.ranking
demo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + (x[2] - 0.5) 2
my_ga = MyGA(func=demo_func, n_dim=3, size_pop=100, max_iter=500, lb=[-1, -10, -5], ub=[2, 10, 2],
precision=[1e-7, 1e-7, 1])
best_x, best_y = my_ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
feature2: continue to run
(New in version 0.3.6)
Run an algorithm for 10 iterations, and then run another 20 iterations base on the 10 iterations before:
from sko.GA import GAfunc = lambda x: x[0] 2
ga = GA(func=func, n_dim=1)
ga.run(10)
ga.run(20)
feature3: 4-ways to accelerate
- vectorization
- multithreading
- multiprocessing
- cached
see https://github.com/guofei9987/scikit-opt/blob/master/examples/example_function_modes.py
feature4: GPU computation
We are developing GPU computation, which will be stable on version 1.0.0
An example is already available: https://github.com/guofei9987/scikit-opt/blob/master/examples/demo_ga_gpu.py
Quick start
1. Differential Evolution
Step1：define your problem
-> Demo code: examples/demo_de.py#s1
'''
min f(x1, x2, x3) = x1^2 + x2^2 + x3^2
s.t.
x1*x2 >= 1
x1*x2 <= 5
x2 + x3 = 1
0 <= x1, x2, x3 <= 5
'''
def obj_func(p):
x1, x2, x3 = p
return x1 2 + x2 2 + x3 2
constraint_eq = [
lambda x: 1 - x[1] - x[2]
]
constraint_ueq = [
lambda x: 1 - x[0] * x[1],
lambda x: x[0] * x[1] - 5
]
Step2: do Differential Evolution
-> Demo code: examples/demo_de.py#s2
from sko.DE import DEde = DE(func=obj_func, n_dim=3, size_pop=50, max_iter=800, lb=[0, 0, 0], ub=[5, 5, 5],
constraint_eq=constraint_eq, constraint_ueq=constraint_ueq)
best_x, best_y = de.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
2. Genetic Algorithm
Step1：define your problem
-> Demo code: examples/demo_ga.py#s1
import numpy as np
def schaffer(p):
'''
This function has plenty of local minimum, with strong shocks
global minimum at (0,0) with value 0
https://en.wikipedia.org/wiki/Test_functions_for_optimization
'''
x1, x2 = p
part1 = np.square(x1) - np.square(x2)
part2 = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(part1)) - 0.5) / np.square(1 + 0.001 * part2)
Step2: do Genetic Algorithm
-> Demo code: examples/demo_ga.py#s2
from sko.GA import GAga = GA(func=schaffer, n_dim=2, size_pop=50, max_iter=800, prob_mut=0.001, lb=[-1, -1], ub=[1, 1], precision=1e-7)
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
-> Demo code: examples/demo_ga.py#s3
import pandas as pd
import matplotlib.pyplot as pltY_history = pd.DataFrame(ga.all_history_Y)
fig, ax = plt.subplots(2, 1)
ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
Y_history.min(axis=1).cummin().plot(kind='line')
plt.show()
2.2 Genetic Algorithm for TSP(Travelling Salesman Problem)
Just import the
GA_TSP, it overloads the crossover, mutation to solve the TSPStep1: define your problem. Prepare your points coordinate and the distance matrix.
Here I generate the data randomly as a demo:
-> Demo code: examples/demo_ga_tsp.py#s1
import numpy as np
from scipy import spatial
import matplotlib.pyplot as pltnum_points = 50
points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
num_points, = routine.shape
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
Step2: do GA
-> Demo code: examples/demo_ga_tsp.py#s2
from sko.GA import GA_TSPga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
best_points, best_distance = ga_tsp.run()
Step3: Plot the result:
-> Demo code: examples/demo_ga_tsp.py#s3
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
ax[1].plot(ga_tsp.generation_best_Y)
plt.show()3. PSO(Particle swarm optimization)
3.1 PSO
Step1: define your problem:
-> Demo code: examples/demo_pso.py#s1
def demo_func(x):
x1, x2, x3 = x
return x1 2 + (x2 - 0.05) 2 + x3 2Step2: do PSO
-> Demo code: examples/demo_pso.py#s2
from sko.PSO import PSOpso = PSO(func=demo_func, n_dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
Step3: Plot the result
-> Demo code: examples/demo_pso.py#s3
import matplotlib.pyplot as pltplt.plot(pso.gbest_y_hist)
plt.show()
3.2 PSO with nonlinear constraint
If you need nolinear constraint like (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2<=0
Codes are like this:
constraint_ueq = (
lambda x: (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2
,
)
pso = PSO(func=demo_func, n_dim=2, pop=40, max_iter=max_iter, lb=[-2, -2], ub=[2, 2]
, constraint_ueq=constraint_ueq)Note that, you can add more then one nonlinear constraint. Just add it to constraint_ueq
More over, we have an animation:
↑see examples/demo_pso_ani.py
4. SA(Simulated Annealing)
4.1 SA for multiple function
Step1: define your problem
-> Demo code: examples/demo_sa.py#s1
demo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + x[2] 2Step2: do SA
-> Demo code: examples/demo_sa.py#s2
from sko.SA import SAsa = SA(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150)
best_x, best_y = sa.run()
print('best_x:', best_x, 'best_y', best_y)
Step3: Plot the result
-> Demo code: examples/demo_sa.py#s3
import matplotlib.pyplot as plt
import pandas as pdplt.plot(pd.DataFrame(sa.best_y_history).cummin(axis=0))
plt.show()
Moreover, scikit-opt provide 3 types of Simulated Annealing: Fast, Boltzmann, Cauchy. See more sa
4.2 SA for TSP
Step1: oh, yes, define your problems. To boring to copy this step.
Step2: DO SA for TSP
-> Demo code: examples/demo_sa_tsp.py#s2
from sko.SA import SA_TSPsa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points)
best_points, best_distance = sa_tsp.run()
print(best_points, best_distance, cal_total_distance(best_points))
Step3: plot the result
-> Demo code: examples/demo_sa_tsp.py#s3
from matplotlib.ticker import FormatStrFormatterfig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(sa_tsp.best_y_history)
ax[0].set_xlabel("Iteration")
ax[0].set_ylabel("Distance")
ax[1].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1],
marker='o', markerfacecolor='b', color='c', linestyle='-')
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].set_xlabel("Longitude")
ax[1].set_ylabel("Latitude")
plt.show()
More: Plot the animation:
5. ACA (Ant Colony Algorithm) for tsp
-> Demo code: examples/demo_aca_tsp.py#s2
from sko.ACA import ACA_TSPaca = ACA_TSP(func=cal_total_distance, n_dim=num_points,
size_pop=50, max_iter=200,
distance_matrix=distance_matrix)
best_x, best_y = aca.run()
6. immune algorithm (IA)
-> Demo code: examples/demo_ia.py#s2
from sko.IA import IA_TSPia_tsp = IA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=500, max_iter=800, prob_mut=0.2,
T=0.7, alpha=0.95)
best_points, best_distance = ia_tsp.run()
print('best routine:', best_points, 'best_distance:', best_distance)
7. Artificial Fish Swarm Algorithm (AFSA)
-> Demo code: examples/demo_afsa.py#s1
def func(x):
x1, x2 = x
return 1 / x1 2 + x1 2 + 1 / x2 2 + x2 2
from sko.AFSA import AFSA
afsa = AFSA(func, n_dim=2, size_pop=50, max_iter=300,
max_try_num=100, step=0.5, visual=0.3,
q=0.98, delta=0.5)
best_x, best_y = afsa.run()
print(best_x, best_y)
Projects using scikit-opt
- Yu, J., He, Y., Yan, Q., & Kang, X. (2021). SpecView: Malware Spectrum Visualization Framework With Singular Spectrum Transformation. IEEE Transactions on Information Forensics and Security, 16, 5093-5107.
- Zhen, H., Zhai, H., Ma, W., Zhao, L., Weng, Y., Xu, Y., ... & He, X. (2021). Design and tests of reinforcement-learning-based optimal power flow solution generator. Energy Reports.
- Heinrich, K., Zschech, P., Janiesch, C., & Bonin, M. (2021). Process data properties matter: Introducing gated convolutional neural networks (GCNN) and key-value-predict attention networks (KVP) for next event prediction with deep learning. Decision Support Systems, 143, 113494.
- Tang, H. K., & Goh, S. K. (2021). A Novel Non-population-based Meta-heuristic Optimizer Inspired by the Philosophy of Yi Jing. arXiv preprint arXiv:2104.08564.
- Wu, G., Li, L., Li, X., Chen, Y., Chen, Z., Qiao, B., ... & Xia, L. (2021). Graph embedding based real-time social event matching for EBSNs recommendation. World Wide Web, 1-22.
- Pan, X., Zhang, Z., Zhang, H., Wen, Z., Ye, W., Yang, Y., ... & Zhao, X. (2021). A fast and robust mixture gases identification and concentration detection algorithm based on attention mechanism equipped recurrent neural network with double loss function. Sensors and Actuators B: Chemical, 342, 129982.
- Castella Balcell, M. (2021). Optimization of the station keeping system for the WindCrete floating offshore wind turbine.
- Zhai, B., Wang, Y., Wang, W., & Wu, B. (2021). Optimal Variable Speed Limit Control Strategy on Freeway Segments under Fog Conditions. arXiv preprint arXiv:2107.14406.
- Yap, X. H. (2021). Multi-label classification on locally-linear data: Application to chemical toxicity prediction.
- Gebhard, L. (2021). Expansion Planning of Low-Voltage Grids Using Ant Colony Optimization Ausbauplanung von Niederspannungsnetzen mithilfe eines Ameisenalgorithmus.
- Ma, X., Zhou, H., & Li, Z. (2021). Optimal Design for Interdependencies between Hydrogen and Power Systems. IEEE Transactions on Industry Applications.
- de Curso, T. D. C. (2021). Estudo do modelo Johansen-Ledoit-Sornette de bolhas financeiras.
- Wu, T., Liu, J., Liu, J., Huang, Z., Wu, H., Zhang, C., ... & Zhang, G. (2021). A Novel AI-based Framework for AoI-optimal Trajectory Planning in UAV-assisted Wireless Sensor Networks. IEEE Transactions on Wireless Communications.
- Liu, H., Wen, Z., & Cai, W. (2021, August). FastPSO: Towards Efficient Swarm Intelligence Algorithm on GPUs. In 50th International Conference on Parallel Processing (pp. 1-10).
- Mahbub, R. (2020). Algorithms and Optimization Techniques for Solving TSP.
- Li, J., Chen, T., Lim, K., Chen, L., Khan, S. A., Xie, J., & Wang, X. (2019). Deep learning accelerated gold nanocluster synthesis. Advanced Intelligent Systems, 1(3), 1900029.
---
En/Speed Up
speed up objective function
Codes in this section see example_function_modes.py, example_method_modes.py
To boost speed performance, scikit-opt supports 3 ways to speed up the objective function: vectorization, parallel, cached
- Vectorization requires that the objective function support vectorization. If so, the vectorization will gain extreme performance.
- multithreading requires nothing. It is usually faster than the common way, better than multiprocessing in io-intensive function
- multiprocessing requires nothing. It is usually faster than the common way, better than multithreading in io-intensive function
- Cached cache all the inputs and outputs. If the input in the next call is already in the cache, this method will fetch out the corresponding output from the cache, instead of calling the function once again. Cached gains extreme performance if the number of possible inputs is not big, such as integer programming or TSP
Totally speaking, vectorization is much faster then parallel, which is faster then common. If the number of input is not big, cached is much better then other ways.
To compare the speed of common, vectorization, parallel:
see /examples/example_function_modes.py
/ Detailed source-code truncated for AI context efficiency. /output:
>on io_costly task,use common mode, costs 5.116588s
on io_costly task,use multithreading mode, costs 3.113499s
on io_costly task,use multiprocessing mode, costs 3.119855s
on io_costly task,use vectorization mode, costs 0.604762s
on cpu_costly task,use common mode, costs 1.625032s
on cpu_costly task,use multithreading mode, costs 1.60131s
on cpu_costly task,use multiprocessing mode, costs 1.673792s
on cpu_costly task,use vectorization mode, costs 0.192595s
To compare the speed of common and cached:
def obj_func4_1(p):
time.sleep(0.1) # say that this function is very complicated and cost 0.1 seconds to run
x1, x2 = p
x = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(x)) - 0.5) / np.square(1 + 0.001 * x)
def obj_func4_2(p):
time.sleep(0.1) # say that this function is very complicated and cost 0.1 seconds to run
x1, x2 = p
x = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(x)) - 0.5) / np.square(1 + 0.001 * x)
set_run_mode(obj_func4_2, 'cached')
ga4_1 = GA(func=obj_func4_1, n_dim=2, size_pop=6, max_iter=10, lb=[-2, -2], ub=[2, 2], precision=1)
ga4_2 = GA(func=obj_func4_2, n_dim=2, size_pop=6, max_iter=10, lb=[-2, -2], ub=[2, 2], precision=1)
start_time = datetime.datetime.now()
best_x, best_y = ga4_1.run()
print('common mode, time costs: ', (datetime.datetime.now() - start_time).total_seconds())
start_time = datetime.datetime.now()
best_x, best_y = ga4_2.run()
print('cache mode, time costs: ', (datetime.datetime.now() - start_time).total_seconds())
print('cache mode, time costs: ', (datetime.datetime.now() - start_time).total_seconds())
output:
>on io_costly task,use common mode, costs 6.120317s
on io_costly task,use cached mode, costs 1.106842s
on cpu_costly task,use common mode, costs 1.914744s
on cpu_costly task,use cached mode, costs 0.222713s
---
Zh/ Coverpage
scikit-opt
群体智能算法
* [](https://pypi.org/project/scikit-opt/)
[](https://travis-ci.com/guofei9987/scikit-opt)
[](https://codecov.io/gh/guofei9987/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/blob/master/LICENSE)
[](https://pepy.tech/project/scikit-opt)
* 差分进化算法，遗传算法
* 粒子群算法，模拟退火算法
* 蚁群算法，免疫优化算法
* 鱼群算法
---
Zh/ Sidebar
* 文档
* 参数说明
* 更多遗传算法
* 更多粒子群算法
* 更多模拟退火算法
* 遗传算法做曲线拟合
* 提升速度
---
Zh/Args
入参一览
可以使用类似 help(GA), GA? 查看详细介绍，例如：
import skohelp(sko.GA.GA)
help(sko.GA.GA_TSP)
help(sko.PSO.PSO)
help(sko.DE.DE)
help(sko.SA.SA)
help(sko.SA.SA_TSP)
help(sko.ACA.ACA_TSP)
help(sko.IA.IA_TSP)
help(sko.AFSA.AFSA)
GA
| 入参 | 默认值 | 意义 |
|-----------------|--------|------------------------|
| func | \- | 目标函数 |
| n\_dim | \- | 目标函数的维度 |
| size\_pop | 50 | 种群规模 |
| max\_iter | 200 | 最大迭代次数 |
| prob\_mut | 0\.001 | 变异概率 |
| lb | \-1 | 每个自变量的最小值 |
| ub | 1 | 每个自变量的最大值 |
| constraint\_eq | 空元组 | 等式约束 |
| constraint\_ueq | 空元组 | 不等式约束 |
| precision | 1e\-7 | 精准度，int/float或者它们组成的列表 |
| n_processes | 0 | 使用的核数, 0表示使用全部的核 |
GA_TSP
| 入参 | 默认值 | 意义 |
|-----------|--------|--------|
| func | \- | 目标函数 |
| n\_dim | \- | 城市个数 |
| size\_pop | 50 | 种群规模 |
| max\_iter | 200 | 最大迭代次数 |
| prob\_mut | 0\.001 | 变异概率 |
PSO
| 入参 | 默认值 | 意义 |
|-----------|------|----------|
| func | \- | 目标函数 |
| n\_dim | \- | 目标函数的维度 |
| size\_pop | 50 | 种群规模 |
| max\_iter | 200 | 最大迭代次数 |
| lb | None | 每个参数的最小值 |
| ub | None | 每个参数的最大值 |
| w | 0\.8 | 惯性权重 |
| c1 | 0\.5 | 个体记忆 |
| c2 | 0\.5 | 集体记忆 |
| constraint\_ueq | 空元组 | 不等式约束 |
| n_processes | 0 | 使用的核数, 0表示使用全部的核 |
DE
| 入参 | 默认值 | 意义 |
|-----------|--------|----------|
| func | \- | 目标函数 |
| n\_dim | \- | 目标函数的维度 |
| size\_pop | 50 | 种群规模 |
| max\_iter | 200 | 最大迭代次数 |
| prob\_mut | 0\.001 | 变异概率 |
| F | 0\.5 | 变异系数 |
| lb | \-1 | 每个自变量的最小值 |
| ub | 1 | 每个自变量的最大值 |
| constraint\_eq | 空元组 | 等式约束 |
| constraint\_ueq | 空元组 | 不等式约束 |
| n_processes | 0 | 使用的核数, 0表示使用全部的核 |
SA
SA and SAFast and SABoltzmann and SACauchy
| 入参 | 默认值 | 意义 |
|--------------------|-------|-------|
| func | \- | 目标函数 |
| x0 | \- | 迭代初始点 |
| T\_max | 100 | 最大温度 |
| T\_min | 1e\-7 | 最小温度 |
| L | 300 | 链长 |
| max\_stay\_counter | 150 | 冷却耗时 |
| lb | | 每个自变量的最小值 |
| ub | | 每个自变量的最大值 |
SA_TSP
| 入参 | 默认值 | 意义 |
|--------------------|-------|-------|
| func | \- | 目标函数 |
| x0 | \- | 迭代初始点 |
| T\_max | 100 | 最大温度 |
| T\_min | 1e\-7 | 最小温度 |
| L | 300 | 链长 |
| max\_stay\_counter | 150 | 冷却耗时 |
ACA_TSP
| 入参 | 默认值 | 意义 |
|------------------|------|----------------------|
| func | \- | 目标函数 |
| n\_dim | \- | 城市个数 |
| size\_pop | 10 | 蚂蚁数量 |
| max\_iter | 20 | 最大迭代次数 |
| distance\_matrix | \- | 城市之间的距离矩阵，用于计算信息素的挥发 |
| alpha | 1 | 信息素重要程度 |
| beta | 2 | 适应度的重要程度 |
| rho | 0\.1 | 信息素挥发速度 |
IA_TSP
| 入参 | 默认值 | 意义 |
|-----------|--------|----------------------------------|
| func | \- | 目标函数 |
| n\_dim | \- | 城市个数 |
| size\_pop | 50 | 种群规模 |
| max\_iter | 200 | 最大迭代次数 |
| prob\_mut | 0\.001 | 变异概率 |
| T | 0\.7 | 抗体与抗体之间的亲和度阈值，大于这个阈值认为亲和，否则认为不亲和 |
| alpha | 0\.95 | 多样性评价指数，也就是抗体和抗原的重要性/抗体浓度重要性 |
AFSA
| 入参 | 默认值 | 意义 |
|---------------|-------|------------------|
| func | \- | 目标函数 |
| n\_dim | \- | 目标函数的维度 |
| size\_pop | 50 | 种群规模 |
| max\_iter | 300 | 最大迭代次数 |
| max\_try\_num | 100 | 最大尝试捕食次数 |
| step | 0\.5 | 每一步的最大位移比例 |
| visual | 0\.3 | 鱼的最大感知范围 |
| q | 0\.98 | 鱼的感知范围衰减系数 |
| delta | 0\.5 | 拥挤度阈值，越大越容易聚群和追尾 |
输出一览
GA&GA_TSP
- ga.generation_best_Y 每一代的最优函数值
- ga.generation_best_X 每一代的最优函数值对应的输入值
- ga.all_history_FitV 每一代的每个个体的适应度
- ga.all_history_Y 每一代每个个体的函数值
- ga.best_y 最优函数值
- ga.best_x 最优函数值对应的输入值
DE
- de.generation_best_Y 每一代的最优函数值
- de.generation_best_X 每一代的最优函数值对应的输入值
- de.all_history_Y 每一代每个个体的函数值
- de.best_y 最优函数值
- de.best_x 最优函数值对应的输入值
PSO
-
pso.record_value 每一代的粒子位置、粒子速度、对应的函数值。pso.record_mode = True 才开启记录-
pso.gbest_y_hist 历史最优函数值-
pso.best_y 最优函数值 （迭代中使用的是 pso.gbest_x, pso.gbest_y）-
pso.best_x 最优函数值对应的输入值SA
- de.generation_best_Y 每一代的最优函数值
- de.generation_best_X 每一代的最优函数值对应的输入值
- sa.best_x 最优函数值
- sa.best_y 最优函数值对应的输入值
ACA
-
de.generation_best_Y 每一代的最优函数值-
de.generation_best_X 每一代的最优函数值对应的输入值-
aca.best_y 最优函数值-
aca.best_x 最优函数值对应的输入值AFSA
-
afsa.best_x 最优函数值-
afsa.best_y 最优函数值对应的输入值IA
- ia.generation_best_Y 每一代的最优函数值
- ia.generation_best_X 每一代的最优函数值对应的输入值
- ia.all_history_FitV 每一代的每个个体的适应度
- ia.all_history_Y 每一代每个个体的函数值
- ia.best_y 最优函数值
- ia.best_x 最优函数值对应的输入值
---
Zh/Curve Fitting
使用遗传算法进行曲线拟合
随机生成训练数据
import numpy as np
import matplotlib.pyplot as plt
from sko.GA import GAx_true = np.linspace(-1.2, 1.2, 30)
y_true = x_true 3 - x_true + 0.4 * np.random.rand(30)
plt.plot(x_true, y_true, 'o')
构造残差
def f_fun(x, a, b, c, d):
return a x 3 + b x 2 + c * x + d
def obj_fun(p):
a, b, c, d = p
residuals = np.square(f_fun(x_true, a, b, c, d) - y_true).sum()
return residuals
使用 scikit-opt 做最优化
ga = GA(func=obj_fun, n_dim=4, size_pop=100, max_iter=500,
lb=[-2] 4, ub=[2] 4)best_params, residuals = ga.run()
print('best_x:', best_params, '\n', 'best_y:', residuals)
画出拟合效果图
y_predict = f_fun(x_true, *best_params)fig, ax = plt.subplots()
ax.plot(x_true, y_true, 'o')
ax.plot(x_true, y_predict, '-')
plt.show()
---
Zh/More Ga
遗传算法进行整数规划
在多维优化时，想让哪个变量限制为整数，就设定 precision 为 整数 即可。
例如，我想让我的自定义函数 demo_func 的某些变量限制为整数+浮点数（分别是隔2个，隔1个，浮点数），那么就设定 precision=[2, 1, 1e-7]
例子如下：
from sko.GA import GAdemo_func = lambda x: (x[0] - 1) 2 + (x[1] - 0.05) 2 + x[2] 2
ga = GA(func=demo_func, n_dim=3, max_iter=500, lb=[-1, -1, -1], ub=[5, 1, 1], precision=[2, 1, 1e-7])
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
说明：
- 当 precision 为整数时，对应的自变量会启用整数规划模式。
- 在整数规划模式下，变量的取值可能个数最好是 $2^n$，这样收敛速度快，效果好。
- 如果 precision 不是整数（例如是0.5）,则不会进入整数规划模式，如果还想用这个模式，那么把对应自变量乘以2，这样 precision 就是整数了。
遗传TSP问题如何固定起点和终点？
固定起点和终点要求路径不闭合（因为如果路径是闭合的，固定与不固定结果实际上是一样的）
假设你的起点和终点坐标指定为(0, 0) 和 (1, 1)，这样构建目标函数
- 起点和终点不参与优化。假设共有n+2个点，优化对象是中间n个点
- 目标函数（总距离）按实际去写。
import numpy as np
from scipy import spatial
import matplotlib.pyplot as pltnum_points = 20
points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points
start_point=[[0,0]]
end_point=[[1,1]]
points_coordinate=np.concatenate([points_coordinate,start_point,end_point])
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
num_points, = routine.shape
# start_point,end_point 本身不参与优化。给一个固定的值，参与计算总路径
routine = np.concatenate([[num_points], routine, [num_points+1]])
return sum([distance_matrix[routine[i], routine[i + 1]] for i in range(num_points+2-1)])
正常运行并画图：
from sko.GA import GA_TSPga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
best_points, best_distance = ga_tsp.run()
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([[num_points],best_points, [num_points+1]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
ax[1].plot(ga_tsp.generation_best_Y)
plt.show()
更多说明，这里
如何设定初始点或初始种群
- 对于遗传算法 GA, 运行 ga=GA(params) 生成模型后，赋值设定初始种群，例如 ga.Chrom = np.random.randint(0,2,size=(80,20))
- 对于差分进化算法 DE，设定 de.X 为初始 X.
- 对于模拟退火算法 SA，入参 x0 就是初始点.
- 对于粒子群算法 PSO，手动赋值 pso.X 为初始 X, 然后执行 pso.cal_y(); pso.update_gbest(); pso.update_pbest() 来更新历史最优点
---
Zh/More Pso
粒子群算法的动画展示
step1:做pso
-> Demo code: examples/demo_pso_ani.py#s1
import numpy as np
from sko.PSO import PSO
def demo_func(x):
x1, x2 = x
return -20 np.exp(-0.2 np.sqrt(0.5 * (x1 2 + x2 2))) - np.exp(
0.5 (np.cos(2 np.pi x1) + np.cos(2 np.pi * x2))) + 20 + np.e
constraint_ueq = (
lambda x: (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2
,
)
max_iter = 50
pso = PSO(func=demo_func, n_dim=2, pop=40, max_iter=max_iter, lb=[-2, -2], ub=[2, 2]
, constraint_ueq=constraint_ueq)
pso.record_mode = True
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
step2:画图
-> Demo code: examples/demo_pso_ani.py#s2
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimationrecord_value = pso.record_value
X_list, V_list = record_value['X'], record_value['V']
fig, ax = plt.subplots(1, 1)
ax.set_title('title', loc='center')
line = ax.plot([], [], 'b.')
X_grid, Y_grid = np.meshgrid(np.linspace(-2.0, 2.0, 40), np.linspace(-2.0, 2.0, 40))
Z_grid = demo_func((X_grid, Y_grid))
ax.contour(X_grid, Y_grid, Z_grid, 30)
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
t = np.linspace(0, 2 * np.pi, 40)
ax.plot(0.5 np.cos(t) + 1, 0.5 np.sin(t), color='r')
plt.ion()
p = plt.show()
def update_scatter(frame):
i, j = frame // 10, frame % 10
ax.set_title('iter = ' + str(i))
X_tmp = X_list[i] + V_list[i] * j / 10.0
plt.setp(line, 'xdata', X_tmp[:, 0], 'ydata', X_tmp[:, 1])
return line
ani = FuncAnimation(fig, update_scatter, blit=True, interval=25, frames=max_iter * 10)
plt.show()
ani.save('pso.gif', writer='pillow')
---
Zh/More Sa
3 types of Simulated Annealing
模拟退火有三种具体形式
fast:
u ~ Uniform(0, 1, size = d)
y = sgn(u - 0.5) T ((1 + 1/T)abs(2*u - 1) - 1.0)xc = y * (upper - lower)
x_new = x_old + xc
c = n exp(-n quench)
T_new = T0 exp(-c kquench)
cauchy:
u ~ Uniform(-pi/2, pi/2, size=d)
xc = learn_rate T tan(u)
x_new = x_old + xcT_new = T0 / (1 + k)
boltzmann:
std = minimum(sqrt(T) ones(d), (upper - lower) / (3learn_rate))
y ~ Normal(0, std, size = d)
x_new = x_old + learn_rate * yT_new = T0 / log(1 + k)
代码示例
#### 1.1 Fast Simulated Annealing
-> Demo code: examples/demo_sa.py#s4
from sko.SA import SAFastsa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_fast.run()
print('Fast Simulated Annealing: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)
#### 1.2 Fast Simulated Annealing with bounds
-> Demo code: examples/demo_sa.py#s5
from sko.SA import SAFastsa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
lb=[-1, 1, -1], ub=[2, 3, 4])
sa_fast.run()
print('Fast Simulated Annealing with bounds: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)
#### 2.1 Boltzmann Simulated Annealing
-> Demo code: examples/demo_sa.py#s6
from sko.SA import SABoltzmannsa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_boltzmann.run()
print('Boltzmann Simulated Annealing: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_boltzmann.best_y)
#### 2.2 Boltzmann Simulated Annealing with bounds
-> Demo code: examples/demo_sa.py#s7
from sko.SA import SABoltzmannsa_boltzmann = SABoltzmann(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
lb=-1, ub=[2, 3, 4])
sa_boltzmann.run()
print('Boltzmann Simulated Annealing with bounds: best_x is ', sa_boltzmann.best_x, 'best_y is ', sa_boltzmann.best_y)
#### 3. Cauchy Simulated Annealing
-> Demo code: examples/demo_sa.py#s8
from sko.SA import SACauchysa_cauchy = SACauchy(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150)
sa_cauchy.run()
print('Cauchy Simulated Annealing: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)
#### 3. Cauchy Simulated Annealing with bounds
-> Demo code: examples/demo_sa.py#s9
from sko.SA import SACauchysa_cauchy = SACauchy(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
lb=[-1, 1, -1], ub=[2, 3, 4])
sa_cauchy.run()
print('Cauchy Simulated Annealing with bounds: best_x is ', sa_cauchy.best_x, 'best_y is ', sa_cauchy.best_y)
---
Zh/README
scikit-opt
[](https://pypi.org/project/scikit-opt/)
[](https://github.com/guofei9987/scikit-opt)
[](https://travis-ci.com/guofei9987/scikit-opt)
[](https://codecov.io/gh/guofei9987/scikit-opt)
[](https://pepy.tech/project/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/stargazers)
[](https://github.com/guofei9987/scikit-opt/network/members)
[](https://gitter.im/guofei9987/scikit-opt?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)
一个封装了7种启发式算法的 Python 代码库
（差分进化算法、遗传算法、粒子群算法、模拟退火算法、蚁群算法、鱼群算法、免疫优化算法）
安装
pip install scikit-opt或者直接把源代码中的
sko 文件夹下载下来放本地也调用可以特性
特性1：UDF（用户自定义算子）
举例来说，你想出一种新的“选择算子”，如下
-> Demo code: examples/demo_ga_udf.py#s1
step1: define your own operator:
def selection_tournament(algorithm, tourn_size):
FitV = algorithm.FitV
sel_index = []
for i in range(algorithm.size_pop):
aspirants_index = np.random.choice(range(algorithm.size_pop), size=tourn_size)
sel_index.append(max(aspirants_index, key=lambda i: FitV[i]))
algorithm.Chrom = algorithm.Chrom[sel_index, :] # next generation
return algorithm.Chrom导入包，并且创建遗传算法实例
-> Demo code: examples/demo_ga_udf.py#s2
import numpy as np
from sko.GA import GA, GA_TSPdemo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + (x[2] - 0.5) 2
ga = GA(func=demo_func, n_dim=3, size_pop=100, max_iter=500, prob_mut=0.001,
lb=[-1, -10, -5], ub=[2, 10, 2], precision=[1e-7, 1e-7, 1])
把你的算子注册到你创建好的遗传算法实例上
-> Demo code: examples/demo_ga_udf.py#s3
ga.register(operator_name='selection', operator=selection_tournament, tourn_size=3)scikit-opt 也提供了十几个算子供你调用
-> Demo code: examples/demo_ga_udf.py#s4
from sko.operators import ranking, selection, crossover, mutationga.register(operator_name='ranking', operator=ranking.ranking). \
register(operator_name='crossover', operator=crossover.crossover_2point). \
register(operator_name='mutation', operator=mutation.mutation)
做遗传算法运算
-> Demo code: examples/demo_ga_udf.py#s5
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)现在 udf 支持遗传算法的这几个算子：crossover,mutation,selection,ranking
Scikit-opt 也提供了十来个算子，参考这里
提供一个面向对象风格的自定义算子的方法，供进阶用户使用:
-> Demo code: examples/demo_ga_udf.py#s6
class MyGA(GA):
def selection(self, tourn_size=3):
FitV = self.FitV
sel_index = []
for i in range(self.size_pop):
aspirants_index = np.random.choice(range(self.size_pop), size=tourn_size)
sel_index.append(max(aspirants_index, key=lambda i: FitV[i]))
self.Chrom = self.Chrom[sel_index, :] # next generation
return self.Chrom ranking = ranking.ranking
demo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + (x[2] - 0.5) 2
my_ga = MyGA(func=demo_func, n_dim=3, size_pop=100, max_iter=500, lb=[-1, -10, -5], ub=[2, 10, 2],
precision=[1e-7, 1e-7, 1])
best_x, best_y = my_ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
特性2：断点继续运行
例如，先跑10代，然后在此基础上再跑20代，可以这么写：
from sko.GA import GAfunc = lambda x: x[0] 2
ga = GA(func=func, n_dim=1)
ga.run(10)
ga.run(20)
特性3：4种加速方法
- [x] 矢量化计算：vectorization
- [x] 多线程计算：multithreading，适用于 IO 密集型目标函数
- [x] 多进程计算：multiprocessing，适用于 CPU 密集型目标函数
- [x] 缓存化计算：cached，适用于目标函数的每次输入有大量重复
see https://github.com/guofei9987/scikit-opt/blob/master/examples/example_function_modes.py
特性4: GPU 加速
GPU加速功能还比较简单，将会在 1.0.0 版本大大完善。
有个 demo 已经可以在现版本运行了: https://github.com/guofei9987/scikit-opt/blob/master/examples/demo_ga_gpu.py
快速开始
1. 差分进化算法
Step1：定义你的问题，这个demo定义了有约束优化问题
-> Demo code: examples/demo_de.py#s1
'''
min f(x1, x2, x3) = x1^2 + x2^2 + x3^2
s.t.
x1*x2 >= 1
x1*x2 <= 5
x2 + x3 = 1
0 <= x1, x2, x3 <= 5
'''
def obj_func(p):
x1, x2, x3 = p
return x1 2 + x2 2 + x3 2
constraint_eq = [
lambda x: 1 - x[1] - x[2]
]
constraint_ueq = [
lambda x: 1 - x[0] * x[1],
lambda x: x[0] * x[1] - 5
]
Step2: 做差分进化算法
-> Demo code: examples/demo_de.py#s2
from sko.DE import DEde = DE(func=obj_func, n_dim=3, size_pop=50, max_iter=800, lb=[0, 0, 0], ub=[5, 5, 5],
constraint_eq=constraint_eq, constraint_ueq=constraint_ueq)
best_x, best_y = de.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
2. 遗传算法
第一步：定义你的问题
-> Demo code: examples/demo_ga.py#s1
import numpy as np
def schaffer(p):
'''
This function has plenty of local minimum, with strong shocks
global minimum at (0,0) with value 0
https://en.wikipedia.org/wiki/Test_functions_for_optimization
'''
x1, x2 = p
part1 = np.square(x1) - np.square(x2)
part2 = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(part1)) - 0.5) / np.square(1 + 0.001 * part2)
第二步：运行遗传算法
-> Demo code: examples/demo_ga.py#s2
from sko.GA import GAga = GA(func=schaffer, n_dim=2, size_pop=50, max_iter=800, prob_mut=0.001, lb=[-1, -1], ub=[1, 1], precision=1e-7)
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
第三步：用 matplotlib 画出结果
-> Demo code: examples/demo_ga.py#s3
import pandas as pd
import matplotlib.pyplot as pltY_history = pd.DataFrame(ga.all_history_Y)
fig, ax = plt.subplots(2, 1)
ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
Y_history.min(axis=1).cummin().plot(kind='line')
plt.show()
2.2 遗传算法用于旅行商问题
GA_TSP 针对TSP问题重载了 交叉(crossover)、变异(mutation) 两个算子第一步，定义问题。
这里作为demo，随机生成距离矩阵. 实战中从真实数据源中读取。
-> Demo code: examples/demo_ga_tsp.py#s1
import numpy as np
from scipy import spatial
import matplotlib.pyplot as pltnum_points = 50
points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
num_points, = routine.shape
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
第二步，调用遗传算法进行求解
-> Demo code: examples/demo_ga_tsp.py#s2
from sko.GA import GA_TSPga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
best_points, best_distance = ga_tsp.run()
第三步，画出结果：
-> Demo code: examples/demo_ga_tsp.py#s3
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
ax[1].plot(ga_tsp.generation_best_Y)
plt.show()3. 粒子群算法
(PSO, Particle swarm optimization)
3.1 粒子群算法
第一步，定义问题
-> Demo code: examples/demo_pso.py#s1
def demo_func(x):
x1, x2, x3 = x
return x1 2 + (x2 - 0.05) 2 + x3 2第二步，做粒子群算法
-> Demo code: examples/demo_pso.py#s2
from sko.PSO import PSOpso = PSO(func=demo_func, n_dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
第三步，画出结果
-> Demo code: examples/demo_pso.py#s3
import matplotlib.pyplot as pltplt.plot(pso.gbest_y_hist)
plt.show()
↑see examples/demo_pso.py
3.2 带非线性约束的粒子群算法
3.2 PSO with nonlinear constraint
假如你的非线性约束是个圆内的面积 (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2<=0
这样写代码:
constraint_ueq = (
lambda x: (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2
,
)
pso = PSO(func=demo_func, n_dim=2, pop=40, max_iter=max_iter, lb=[-2, -2], ub=[2, 2]
, constraint_ueq=constraint_ueq)可以有多个非线性约束，向
constraint_ueq 加就行了。4. 模拟退火算法
(SA, Simulated Annealing)
4.1 模拟退火算法用于多元函数优化
第一步：定义问题
-> Demo code: examples/demo_sa.py#s1
demo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + x[2] 2第二步，运行模拟退火算法
-> Demo code: examples/demo_sa.py#s2
from sko.SA import SAsa = SA(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150)
best_x, best_y = sa.run()
print('best_x:', best_x, 'best_y', best_y)
第三步，画出结果
-> Demo code: examples/demo_sa.py#s3
import matplotlib.pyplot as plt
import pandas as pdplt.plot(pd.DataFrame(sa.best_y_history).cummin(axis=0))
plt.show()
另外，scikit-opt 还提供了三种模拟退火流派: Fast, Boltzmann, Cauchy. 更多参见 more sa
4.2 模拟退火算法解决TSP问题（旅行商问题）
第一步，定义问题。（我猜你已经无聊了，所以不黏贴这一步了）
第二步，调用模拟退火算法
-> Demo code: examples/demo_sa_tsp.py#s2
from sko.SA import SA_TSPsa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points)
best_points, best_distance = sa_tsp.run()
print(best_points, best_distance, cal_total_distance(best_points))
第三步，画出结果
-> Demo code: examples/demo_sa_tsp.py#s3
from matplotlib.ticker import FormatStrFormatterfig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(sa_tsp.best_y_history)
ax[0].set_xlabel("Iteration")
ax[0].set_ylabel("Distance")
ax[1].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1],
marker='o', markerfacecolor='b', color='c', linestyle='-')
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].set_xlabel("Longitude")
ax[1].set_ylabel("Latitude")
plt.show()
咱还有个动画
↑参考代码 examples/demo_sa_tsp.py
5. 蚁群算法
蚁群算法(ACA, Ant Colony Algorithm)解决TSP问题
-> Demo code: examples/demo_aca_tsp.py#s2
from sko.ACA import ACA_TSPaca = ACA_TSP(func=cal_total_distance, n_dim=num_points,
size_pop=50, max_iter=200,
distance_matrix=distance_matrix)
best_x, best_y = aca.run()
6. 免疫优化算法
(immune algorithm, IA)
-> Demo code: examples/demo_ia.py#s2
from sko.IA import IA_TSPia_tsp = IA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=500, max_iter=800, prob_mut=0.2,
T=0.7, alpha=0.95)
best_points, best_distance = ia_tsp.run()
print('best routine:', best_points, 'best_distance:', best_distance)
7. 人工鱼群算法
人工鱼群算法(artificial fish swarm algorithm, AFSA)
-> Demo code: examples/demo_afsa.py#s1
def func(x):
x1, x2 = x
return 1 / x1 2 + x1 2 + 1 / x2 2 + x2 2
from sko.AFSA import AFSA
afsa = AFSA(func, n_dim=2, size_pop=50, max_iter=300,
max_try_num=100, step=0.5, visual=0.3,
q=0.98, delta=0.5)
best_x, best_y = afsa.run()
print(best_x, best_y)
使用本项目的论文
- Yu, J., He, Y., Yan, Q., & Kang, X. (2021). SpecView: Malware Spectrum Visualization Framework With Singular Spectrum Transformation. IEEE Transactions on Information Forensics and Security, 16, 5093-5107.
- Zhen, H., Zhai, H., Ma, W., Zhao, L., Weng, Y., Xu, Y., ... & He, X. (2021). Design and tests of reinforcement-learning-based optimal power flow solution generator. Energy Reports.
- Heinrich, K., Zschech, P., Janiesch, C., & Bonin, M. (2021). Process data properties matter: Introducing gated convolutional neural networks (GCNN) and key-value-predict attention networks (KVP) for next event prediction with deep learning. Decision Support Systems, 143, 113494.
- Tang, H. K., & Goh, S. K. (2021). A Novel Non-population-based Meta-heuristic Optimizer Inspired by the Philosophy of Yi Jing. arXiv preprint arXiv:2104.08564.
- Wu, G., Li, L., Li, X., Chen, Y., Chen, Z., Qiao, B., ... & Xia, L. (2021). Graph embedding based real-time social event matching for EBSNs recommendation. World Wide Web, 1-22.
- Pan, X., Zhang, Z., Zhang, H., Wen, Z., Ye, W., Yang, Y., ... & Zhao, X. (2021). A fast and robust mixture gases identification and concentration detection algorithm based on attention mechanism equipped recurrent neural network with double loss function. Sensors and Actuators B: Chemical, 342, 129982.
- Castella Balcell, M. (2021). Optimization of the station keeping system for the WindCrete floating offshore wind turbine.
- Zhai, B., Wang, Y., Wang, W., & Wu, B. (2021). Optimal Variable Speed Limit Control Strategy on Freeway Segments under Fog Conditions. arXiv preprint arXiv:2107.14406.
- Yap, X. H. (2021). Multi-label classification on locally-linear data: Application to chemical toxicity prediction.
- Gebhard, L. (2021). Expansion Planning of Low-Voltage Grids Using Ant Colony Optimization Ausbauplanung von Niederspannungsnetzen mithilfe eines Ameisenalgorithmus.
- Ma, X., Zhou, H., & Li, Z. (2021). Optimal Design for Interdependencies between Hydrogen and Power Systems. IEEE Transactions on Industry Applications.
- de Curso, T. D. C. (2021). Estudo do modelo Johansen-Ledoit-Sornette de bolhas financeiras.
- Wu, T., Liu, J., Liu, J., Huang, Z., Wu, H., Zhang, C., ... & Zhang, G. (2021). A Novel AI-based Framework for AoI-optimal Trajectory Planning in UAV-assisted Wireless Sensor Networks. IEEE Transactions on Wireless Communications.
- Liu, H., Wen, Z., & Cai, W. (2021, August). FastPSO: Towards Efficient Swarm Intelligence Algorithm on GPUs. In 50th International Conference on Parallel Processing (pp. 1-10).
- Mahbub, R. (2020). Algorithms and Optimization Techniques for Solving TSP.
- Li, J., Chen, T., Lim, K., Chen, L., Khan, S. A., Xie, J., & Wang, X. (2019). Deep learning accelerated gold nanocluster synthesis. Advanced Intelligent Systems, 1(3), 1900029.
---
Zh/Speed Up
目标函数加速
本章节代码见于 example_function_modes.py, example_method_modes.py
为了提升速度，scikit-opt 支持3种提升速度的方案：矢量化，并行化，缓存化
- 矢量化：要求目标函数本身支持矢量化运算（详见代码）。矢量化运算拥有极高的性能，通常比并行化运算要快。算法中，每代对应1次矢量化运算
- 多线程：对目标函数没什么要求，通常比一般运算要快。如果目标函数是 IO 密集型，能达到更优的性能
- 多进程：对目标函数没什么要求，通常比一般运算要快。如果目标函数是 CPU 密集型，能达到更优的性能
- 缓存化：把每次计算的输入和输出缓存下来，下次调用时，如果已经缓存中已经存在，那么直接取出结果，而不再调用。缓存化特别适用于输入值有限的情况，例如纯整数规划、迭代到后期的TSP问题等。
总的来说，性能上，矢量化 远远大于 多线程/多进程 大于 不加速，如果是输入值得可能个数有限，缓存化 远大于其他方案。
下面比较 不加速、矢量化、多线程、多进程 的性能：
see /examples/example_function_modes.py
/ Detailed source-code truncated for AI context efficiency. /output:
>>on io_costly task,use common mode, costs 5.116588s
on io_costly task,use multithreading mode, costs 3.113499s
on io_costly task,use multiprocessing mode, costs 3.119855s
on io_costly task,use vectorization mode, costs 0.604762s
on cpu_costly task,use common mode, costs 1.625032s
on cpu_costly task,use multithreading mode, costs 1.60131s
on cpu_costly task,use multiprocessing mode, costs 1.673792s
on cpu_costly task,use vectorization mode, costs 0.192595s
下面比较 不加速 和 缓存化 的性能
def obj_func4_1(p):
time.sleep(0.1) # say that this function is very complicated and cost 0.1 seconds to run
x1, x2 = p
x = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(x)) - 0.5) / np.square(1 + 0.001 * x)
def obj_func4_2(p):
time.sleep(0.1) # say that this function is very complicated and cost 0.1 seconds to run
x1, x2 = p
x = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(x)) - 0.5) / np.square(1 + 0.001 * x)
set_run_mode(obj_func4_2, 'cached')
ga4_1 = GA(func=obj_func4_1, n_dim=2, size_pop=6, max_iter=10, lb=[-2, -2], ub=[2, 2], precision=1)
ga4_2 = GA(func=obj_func4_2, n_dim=2, size_pop=6, max_iter=10, lb=[-2, -2], ub=[2, 2], precision=1)
start_time = datetime.datetime.now()
best_x, best_y = ga4_1.run()
print('common mode, time costs: ', (datetime.datetime.now() - start_time).total_seconds())
start_time = datetime.datetime.now()
best_x, best_y = ga4_2.run()
print('cache mode, time costs: ', (datetime.datetime.now() - start_time).total_seconds())
print('cache mode, time costs: ', (datetime.datetime.now() - start_time).total_seconds())
output:
>on io_costly task,use common mode, costs 6.120317s
on io_costly task,use cached mode, costs 1.106842s
on cpu_costly task,use common mode, costs 1.914744s
on cpu_costly task,use cached mode, costs 0.222713s
算子优化加速
主要手段是 矢量化 和 逻辑化
对于可以矢量化的算子，scikit-opt 都尽量做了矢量化，并且默认调用矢量化的算子，且 无须用户额外操作。
另外，考虑到有些算子矢量化后，代码可读性下降，因此矢量化前的算子也会保留，为用户进阶学习提供方便。
0/1 基因的mutation
做一个mask，是一个与
Chrom 大小一致的0/1矩阵，如果值为1，那么对应位置进行变异（0变1或1变0） 自然想到用整除2的方式进行
def mutation(self):
# mutation of 0/1 type chromosome
mask = (np.random.rand(self.size_pop, self.len_chrom) < self.prob_mut) * 1
self.Chrom = (mask + self.Chrom) % 2
return self.Chrom如此就实现了一次性对整个种群所有基因变异的矢量化运算。用pycharm的profile功能试了一下，效果良好
再次改进。我还嫌求余数这一步速度慢，画一个真值表
|A|mask：是否变异|A变异后|
|--|--|--|
|1|0|1|
|0|0|0|
|1|1|0|
|0|1|1|
发现这就是一个 异或
def mutation2(self):
mask = (np.random.rand(self.size_pop, self.len_chrom) < self.prob_mut)
self.Chrom ^= mask
return self.Chrom测试发现运行速度又快了1~3倍，与最原始的双层循环相比，快了约20倍。
0/1基因的crossover
同样思路，试试crossover.
- mask同样，1表示对应点交叉，0表示对应点不交叉
做一个真值表，总共8种可能，发现其中只有2种可能基因有变化（等位基因一样时，交叉后的结果与交叉前一样）
|A基因|B基因|是否交叉|交叉后的A基因|交叉后的B基因|
|--|--|--|--|--|
|1|0|1|0|1|
|0|1|1|1|0|
可以用 异或 和 且 来表示是否变化的表达式: mask = (A^B)&C，然后可以计算了A^=mask, B^=mask
代码实现
def crossover_2point_bit(self):
Chrom, size_pop, len_chrom = self.Chrom, self.size_pop, self.len_chrom
Chrom1, Chrom2 = Chrom[::2], Chrom[1::2]
mask = np.zeros(shape=(int(size_pop / 2),len_chrom),dtype=int)
for i in range(int(size_pop / 2)):
n1, n2 = np.random.randint(0, self.len_chrom, 2)
if n1 > n2:
n1, n2 = n2, n1
mask[i, n1:n2] = 1
mask2 = (Chrom1 ^ Chrom2) & mask
Chrom1^=mask2
Chrom2^=mask2
Chrom[::2], Chrom[1::2]=Chrom1,Chrom2
self.Chrom=Chrom
return self.Chrom测试结果，效率提升约1倍。
锦标赛选择算子selection_tournament
实战发现，selection_tournament 往往是最耗时的，几乎占用一半时间，因此需要优化。
优化前的算法是遍历，每次选择一组进行锦标赛。但可以在二维array上一次性操作。
def selection_tournament_faster(self, tourn_size=3):
'''
Select the best individual among tournsize randomly chosen
Same with selection_tournament but much faster using numpy
individuals,
:param self:
:param tourn_size:
:return:
'''
aspirants_idx = np.random.randint(self.size_pop, size=(self.size_pop, tourn_size))
aspirants_values = self.FitV[aspirants_idx]
winner = aspirants_values.argmax(axis=1) # winner index in every team
sel_index = [aspirants_idx[i, j] for i, j in enumerate(winner)]
self.Chrom = self.Chrom[sel_index, :]
return self.Chrom发现own time 和time 都降为原来的10%~15%，效率提升了约9倍
---
Coverpage
scikit-opt
Powerful Python module for Heuristic Algorithms
* Genetic Algorithm
* Particle Swarm Optimization
* Simulated Annealing
* Ant Colony Algorithm
* Immune Algorithm
* Artificial Fish Swarm Algorithm
---
CONTRIBUTING
Contributing guidelines
This page explains how you can contribute to the development of
scikit-opt by submitting patches, tests, new models, or examples.
scikit-opt is developed on
Github using the
Git version control system.
Submitting a Bug Report
- Include a short, self-contained code snippet that reproduces the
problem
- Ensure that the bug still exists on latest version.
Making Changes to the Code
For a pull request to be accepted, you must meet the below requirements.
This greatly helps in keeping the job of maintaining and releasing the
software a shared effort.
- One branch. One feature. Branches are cheap and github makes it
easy to merge and delete branches with a few clicks. Avoid the
temptation to lump in a bunch of unrelated changes when working on a
feature, if possible. This helps us keep track of what has changed
when preparing a release.
- Commit messages should be clear and concise. If your commit references or
closes a specific issue, you can close it by mentioning it in the
commit
message.
(For maintainers: These suggestions go for Merge commit comments
too. These are partially the record for release notes.)
- Each function, class, method, and attribute needs to be documented.
- If you are adding new functionality, you need to add it to the
documentation by editing (or creating) the appropriate file in
docs/.
How to Submit a Pull Request
So you want to submit a patch to scikit-opt but are not too familiar
with github? Here are the steps you need to take.
1. Fork the
scikit-opt repository
on Github.
2. Create a new feature
branch.
Each branch must be self-contained, with a single new feature or
bugfix.
3. Make sure the test suite passes. This includes testing on Python 3.
The easiest way to do this is to either enable
Travis-CI on your fork, or to make a pull
request and check there.
4. If it is a big, new feature please submit an example to path /examples/ or /tests/ and add it to .travis.yml.
5. Submit a pull
request
License
scikit-opt is released under the MIT license.
---
README
scikit-opt
[](https://pypi.org/project/scikit-opt/)
[](https://travis-ci.com/guofei9987/scikit-opt)
[](https://codecov.io/gh/guofei9987/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/blob/master/LICENSE)
[](https://github.com/guofei9987/scikit-opt/fork)
[](https://pepy.tech/project/scikit-opt)
[](https://github.com/guofei9987/scikit-opt/discussions)
<a href="https://hellogithub.com/repository/guofei9987/scikit-opt" target="_blank"><img src="https://abroad.hellogithub.com/v1/widgets/recommend.svg?rid=6763d615842e4449a02f024f3e2e345c&claim_uid=se0WHo8cbiLv2w1&theme=small" alt="Featured｜HelloGitHub" /></a>
Swarm Intelligence in Python
(Genetic Algorithm, Particle Swarm Optimization, Simulated Annealing, Ant Colony Algorithm, Immune Algorithm, Artificial Fish Swarm Algorithm in Python)
- Documentation: https://scikit-opt.github.io/scikit-opt/#/en/
- 文档： https://scikit-opt.github.io/scikit-opt/#/zh/
- Source code: https://github.com/guofei9987/scikit-opt
- Help us improve scikit-opt https://www.wjx.cn/jq/50964691.aspx
install
pip install scikit-optFor the current developer version:
git clone [email protected]:guofei9987/scikit-opt.git
cd scikit-opt
pip install .Features
Feature1: UDF
UDF (user defined function) is available now!
For example, you just worked out a new type of selection function.
Now, your selection function is like this:
-> Demo code: examples/demo_ga_udf.py#s1
step1: define your own operator:
def selection_tournament(algorithm, tourn_size):
FitV = algorithm.FitV
sel_index = []
for i in range(algorithm.size_pop):
aspirants_index = np.random.choice(range(algorithm.size_pop), size=tourn_size)
sel_index.append(max(aspirants_index, key=lambda i: FitV[i]))
algorithm.Chrom = algorithm.Chrom[sel_index, :] # next generation
return algorithm.ChromImport and build ga
-> Demo code: examples/demo_ga_udf.py#s2
import numpy as np
from sko.GA import GA, GA_TSPdemo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + (x[2] - 0.5) 2
ga = GA(func=demo_func, n_dim=3, size_pop=100, max_iter=500, prob_mut=0.001,
lb=[-1, -10, -5], ub=[2, 10, 2], precision=[1e-7, 1e-7, 1])
Regist your udf to GA
-> Demo code: examples/demo_ga_udf.py#s3
ga.register(operator_name='selection', operator=selection_tournament, tourn_size=3)scikit-opt also provide some operators
-> Demo code: examples/demo_ga_udf.py#s4
from sko.operators import ranking, selection, crossover, mutationga.register(operator_name='ranking', operator=ranking.ranking). \
register(operator_name='crossover', operator=crossover.crossover_2point). \
register(operator_name='mutation', operator=mutation.mutation)
Now do GA as usual
-> Demo code: examples/demo_ga_udf.py#s5
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)Until Now, the udf surportcrossover,mutation,selection,rankingof GA
scikit-opt provide a dozen of operators, see here
For advanced users:
-> Demo code: examples/demo_ga_udf.py#s6
class MyGA(GA):
def selection(self, tourn_size=3):
FitV = self.FitV
sel_index = []
for i in range(self.size_pop):
aspirants_index = np.random.choice(range(self.size_pop), size=tourn_size)
sel_index.append(max(aspirants_index, key=lambda i: FitV[i]))
self.Chrom = self.Chrom[sel_index, :] # next generation
return self.Chrom ranking = ranking.ranking
demo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + (x[2] - 0.5) 2
my_ga = MyGA(func=demo_func, n_dim=3, size_pop=100, max_iter=500, lb=[-1, -10, -5], ub=[2, 10, 2],
precision=[1e-7, 1e-7, 1])
best_x, best_y = my_ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
feature2: continue to run
(New in version 0.3.6)
Run an algorithm for 10 iterations, and then run another 20 iterations base on the 10 iterations before:
from sko.GA import GAfunc = lambda x: x[0] 2
ga = GA(func=func, n_dim=1)
ga.run(10)
ga.run(20)
feature3: 4-ways to accelerate
- vectorization
- multithreading
- multiprocessing
- cached
see https://github.com/guofei9987/scikit-opt/blob/master/examples/example_function_modes.py
feature4: GPU computation
We are developing GPU computation, which will be stable on version 1.0.0
An example is already available: https://github.com/guofei9987/scikit-opt/blob/master/examples/demo_ga_gpu.py
Quick start
1. Differential Evolution
Step1：define your problem
-> Demo code: examples/demo_de.py#s1
'''
min f(x1, x2, x3) = x1^2 + x2^2 + x3^2
s.t.
x1*x2 >= 1
x1*x2 <= 5
x2 + x3 = 1
0 <= x1, x2, x3 <= 5
'''
def obj_func(p):
x1, x2, x3 = p
return x1 2 + x2 2 + x3 2
constraint_eq = [
lambda x: 1 - x[1] - x[2]
]
constraint_ueq = [
lambda x: 1 - x[0] * x[1],
lambda x: x[0] * x[1] - 5
]
Step2: do Differential Evolution
-> Demo code: examples/demo_de.py#s2
from sko.DE import DEde = DE(func=obj_func, n_dim=3, size_pop=50, max_iter=800, lb=[0, 0, 0], ub=[5, 5, 5],
constraint_eq=constraint_eq, constraint_ueq=constraint_ueq)
best_x, best_y = de.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
2. Genetic Algorithm
Step1：define your problem
-> Demo code: examples/demo_ga.py#s1
import numpy as np
def schaffer(p):
'''
This function has plenty of local minimum, with strong shocks
global minimum at (0,0) with value 0
https://en.wikipedia.org/wiki/Test_functions_for_optimization
'''
x1, x2 = p
part1 = np.square(x1) - np.square(x2)
part2 = np.square(x1) + np.square(x2)
return 0.5 + (np.square(np.sin(part1)) - 0.5) / np.square(1 + 0.001 * part2)
Step2: do Genetic Algorithm
-> Demo code: examples/demo_ga.py#s2
from sko.GA import GAga = GA(func=schaffer, n_dim=2, size_pop=50, max_iter=800, prob_mut=0.001, lb=[-1, -1], ub=[1, 1], precision=1e-7)
best_x, best_y = ga.run()
print('best_x:', best_x, '\n', 'best_y:', best_y)
-> Demo code: examples/demo_ga.py#s3
import pandas as pd
import matplotlib.pyplot as pltY_history = pd.DataFrame(ga.all_history_Y)
fig, ax = plt.subplots(2, 1)
ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
Y_history.min(axis=1).cummin().plot(kind='line')
plt.show()
2.2 Genetic Algorithm for TSP(Travelling Salesman Problem)
Just import the
GA_TSP, it overloads the crossover, mutation to solve the TSPStep1: define your problem. Prepare your points coordinate and the distance matrix.
Here I generate the data randomly as a demo:
-> Demo code: examples/demo_ga_tsp.py#s1
import numpy as np
from scipy import spatial
import matplotlib.pyplot as pltnum_points = 50
points_coordinate = np.random.rand(num_points, 2) # generate coordinate of points
distance_matrix = spatial.distance.cdist(points_coordinate, points_coordinate, metric='euclidean')
def cal_total_distance(routine):
'''The objective function. input routine, return total distance.
cal_total_distance(np.arange(num_points))
'''
num_points, = routine.shape
return sum([distance_matrix[routine[i % num_points], routine[(i + 1) % num_points]] for i in range(num_points)])
Step2: do GA
-> Demo code: examples/demo_ga_tsp.py#s2
from sko.GA import GA_TSPga_tsp = GA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=50, max_iter=500, prob_mut=1)
best_points, best_distance = ga_tsp.run()
Step3: Plot the result:
-> Demo code: examples/demo_ga_tsp.py#s3
fig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1], 'o-r')
ax[1].plot(ga_tsp.generation_best_Y)
plt.show()3. PSO(Particle swarm optimization)
3.1 PSO
Step1: define your problem:
-> Demo code: examples/demo_pso.py#s1
def demo_func(x):
x1, x2, x3 = x
return x1 2 + (x2 - 0.05) 2 + x3 2Step2: do PSO
-> Demo code: examples/demo_pso.py#s2
from sko.PSO import PSOpso = PSO(func=demo_func, n_dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
Step3: Plot the result
-> Demo code: examples/demo_pso.py#s3
import matplotlib.pyplot as pltplt.plot(pso.gbest_y_hist)
plt.show()
3.2 PSO with nonlinear constraint
If you need nolinear constraint like (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2<=0
Codes are like this:
constraint_ueq = (
lambda x: (x[0] - 1) 2 + (x[1] - 0) 2 - 0.5 2
,
)
pso = PSO(func=demo_func, n_dim=2, pop=40, max_iter=max_iter, lb=[-2, -2], ub=[2, 2]
, constraint_ueq=constraint_ueq)Note that, you can add more then one nonlinear constraint. Just add it to constraint_ueq
More over, we have an animation:
↑see examples/demo_pso_ani.py
4. SA(Simulated Annealing)
4.1 SA for multiple function
Step1: define your problem
-> Demo code: examples/demo_sa.py#s1
demo_func = lambda x: x[0] 2 + (x[1] - 0.05) 2 + x[2] 2Step2: do SA
-> Demo code: examples/demo_sa.py#s2
from sko.SA import SAsa = SA(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150)
best_x, best_y = sa.run()
print('best_x:', best_x, 'best_y', best_y)
Step3: Plot the result
-> Demo code: examples/demo_sa.py#s3
import matplotlib.pyplot as plt
import pandas as pdplt.plot(pd.DataFrame(sa.best_y_history).cummin(axis=0))
plt.show()
Moreover, scikit-opt provide 3 types of Simulated Annealing: Fast, Boltzmann, Cauchy. See more sa
4.2 SA for TSP
Step1: oh, yes, define your problems. To boring to copy this step.
Step2: DO SA for TSP
-> Demo code: examples/demo_sa_tsp.py#s2
from sko.SA import SA_TSPsa_tsp = SA_TSP(func=cal_total_distance, x0=range(num_points), T_max=100, T_min=1, L=10 * num_points)
best_points, best_distance = sa_tsp.run()
print(best_points, best_distance, cal_total_distance(best_points))
Step3: plot the result
-> Demo code: examples/demo_sa_tsp.py#s3
from matplotlib.ticker import FormatStrFormatterfig, ax = plt.subplots(1, 2)
best_points_ = np.concatenate([best_points, [best_points[0]]])
best_points_coordinate = points_coordinate[best_points_, :]
ax[0].plot(sa_tsp.best_y_history)
ax[0].set_xlabel("Iteration")
ax[0].set_ylabel("Distance")
ax[1].plot(best_points_coordinate[:, 0], best_points_coordinate[:, 1],
marker='o', markerfacecolor='b', color='c', linestyle='-')
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax[1].set_xlabel("Longitude")
ax[1].set_ylabel("Latitude")
plt.show()
More: Plot the animation:
5. ACA (Ant Colony Algorithm) for tsp
-> Demo code: examples/demo_aca_tsp.py#s2
from sko.ACA import ACA_TSPaca = ACA_TSP(func=cal_total_distance, n_dim=num_points,
size_pop=50, max_iter=200,
distance_matrix=distance_matrix)
best_x, best_y = aca.run()
6. immune algorithm (IA)
-> Demo code: examples/demo_ia.py#s2
from sko.IA import IA_TSPia_tsp = IA_TSP(func=cal_total_distance, n_dim=num_points, size_pop=500, max_iter=800, prob_mut=0.2,
T=0.7, alpha=0.95)
best_points, best_distance = ia_tsp.run()
print('best routine:', best_points, 'best_distance:', best_distance)
7. Artificial Fish Swarm Algorithm (AFSA)
-> Demo code: examples/demo_afsa.py#s1
def func(x):
x1, x2 = x
return 1 / x1 2 + x1 2 + 1 / x2 2 + x2 2
from sko.AFSA import AFSA
afsa = AFSA(func, n_dim=2, size_pop=50, max_iter=300,
max_try_num=100, step=0.5, visual=0.3,
q=0.98, delta=0.5)
best_x, best_y = afsa.run()
print(best_x, best_y)
Projects using scikit-opt
- Yu, J., He, Y., Yan, Q., & Kang, X. (2021). SpecView: Malware Spectrum Visualization Framework With Singular Spectrum Transformation. IEEE Transactions on Information Forensics and Security, 16, 5093-5107.
- Zhen, H., Zhai, H., Ma, W., Zhao, L., Weng, Y., Xu, Y., ... & He, X. (2021). Design and tests of reinforcement-learning-based optimal power flow solution generator. Energy Reports.
- Heinrich, K., Zschech, P., Janiesch, C., & Bonin, M. (2021). Process data properties matter: Introducing gated convolutional neural networks (GCNN) and key-value-predict attention networks (KVP) for next event prediction with deep learning. Decision Support Systems, 143, 113494.
- Tang, H. K., & Goh, S. K. (2021). A Novel Non-population-based Meta-heuristic Optimizer Inspired by the Philosophy of Yi Jing. arXiv preprint arXiv:2104.08564.
- Wu, G., Li, L., Li, X., Chen, Y., Chen, Z., Qiao, B., ... & Xia, L. (2021). Graph embedding based real-time social event matching for EBSNs recommendation. World Wide Web, 1-22.
- Pan, X., Zhang, Z., Zhang, H., Wen, Z., Ye, W., Yang, Y., ... & Zhao, X. (2021). A fast and robust mixture gases identification and concentration detection algorithm based on attention mechanism equipped recurrent neural network with double loss function. Sensors and Actuators B: Chemical, 342, 129982.
- Castella Balcell, M. (2021). Optimization of the station keeping system for the WindCrete floating offshore wind turbine.
- Zhai, B., Wang, Y., Wang, W., & Wu, B. (2021). Optimal Variable Speed Limit Control Strategy on Freeway Segments under Fog Conditions. arXiv preprint arXiv:2107.14406.
- Yap, X. H. (2021). Multi-label classification on locally-linear data: Application to chemical toxicity prediction.
- Gebhard, L. (2021). Expansion Planning of Low-Voltage Grids Using Ant Colony Optimization Ausbauplanung von Niederspannungsnetzen mithilfe eines Ameisenalgorithmus.
- Ma, X., Zhou, H., & Li, Z. (2021). Optimal Design for Interdependencies between Hydrogen and Power Systems. IEEE Transactions on Industry Applications.
- de Curso, T. D. C. (2021). Estudo do modelo Johansen-Ledoit-Sornette de bolhas financeiras.
- Wu, T., Liu, J., Liu, J., Huang, Z., Wu, H., Zhang, C., ... & Zhang, G. (2021). A Novel AI-based Framework for AoI-optimal Trajectory Planning in UAV-assisted Wireless Sensor Networks. IEEE Transactions on Wireless Communications.
- Liu, H., Wen, Z., & Cai, W. (2021, August). FastPSO: Towards Efficient Swarm Intelligence Algorithm on GPUs. In 50th International Conference on Parallel Processing (pp. 1-10).
- Mahbub, R. (2020). Algorithms and Optimization Techniques for Solving TSP.
- Li, J., Chen, T., Lim, K., Chen, L., Khan, S. A., Xie, J., & Wang, X. (2019). Deep learning accelerated gold nanocluster synthesis. Advanced Intelligent Systems, 1(3), 1900029.
---