solution[1:]] sol_dists = numpy.append(sol_dists, (distances[solution[0], solution[-1]],)) sol_pens = penalties[solution[0:-1], solution[1:]] sol_pens = numpy.append(sol_pens, (penalties[solution[0], solution[-1]],)) utility = sol_dists/(sol_pens+1) index_a = numpy.argmax(utility) index_b = index_a+1 if index_a != count-1 else 0 a = solution[index_a] b = solution[index_b] penalties[a,b] += 1.0 penalties[b,a] += 1.0