Pancake's Personal Website – 第 4 页

NSGA-ii 双目标遗传算法

论文原文《A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II》，可在IEEE上搜索。

简要叙述：编码交叉变异依然不变，在在此基础上加上了

快速非支配排序
拥挤度
精英策略

非支配排序

支配简述：对于一个点a，他的f1(a), f2(a)……fn(a)都比另一个点b，他的f1(b), f2(b)……fn(b)大，那么a支配b。
非支配即a他的所有值至少有一个小于b，且至少有一个大于b，那么，那么a和b非支配。

NSGA-ii 使用快速非支配排序，具体算法可见代码，根据代码思路自己去画一画。

非支配排序会最终将很多点分成不同rank的类，根据目标所需最大最小不同后续选择对rank不同的排序。

例子代码

import matplotlib.pyplot as plt
import numpy as np

P = np.array([[1, 9], [2, 7], [3, 3],
              [2.5, 10], [4, 4.5], [5, 1],
              [4.5, 7.5], [6, 4], [7, 2], [2, 4], [4, 2]])


def fast_non_dominated_sort(values1, values2):
    S = [[] for i in range(0, len(values1))]
    front = [[]]
    n = [0 for i in range(0, len(values1))]
    rank = [0 for i in range(0, len(values1))]

    for p in range(0, len(values1)):
        S[p] = []
        n[p] = 0
        for q in range(0, len(values1)):
            if (values1[p] > values1[q] and values2[p] > values2[q]) or (values1[p] >= values1[q] and values2[p] > values2[q]) or (
                    values1[p] > values1[q] and values2[p] >= values2[q]):
                if q not in S[p]:
                    S[p].append(q)
            elif (values1[q] > values1[p] and values2[q] > values2[p]) or (values1[q] >= values1[p] and values2[q] > values2[p]) or (
                    values1[q] > values1[p] and values2[q] >= values2[p]):
                n[p] = n[p] + 1
        if n[p] == 0:
            rank[p] = 0
            if p not in front[0]:
                front[0].append(p)
    print(S, n)
    i = 0
    while front[i] != []:
        Q = []
        for p in front[i]:
            for q in S[p]:
                n[q] = n[q] - 1
                if n[q] == 0:
                    rank[q] = i + 1
                    if q not in Q:
                        Q.append(q)
        i = i + 1
        front.append(Q)

    del front[len(front) - 1]
    return front


ans = fast_non_dominated_sort(P[:, 0], P[:, 1])

plt.plot([P[i, 0] for i in ans[0]], [P[i, 1] for i in ans[0]], "b.", markersize=8, label="1")
plt.plot([P[i, 0] for i in ans[1]], [P[i, 1] for i in ans[1]], "y.", markersize=8, label="2")
plt.plot([P[i, 0] for i in ans[2]], [P[i, 1] for i in ans[2]], "k.", markersize=8, label="3")
plt.legend()
plt.show()

import matplotlib.pyplot as plt

import numpy as np

P = np.array([[1, 9], [2, 7], [3, 3],

[2.5, 10], [4, 4.5], [5, 1],

[4.5, 7.5], [6, 4], [7, 2], [2, 4], [4, 2]])

def fast_non_dominated_sort(values1, values2):

S = [[] for i in range(0, len(values1))]

front = [[]]

n = [0 for i in range(0, len(values1))]

rank = [0 for i in range(0, len(values1))]

for p in range(0, len(values1)):

S[p] = []

n[p] = 0

for q in range(0, len(values1)):

if (values1[p] > values1[q] and values2[p] > values2[q]) or (values1[p] >= values1[q] and values2[p] > values2[q]) or (

values1[p] > values1[q] and values2[p] >= values2[q]):

if q not in S[p]:

S[p].append(q)

elif (values1[q] > values1[p] and values2[q] > values2[p]) or (values1[q] >= values1[p] and values2[q] > values2[p]) or (

values1[q] > values1[p] and values2[q] >= values2[p]):

n[p] = n[p] + 1

if n[p] == 0:

rank[p] = 0

if p not in front[0]:

front[0].append(p)

print(S, n)

i = 0

while front[i] != []:

Q = []

for p in front[i]:

for q in S[p]:

n[q] = n[q] - 1

if n[q] == 0:

rank[q] = i + 1

if q not in Q:

Q.append(q)

i = i + 1

front.append(Q)

del front[len(front) - 1]

return front

ans = fast_non_dominated_sort(P[:, 0], P[:, 1])

plt.plot([P[i, 0] for i in ans[0]], [P[i, 1] for i in ans[0]], "b.", markersize=8, label="1")

plt.plot([P[i, 0] for i in ans[1]], [P[i, 1] for i in ans[1]], "y.", markersize=8, label="2")

plt.plot([P[i, 0] for i in ans[2]], [P[i, 1] for i in ans[2]], "k.", markersize=8, label="3")

plt.legend()

plt.show()

拥挤度计算

个人认为很多时候可以用用拥挤距离替代，每个类都分别计算拥挤度，但最终很多时候在算法中只需要计算一个rank的拥挤度。

一个rank中的所有点，其中第一个点和最后一个点的拥挤距离为∞，中间点的距离即delta_y1 + delta_y2。拥挤距离可用delta_y1 / (max(y1) – min(y1)) + delta_y2 / (max(y2) – min(y2))

精英策略

每次从n个种群交叉变异后会生成2n的种群，通过精英策略选取其中前n个进行下一次迭代。策略如下：

rank大/小的优先，看具体题目
如果有个rank加上后在超过了n个个体的种群，则进行拥挤度排序，选取最大/小的直到种群为n个个体。

Github上大佬的代码

# Program Name: NSGA-II.py
# Description: This is a python implementation of Prof. Kalyanmoy Deb's popular NSGA-II algorithm
# Author: Haris Ali Khan 
# Supervisor: Prof. Manoj Kumar Tiwari

# Importing required modules
import math
import random
import matplotlib.pyplot as plt


# First function to optimize
def function1(x):
    # value = -x ** 2
    value = x * math.sin(10 * math.pi * x) + 2
    return value


