scikit-opt

To enforce integer constraints on specific variables in a Genetic Algorithm, set the corresponding values in the precision parameter to an integer. For example, if you want the first variable to be an integer with interval 2, the second to be an integer with interval 1, and the third to be a float, set precision=[2, 1, 1e-7].

Note: Performance is optimized when the number of possible values for an integer precision is a power of 2 ($2^n$). If you need non-integer precision (e.g., 0.5), it is recommended to manually transform the variable (e.g., create a new variable and multiply it by 2).

from sko.GA import GA

demo_func = lambda x: (x[0] - 1) ** 2 + (x[1] - 0.05) ** 2 + x[2] ** 2
# precision=[2, 1, 1e-7] makes x[0] and x[1] integers
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)

You can accelerate the execution of your objective functions in scikit-opt by using sko.tools.set_run_mode. There are four primary modes available:

1. common: The default execution mode.
2. multithreading: Uses multiple threads. Generally faster than common and often better than multiprocessing for I/O-intensive functions.
3. multiprocessing: Uses multiple processes. Generally faster than common and often better than multithreading for I/O-intensive functions.
4. vectorization: Requires the objective function itself to be written to support vectorized inputs (e.g., accepting a 2D array of parameters). This provides extreme performance gains.
5. cached: Caches all inputs and outputs. If a subsequent call uses an input already in the cache, the cached output is returned instead of re-running the function. This is highly effective for problems with a limited number of possible inputs, such as integer programming or the Travelling Salesman Problem (TSP).

When solving the Travelling Salesman Problem (TSP) where the path must start and end at specific coordinates (and is not a simple cycle), you must incorporate these points into your objective function.

1. Include the start and end points in your coordinate set.
2. The objective function should accept the sequence of intermediate points as input.
3. Inside the objective function, concatenate the start point index and end point index to the input routine to calculate the total distance including the fixed points.

import numpy as np
from scipy import spatial
from sko.GA import GA_TSP

num_points = 20
points_coordinate = np.random.rand(num_points, 2)
start_point=[[0,0]]
end_point=[[1,1]]
# Combine points with start and end
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.'''
    num_points, = routine.shape
    # Concatenate start point index (num_points) and end point index (num_points+1)
    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)])

ga_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()

To visualize the optimization process, you can use the PSO class from sko.PSO and set record_mode = True. This enables the storage of particle positions (X) and velocities (V) at each iteration in the record_value attribute.

Key parameters for PSO:

• func: The objective function to minimize.
• n_dim: Number of dimensions.
• pop: Population size.
• max_iter: Maximum number of iterations.
• lb: Lower bounds.
• ub: Upper bounds.
• constraint_ueq: Equality constraints (lambda function).

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)

If pso.record_mode was set to True during the run() method, you can access pso.record_value['X'] (positions) and pso.record_value['V'] (velocities) to create an animation of the particles moving through the search space using matplotlib.animation.FuncAnimation.

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Assuming 'pso' was run with record_mode = True
record_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')