titanic 神经网络 神经网络解决tsp问题

目录

• 基于CHNN求解TSP问题

• CHNN

• CHNN的网络结构
• CHNN状态方程
• CHNN的能量函数

• TSP问题
• 基于CHNN求解TSP问题的思路

• 1. 分析问题，将TSP的状态用数学符号描述出来

• 1.1 TSP问题描述
• 1.2 TSP问题的约束条件
• 1.3 TSP问题的目标函数

• 2. 构造网络能量函数
• 3. 优化网络能量函数
• 4. CHNN的动态方程
• 5. 状态更新方程

• python代码实现

基于CHNN求解TSP问题

CHNN

CHNN是连续型Hopfield神经网络。CHNN是以模拟量作为网络的输入和输出，而且各个神经元之间的工作方式是并行的。因此相对于DHNN，CHNN更接近于生物神经网络。CHNN更适合于求解组合优化问题。

CHNN的网络结构

CHNN的神经网络结构可以等效为放大电子电路。每一个神经元可以等效为一个电路放大器元件。CHNN的等效电路拓扑图如下:

其中：

$U_i$ 是放大电子元件的输入电压，对应于CHNN神经网络中神经元的输入，包括恒定的外部电流和其他电子元件的反馈连接。

$V_i$是放大电子元件的输出电压，对应于CHNN神经网络中神经元的输出，其输出有正负值，正值代表兴奋，负值代表抑制。

运算放大器$i$表示第$i$个神经元。

CHNN状态方程

设电容C的两端电压为$U_c$，存储的电荷量为Q，则有$C=\frac{Q}{U_c}$
则经过电容C的电流为：$i = \frac{dQ}{dt}=C\frac{dU_c}{dt}$
根据基尔霍夫电流定律，CHNN等效电路的电流关系为：
$C_i\frac{du_i}{dt}+\frac{u_i}{R_{i0}}=\sum_{j=1}^{n}\frac{1}{R_{ij}}(v_j-u_i)+I_i$
令$T_{ij}$表示神经网络中神经元之间的连接权重：$T_{ij}=\frac{1}{R_{ij}}$
则电流关系式可以化简为
$C_i\frac{du_i}{dt}=\sum_{j=1}^{n}T_{ij}(v_j-u_i)+I_i-T_{i0}u_i \\ \Rightarrow C_i\frac{du_i}{dt}=\sum_{j=1}^{n}T_{ij}v_j+I_i-T_{i}u_i$
上式就是CHNN神经网络的状态方程，其中$v_i=f_i(u_i)$，即输入电压是输出电压的非线性映射。$f_i$是S型激励函数

CHNN的能量函数

由于CHNN的神经网络模型中的权重是不变的，并且不需要学习，因此我们采用能量函数的方式来衡量网路的稳定性。
CHNN网络能量函数公式如下：
$E=-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}T_{ij}v_iv_j-\sum_{i=1}^{n}v_iI_i+\sum_{i=1}^{n}\frac{1}{R_i}\int_0^{u_i}f^{-1}(v_i)dv_i$

TSP问题

旅行商(TSP)问题是人工智能领域一个组合优化问题。

问题描述：有一个推销员，要到n个城市推销商品，他要找出一个包含所有n个城市的具有最短路程的环路。

基于CHNN求解TSP问题的思路

1. 分析问题，将TSP的状态用数学符号描述出来

1.1 TSP问题描述

我们假设一共有5个城市，A、B、C、D、E。我们使用一个$5\times5$矩阵表示每一种走法(状态)。

城市	第一步	第二步	第三步	第四步	第五步
A	0	0	1	0	0
B	1	0	0	0	0
C	0	1	0	0	0
D	0	0	0	1	0
E	0	0	0	0	1

如上表所示，$v_{xi}=1$ 表示第i步在第x城市, $v_{xi}=0$ 表示第i步不在第x城市 。
故，每个元素的取值为0或1
注：x,y表示城市，i，j表示的第几步

1.2 TSP问题的约束条件

TSP问题有三个约束条件

2. 一个城市只去一次 =====> 每一行只有一个‘1’
   $第x行的所有元素V_{xi}按顺序两两相乘之和为0 \\ \sum_{i=1}^{n-1}\sum_{j=i+1}^{n}V_{x_i}V_{x_j} = 0$