# Second function to optimize
def function2(x):
    # value = -(x - 2) ** 2
    value = 1
    return value


# Function to find index of list
def index_of(a, list):
    for i in range(0, len(list)):
        if list[i] == a:
            return i
    return -1


# Function to sort by values
def sort_by_values(list1, values):
    sorted_list = []
    while len(sorted_list) != len(list1):
        if index_of(min(values), values) in list1:
            sorted_list.append(index_of(min(values), values))
        values[index_of(min(values), values)] = math.inf
    return sorted_list


# Function to carry out NSGA-II's fast non dominated sort
def fast_non_dominated_sort(values1, values2):
    S = [[] for i in range(0, len(values1))]
    front = [[]]
    n = [0 for i in range(0, len(values1))]
    rank = [0 for i in range(0, len(values1))]

    for p in range(0, len(values1)):
        S[p] = []
        n[p] = 0
        for q in range(0, len(values1)):
            if (values1[p] > values1[q] and values2[p] > values2[q]) or (values1[p] >= values1[q] and values2[p] > values2[q]) or (
                    values1[p] > values1[q] and values2[p] >= values2[q]):
                if q not in S[p]:
                    S[p].append(q)
            elif (values1[q] > values1[p] and values2[q] > values2[p]) or (values1[q] >= values1[p] and values2[q] > values2[p]) or (
                    values1[q] > values1[p] and values2[q] >= values2[p]):
                n[p] = n[p] + 1
        if n[p] == 0:
            rank[p] = 0
            if p not in front[0]:
                front[0].append(p)

    i = 0
    while front[i] != []:
        Q = []
        for p in front[i]:
            for q in S[p]:
                n[q] = n[q] - 1
                if n[q] == 0:
                    rank[q] = i + 1
                    if q not in Q:
                        Q.append(q)
        i = i + 1
        front.append(Q)

    del front[len(front) - 1]
    return front


# Function to calculate crowding distance
def crowding_distance(values1, values2, front):
    distance = [0 for i in range(0, len(front))]
    sorted1 = sort_by_values(front, values1[:])
    sorted2 = sort_by_values(front, values2[:])
    distance[0] = 0x6f6f6f6f
    distance[len(front) - 1] = 0x6f6f6f6f
    for k in range(1, len(front) - 1):
        distance[k] = distance[k] + (values1[sorted1[k + 1]] - values1[sorted1[k - 1]]) / (max(values1) - min(values1))
    for k in range(1, len(front) - 1):
        distance[k] = distance[k] + (values2[sorted2[k + 1]] - values2[sorted2[k - 1]]) / (max(values2) - min(values2))
    return distance


# Function to carry out the crossover
def crossover(a, b):
    r = random.random()
    if r > 0.5:
        return mutation((a + b) / 2)
    else:
        return mutation((a - b) / 2)


# Function to carry out the mutation operator
def mutation(solution):
    mutation_prob = random.random()
    if mutation_prob < 1:
        solution = min_x + (max_x - min_x) * random.random()
    return solution


# Main program starts here
pop_size = 20
max_gen = 100

# Initialization
min_x = -1
max_x = 2
solution = [min_x + (max_x - min_x) * random.random() for i in range(0, pop_size)]
gen_no = 0
while gen_no < max_gen:
    function1_values = [function1(solution[i]) for i in range(0, pop_size)]
    function2_values = [function2(solution[i]) for i in range(0, pop_size)]
    non_dominated_sorted_solution = fast_non_dominated_sort(function1_values[:], function2_values[:])
    print("The best front for Generation number ", gen_no, " is")
    for valuez in non_dominated_sorted_solution[0]:
        print(round(solution[valuez], 3), end=" ")
    print("\n")
    crowding_distance_values = []
    for i in range(0, len(non_dominated_sorted_solution)):
        crowding_distance_values.append(crowding_distance(function1_values[:], function2_values[:], non_dominated_sorted_solution[i][:]))
    solution2 = solution[:]
    # Generating offsprings
    while len(solution2) != 2 * pop_size:
        a1 = random.randint(0, pop_size - 1)
        b1 = random.randint(0, pop_size - 1)
        solution2.append(crossover(solution[a1], solution[b1]))
    function1_values2 = [function1(solution2[i]) for i in range(0, 2 * pop_size)]
    function2_values2 = [function2(solution2[i]) for i in range(0, 2 * pop_size)]
    non_dominated_sorted_solution2 = fast_non_dominated_sort(function1_values2[:], function2_values2[:])
    crowding_distance_values2 = []
    for i in range(0, len(non_dominated_sorted_solution2)):
        crowding_distance_values2.append(crowding_distance(function1_values2[:], function2_values2[:], non_dominated_sorted_solution2[i][:]))
    new_solution = []
    for i in range(0, len(non_dominated_sorted_solution2)):
        non_dominated_sorted_solution2_1 = [index_of(non_dominated_sorted_solution2[i][j], non_dominated_sorted_solution2[i]) for j in
                                            range(0, len(non_dominated_sorted_solution2[i]))]
        front22 = sort_by_values(non_dominated_sorted_solution2_1[:], crowding_distance_values2[i][:])
        front = [non_dominated_sorted_solution2[i][front22[j]] for j in range(0, len(non_dominated_sorted_solution2[i]))]
        front.reverse()
        for value in front:
            new_solution.append(value)
            if len(new_solution) == pop_size:
                break
        if len(new_solution) == pop_size:
            break
    solution = [solution2[i] for i in new_solution]
    gen_no = gen_no + 1

# Lets plot the final front now
print(solution[0])
print(function1(solution[0]))
function1 = [i * -1 for i in function1_values]
function2 = [j * -1 for j in function2_values]
plt.xlabel('Function 1', fontsize=15)
plt.ylabel('Function 2', fontsize=15)
plt.scatter(function1, function2)
plt.show()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

# Program Name: NSGA-II.py

# Description: This is a python implementation of Prof. Kalyanmoy Deb's popular NSGA-II algorithm

# Author: Haris Ali Khan

# Supervisor: Prof. Manoj Kumar Tiwari

