archieve: homework8

homework8/main.py

@@ -0,0 +1,179 @@
+"""
+连续型——Hopfield神经网络求解TSP
+1、初始化权值（A, D, U0）
+2、计算N个城市的距离矩阵dxy
+3、初始化神经网络的输入Uxi和输出Vxi
+4、利用动力微分方程计算：dUxi / dt
+5、由一阶欧拉方法更新计算：Uxi(t + 1) = Uxi(t) + dUxi / dt * step
+6、由非线性函数sigmoid更新计算：Vxi(t) = 0.5 * (1 + th(Uxi / U0))
+7、计算能量函数E
+8、检查路径是否合法
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+
+
+# 代价函数（具有三角不等式性质）
+def price_cn(vec1, vec2):
+    # 元素的平方和再开根号
+    return np.linalg.norm(np.array(vec1) - np.array(vec2))
+
+
+# 计算该方案下，总的路径长度
+def calc_distance(path):
+    dis = 0.0
+    for i in range(len(path) - 1):
+        dis += distance[path[i]][path[i + 1]]
+    return dis
+
+
+# 得到城市之间的距离矩阵
+def get_distance(citys):
+    N = len(citys)
+    # 构造一个N*N的零矩阵
+    distance = np.zeros((N, N))
+    for i, curr_point in enumerate(citys):
+        line = []
+        # 计算不同城市之间的距离
+        [line.append(price_cn(curr_point, other_point)) if i != j else line.append(0.0) for j, other_point in
+         enumerate(citys)]
+        # 把距离添加到矩阵相应的位置上
+        distance[i] = line
+    return distance
+
+
+# 动态方程计算微分方程du
+def calc_du(V, distance):
+    a = np.sum(V, axis=0) - 1  # 按列相加 - 1
+    b = np.sum(V, axis=1) - 1  # 按行相加 - 1
+    t1 = np.zeros((N, N))
+    t2 = np.zeros((N, N))
+    for i in range(N):
+        for j in range(N):
+            t1[i, j] = a[j]
+    for i in range(N):
+        for j in range(N):
+            t2[j, i] = b[j]
+    # 将第一列移动到最后一列
+    c_1 = V[:, 1:N]
+    # 构造一个N行1列的零举证
+    c_0 = np.zeros((N, 1))
+    c_0[:, 0] = V[:, 0]
+    # 把c_1和c_0在行方向上连接
+    c = np.concatenate((c_1, c_0), axis=1)
+    c = np.dot(distance, c)
+    return -A * (t1 + t2) - D * c
+
+
+# 更新神经网络的输入U
+# Uxi(t+1) = Uxi(t) + dUxi/dt * step
+def calc_U(U, du, step):
+    return U + du * step
+
+
+# 更新神经网络的输出V
+# Vxi(t) = 0.5 * (1 + th(Uxi/U0))
+def calc_V(U, U0):
+    return 1 / 2 * (1 + np.tanh(U / U0))
+
+
+# 计算当前网络的能量
+def calc_energy(V, distance):
+    t1 = np.sum(np.power(np.sum(V, axis=0) - 1, 2))
+    t2 = np.sum(np.power(np.sum(V, axis=1) - 1, 2))
+    idx = [i for i in range(1, N)]
+    idx = idx + [0]
+    Vt = V[:, idx]
+    t3 = distance * Vt
+    t3 = np.sum(np.sum(np.multiply(V, t3)))
+    e = 0.5 * (A * (t1 + t2) + D * t3)
+    return e
+
+
+# 检查路径的正确性
+def check_path(V):
+    newV = np.zeros([N, N])
+    route = []
+    for i in range(N):
+        mm = np.max(V[:, i])
+        for j in range(N):
+            if V[j, i] == mm:
+                newV[j, i] = 1
+                route += [j]
+                break
+    return route, newV
+
+
+# 可视化画出哈密顿回路和能量趋势
+def draw_H_and_E(citys, H_path, energys):
+    fig = plt.figure()
+    # 绘制哈密顿回路
+    ax1 = fig.add_subplot(121)
+    # 设置x轴的数值显示范围
+    plt.xlim(0, 7)
+    # 设置y轴的数值显示范围
+    plt.ylim(0, 7)
+    for (from_, to_) in H_path:
+        # 绘制城市点,大小为0.2，颜色为红色
+        p1 = plt.Circle(citys[from_], 0.2, color='red')
+        p2 = plt.Circle(citys[to_], 0.2, color='red')
+        ax1.add_patch(p1)
+        ax1.add_patch(p2)
+        ax1.plot((citys[from_][0], citys[to_][0]), (citys[from_][1], citys[to_][1]), color='red')
+        ax1.annotate(text=chr(97 + to_), xy=citys[to_], xytext=(-8, -4), textcoords='offset points', fontsize=20)
+    ax1.axis('equal')
+    ax1.grid()
+    # 绘制能量趋势图
+    ax2 = fig.add_subplot(122)
+    ax2.plot(np.arange(0, len(energys), 1), energys, color='red')
+    plt.show()
+
+
+# 定义城市坐标
+citys = np.array([[2, 6], [2, 4], [1, 3], [4, 6], [5, 5], [4, 4], [6, 4], [3, 2], [7, 7], [9, 3]])
+# 定义城市与城市之间的距离矩阵
+distance = get_distance(citys)
+# 计算城市个数
+N = len(citys)
+# 初始化参数A和D
+A = N * N
+D = N / 2
+U0 = 0.0009  # 初始输入
+step = 0.0001  # 步长
+num_iter = 10000  # 迭代次数
+# 初始化神经网络的输入状态（电路的输入U）
+U = 1 / 2 * U0 * np.log(N - 1) + (2 * (np.random.random((N, N))) - 1)
+# 初始化神经网络的输出状态（电路的输出V）
+V = calc_V(U, U0)
+energys = np.array([0.0 for x in range(num_iter)])  # 每次迭代的能量
+best_distance = np.inf  # 最优距离
+best_route = []  # 最优路线
+H_path = []  # 哈密顿回路
+# 开始迭代训练网络
+for n in range(num_iter):
+    # 利用动态方程计算du
+    du = calc_du(V, distance)
+    # 由一阶欧拉法更新下一个时间的输入状态（电路的输入U）
+    U = calc_U(U, du, step)
+    # 由sigmoid函数更新下一个时间的输出状态（电路的输出V）
+    V = calc_V(U, U0)
+    # 计算当前网络的能量E
+    energys[n] = calc_energy(V, distance)
+    # 检查路径的合法性
+    route, newV = check_path(V)
+    if len(np.unique(route)) == N:
+        route.append(route[0])
+        dis = calc_distance(route)
+        if dis < best_distance:  # 如果dis小于现有最好的best_distance则把best_distance替换为dis
+            H_path = []
+            # 更新dis
+            best_distance = dis
+            # 跟新route
+            best_route = route
+            [H_path.append((route[i], route[i + 1])) for i in range(len(route) - 1)]
+            print('第{}次迭代找到的次优解距离为：{}，能量为：{}，路径为：'.format(n, best_distance, energys[n]))
+            [print(chr(97 + v), end=',' if i < len(best_route) - 1 else '\n') for i, v in enumerate(best_route)]
+if len(H_path) > 0:
+    draw_H_and_E(citys, H_path, energys)
+else:
+    print('没有找到最优解')