3. $那么所有行的所有元素V_{xi}按顺序两两相乘之和为0 \\ \sum_{x=1}^{n}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}V_{x_i}V_{x_j} = 0$
4. 一次只去一个城市 =====> 每一列只有一个‘1’
   $第i列的所有元素V_{xi}按顺序两两相乘之和为0 \\ \sum_{x=1}^{n-1}\sum_{y=x+1}^{n}V_{x_i}V_{y_i} = 0$
   $那么所有列的所有元素V_{xi}按顺序两两相乘之和为0 \\ \sum_{i=1}^{n}\sum_{x=1}^{n-1}\sum_{y=x+1}^{n}V_{x_i}V_{y_i} = 0$
5. 一共有n个城市 =====> 矩阵之和为n $\sum_{x=1}^{n}\sum_{i=1}^{n}V_{x_i}=n$

1.3 TSP问题的目标函数

我们想要得到的是访问n个城市的总距离最小。我们使用$d_{xy}$表示$x$与$y$两个城市的距离，那么顺序访问城市$x$和城市$y$，有$d_{xy}V_{xi}V_{y, i+1}和d_{xy}V_{xi}V_{y,i-1}$至少有一个为0。

那么顺序访问$x$和$y$的所有可能的路径有
$D_{xy}=\sum_{i=1}^{n}(d_{xy}V_{xi}V_{y, i+1}+d_{xy}V_{xi}V_{y,i-1}) \\ =\sum_{i=1}^{n}d_{xy}V_{xi}(V_{y, i+1}+V_{y,i-1})$

那么顺序访问所有城市的可能的路径有
$D=\sum_{x=1}^{n}\sum_{y=1}^{n}\sum_{i=1}^{n}d_{xy}V_{xi}(V_{y, i+1}+V_{y,i-1})$
我们想要得到总距离最小，即上式取最小值。

2. 构造网络能量函数

现在我们就把求解TSP问题转化为求解有三个约束条件下求解最小值的最优化问题。
我们根据三个约束条件和一个优化目标函数来构造网络能量函数
$E=\frac{A}{2}\sum_{x=1}^{n}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}V_{x_i}V_{x_i} \\ +\frac{B}{2} \sum_{i=1}^{n}\sum_{x=1}^{n-1}\sum_{y=x+1}^{n}V_{x_i}V_{y_i} \\ +\frac{C}{2}(\sum_{x=1}^{n}\sum_{i=1}^{n}V_{x_i}-n)^2 \\ +\frac{D}{2}\sum_{x=1}^{n}\sum_{y=1}^{n}\sum_{i=1}^{n}d_{xy}V_{xi}(V_{y, i+1}+V_{y,i-1})$
其中A、B、C、D是常数。
当E达到最小值的时候，就选择了最优路径。

3. 优化网络能量函数

在’Theories on the Hopfield neural networks’这篇论文中对上式进行了改进，提高了收敛速度。改进后的公式如下：
$E=\frac{A}{2}\sum_{x=1}^{n}(\sum_{i=1}^{n}V_{xi}-1)^2+\frac{A}{2}\sum_{i=1}^{n}(\sum_{x=1}^{n}V_{xi}-1)^2+\frac{D}{2}\sum_{x=1}^{n}\sum_{y=1}^{n}\sum_{i=1}^{n}V_{xi}d_{xy}V_{y, i+1}$

4. CHNN的动态方程

根据优化后的能量函数可以得到动态方程如下：
$\frac{dU_{ij}}{dt}=-\frac{\partial E}{\partial v_{xi}}=-A(\sum_{i=1}^{n}V_{xi}-1)-A(\sum_{x=1}^{n}V_{xi}-1)-D\sum_{y=1}^{n}d_{xy}v_{y,i+1}$

5. 状态更新方程

输入状态的更新方程如下：
$U_{ij}(t+1)=U_{ij}(t)+\frac{dU_{ij}}{dt}\Delta t$
输出状态的更新方程为(激活函数使用sigmoid函数)：
$V_{ij}=\frac{1}{2}(1+tanh(\frac{U_{ij}}{U_0}))$

python代码实现

下面我们就可以根据以上过程来实现代码