# Importing required modules

import math

import random

import matplotlib.pyplot as plt

# First function to optimize

def function1(x):

# value = -x ** 2

value = x * math.sin(10 * math.pi * x) + 2

return value

# Second function to optimize

def function2(x):

# value = -(x - 2) ** 2

value = 1

return value

# Function to find index of list

def index_of(a, list):

for i in range(0, len(list)):

if list[i] == a:

return i

return -1

# Function to sort by values

def sort_by_values(list1, values):

sorted_list = []

while len(sorted_list) != len(list1):

if index_of(min(values), values) in list1:

sorted_list.append(index_of(min(values), values))

values[index_of(min(values), values)] = math.inf

return sorted_list

# Function to carry out NSGA-II's fast non dominated sort

def fast_non_dominated_sort(values1, values2):

S = [[] for i in range(0, len(values1))]

front = [[]]

n = [0 for i in range(0, len(values1))]

rank = [0 for i in range(0, len(values1))]

for p in range(0, len(values1)):

S[p] = []

n[p] = 0

for q in range(0, len(values1)):

if (values1[p] > values1[q] and values2[p] > values2[q]) or (values1[p] >= values1[q] and values2[p] > values2[q]) or (

values1[p] > values1[q] and values2[p] >= values2[q]):

if q not in S[p]:

S[p].append(q)

elif (values1[q] > values1[p] and values2[q] > values2[p]) or (values1[q] >= values1[p] and values2[q] > values2[p]) or (

values1[q] > values1[p] and values2[q] >= values2[p]):

n[p] = n[p] + 1

if n[p] == 0:

rank[p] = 0

if p not in front[0]:

front[0].append(p)

i = 0

while front[i] != []:

Q = []

for p in front[i]:

for q in S[p]:

n[q] = n[q] - 1

if n[q] == 0:

rank[q] = i + 1

if q not in Q:

Q.append(q)

i = i + 1

front.append(Q)

del front[len(front) - 1]

return front

# Function to calculate crowding distance

def crowding_distance(values1, values2, front):

distance = [0 for i in range(0, len(front))]

sorted1 = sort_by_values(front, values1[:])

sorted2 = sort_by_values(front, values2[:])

distance[0] = 0x6f6f6f6f

distance[len(front) - 1] = 0x6f6f6f6f

for k in range(1, len(front) - 1):

distance[k] = distance[k] + (values1[sorted1[k + 1]] - values1[sorted1[k - 1]]) / (max(values1) - min(values1))

for k in range(1, len(front) - 1):

distance[k] = distance[k] + (values2[sorted2[k + 1]] - values2[sorted2[k - 1]]) / (max(values2) - min(values2))

return distance

# Function to carry out the crossover

def crossover(a, b):

r = random.random()

if r > 0.5:

return mutation((a + b) / 2)

else:

return mutation((a - b) / 2)

# Function to carry out the mutation operator

def mutation(solution):

mutation_prob = random.random()

if mutation_prob < 1:

solution = min_x + (max_x - min_x) * random.random()

return solution

# Main program starts here

pop_size = 20

max_gen = 100

# Initialization

min_x = -1

max_x = 2

solution = [min_x + (max_x - min_x) * random.random() for i in range(0, pop_size)]

gen_no = 0

while gen_no < max_gen:

function1_values = [function1(solution[i]) for i in range(0, pop_size)]

function2_values = [function2(solution[i]) for i in range(0, pop_size)]

non_dominated_sorted_solution = fast_non_dominated_sort(function1_values[:], function2_values[:])

print("The best front for Generation number ", gen_no, " is")

for valuez in non_dominated_sorted_solution[0]:

print(round(solution[valuez], 3), end=" ")

print("\n")

crowding_distance_values = []

for i in range(0, len(non_dominated_sorted_solution)):

crowding_distance_values.append(crowding_distance(function1_values[:], function2_values[:], non_dominated_sorted_solution[i][:]))

solution2 = solution[:]

# Generating offsprings

while len(solution2) != 2 * pop_size:

a1 = random.randint(0, pop_size - 1)

b1 = random.randint(0, pop_size - 1)

solution2.append(crossover(solution[a1], solution[b1]))

function1_values2 = [function1(solution2[i]) for i in range(0, 2 * pop_size)]

function2_values2 = [function2(solution2[i]) for i in range(0, 2 * pop_size)]

non_dominated_sorted_solution2 = fast_non_dominated_sort(function1_values2[:], function2_values2[:])

crowding_distance_values2 = []

for i in range(0, len(non_dominated_sorted_solution2)):

crowding_distance_values2.append(crowding_distance(function1_values2[:], function2_values2[:], non_dominated_sorted_solution2[i][:]))

new_solution = []

for i in range(0, len(non_dominated_sorted_solution2)):

non_dominated_sorted_solution2_1 = [index_of(non_dominated_sorted_solution2[i][j], non_dominated_sorted_solution2[i]) for j in

range(0, len(non_dominated_sorted_solution2[i]))]

front22 = sort_by_values(non_dominated_sorted_solution2_1[:], crowding_distance_values2[i][:])

front = [non_dominated_sorted_solution2[i][front22[j]] for j in range(0, len(non_dominated_sorted_solution2[i]))]

front.reverse()

for value in front:

new_solution.append(value)

if len(new_solution) == pop_size:

break

if len(new_solution) == pop_size:

break

solution = [solution2[i] for i in new_solution]

gen_no = gen_no + 1

# Lets plot the final front now

print(solution[0])

print(function1(solution[0]))

function1 = [i * -1 for i in function1_values]

function2 = [j * -1 for j in function2_values]

plt.xlabel('Function 1', fontsize=15)

plt.ylabel('Function 2', fontsize=15)

plt.scatter(function1, function2)

plt.show()

优化算法的学习

个人关于优化算法的学习历程：SA GA TSP PSO NSGA-ii

下列代码为个人学习及复现的过程，不一定很正确，有一定参考性，NSGA-ii不在此展示，详情可见下一篇。

退火（SA）

退火简述：个人认为有点像穷举，在x附近左右横跳，达到条件降温，没达到升温

