Podręcznik
4. Particle Swarm Optimization
4.4. Przykład PSO
Będziemy korzystali z biblioteki PSO zawartej jako część biblioteki SciKit. Nie jest to jedyna implementacja PSO którą mamy dostępną - ale jest stabilna, powszechnie używana, i ciągle uzupełniana.
By można było z niej korzystać - należy ją najpierw zainstalować. Robimy to, korzystając z menager-a pakietów pythona (jeśi ćwiczenie wykonujecie poza Colab - to te polecenia muszą zostać wykonane w terminalu, ewentualnie możecie zainstalować bibliotekę korzystając z opcji środowiska programistycznego)
!pip install scikit-opt
Po zainstalowaniu bibliotek należy zaimportować potrzebne biblioteki (sko - to SciKit Optimization).
import matplotlib.pyplot as plt
import numpy as np
from sko.PSO import PSO
Pierwszym krokiem optymalizacji jest definicja funkcji celu. W przykładzie wykorzystamy prostą, trójwymiarową funkcję kwadratową
def demo_func(x):
x1, x2, x3 = x
return x1;**2+(x2-0.05)**2+ x3 ** 2
W drugim kroku uruchamiamy algorytm optymalizacyjny.
Parametry przekazywane do PSO:
- func - funkcja celu
- n_dim - liczba wymiarów funkcji celu
- size_pop - wielkość roju
- max_iter - liczba iteracji
- lb - dolne ograniczenie zmienności parametrów
- ub - górne ograniczenie zmienności parametrów
- w - inercja cząstki
- c1 - waga dążenia do najlepszej dotychczasowej pozycji (moja wiedza)
- c2 - waga dążenia do pozycji lidera roju (wiedza najlepszego)
- constraint_ueq - ograniczenia nierównościowe
pso = PSO(func=demo_func, n_dim=3, pop=40, max_iter=300, lb=[-100, -100, -100], ub=[100, 100, 100], w=0.8, c1=0.5, c2=0.5)
pso.run()
print('Znaleziony wektor parametrów ', pso.gbest_x, 'wartość funkcji', pso.gbest_y
W trzecim kroku - rysujemy jak zmieniała się wartość funkcji celu w trakcie obliczeń
import matplotlib.pyplot as plt
plt.plot(pso.gbest_y_hist)
plt.show()
Komplet parametrów
Teraz przygotujemy animację która pokaże jak cząstki zachowywały się w trakcie obliczeń. Wykorzystana zostanie znów standardowa wersja PSO zaimplementowana w SCIKIT.
Wykorzystamy do tego celu dwuwymiarową funkcję celu, i dodamy do niej ograniczenie przeszukiwania do okręgu:
import numpy as np
from sko.PSO import PSO
def demo_func_2(x):
x1, x2 = x
return 2 * x1 ** 2 - 1.05 * x1 ** 4 + (x1 ** 6) / 6 + x1 * x2 + x2 ** 2
# 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,
)
Narysujmy naszą funkcję celu z ograniczeniami
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
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_2((X_grid, Y_grid))
ax.contour(X_grid, Y_grid, Z_grid, 30)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
t = np.linspace(0, 2 * np.pi, 40)
ax.plot(2 * np.cos(t) + 1, 2 * np.sin(t), color='r')
plt.ion()
p = plt.show()
Wykonamy obliczenia - poszukajcie tutaj innych parametrów optymalizacji
max_iter = 50
pso = PSO(func=demo_func_2, 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('Znaleziony wektor parametrów ', pso.gbest_x, 'Wartość funkcji', pso.gbest_y)
W ostatnim kroku przygotujemy animację - jak wyglądały zmiany w roju w trakcie obliczeń
record_value = pso.record_value
X_list, V_list = record_value['X'], record_value['V']
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')