import numpy as np
import random
from math import sqrt, log, tanh, exp
from matplotlib import pyplot as plt
import logging
logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(filename)s[line:%(lineno)d]- %(levelname)s:%(message)s')

# 动态方程中的两个系数
A = 300
D = 200
U0 = 0.1
max_iter = 120  # 最大迭代次数
step = 0.0001  # 步长


def d_u(state_v, distance, n):
    """ 动态方程
    :param state_v: 输出矩阵
    :param distance: 城市之间的距离矩阵
    :param n: 城市的个数
    :return: 返回 动态方程的矩阵
    """
    a = np.sum(state_v, axis=0) - 1  # 得到行向量，把每一列的元素都加到第一行  对前一个下标求和
    b = np.sum(state_v, axis=1) - 1  # 得到列向量(列向量的值，行向量的形式)，把每一行的元素都加到第一列  对后一个下标求和
    # 把上述两个向量扩展为两个矩阵
    x = np.array([[0.0 for i in range(n)] for i in range(n)])   # 用来扩展a向量
    y = np.array([[0.0 for i in range(n)] for i in range(n)])   # 用来扩展b向量
    # 扩展a向量为一个矩阵
    for i in range(n):
        for j in range(n):
            x[i, j] = a[j]  # x的每一行都和a向量相同，相当于x矩阵是n行a
    # 扩展b向量为一个矩阵
    for i in range(n):
        for j in range(n):
            y[j, i] = b[j]  # y的每一列都和b相同，相当于y矩阵是n列b，但是b在python里的形式是行向量

    # 创建C矩阵: 将V矩阵的第一列移到最后一列，并与距离矩阵相乘？？？ 为什么要移动第一列
    c = np.array([[0.0 for i in range(n)] for i in range(n)])
    c[:, 0:n-1] = state_v[:, 1:n]
    c[:, n-1] = state_v[:, 0]  # 将V的第一列放到最后一列
    c = np.dot(distance, c)  # 距离矩阵乘以V
    return -A * (x + y) - D * c


def energy(state_v, distance, n):
    """ 能量函数
    :param state_v: 输出矩阵
    :param distance: 距离矩阵
    :param n: 城市的个数
    :return: 能量函数的值,是一个数
    """
    a = sum([x*x for x in (np.sum(state_v, axis=0)-1)])
    b = sum([x*x for x in (np.sum(state_v, axis=1)-1)])
    c = np.array([[0.0 for i in range(n)] for i in range(n)])
    c[:, 0:n-1] = state_v[:, 1:n]
    c[:, n-1] = state_v[:, 0]  # 将V的第一列放到最后一列
    c = np.sum(np.sum(np.multiply(state_v, np.dot(distance, c)), axis=0))
    return 0.5 * (A * (a + b) + D * c)

def get_v(state_v, n):
    ''' 得到稳定状态下的矩阵信息
    :param state_v: 稳定输出矩阵
    :param n: 城市的个数
    :return:
    '''
    (row, col) = state_v.shape
    V_H = np.array([0.0 for i in range(row)])  # 用来存放每一列的最大值，共row行
    V_W = np.array([0 for i in range(row)])  # 用来存放每一列最大值的行号
    for i in range(n):
        for j in range(n):
            if V[j, i] > V_H[i]:
                V_H[i] = V[j, i]
                V_W[i] = j
    # 创建一个数组将V中每列最大值的位置设为1，其他设为0
    V_1 = np.array([[0 for i in range(row)] for i in range(col)])
    for i in range(col):
        V_1[V_W[i], i] = 1
    return V_H, V_W, V_1

# 归一化处理用到的函数
def normalization(energy_all):
    # return [float(x-energy_all.mean())/energy_all.std() for x in energy_all]  #z-score
    return [log(x)/log(np.max(energy_all)) for x in energy_all]   #log函数转化

def draw_enegry(energy_all):
    ''' 画出能量函数的曲线图
    :param energy_all: 能量函数数组
    :return: 无返回值
    '''
    y = normalization(energy_all)  # 对能量函数值进行归一化处理
    plt.plot(y)
    plt.xlabel('迭代次数', loc='right')
    plt.ylabel('能量函数值', loc='top')
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文字体设置-黑体
    plt.show()