from math import *
from random import random

# 方程或问题的计算
func = lambda x, y: cos(x * y) + x * y + y ** 3 + x ** 2
x_min = -100
y_min = -100
x_max = 100
y_max = 100

# 设置 SA 参数
e = 1e-50
L = 2000000
at = 0.999
T = 1

# 选取 x(0)
x = random() * (x_max - x_min) + x_min
y = random() * (y_max - y_min) + y_min
result = func(x, y)

# 循环
for _ in range(L):
    # 计算 x‘ ，需要选取合适的计算方法，以下为两种计算方法，各有利弊
    # x_ = random() * (x_max - x_min) + x_min
    # y_ = random() * (y_max - y_min) + y_min
    step = 0.1
    x_ = x + (random() * (x_max - x_min) + x_min) * step
    x_ = x_max if x_ > x_max else x_min if x_ < x_min else x_
    y_ = y + (random() * (y_max - y_min) + y_min) * step
    y_ = y_max if y_ > y_max else y_min if y_ < y_min else y_
    result_ = func(x_, y_)

    # 判断是否满足退火条件
    if result_ - result < 0 or exp((result - result_) / T) > random():
        # 更新
        result = result_
        x = x_
        y = y_

    # 温度降低 降低到一定值推出
    T *= at
    if T < e:
        break

print(result, x, y)

from math import *

from random import random

# 方程或问题的计算

func = lambda x, y: cos(x * y) + x * y + y ** 3 + x ** 2

x_min = -100

y_min = -100

x_max = 100

y_max = 100

# 设置 SA 参数

e = 1e-50

L = 2000000

at = 0.999

T = 1

# 选取 x(0)

x = random() * (x_max - x_min) + x_min

y = random() * (y_max - y_min) + y_min

result = func(x, y)

# 循环

for _ in range(L):

# 计算 x‘ ，需要选取合适的计算方法，以下为两种计算方法，各有利弊

# x_ = random() * (x_max - x_min) + x_min

# y_ = random() * (y_max - y_min) + y_min

step = 0.1

x_ = x + (random() * (x_max - x_min) + x_min) * step

x_ = x_max if x_ > x_max else x_min if x_ < x_min else x_

y_ = y + (random() * (y_max - y_min) + y_min) * step

y_ = y_max if y_ > y_max else y_min if y_ < y_min else y_

result_ = func(x_, y_)

# 判断是否满足退火条件

if result_ - result < 0 or exp((result - result_) / T) > random():

# 更新

result = result_

x = x_

y = y_

# 温度降低降低到一定值推出

T *= at

if T < e:

break

print(result, x, y)

遗传（GA）

遗传简述：对问题编码，可以是二进制也可以是十进制或者是排列（各位都不一样），交叉，变异进化，其中会包含精英保留等操作，NSGA在此基础上加上了非支配排序。

from math import sin, pi
from random import randint, random
from operator import itemgetter, attrgetter

# 方程或问题的计算
func = lambda x: x * sin(10 * pi * x) + 2
x_min = -1
x_max = 2

# 设置遗传算法各种参数，分别代表：小数点后精度(epsilon)，编码位数(coding digit)，交叉变换位置(t)，交叉概率(pc)，变异概率(pm)，种群大小(M)，循环次数(L)
epsilon = 6
coding_digit = 22
t = 4
pc = 0.8
pm = 0.2
M = 1000
L = 100


# 初始化种群
class Population:
    x = None
    xFlag = 1
    result = None

    def __init__(self, x):
        self.x = abs(x_max if x > x_max else x_min if x < x_min else x)
        self.xFlag = 1 if x > 0 else -1
        self.result = func(self.x * self.xFlag)


population = []
for i in range(M):
    x = (x_max - x_min) / M * i + x_min
    population.append(Population(x))


def binary(X: float) -> list:
    X = list(bin(int(X * 10 ** epsilon)))[2:]
    return ["0" for i in range(coding_digit - len(X))] + X


def decimal(X: list) -> float:
    return int("".join(X), 2) / 10 ** epsilon


def cross(X: list, Y: list) -> [list, list]:
    X_ = X[:t] + Y[t:]
    Y_ = Y[:t] + X[t:]
    return X_, Y_


def variation(X: list) -> list:
    # 可以固定变异的位置也可以随机，以下设置为随机
    position = randint(0, coding_digit - 1)
    X[position] = "1" if X[position] == "0" else "0"
    return X


for _ in range(L):
    # 对种群进行排序，并选出最小的种群
    population = sorted(population, key=lambda X: X.result, reverse=True)
    population = population[:M]
    # 交叉
    for i in range(0, M, 2):
        if random() < pc:
            x1, x2 = cross(binary(population[i].x), binary(population[i + 1].x))
            x1 = decimal(x1) * population[i].xFlag
            x2 = decimal(x2) * population[i + 1].xFlag
            population.append(Population(x1))
            population.append(Population(x2))
    # 变异
    length = len(population)
    for i in range(length):
        if random() < pm:
            x = decimal(variation(binary(population[i].x))) * population[i].xFlag
            population.append(Population(x))

sorted(population, key=lambda X: X.result)
ans = population[0]
print(ans.x * ans.xFlag, ans.result)

from math import sin, pi

from random import randint, random

from operator import itemgetter, attrgetter

# 方程或问题的计算

func = lambda x: x * sin(10 * pi * x) + 2

x_min = -1

x_max = 2

# 设置遗传算法各种参数，分别代表：小数点后精度(epsilon)，编码位数(coding digit)，交叉变换位置(t)，交叉概率(pc)，变异概率(pm)，种群大小(M)，循环次数(L)

epsilon = 6

coding_digit = 22

t = 4

pc = 0.8

pm = 0.2

M = 1000

L = 100

# 初始化种群

class Population:

x = None

xFlag = 1

result = None

def __init__(self, x):

self.x = abs(x_max if x > x_max else x_min if x < x_min else x)

self.xFlag = 1 if x > 0 else -1

self.result = func(self.x * self.xFlag)

population = []

for i in range(M):

x = (x_max - x_min) / M * i + x_min

population.append(Population(x))

def binary(X: float) -> list:

X = list(bin(int(X * 10 ** epsilon)))[2:]

return ["0" for i in range(coding_digit - len(X))] + X

def decimal(X: list) -> float:

return int("".join(X), 2) / 10 ** epsilon

def cross(X: list, Y: list) -> [list, list]:

X_ = X[:t] + Y[t:]

Y_ = Y[:t] + X[t:]

return X_, Y_

def variation(X: list) -> list:

# 可以固定变异的位置也可以随机，以下设置为随机

position = randint(0, coding_digit - 1)

X[position] = "1" if X[position] == "0" else "0"

return X

for _ in range(L):

# 对种群进行排序，并选出最小的种群

population = sorted(population, key=lambda X: X.result, reverse=True)

population = population[:M]

# 交叉

for i in range(0, M, 2):

if random() < pc:

x1, x2 = cross(binary(population[i].x), binary(population[i + 1].x))

x1 = decimal(x1) * population[i].xFlag

x2 = decimal(x2) * population[i + 1].xFlag

population.append(Population(x1))

population.append(Population(x2))

# 变异

length = len(population)

for i in range(length):

if random() < pm:

x = decimal(variation(binary(population[i].x))) * population[i].xFlag

population.append(Population(x))

sorted(population, key=lambda X: X.result)

ans = population[0]

print(ans.x * ans.xFlag, ans.result)

蚁群（TSP）

蚁群简述：信息素的更新与获取，信息素可以不用完全按照公式，根据具体题目可以有所不同。

from math import sin, pi
from random import random, choice
import numpy as np
import matplotlib.pyplot as plt

func = lambda x: x * sin(10 * pi * x) + 2
x_min = -1
x_max = 2

# alpha = 2.5
# beta = 1
rou = 0.2
M = 100
L = 700
step = 1
Q = 1

population = []
t = []
possibility = []
for i in range(M):
    population.append(random() * (x_max - x_min) + x_min)
    t.append(func(population[i]))
    possibility.append(1 / t[i])

t_max = []
for loop in range(1, L + 1):
    step = 1 / loop
    tBest = max(t)
    for i in range(M):
        possibility[i] = (tBest - t[i]) / tBest
        if possibility[i] < random():
            population[i] += step * choice([-1, 1])
        else:
            population[i] += choice([-1, 1]) * random() * (x_max - x_min) / 2
        population[i] = x_min if population[i] < x_min else x_max if population[i] > x_max else population[i]
        t[i] = rou * t[i] + Q * func(population[i])
    t_max.append(func(population[t.index(max(t))]))

I = t.index(max(t))
print(population[I], func(population[I]))
plt.plot(np.linspace(1, L, L), t_max)
plt.show()

from math import sin, pi

from random import random, choice

import numpy as np

import matplotlib.pyplot as plt

func = lambda x: x * sin(10 * pi * x) + 2

x_min = -1

x_max = 2

# alpha = 2.5

# beta = 1

rou = 0.2

M = 100

L = 700

step = 1

Q = 1

population = []

t = []

possibility = []

for i in range(M):

population.append(random() * (x_max - x_min) + x_min)

t.append(func(population[i]))

possibility.append(1 / t[i])

t_max = []

for loop in range(1, L + 1):

step = 1 / loop

tBest = max(t)

for i in range(M):

possibility[i] = (tBest - t[i]) / tBest

if possibility[i] < random():

population[i] += step * choice([-1, 1])

else:

population[i] += choice([-1, 1]) * random() * (x_max - x_min) / 2

population[i] = x_min if population[i] < x_min else x_max if population[i] > x_max else population[i]

t[i] = rou * t[i] + Q * func(population[i])

t_max.append(func(population[t.index(max(t))]))

I = t.index(max(t))

print(population[I], func(population[I]))

plt.plot(np.linspace(1, L, L), t_max)

plt.show()

粒子群（PSO）

粒子群简述：速度的更新，包含个人学习速率和社会学习速率，不同点的速度不同，更新的速度不同。

from math import cos
from random import random
import numpy as np

# 方程或问题的计算
func = lambda x, y: -(cos(x * y) + x * y + y ** 3 + x ** 2)
x_min = -100
y_min = -100
x_max = 100
y_max = 100

# 设置参数
L = 100
omega = 1
groupLR = 0.5
socialLR = 1.5
M = 20
v_min = -1
v_max = 1

# 初始化种群
population = np.zeros(M * 2).reshape(M, 2)
v = np.zeros(M * 2).reshape(M, 2)
fitness = np.zeros(M)
for i in range(M):
    population[i][0] = random() * (x_max - x_min) + x_min
    population[i][1] = random() * (x_max - x_min) + x_min
    v[i][0] = random() * (x_max - x_min) + x_min
    v[i][1] = random() * (x_max - x_min) + x_min
    fitness[i] = func(population[i][0], population[i][1])

# 初始化个人与种群最优解
# M * 2
personalBestFitness = fitness.copy()
personalPopulation = population.copy()
# 2 * 2
groupBest = population[fitness.argmax()].copy()
groupBestFitness = fitness.max()

for _ in range(L):
    for i in range(len(fitness)):
        # 更新速度
        # 随时间适应度 + 个人学习速率 * （个人最好 - 当前） + 社会学习速率 * （团队最好 - 当前）
        v[i] = omega * v[i] + groupLR * random() * (personalPopulation[i] - population[i]) + socialLR * random() * (groupBest - population[i])
    v[v < v_min] = v_min
    v[v > v_max] = v_max

    for i in range(len(fitness)):
        # 更新位置
        population[i] += v[i]
        population[i][0] = x_min if population[i][0] < x_min else x_max if population[i][0] > x_max else population[i][0]
        population[i][1] = y_min if population[i][1] < y_min else y_max if population[i][1] > y_max else population[i][1]

        # 适应度更新
        fitness[i] = func(population[i][0], population[i][1])

        # 个人最好更新
        if fitness[i] > personalBestFitness[i]:
            personalBestFitness[i] = fitness[i]
            personalPopulation[i] = population[i].copy()

    # 团队最好更新
    if personalBestFitness.max() > groupBestFitness:
        groupBestFitness = personalBestFitness.max()
        groupBest = population[personalBestFitness.argmax()].copy()