def draw_rode(state_v, V_1, V_W, n):
    ''' 画出路线图
    :param state_v: 稳定输出矩阵
    :param V_1: 0-1状态矩阵
    :param V_W: 每列最大值所在的行向量
    :param n: 城市个数
    :return:
    '''
    (row, col) = state_v.shape
    V_1H = np.array([0.0 for i in range(n)])
    V_1W = np.array([0 for i in range(n)])
    for i in range(col):
        V_1[V_W[i], i] = 1
    for i in range(n):
        for j in range(n):
            if V_1[j, i] > V_1H[i]:
                V_1H[i] = V_1[j, i]
                V_1W[i] = j
    # 随机一种路径和优化后的路径进行比较
    city_i = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    V_0W = random.sample(list(city_i), 10)
    # 创建两个城市矩阵，将他们按照我们所给的初始顺序和最终顺序进行排序
    city_b = np.array(
        [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0],
         [0.0, 0.0]])
    city_f = np.array(
        [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0],
         [0.0, 0.0]])
    k = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    k = [city[i] for i in V_0W]
    for i in range(n):
        city_b[i] = k[i]
    j = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    j = [city[i] for i in V_1W]
    for i in range(n):
        city_f[i] = j[i]
    # 下面将坐标点的X，Y轴数值分别提取进行画图
    X_b = []
    Y_b = []
    X_f = []
    Y_f = []
    for i in range(n):
        X_b.append(city_b[i][0])
        Y_b.append(city_b[i][1])
        X_f.append(city_f[i][0])
        Y_f.append(city_f[i][1])
    # 将第一个点加入，因为要回到初始点
    X_b.append(city_b[0][0])
    Y_b.append(city_b[0][1])
    X_f.append(city_f[0][0])
    Y_f.append(city_f[0][1])
    plt.plot(X_b, Y_b, 'r', X_f, Y_f, '-.b')   # 红线是随机，蓝线是优化后的结果
    plt.scatter(X_f, Y_f)
    plt.legend(labels=['随机路线', '优化后路线'], loc='lower right', fontsize=6)    # 图例
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文字体设置-黑体
    plt.show()


if __name__ == '__main__':
    city = np.array([
        [0.7000, 0.2000], [0.4000, 0.3000], [0.5000, 0.8000],
        [0.3000, 0.4000], [0.1000, 0.9000], [0.9000, 0.4000],
        [0.8000, 0.6000], [0.6000, 0.9000], [0.3000, 0.6000],
        [0.2000, 0.8000]])  # 10个城市的坐标
    n = len(city)  # 城市的个数

    # 计算距离矩阵
    distance = np.array([[0.0 for i in range(n)] for i in range(n)])
    for i in range(n):
        for j in range(n):
            a = sqrt((city[i, 0]-city[j, 0]) ** 2 + (city[i, 1]-city[j, 1]) ** 2)
            distance[i, j] = a
            distance[j, i] = a

    # 随机给定网络初始状态矩阵
    delta = 2 * (np.random.random((n, n))) - 1   # 随机产生一个n*n的矩阵，用来随机初始化状态矩阵
    U = U0 * log(n-1) + delta  # 随机定义一个初始输入矩阵
    V = np.array([[0.0 for i in range(n)] for i in range(n)])
    for i in range(n):
        for j in range(n):
            # V[i, j] = (1+2/(1+exp(-2*(U[i, j]/U0)))-1)/2
            V[i, j] = 0.5 * (1 + tanh(U[i, j] / U0))   # 输出矩阵、状态矩阵

    # 创建一个向量，用来存储每一步的能量值
    energy_all = np.array([0.0 for i in range(max_iter)])

    # 迭代
    for k in range(max_iter):
        du = d_u(V, distance, n)  # 动态方程、梯度
        U = U + du * step  # 更新输入矩阵
        for i in range(n):
            for j in range(n):
                V[i, j] = 0.5 * (1 + tanh(U[i, j] / U0))  # 输出矩阵、状态矩阵
                # V[i, j] = (1 + 2 / (1 + exp(-2 * (U[i, j] / U0))) - 1) / 2
        energy_all[k] = energy(V, distance, n)  # 求能量函数的值，放入向量中
    Vs = get_v(V, n)
    draw_enegry(energy_all)  # 画出能量函数图
    draw_rode(V, V_W=Vs[1], V_1=Vs[2], n=n)  # 画出路线图