print(groupBestFitness, groupBest[0], groupBest[1])

from math import cos

from random import random

import numpy as np

# 方程或问题的计算

func = lambda x, y: -(cos(x * y) + x * y + y ** 3 + x ** 2)

x_min = -100

y_min = -100

x_max = 100

y_max = 100

# 设置参数

L = 100

omega = 1

groupLR = 0.5

socialLR = 1.5

M = 20

v_min = -1

v_max = 1

# 初始化种群

population = np.zeros(M * 2).reshape(M, 2)

v = np.zeros(M * 2).reshape(M, 2)

fitness = np.zeros(M)

for i in range(M):

population[i][0] = random() * (x_max - x_min) + x_min

population[i][1] = random() * (x_max - x_min) + x_min

v[i][0] = random() * (x_max - x_min) + x_min

v[i][1] = random() * (x_max - x_min) + x_min

fitness[i] = func(population[i][0], population[i][1])

# 初始化个人与种群最优解

# M * 2

personalBestFitness = fitness.copy()

personalPopulation = population.copy()

# 2 * 2

groupBest = population[fitness.argmax()].copy()

groupBestFitness = fitness.max()

for _ in range(L):

for i in range(len(fitness)):

# 更新速度

# 随时间适应度 + 个人学习速率 * （个人最好 - 当前） + 社会学习速率 * （团队最好 - 当前）

v[i] = omega * v[i] + groupLR * random() * (personalPopulation[i] - population[i]) + socialLR * random() * (groupBest - population[i])

v[v < v_min] = v_min

v[v > v_max] = v_max

for i in range(len(fitness)):

# 更新位置

population[i] += v[i]

population[i][0] = x_min if population[i][0] < x_min else x_max if population[i][0] > x_max else population[i][0]

population[i][1] = y_min if population[i][1] < y_min else y_max if population[i][1] > y_max else population[i][1]

# 适应度更新

fitness[i] = func(population[i][0], population[i][1])

# 个人最好更新

if fitness[i] > personalBestFitness[i]:

personalBestFitness[i] = fitness[i]

personalPopulation[i] = population[i].copy()

# 团队最好更新

if personalBestFitness.max() > groupBestFitness:

groupBestFitness = personalBestFitness.max()

groupBest = population[personalBestFitness.argmax()].copy()

print(groupBestFitness, groupBest[0], groupBest[1])

Pytorch 神经网络基础

请先看该文章。

当要运用到自己的图片文件的时候可在上文中加上以下代码，并相应进行改变神经网络中data运算。参考可见。

from torch.utils.data import Dataset, DataLoader
class picture(Dataset):
    def __init__(self):
        self.x = data
        self.y = label
        self.len = len(self.y)

    def __len__(self):
        return self.len

    def __getitem__(self, item):
        return self.x[item], self.y[item]

trainset = word()
trainsetloader = torch.utils.data.DataLoader(trainset, batch_size=30, shuffle=False)

from torch.utils.data import Dataset, DataLoader

class picture(Dataset):

def __init__(self):

self.x = data

self.y = label

self.len = len(self.y)

def __len__(self):

return self.len

def __getitem__(self, item):

return self.x[item], self.y[item]

trainset = word()

trainsetloader = torch.utils.data.DataLoader(trainset, batch_size=30, shuffle=False)

准确率计算参考代码：

import torch
from torch import nn
from torchvision import datasets, transforms

model = torch.load("my_handwrite_recognize_model.pt")
transform = transforms.Compose(
    [transforms.ToTensor(),
     transforms.Normalize((0.5,), (0.5,))])
testset = datasets.MNIST(root='./data', train=True,
                         download=True, transform=transform)
testdataiter = iter(testset)

sum = 0
correct = 0
for testimages, testlabels in testdataiter:
    sum += 1
    img_vector = testimages[0].squeeze().view(1, -1)
    # 模型返回的是一个1x10的矩阵，第几个元素值最大那就是表示模型认为当前图片是数字几
    result_digit = model(img_vector)
    if int(result_digit.max(1)[1]) == testlabels:
        correct += 1

print(correct / sum)

import torch

from torch import nn

from torchvision import datasets, transforms

model = torch.load("my_handwrite_recognize_model.pt")

transform = transforms.Compose(

[transforms.ToTensor(),

transforms.Normalize((0.5,), (0.5,))])

testset = datasets.MNIST(root='./data', train=True,

download=True, transform=transform)

testdataiter = iter(testset)

sum = 0

correct = 0

for testimages, testlabels in testdataiter:

sum += 1

img_vector = testimages[0].squeeze().view(1, -1)

# 模型返回的是一个1x10的矩阵，第几个元素值最大那就是表示模型认为当前图片是数字几

result_digit = model(img_vector)

if int(result_digit.max(1)[1]) == testlabels:

correct += 1

print(correct / sum)

下面用分类损失与回归损失代码对比。可供测试data文件。测试文件说明：最后一列为label，除最后一列外为data

import torch
import numpy as np

# 打开文件
with open('data.csv') as f:
    data = np.loadtxt(f, delimiter=",", skiprows=0)
    # 注意此时 label 是一维，不能是二维
    label = data[:, -1]
    data = data[:, :-1]

# numpy2torch
data = torch.from_numpy(data).float()
label = torch.from_numpy(label).long()

# 网络结构
hidden = [32, 32, 32, 32]
dim_in, dim_out = 7, 28
model = torch.nn.Sequential(
    torch.nn.Linear(dim_in, hidden[0]),
    torch.nn.ReLU(),
    torch.nn.Linear(hidden[0], hidden[1]),
    torch.nn.ReLU(),
    torch.nn.Linear(hidden[1], hidden[2]),
    torch.nn.ReLU(),
    torch.nn.Linear(hidden[2], dim_out),
)

# 损失函数
loss_fn = torch.nn.CrossEntropyLoss()

# 优化器
learning_rate = 1e-4
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

# 训练神经网络
for t in range(100000):
    # 让神经网络根据现有的参数，根据当前的输入计算一个输出
    model_output = model(data)
    # 用所设计算损失(误差)函数计算误差
    loss = loss_fn(model_output, label)
    print(loss)
    # 每次训练前清零之前计算的梯度(导数)
    optimizer.zero_grad()
    # 根据误差反向传播计算误差对各个权重的导数
    loss.backward()
    # 根据优化器里面的算法自动调整神经网络权重
    optimizer.step()
total = len(data)
corr = 0
for i in range(len(data)):
    label_ = model(data[i]).argmax()
    if label_ == label[i]:
        corr += 1
print(corr / total)
# 保存下训练好的模型,省得下次再重新训练
torch.save(model, './word.pt')

import torch

import numpy as np

# 打开文件

with open('data.csv') as f:

data = np.loadtxt(f, delimiter=",", skiprows=0)

# 注意此时 label 是一维，不能是二维

label = data[:, -1]

data = data[:, :-1]

# numpy2torch

data = torch.from_numpy(data).float()

label = torch.from_numpy(label).long()

# 网络结构

hidden = [32, 32, 32, 32]

dim_in, dim_out = 7, 28

model = torch.nn.Sequential(

torch.nn.Linear(dim_in, hidden[0]),

torch.nn.ReLU(),

torch.nn.Linear(hidden[0], hidden[1]),

torch.nn.ReLU(),

torch.nn.Linear(hidden[1], hidden[2]),

torch.nn.ReLU(),

torch.nn.Linear(hidden[2], dim_out),

)

# 损失函数

loss_fn = torch.nn.CrossEntropyLoss()

# 优化器

learning_rate = 1e-4

optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

# 训练神经网络

for t in range(100000):

# 让神经网络根据现有的参数，根据当前的输入计算一个输出

model_output = model(data)

# 用所设计算损失(误差)函数计算误差

loss = loss_fn(model_output, label)

print(loss)

# 每次训练前清零之前计算的梯度(导数)

optimizer.zero_grad()

# 根据误差反向传播计算误差对各个权重的导数

loss.backward()

# 根据优化器里面的算法自动调整神经网络权重

optimizer.step()

total = len(data)

corr = 0

for i in range(len(data)):

label_ = model(data[i]).argmax()

if label_ == label[i]:

corr += 1

print(corr / total)

# 保存下训练好的模型,省得下次再重新训练

torch.save(model, './word.pt')

以上运用分类损失，效果一般，只有46%的分类准确率。

import torch
import numpy as np

# 打开文件
with open('data.csv') as f:
    data = np.loadtxt(f, delimiter=",", skiprows=0)
    # 注意此时 label 是二维，不能是一维
    label = data[:, -1:]
    data = data[:, :-1]

# numpy2torch
data = torch.from_numpy(data).float()
label = torch.from_numpy(label).float()

# 网络结构
hidden = [32, 32, 32, 32]
dim_in, dim_out = 7, 1
model = torch.nn.Sequential(
    torch.nn.Linear(dim_in, hidden[0]),
    torch.nn.ReLU(),
    torch.nn.Linear(hidden[0], hidden[1]),
    torch.nn.ReLU(),
    torch.nn.Linear(hidden[1], hidden[2]),
    torch.nn.ReLU(),
    torch.nn.Linear(hidden[2], dim_out),
)

# 损失函数
loss_fn = torch.nn.MSELoss()


def r2_loss(y_test, y):
    return 1 - ((y_test - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()


# 优化器
learning_rate = 1e-4
optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

# 训练神经网络
for t in range(100000):
    # 让神经网络根据现有的参数，根据当前的输入计算一个输出
    model_output = model(data)
    # 用所设计算损失(误差)函数计算误差
    loss = loss_fn(model_output, label)
    print(loss)
    # 每次训练前清零之前计算的梯度(导数)
    optimizer.zero_grad()
    # 根据误差反向传播计算误差对各个权重的导数
    loss.backward()
    # 根据优化器里面的算法自动调整神经网络权重
    optimizer.step()

print(r2_loss(model(data), label))
# 保存下训练好的模型,省得下次再重新训练
torch.save(model, './word2.pt')

import torch

import numpy as np

# 打开文件

with open('data.csv') as f:

data = np.loadtxt(f, delimiter=",", skiprows=0)

# 注意此时 label 是二维，不能是一维

label = data[:, -1:]

data = data[:, :-1]

# numpy2torch

data = torch.from_numpy(data).float()

label = torch.from_numpy(label).float()

# 网络结构

hidden = [32, 32, 32, 32]

dim_in, dim_out = 7, 1

model = torch.nn.Sequential(

torch.nn.Linear(dim_in, hidden[0]),

torch.nn.ReLU(),

torch.nn.Linear(hidden[0], hidden[1]),

torch.nn.ReLU(),

torch.nn.Linear(hidden[1], hidden[2]),

torch.nn.ReLU(),

torch.nn.Linear(hidden[2], dim_out),

)

# 损失函数

loss_fn = torch.nn.MSELoss()

def r2_loss(y_test, y):

return 1 - ((y_test - y) ** 2).sum() / ((y.mean() - y) ** 2).sum()

# 优化器

learning_rate = 1e-4

optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)

# 训练神经网络

for t in range(100000):

# 让神经网络根据现有的参数，根据当前的输入计算一个输出

model_output = model(data)

# 用所设计算损失(误差)函数计算误差

loss = loss_fn(model_output, label)

print(loss)

# 每次训练前清零之前计算的梯度(导数)

optimizer.zero_grad()

# 根据误差反向传播计算误差对各个权重的导数

loss.backward()

# 根据优化器里面的算法自动调整神经网络权重

optimizer.step()

print(r2_loss(model(data), label))

# 保存下训练好的模型,省得下次再重新训练

torch.save(model, './word2.pt')

以上运用回归损失，效果显著，R²=0.96（评价函数）

二重积分参考代码

class doubleIntegral:
    def __init__(self, xMin, xMax, yMin, yMax, func):
        self.xMin = xMin
        self.xMax = xMax
        self.yMin = yMin
        self.yMax = yMax
        self.func = func

    def trapezoid(self, nX, nY):
        result = 0
        betweenX = (self.xMax - self.xMin) / nX
        betweenY = (self.yMax - self.yMin) / nY
        sums = []
        for i in range(nY + 1):
            temp = 0
            for j in range(nX):
                temp += betweenX / 2 * (self.func(self.xMin + j * betweenX, self.yMin + i * betweenY) +
                                        self.func(self.xMin + (j + 1) * betweenX, self.yMin + i * betweenY))
            sums.append(temp)
        ans = 0
        for i in range(nY):
            ans += betweenY / 2 * (sums[i] + sums[i + 1])
        return ans

    def simpsons(self, nX, nY):
        result = 0
        betweenX = (self.xMax - self.xMin) / nX
        betweenY = (self.yMax - self.yMin) / nY / 2
        sums = []
        for i in range(2 * nY + 1):
            temp = 0
            for j in range(nX):
                a = self.xMin + j * betweenX
                b = self.xMin + (j + 1) * betweenX
                c = (a + b) / 2
                temp += betweenX / 6 * (self.func(a, self.yMin + i * betweenY) + self.func(
                    b, self.yMin + i * betweenY) + 4 * self.func(c, self.yMin + i * betweenY))
            sums.append(temp)
        ans = 0
        for i in range(nY):
            ans += betweenY / 3 * \
                (sums[i * 2] + sums[(i + 1) * 2] + 4 * sums[i * 2 + 1])
        return ans


if __name__ == "__main__":
    # 个人不知道为什么梯形法求出来不对，有兴趣的可以自己去看看
    # 输入格式：xMin, xMax, yMin, yMax, func
    # 输入格式'： nx, ny
    print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).trapezoid(4, 8))
    print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).simpsons(4, 8))

class doubleIntegral:

def __init__(self, xMin, xMax, yMin, yMax, func):

self.xMin = xMin

self.xMax = xMax

self.yMin = yMin

self.yMax = yMax

self.func = func

def trapezoid(self, nX, nY):

result = 0

betweenX = (self.xMax - self.xMin) / nX

betweenY = (self.yMax - self.yMin) / nY

sums = []

for i in range(nY + 1):

temp = 0

for j in range(nX):

temp += betweenX / 2 * (self.func(self.xMin + j * betweenX, self.yMin + i * betweenY) +

self.func(self.xMin + (j + 1) * betweenX, self.yMin + i * betweenY))

sums.append(temp)

ans = 0

for i in range(nY):

ans += betweenY / 2 * (sums[i] + sums[i + 1])

return ans

def simpsons(self, nX, nY):

result = 0

betweenX = (self.xMax - self.xMin) / nX

betweenY = (self.yMax - self.yMin) / nY / 2

sums = []

for i in range(2 * nY + 1):

temp = 0

for j in range(nX):

a = self.xMin + j * betweenX

b = self.xMin + (j + 1) * betweenX

c = (a + b) / 2

temp += betweenX / 6 * (self.func(a, self.yMin + i * betweenY) + self.func(

b, self.yMin + i * betweenY) + 4 * self.func(c, self.yMin + i * betweenY))

sums.append(temp)

ans = 0

for i in range(nY):

ans += betweenY / 3 * \

(sums[i * 2] + sums[(i + 1) * 2] + 4 * sums[i * 2 + 1])

return ans

if __name__ == "__main__":

# 个人不知道为什么梯形法求出来不对，有兴趣的可以自己去看看

# 输入格式：xMin, xMax, yMin, yMax, func

# 输入格式'： nx, ny

print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).trapezoid(4, 8))

print(doubleIntegral(0, 1, 0, 2, lambda x, y: x**2+2*y).simpsons(4, 8))

树莓派4b jupyter frp端口映射

树莓派 frpc 设置：
从 github 中找到 arm 版本 wget 下载
frpc.ini 参考配置

[common]
server_addr = you server ip
server_port = 7000
token = your password

[jupyter]
type = http
local_ip = 127.0.0.1
local_port = 8888
custom_domains = your domain

[common]

server_addr = you server ip

server_port = 7000

token = your password

[jupyter]

type = http

local_ip = 127.0.0.1

local_port = 8888

custom_domains = your domain

jupyter 设置跟之前文章相同，即 ip, allow_remote_access等设置

服务器端 frps 配置: 云端开启7000~7003端口

[common]
bind_port = 7000
vhost_http_port = 7002
token = your password
dashboard_port = 7001
dashboard_user = Your dashboard usr name
dashboard_pwd = Your dashboard password
vhost_http_port = 7003

[common]

bind_port = 7000

vhost_http_port = 7002

token = your password

dashboard_port = 7001

dashboard_user = Your dashboard usr name

dashboard_pwd = Your dashboard password

vhost_http_port = 7003

云端 nginx 参考配置：

server{
    listen 80;
    listen [::]:80;
    server_name yourdomain; 
    location / {
        proxy_pass http://127.0.0.1:7003;
        proxy_set_header X-Real-IP $remote_addr;
        proxy_set_header Host $host; 
        proxy_set_header X-Forwarded-For $proxy_add_x_forwarded_for;
        proxy_http_version 1.1;
        proxy_set_header Upgrade $http_upgrade;
        proxy_set_header Connection "upgrade";
        proxy_read_timeout 120s;
        proxy_next_upstream error;
    }
}

server{

listen 80;

listen [::]:80;

server_name yourdomain;

location / {

proxy_pass http://127.0.0.1:7003;

proxy_set_header X-Real-IP $remote_addr;

proxy_set_header Host $host;

proxy_set_header X-Forwarded-For $proxy_add_x_forwarded_for;

proxy_http_version 1.1;

proxy_set_header Upgrade $http_upgrade;

proxy_set_header Connection "upgrade";

proxy_read_timeout 120s;

proxy_next_upstream error;

}

最终效果即通过域名 python.pancake2021.work 访问树莓派 jupyter