关于最优化方法理论的算法有很多,但是每种算法都有一个适用的场景。所以根据每种算法的特性来合理地使用和升级算法才是最重要的。学习最优化方法的知识背景如下(并不是所有都需要):

必需知识基础:

  1. 高等数学、微积分
  2. 线性代数
  3. 概率论与数理统计

知识技能提升:

  1. 集合论与图论
  2. 运筹学
  3. 数据结构

知乎和CSDN上有很多大神总结的笔记非常好,通俗易懂,还有一些相关公众号也不错,可以自行搜索收藏。其实关于最优化方法的许多方法还是用MATLAB实现起来更方便,Python也有很多相关的库可以直接调用。这里只是通过Python手撕最优化方法的算法的方式来进一步加深对这些理论的理解。

手撕单纯形法(Simplex Method)

参考博客1:https://blog.csdn.net/weixin_43902708/article/details/89253943?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161570979116780261967484%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=161570979116780261967484&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2alltop_click~default-1-89253943.first_rank_v2_pc_rank_v29&utm_term=%E5%8D%95%E7%BA%AF%E5%BD%A2%E6%B3%95

参考博客2:https://blog.csdn.net/weixin_40955254/article/details/85073962?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161573480016780264052316%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=161573480016780264052316&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allsobaiduend~default-4-85073962.first_rank_v2_pc_rank_v29&utm_term=%E5%8D%95%E7%BA%AF%E5%BD%A2%E6%B3%95

算法实现

单纯形法是一种线性规划算法,在基本单纯形法的基础上还衍生出了优化单纯形法和对偶单纯形法。

单纯形法的代码实现如下(最基本的):

  1. import numpy as np
  2. EQUATION_NUM = 3
  3. A_MATRIX = np.array([[0.25,0.50,1,0,0,120],[0.50,0.50,0,1,0,150],[0.25,0,0,0,1,50],[-12,-15,0,0,0,0]])
  4. IN_VARIABLE = [3,4,5]
  5. #基本步骤:检查系数小于0——找到f(z)中各个表达变量并锁定各个表达变量
  6. #对应的方程——进行矩阵变换——循环知道系数全部大于0
  7. #设置零值为非常小的正数
  8. def add_tiny(x):
  9. if x == 0:
  10. return x + 0.002
  11. else:
  12. return x
  13. process = A_MATRIX.copy()
  14. while(True):
  15. weightlist = process[EQUATION_NUM].copy()[:-1]
  16. piece = process[:EQUATION_NUM].copy()
  17. #检查是不是f(z)的所有表示项前的系数都大于0,都大于0则循环中止
  18. if (weightlist >= 0).all():
  19. break
  20. #得到b向量
  21. b = [i[5] for i in piece]
  22. b = np.array(b)
  23. now_variable = np.argmin(weightlist)
  24. temp = [add_tiny(i[now_variable]) for i in piece]
  25. temp = np.array(temp)
  26. temp2 = b/temp
  27. #找到锁定行
  28. now_row = np.argmin(temp2)
  29. #修改换入变量
  30. IN_VARIABLE[now_row] = now_variable + 1
  31. #矩阵变换
  32. for i in range(EQUATION_NUM+1):
  33. if i != now_row:
  34. ratio = process[i][now_variable]/process[now_row][now_variable]
  35. process[i] -= ratio * process[now_row]
  36. else:
  37. ratio = 1.0/process[i][now_variable]
  38. process[i] *= ratio
  39. print(process)

手撕遗传算法(Genetic Algorithm)

参考博客1:https://blog.csdn.net/ha_ha_ha233/article/details/91364937

参考博客2:https://blog.csdn.net/weixin_41503009/article/details/104782088

思维框架

算法实现

遗传算法是现代智能优化算法中应用最广的一类算法了,很多最优化方法都是从遗传算法/进化算法衍生而来。

遗传算法实现(截至选择部分):

  1. import numpy as np
  2. import math
  3. import matplotlib.pyplot as plt
  4. from matplotlib import cm
  5. from mpl_toolkits.mplot3d import Axes3D
  6. N_GENERATIONS = 50
  7. DNA_LEN=10
  8. POP_SIZE=100
  9. XBound=[-1,2]
  10. #f(x)为需要求解的函数
  11. def f(pop):
  12. flist=[]
  13. for one in pop:
  14. flist.append(one * np.sin(10*math.pi*one) + 2)
  15. return flist
  16. #prese(x)为把给定的数值矩阵压缩为指定范围的列表
  17. def press(fx,Xbound):
  18. coeff = np.array([2**i for i in range(len(fx[0]))])
  19. ans = np.dot(fx,coeff)/(2**(len(fx[0]))-1)
  20. mappx = [i * (Xbound[1]-Xbound[0])+ Xbound[0] for i in ans]
  21. return mappx
  22. #这是根据每个点的值得到每个点的适应度
  23. def get_fitness(fpred):
  24. return fpred-np.min(fpred)+1e-3
  25. #选择函数,得到新的种群矩阵
  26. def select(pop, fitness): # nature selection wrt pop's fitness
  27. temp=np.array(fitness)
  28. idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,
  29. p=(temp)/(temp.sum()) )
  30. return pop[idx]
  31. pop=np.random.randint(low=0,high=2,size=(POP_SIZE,DNA_LEN))
  32. save_1 = f(press(pop,XBound))
  33. fitness=get_fitness(save_1)
  34. select=select(pop,save_1)
  35. print(select)

至此我们已经基本完成了遗传算法的骨干,后面还要实现交叉和变异。需要交叉和变异的原因是:我们在利用np.random.choice函数进行样本点的重新筛选时有可能没有选到真正的最大点,从而进入局部最优解而不是全局最优解,所以需要此步骤提供其他点的可能性。所以下面我们来实现下交叉和变异。

遗传算法实现(包括选择、交叉和变异):

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. import seaborn as sns
  4. sns.set()
  5. '''
  6. 下面是几个城市基本参数的设置,N_GENERATIONS是迭代的次数
  7. 禁忌长度为2,候选集合长度为6
  8. '''
  9. CHENGDU = [1,1]
  10. SHANGHAI = [5,2]
  11. XIAN = [2,3]
  12. WUHAN = [3,2]
  13. BEIJING = [4,3.5]
  14. SHENZHEN = [3.5,0]
  15. FUZHOU = [4.5,0.5]
  16. CITY_LEN = 7
  17. CHENGSHI = [CHENGDU,SHANGHAI,XIAN,WUHAN,BEIJING,SHENZHEN,FUZHOU]
  18. N_GENERATIONS = 15
  19. Tabu_list = [10,10]
  20. #计算两个城市间的距离
  21. def distance(city1,city2):
  22. temp = [i-j for i,j in zip(city1,city2)]
  23. return (sum([x**2 for x in temp]))**(1/2)
  24. #计算整个周游下来的总路程,即适应度函数
  25. def total_distance(cityset):
  26. citycopy = cityset[:]
  27. citycopy.append(cityset[0])
  28. dist = 0
  29. for i in range(len(cityset)):
  30. dist += distance(citycopy[i],citycopy[i+1])
  31. return dist
  32. #产生候选集合,长度为6
  33. def alter_cityset(cityset):
  34. alterset=[]
  35. for i in range(CITY_LEN - 1):
  36. citys = cityset[:]
  37. if i != CITY_LEN - 2:
  38. citys[1+i],citys[2+i] = citys[2+i],citys[1+i]
  39. else:
  40. citys[1+i],citys[1] = citys[1],citys[1+i]
  41. alterset.append(citys)
  42. return alterset
  43. #计算长度为6的候选集合中各个路径对应的总路程,即计算候选的适应度函数值
  44. def set_to_distance(altercities):
  45. dist = []
  46. for altercity in altercities:
  47. dist.append(total_distance(altercity))
  48. return dist
  49. #打印当前的路程情况
  50. def print_plot(cityset):
  51. x_list=[];y_list=[]
  52. for city in cityset:
  53. x_list.append(city[0])
  54. y_list.append(city[1])
  55. x_list.append(cityset[0][0])
  56. y_list.append(cityset[0][1])
  57. plt.plot(x_list,y_list)
  58. plt.scatter(x_list,y_list, linewidth=1, alpha=0.6)
  59. plt.scatter(x_list[0],y_list[0],color="red",marker="o")
  60. plt.scatter(x_list[1],y_list[1],color="lawngreen",marker="o")
  61. print('最优适应度为:',total_distance(cityset))
  62. plt.show()
  63. #完成禁忌表Tabu_list的变换
  64. def swtich_tabu(number=0):
  65. Tabu_list.pop(0)
  66. Tabu_list.append(number)
  67. #找到列表中第n小元素的索引
  68. def find_min(list1=None,n=1):
  69. temp = list1[:];count = 0
  70. while count < n-1:
  71. temp[np.argmin(temp)]=np.max(temp)
  72. count += 1
  73. return np.argmin(temp)
  74. #这就是循环主体部分,不断更新目前最优路径和最优适应度
  75. def explore(chengshi,optimization,route_best):
  76. alter = alter_cityset(chengshi)
  77. dis = set_to_distance(alter)
  78. min_index = np.argmin(dis)
  79. route_be = []
  80. if dis[min_index] < optimization:
  81. opti = dis[min_index]
  82. route_be = city_now = alter[min_index]
  83. swtich_tabu(min_index)
  84. elif min_index not in Tabu_list:
  85. [opti,route_be] = [optimization,route_best]
  86. city_now = alter[min_index]
  87. swtich_tabu(min_index)
  88. else:
  89. min2_index = find_min(dis,2)
  90. if min2_index not in Tabu_list:
  91. [opti,route_be] = [optimization,route_best]
  92. city_now = alter[min2_index]
  93. swtich_tabu(min2_index)
  94. else:
  95. min3_index = find_min(dis,3)
  96. [opti,route_be] = [optimization,route_best]
  97. city_now = alter[min3_index]
  98. swtich_tabu(min3_index)
  99. return [city_now,opti,route_be]
  100. if __name__ == "__main__":
  101. route_now = CHENGSHI
  102. route_best = []
  103. np.random.shuffle(route_now)
  104. optimization = total_distance(CHENGSHI)
  105. print('初始适应度:',optimization)
  106. #在迭代次数内不断迭代
  107. for count in range(N_GENERATIONS):
  108. [route_now,optimization,route_best] = explore(route_now,optimization,route_best)
  109. print('当前最优路径是:',route_now)
  110. print('当前最优适应度是:',optimization)
  111. [route_now,optimization,route_best] = explore(route_now,optimization,route_best)
  112. print_plot(route_best)

手撕模拟退火算法(SA Algorithm)

参考博客:https://blog.csdn.net/WFRainn/article/details/80303138?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161451812216780271524683%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=161451812216780271524683&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allsobaiduend~default-2-80303138.first_rank_v2_pc_rank_v29&utm_term=%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB%E7%AE%97%E6%B3%95Python

思维框架

算法实现

模拟退火算法——最经典的优化算法,收敛较慢但是确实容易找到全局最优解。

模拟退火算法实现:

  1. import numpy as np
  2. import random
  3. import math
  4. '''
  5. 下面几个初始化参数聪上往下分别为:求解范围的上、下界、温度降低的速率以及初始温度
  6. '''
  7. HIGH = 3
  8. LOW = -3
  9. ALF = 0.95
  10. Temper = 50
  11. #目标函数
  12. def f(x):
  13. return 11*np.sin(x) + 7 * np.cos(5*x)
  14. #概率函数p,这里只考虑下山的情况
  15. def p(x,x1):
  16. return math.exp(-abs(f(x)-f(x1))/Temper)
  17. #新解的随机产生
  18. def new(x):
  19. x1 = x + Temper*random.uniform(-1,1)
  20. if (x1 <= HIGH) & (x1 >= LOW):
  21. return x1
  22. elif x1 < LOW:
  23. rand = random.uniform(-1,1)
  24. return rand * LOW + (1 - rand) * x
  25. else:
  26. rand = random.uniform(-1,1)
  27. return rand * HIGH + (1 - rand) * x
  28. def main():
  29. global x,Temper
  30. x=random.uniform(LOW,HIGH)
  31. while Temper > 0.001:
  32. for i in range(500):
  33. x1 =new(x)
  34. if f(x1) <= f(x):
  35. x = x1
  36. elif random.random() <= p(x,x1):
  37. x = x1
  38. else:
  39. continue
  40. Temper = Temper*ALF
  41. print('最小值为:{}'.format(f(x)))
  42. if __name__=="__main__":
  43. main()

手撕禁忌搜索算法(Tabu Search)

参考博客1:https://blog.csdn.net/adkjb/article/details/81712969

参考博客2:https://blog.csdn.net/ChinaJane163/article/details/49095085

参考博客3:https://blog.csdn.net/tyhj_sf/article/details/54235550

参考博客4:https://blog.csdn.net/weixin_42528077/article/details/86762878

思维框架

算法实现

禁忌搜索算法是一类路径搜索类算法,用于确定多点间的最小闭路径。

禁忌搜索算法实现:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

'''
下面是几个城市基本参数的设置,N_GENERATIONS是迭代的次数
禁忌长度为2,候选集合长度为6
'''
CHENGDU = [1,1]
SHANGHAI = [5,2]
XIAN = [2,3]
WUHAN = [3,2]
BEIJING = [4,3.5]
SHENZHEN = [3.5,0]
FUZHOU = [4.5,0.5]
CITY_LEN = 7
CHENGSHI = [CHENGDU,SHANGHAI,XIAN,WUHAN,BEIJING,SHENZHEN,FUZHOU]
N_GENERATIONS = 15
Tabu_list = [10,10]

#计算两个城市间的距离
def distance(city1,city2):
    temp = [i-j for i,j in zip(city1,city2)]
    return (sum([x**2 for x in temp]))**(1/2)

#计算整个周游下来的总路程,即适应度函数
def total_distance(cityset):
    citycopy = cityset[:]
    citycopy.append(cityset[0])
    dist = 0
    for i in range(len(cityset)):
        dist += distance(citycopy[i],citycopy[i+1])
    return dist

#产生候选集合,长度为6
def alter_cityset(cityset):
    alterset=[]
    for i in range(CITY_LEN - 1):
        citys = cityset[:]
        if i != CITY_LEN - 2:
            citys[1+i],citys[2+i] = citys[2+i],citys[1+i]
        else:
            citys[1+i],citys[1] = citys[1],citys[1+i]
        alterset.append(citys)
    return alterset

#计算长度为6的候选集合中各个路径对应的总路程,即计算候选的适应度函数值
def set_to_distance(altercities):
    dist = []
    for altercity in altercities:
        dist.append(total_distance(altercity))
    return dist

#打印当前的路程情况
def print_plot(cityset):
    x_list=[];y_list=[]
    for city in cityset:
        x_list.append(city[0])
        y_list.append(city[1])
    x_list.append(cityset[0][0])
    y_list.append(cityset[0][1])
    plt.plot(x_list,y_list)
    plt.scatter(x_list,y_list, linewidth=1, alpha=0.6)
    plt.scatter(x_list[0],y_list[0],color="red",marker="o")
    plt.scatter(x_list[1],y_list[1],color="lawngreen",marker="o")
    print('最优适应度为:',total_distance(cityset))  

#完成禁忌表Tabu_list的变换
def swtich_tabu(number=0):
    Tabu_list.pop(0)
    Tabu_list.append(number)

#找到列表中第n小元素的索引
def find_min(list1=None,n=1):
    temp = list1[:];count = 0
    while count < n-1:
        temp[np.argmin(temp)]=np.max(temp)
        count += 1
    return np.argmin(temp)

#这就是循环主体部分,不断更新目前最优路径和最优适应度
def explore(chengshi,optimization,route_best):
    alter = alter_cityset(chengshi)
    dis = set_to_distance(alter)
    min_index = np.argmin(dis)
    route_be = []
    if dis[min_index] < optimization:
        opti = dis[min_index]
        route_be = city_now = alter[min_index]
        swtich_tabu(min_index)
    elif min_index not in Tabu_list:
        [opti,route_be] = [optimization,route_best]
        city_now = alter[min_index]
        swtich_tabu(min_index)    
    else:
        min2_index = find_min(dis,2)
        if min2_index not in Tabu_list:
            [opti,route_be] = [optimization,route_best]
            city_now = alter[min2_index]
            swtich_tabu(min2_index)
        else:
            min3_index = find_min(dis,3)
            [opti,route_be] = [optimization,route_best]
            city_now = alter[min3_index]
            swtich_tabu(min3_index)
    return [city_now,opti,route_be]

if __name__ == "__main__":
    route_now = CHENGSHI
    route_best = []
    np.random.shuffle(route_now)
    optimization = total_distance(CHENGSHI)
    print('初始适应度:',optimization)
    for count in range(N_GENERATIONS):
        [route_now,optimization,route_best] = explore(route_now,optimization,route_best)
        print('当前最优路径是:',route_now)
        print('当前最优适应度是:',optimization)
    [route_now,optimization,route_best] = explore(route_now,optimization,route_best)
    print_plot(route_best)

手撕蚁群算法(Antcolony Algorithm)

参考博客1:https://blog.csdn.net/ACA_JSP/article/details/88362006?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161466651216780255259333%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=161466651216780255259333&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allfirst_rank_v2~rank_v29-8-88362006.first_rank_v2_pc_rank_v29&utm_term=%E8%9A%81%E7%BE%A4%E7%AE%97%E6%B3%95

参考博客2:https://blog.csdn.net/fanxin_i/article/details/80380733?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161461094416780266251736%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=161461094416780266251736&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allsobaiduend~default-2-80380733.first_rank_v2_pc_rank_v29&utm_term=%E8%9A%81%E7%BE%A4%E7%AE%97%E6%B3%95Python

参考博客3:https://blog.csdn.net/qq_38526623/article/details/105318909?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161461094416780266251736%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=161461094416780266251736&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2alltop_click~default-1-105318909.first_rank_v2_pc_rank_v29&utm_term=%E8%9A%81%E7%BE%A4%E7%AE%97%E6%B3%95Python

算法实现

蚁群算法是应用极其普遍的一类算法,可以不断优化,且在数模竞赛中使用极广。

蚁群算法如下(下面只是把基本的蚁群算法实现了下,参数还需要调整,模型也可以继续优化):

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
'''
下面是几个城市基本参数的设置
蚁群蚂蚁数量为10
'''
CHENGDU = [1,1]
SHANGHAI = [5,2]
XIAN = [2,3]
WUHAN = [3,2]
BEIJING = [4,3.5]
SHENZHEN = [3.5,0]
FUZHOU = [4.5,0.5]
CITY_LEN = 7
ANT_NUMBER = 35
MOVE_TIMES = 20
CHENGSHI = [CHENGDU,SHANGHAI,XIAN,WUHAN,BEIJING,SHENZHEN,FUZHOU]

#计算两个城市间的距离
def distance(city1,city2):
    temp = [i-j for i,j in zip(city1,city2)]
    return (sum([x**2 for x in temp]))**(1/2)

#计算周游一圈总的长度
def total_distance(cityset):
    citycopy = cityset[:]
    citycopy.append(cityset[0])
    dist = 0
    for i in range(len(cityset)):
        dist += distance(citycopy[i],citycopy[i+1])
    return dist

#根据城市列表生成距离矩阵
def dist_matrix(cityset):
    dist_matrix = np.zeros((CITY_LEN,CITY_LEN))
    for i in range(CITY_LEN):
        for j in range(CITY_LEN):
            dist_matrix[i][j] = distance(cityset[i],cityset[j])
    return dist_matrix

#生成初始信息素矩阵
def phe_matrix():
    phe_matrix = np.diag([-5,-5,-5,-5,-5,-5,-5])
    phe_matrix += 5
    return phe_matrix

#求解最后的预测路径
def find_route(mat,cityset):
    list1 = []
    list1.append(0)
    count = 0
    while(True):
        data = mat[count]
        for item in list1:
            data[item] = 0
        count = max_index = np.argmax(data)
        list1.append(max_index)
        if(len(list1) == 7):
            break
    route = []
    for number in list1:
        route.append(cityset[number])
    return route

#打印当前的路程情况
def print_plot(cityset):
    x_list=[];y_list=[]
    for city in cityset:
        x_list.append(city[0])
        y_list.append(city[1])
    x_list.append(cityset[0][0])
    y_list.append(cityset[0][1])
    plt.plot(x_list,y_list)
    plt.scatter(x_list,y_list, linewidth=1, alpha=0.6)
    plt.scatter(x_list[0],y_list[0],color="red",marker="o")
    plt.scatter(x_list[1],y_list[1],color="lawngreen",marker="o")
    print('最短路径长度为:',total_distance(cityset))

#蚂蚁类,有当前位置和移动禁忌表两个元素
class ant:
    def __init__(self):
        self.position = np.random.randint(low=0,high=CITY_LEN)
        self.Tabu_list = []

    def update_Tabu(self):
        if len(self.Tabu_list) >= CITY_LEN - 1:
            self.Tabu_list = []
        else:
            self.Tabu_list.append(self.position)

    def info(self):
        print('当前此蚂蚁所在位置为:',self.position+1)
        if len(self.Tabu_list)==0:
             print('当前此蚂蚁移动禁忌表为:',self.Tabu_list)
        else:
           print('当前此蚂蚁移动禁忌表为:',[i+1 for i in self.Tabu_list]) 

    def move(self,new_position):
        self.update_Tabu()
        self.position = new_position

#蚁群类,有蚁群蚂蚁数量和十个蚂蚁两个元素
class ants:
    def __init__(self,ant_number):
        self.number = ant_number
        self.antss = []
        for i in range(ant_number):
            self.antss.append(ant())

    def info(self):
        for i in range(self.number):
            print('蚂蚁',i+1)
            self.antss[i].info()

#路径类,有距离矩阵、信息素矩阵和比值矩阵三个元素
#距离矩阵一旦生成不再改变,信息素矩阵会更新
class route_matrix:
    def __init__(self,cityset):
        self.dist_matrix = dist_matrix(cityset)
        self.phe_matrix = phe_matrix()
        for i in range(CITY_LEN):
            self.dist_matrix[i][i] = 1
        self.index_matrix = self.phe_matrix/self.dist_matrix

    def update_index(self,pos1,pos2):
        self.index_matrix[pos1][pos2] += 3.0/self.dist_matrix[pos1][pos2]

    def info(self):
        #print('当前路径的距离矩阵为:')
        #print(self.dist_matrix)
        print('当前路径的信息素矩阵为:')
        print(self.phe_matrix)
        print('当前路径的比值矩阵为:')
        print(self.index_matrix)

#蚁群类,包含了蚂蚁和路程矩阵两个元素
class antcolony:
     def __init__(self,ants,route_matrix):
         self.ants = ants
         self.matrix = route_matrix

     def info(self):
         self.ants.info()
         self.matrix.info()

     def move_prob(self):
         prob_lists = []
         for ant in self.ants.antss:
             temps = self.matrix.index_matrix[ant.position][:]
             prob_list = []
             for temp in temps:
                 prob_list.append(temp)
             for Tabu in ant.Tabu_list:
                 prob_list[Tabu] = 0
             new_list = [x/sum(prob_list) for x in prob_list]
             prob_lists.append(new_list)
         return prob_lists

     def update_matrix(self,pos1,pos2):
         self.matrix.phe_matrix[pos1][pos2] += 3.0/self.matrix.dist_matrix[pos1][pos2]
         self.matrix.update_index(pos1,pos2)

     def wholemove(self):
         if len(self.ants.antss[0].Tabu_list) < 6:
             prob_lists = self.move_prob()
         for i in range(ANT_NUMBER):
             if len(self.ants.antss[i].Tabu_list) < 6:
                 idx = np.random.choice(np.arange(CITY_LEN),p=prob_lists[i])
                 #print(i,"的当前位置为:",self.ants.antss[i].position+1)
                 #print(i,'选择的下一位置为:',idx+1)
                 last_pos = self.ants.antss[i].position
                 self.ants.antss[i].move(idx)
                 self.update_matrix(last_pos,idx)
             else:
                 last_pos = self.ants.antss[i].position
                 origin = self.ants.antss[i].Tabu_list[0]
                 self.ants.antss[i].move(origin)
                 self.update_matrix(last_pos,origin)

if __name__ == "__main__":
    np.random.shuffle(CHENGSHI)
    route = route_matrix(CHENGSHI)
    ants = ants(ANT_NUMBER)
    antcolony = antcolony(ants,route)
    for count in range(MOVE_TIMES):
        antcolony.wholemove()
    matrix = antcolony.matrix.index_matrix
    answer = find_route(matrix,CHENGSHI)
    antcolony.info()
    print_plot(answer)

蚁群算法最大的优点就是:不容易陷入局部最优解,但是与此同时其找寻全局最优解的速度也变慢了。也即不容易早熟收敛,但是收敛也忒慢了点!所以下面优化后的代码是采取了最优-最差蚂蚁策略的蚂蚁系统:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
'''
下面是几个城市基本参数的设置
蚁群蚂蚁数量为10
'''
CHENGDU = [1,1]
SHANGHAI = [5,2]
XIAN = [2,3]
WUHAN = [3,2]
BEIJING = [4,3.5]
SHENZHEN = [3.5,0]
FUZHOU = [4.5,0.5]
CITY_LEN = 7
ANT_NUMBER = 35
MOVE_TIMES = 20
CHENGSHI = [CHENGDU,SHANGHAI,XIAN,WUHAN,BEIJING,SHENZHEN,FUZHOU]

#计算两个城市间的距离
def distance(city1,city2):
    temp = [i-j for i,j in zip(city1,city2)]
    return (sum([x**2 for x in temp]))**(1/2)

#计算周游一圈总的长度
def total_distance(cityset):
    citycopy = cityset[:]
    citycopy.append(cityset[0])
    dist = 0
    for i in range(len(cityset)):
        dist += distance(citycopy[i],citycopy[i+1])
    return dist

#根据城市列表生成距离矩阵
def dist_matrix(cityset):
    dist_matrix = np.zeros((CITY_LEN,CITY_LEN))
    for i in range(CITY_LEN):
        for j in range(CITY_LEN):
            dist_matrix[i][j] = distance(cityset[i],cityset[j])
    return dist_matrix

#根据数字列表生成的路径
def pro_route(index_list,cityset=CHENGSHI):
    alter = []
    for index in index_list:
        alter.append(cityset[index])
    return alter

#生成初始信息素矩阵
def phe_matrix():
    phe_matrix = np.diag([-5,-5,-5,-5,-5,-5,-5])
    phe_matrix += 5
    return phe_matrix

#将路程表分解为所有相邻两个元素组成的列表
def decompose(index_list):
    indexcopy = index_list[:]
    indexcopy.append(index_list[0])
    decom = []
    for i in range(len(index_list)):
        part = indexcopy[i:i+2]
        decom.append(part)
    return decom

#求解最后的预测路径
def find_route(mat,cityset):
    list1 = []
    list1.append(0)
    count = 0
    while(True):
        data = mat[count]
        for item in list1:
            data[item] = 0
        count = max_index = np.argmax(data)
        list1.append(max_index)
        if(len(list1) == 7):
            break
    route = []
    for number in list1:
        route.append(cityset[number])
    return route

#打印当前的路程情况
def print_plot(cityset):
    x_list=[];y_list=[]
    for city in cityset:
        x_list.append(city[0])
        y_list.append(city[1])
    x_list.append(cityset[0][0])
    y_list.append(cityset[0][1])
    plt.plot(x_list,y_list)
    plt.scatter(x_list,y_list, linewidth=1, alpha=0.6)
    plt.scatter(x_list[0],y_list[0],color="red",marker="o")
    plt.scatter(x_list[1],y_list[1],color="lawngreen",marker="o")
    print('最短路径长度为:',total_distance(cityset))

#蚂蚁类,有当前位置和移动禁忌表两个元素
class ant:
    def __init__(self):
        self.position = np.random.randint(low=0,high=CITY_LEN)
        self.Tabu_list = []

    def update_Tabu(self):
        if len(self.Tabu_list) >= CITY_LEN - 1:
            self.Tabu_list = []
        else:
            self.Tabu_list.append(self.position)

    def info(self):
        print('当前此蚂蚁所在位置为:',self.position+1)
        if len(self.Tabu_list)==0:
             print('当前此蚂蚁移动禁忌表为:',self.Tabu_list)
        else:
           print('当前此蚂蚁移动禁忌表为:',[i+1 for i in self.Tabu_list]) 

    def move(self,new_position):
        self.update_Tabu()
        self.position = new_position

#蚁群类,有蚁群蚂蚁数量和十个蚂蚁两个元素
class ants:
    def __init__(self,ant_number):
        self.number = ant_number
        self.antss = []
        for i in range(ant_number):
            self.antss.append(ant())

    def info(self):
        for i in range(self.number):
            print('蚂蚁',i+1)
            self.antss[i].info()

#路径类,有距离矩阵、信息素矩阵和比值矩阵三个元素
#距离矩阵一旦生成不再改变,信息素矩阵会更新
class route_matrix:
    def __init__(self,cityset):
        self.dist_matrix = dist_matrix(cityset)
        self.phe_matrix = phe_matrix()
        for i in range(CITY_LEN):
            self.dist_matrix[i][i] = 1
        self.index_matrix = self.phe_matrix/self.dist_matrix

    def update_index(self,pos1,pos2,num):
        if num == 2:
            up = 5.0
        elif num == 3:
            up = 1.0
        else:
            up = 3.0
        self.index_matrix[pos1][pos2] += up/self.dist_matrix[pos1][pos2]

    def info(self):
        #print('当前路径的距离矩阵为:')
        #print(self.dist_matrix)
        print('当前路径的信息素矩阵为:')
        print(self.phe_matrix)
        print('当前路径的比值矩阵为:')
        print(self.index_matrix)

#蚁群类,包含了蚂蚁和路程矩阵两个元素
class antcolony:
     def __init__(self,ants,route_matrix):
         self.ants = ants
         self.matrix = route_matrix
         self.best = []
         self.worst = []

     def info(self):
         self.ants.info()
         self.matrix.info()

     def move_prob(self):
         prob_lists = []
         for ant in self.ants.antss:
             temps = self.matrix.index_matrix[ant.position][:]
             prob_list = []
             for temp in temps:
                 prob_list.append(temp)
             for Tabu in ant.Tabu_list:
                 prob_list[Tabu] = 0
             new_list = [x/sum(prob_list) for x in prob_list]
             prob_lists.append(new_list)
         return prob_lists

     def update_matrix(self,pos1,pos2):
         if [pos1,pos2] in self.best or [pos2,pos1] in self.best:
             self.matrix.phe_matrix[pos1][pos2] += 5.0/self.matrix.dist_matrix[pos1][pos2]
             self.matrix.update_index(pos1,pos2,2)
         elif [pos1,pos2] in self.worst or [pos2,pos1] in self.worst:
             self.matrix.phe_matrix[pos1][pos2] += 1.0/self.matrix.dist_matrix[pos1][pos2]
             self.matrix.update_index(pos1,pos2,3)
         else:
            self.matrix.phe_matrix[pos1][pos2] += 3.0/self.matrix.dist_matrix[pos1][pos2]
            self.matrix.update_index(pos1,pos2,1)

     def wholemove(self):
         if len(self.ants.antss[0].Tabu_list) < 6:
             prob_lists = self.move_prob()
             turn_on = 0
         else:
             best_list = []
             turn_on = 1
         for i in range(ANT_NUMBER):
             if turn_on == 0:
                 idx = np.random.choice(np.arange(CITY_LEN),p=prob_lists[i])
                 #print(i,"的当前位置为:",self.ants.antss[i].position+1)
                 #print(i,'选择的下一位置为:',idx+1)
                 last_pos = self.ants.antss[i].position
                 self.ants.antss[i].move(idx)
                 self.update_matrix(last_pos,idx)
             else:
                 sums = sum(self.ants.antss[i].Tabu_list)
                 self.ants.antss[i].Tabu_list.append(21-sums)
                 best_list.append(self.ants.antss[i].Tabu_list)
                 last_pos = self.ants.antss[i].position
                 origin = self.ants.antss[i].Tabu_list[0]
                 self.ants.antss[i].move(origin)
                 self.update_matrix(last_pos,origin)
         if turn_on == 1:
             dist_list = [total_distance(pro_route(i)) for i in best_list]
             self.best = decompose(best_list[np.argmin(dist_list)])
             self.worst = decompose(best_list[np.argmax(dist_list)])
             print('当前最优为:',self.best)
             print('当前最不优为:',self.worst)


if __name__ == "__main__":
    np.random.shuffle(CHENGSHI)
    route = route_matrix(CHENGSHI)
    ants = ants(ANT_NUMBER)
    antcolony = antcolony(ants,route)
    for count in range(MOVE_TIMES):
        antcolony.wholemove()
    matrix = antcolony.matrix.index_matrix
    answer = find_route(matrix,CHENGSHI)
    antcolony.info()
    print_plot(answer)

手撕BP神经网络

参考知乎文章1:https://zhuanlan.zhihu.com/p/65472471

参考知乎文章2:https://zhuanlan.zhihu.com/p/59385110

算法实现

这不用说了吧——现代机器学习Machine Learning最好用且还在不断优化的一个热点算法。

BP神经网络的算法实现如下:

import pandas as pd
import numpy as np

INPUTS_LEN = 4

#激活函数,添加非线性成分
def sigmoid(x,method):
    if method == 1:
        return 1/(1+np.exp(-x))
    else:
        return [1/(1 + np.exp(-i)) for i in x]

#激活函数的导数
def deriv(x):
    fx = sigmoid(x,1)
    return fx*(1-fx)

#均方差损失函数
def mse_loss(y_tr, y_pr):
    return ((y_tr - y_pr) ** 2).mean()

#交叉熵损失函数
def log_loss(y_tr, y_pr):
    return (- np.multiply(y_tr,np.log(y_pr))
    - np.multiply(1-y_tr,np.log(1-y_pr))).mean()

#神经网络类
class neural_network():
    def __init__(self):
        self.w1 = np.random.normal(size=INPUTS_LEN*2).reshape(2,INPUTS_LEN)
        self.b1 = np.random.normal(size=2)
        self.w2 = np.random.normal(size=4).reshape(2,2)
        self.b2 = np.random.normal(size=2)
        self.w3 = np.random.normal(size=2)
        self.b3 = np.random.normal()

    def feedforward(self,inputs):
        [hi11,hi12] = np.dot(inputs,self.w1.T) + self.b1
        [ho11,ho12] = sigmoid([hi11,hi12],2)
        [hi21,hi22] = np.dot([ho11,ho12],self.w2.T) + self.b2
        [ho21,ho22] = sigmoid([hi21,hi22],2)
        sum_o = np.dot([ho21,ho22],self.w3.T) + self.b3
        output = 100 * sigmoid(sum_o,1)
        return output

    def train(self,input,answer):
        learm_rate = 0.9
        y_tr = answer
        [hi11,hi12] = np.dot(input,self.w1.T) + self.b1
        [ho11,ho12] = sigmoid([hi11,hi12],2)
        [hi1,ho1] = [[hi11,hi12],[ho11,ho12]]
        [hi21,hi22] = np.dot(ho1,self.w2.T) + self.b2
        [ho21,ho22] = sigmoid([hi21,hi22],2)
        [hi2,ho2] = [[hi21,hi22],[ho21,ho22]]
        sum_o = np.dot(ho2,self.w3.T) + self.b3
        y_pr = sigmoid(sum_o,1)
        print('y_pr为:',y_pr)

        #d_L_d_o = -y_tr/y_pr + (1-y_tr)/(1-y_pr)
        d_L_d_o = -2 * (y_tr - y_pr)

        ratio1 = learm_rate*d_L_d_o*deriv(sum_o)
        dwlist3 = [ratio1* i for i in ho2]

        ratio11 = learm_rate*ratio1*self.w3[0]*deriv(hi21)
        temp1 = [ratio11* i for i in ho1]
        ratio2 = learm_rate*ratio1*self.w3[1]*deriv(hi22)
        temp2 = [ratio2*i for i in ho1]
        dwlist2 = [temp1,temp2]

        ratio3 = learm_rate*(ratio11*self.w2[0][0]+ratio2*self.w2[1][0])*deriv(hi11)
        temp3 = [ratio3*i for i in input]
        ratio4 = learm_rate * (ratio11*self.w2[0][1]+ratio2*self.w2[1][1])*deriv(hi12)
        temp4 = [ratio4* i for i in input]
        dwlist1 = [temp3,temp4]

        self.w3 -= dwlist3
        self.w2 -= dwlist2
        self.w1 -= dwlist1
        self.b3 -= ratio1
        self.b2[0] -= ratio11
        self.b2[1] -= ratio2
        self.b1[0] -= ratio3
        self.b1[1] -= ratio4 

        d_l_d_input = []
        temp = self.w1.copy()
        temp2 = dwlist1.copy()
        data = temp*temp2*10
        for i in range(4):
            d_l_d_input.append((data[0][i]+data[1][i])/(9*input[i]))
        d_l_d_input = np.array(d_l_d_input)
        return d_l_d_input

    def show(self):
        print('w1:',self.w1)
        print('b1:',self.b1)
        print('w2:',self.w2)
        print('b2:',self.b2)
        print('w3:',self.w3)
        print('b3:',self.b3)

手撕CNN神经网络(卷积神经网络)

参考知乎文章1:https://zhuanlan.zhihu.com/p/102119808

参考知乎文章2:https://zhuanlan.zhihu.com/p/42559190

参考知乎文章3:https://www.zhihu.com/question/54504471/answer/332657604

参考知乎文章4:https://zhuanlan.zhihu.com/p/27908027

算法实现

在神经网络上进一步得到,主要用于处理图片像素矩阵。

卷积神经网络算法实现如下(需要导入上一个写的 BP Neural Network模块):

from PIL import Image
import numpy as np
from neural_network import *

#找矩阵中第二大的值
def find_second(np_matrix):
    temp_list = []
    for i in range(len(np_matrix)):
        for j in range(len(np_matrix[0])):
            temp_list.append(np_matrix[i][j])
    temp_list = np.array(temp_list)
    temp_list[np.argmax(temp_list)] = 0
    max2_index = np.argmax(temp_list)
    max_num = np.max(temp_list)
    return (max_num,max2_index)

#卷积神经网络类
class convolnetwork:
    def __init__(self,pixel):
        self.pixel = []
        self.conkernels = np.random.uniform(low = -1.0,high = 1.0,size =9).reshape(3,3)
        self.convol_maps = []
        self.poolings = []
        self.pool_nums = []
        self.softmax = neural_network()

    def show(self):
        print('卷积核如下:')
        print(self.conkernels)
        print('softmax层如下:')
        self.softmax.show()

    #数据处理成8*8的矩阵
    def press(self):
        pixel_matrixs = []
        for i in range(3):
            process = self.pixel[:,:,i]
            pixel_matrix = []
            for j in range(8):
                list_now = []
                for k in range(8):
                    choice = np.array(process[j*130:(j+1)*130,k*195:(k+1)*195])
                    list_now.append(choice.mean())
                pixel_matrix.append(list_now)
            pixel_matrixs.append(pixel_matrix)
        temp = np.array(pixel_matrixs)
        temp = temp / 255 * 2 - 1
        temp[0] = (temp[0] + temp[1] + temp[2])/3
        self.pixel = temp[0]

    #卷积,卷积核与前面处理的数据卷积运算
    def convolve(self):
        temp = []
        for j in range(6):
            list_now = []
            for k in range(6):
                piece = np.array(self.pixel[j:j+3,k:k+3].copy())
                mean_num = (piece*self.conkernels).mean()
                list_now.append(mean_num)
            self.convol_maps.append(list_now)
        self.convol_maps = np.array(self.convol_maps)

    #池化,各个区域取最大值
    def pooling(self):
        #self.last_input = 
        list_now = []
        for i in range(2):
            for j in range(2):
                num = []
                temp = np.array(self.convol_maps[i*3:(i+1)*3,j*3:(j+1)*3].copy())
                max_index = np.argmax(self.convol_maps[i*3:(i+1)*3,j*3:(j+1)*3].copy())
                column_index = max_index%3 + j*3 + 1
                row_index = int((max_index - max_index%3)/3 + i*3) + 1
                max2_num,max2_index = find_second(temp)
                column2_index = max2_index%3 + j*3 + 1
                row2_index = int((max2_index - max2_index%3)/3 + i*3) + 1
                if row_index == -1:
                    row_index += 1
                if row2_index == -1:
                    row2_index += 1
                num.extend([row_index,column_index])
                num.extend([row2_index,column2_index])
                list_now.append((temp.max()+max2_num)/2.0)
                self.pool_nums.append(num)
        self.poolings = np.array(list_now)

    #BP反向传播算法
    def back_propagation(self,answer):
        d_l_d_inputs = self.softmax.train(self.poolings,answer)
        for i in range(4):
            temp = self.pool_nums[i]
            row = temp[0]
            column = temp[1]
            row2 = temp[2]
            column2 = temp[3]
            piece = self.pixel[row-1:row+2,column-1:column+2].copy()
            piece2 = self.pixel[row2-1:row2+2,column2-1:column2+2].copy()
            self.conkernels -= 9.0/2.0*(piece*d_l_d_inputs[i]+piece2*d_l_d_inputs[i]) 

    #采用BP+SGD的方法来训练神经网络
    def train(self,pixel,answer):
        self.pixel = pixel
        self.convol_maps = []
        self.poolings = []
        self.pool_nums = []
        self.press()
        self.convolve()
        self.pooling()
        self.back_propagation(answer)


img = Image.open('xjt.jpg')
pixels = np.array(img)
convol1 = convolnetwork(pixels) 
convol1.show()
convol1.train(pixels,1)
convol1.show()

手撕Dijkstra最短路算法

代码实现

Dijkstra最短路算法是搜索两点间最短路径的一种算法。

Dijkstra最短路算法实现如下:

import numpy as np

NODES_NUM = 9
DISTANCE_MATRIX = np.array([[0,4,0,0,0,0,0,8,0],
                            [4,0,8,0,0,0,0,3,0],
                            [0,8,0,7,0,4,0,0,2],
                            [0,0,7,0,9,14,0,0,0],
                            [0,0,0,9,0,10,0,0,0],
                            [0,0,4,14,10,0,2,0,0],
                            [0,0,0,0,0,2,0,6,6],
                            [8,3,0,0,0,0,6,0,1],
                            [0,0,2,0,0,0,6,1,0]])

VISITED = np.ones(9)
DISTANCE_TOTAL = np.array([0,100,100,100,100,100,100,100,100])

PARENT = np.array([-1,-1,-1,-1,-1,-1,-1,-1,-1])

#先选取一个起点
parent_now = 0
#实现连接点路径覆盖并在未选过的点集合中选取距离最小的点作为新的起点
#直至遍历所有点(所有点都做过起点)
while (sum(VISITED) != 0):
    VISITED[parent_now] = 0
    source = DISTANCE_TOTAL[parent_now]
    for i in range(NODES_NUM):
        if DISTANCE_MATRIX[parent_now][i] != 0 and VISITED[i] != 0:
            if ((source + DISTANCE_MATRIX[parent_now][i]) < DISTANCE_TOTAL[i]):
                DISTANCE_TOTAL[i] = source + DISTANCE_MATRIX[parent_now][i]
                PARENT[i] = parent_now
            else:
                pass
    temp = []
    sequence = []
    for i in range(NODES_NUM):
        if VISITED[i] != 0:
            temp.append(DISTANCE_TOTAL[i])
            sequence.append(i)
    if len(temp) != 0:
        index = np.argmin(temp)
        parent_now = sequence[index]
    print(DISTANCE_TOTAL)

print('所找到的最短路径长度为:',np.max(DISTANCE_TOTAL))

手撕Floyd-Warshall最短路算法

代码实现

Floyd-Warshall最短路算法是数模中应用极广的一种算法,可求出任意两点间最短路径距离,很好编写代码。

Floyd-Warshall最短路算法实现如下:

import numpy as np

NODES_NUM = 9
DISTANCE_MATRIX = np.array([[60,4,60,60,60,60,60,8,60],
                            [4,60,8,60,60,60,60,3,60],
                            [60,8,60,7,60,4,60,60,2],
                            [60,60,7,60,9,14,60,60,60],
                            [60,60,60,9,60,10,60,60,60],
                            [60,60,4,14,10,60,2,60,60],
                            [60,60,60,60,60,2,60,6,6],
                            [8,3,60,60,60,60,6,60,1],
                            [60,60,2,60,60,60,6,1,60]])

#暴力遍历,根据最短路的子路径也为最短路的法则
for k in range(NODES_NUM):
    for i in range(NODES_NUM):
        for j in range(NODES_NUM):
            DISTANCE_MATRIX[i][j] = min(DISTANCE_MATRIX[i][j],DISTANCE_MATRIX[i][k]+DISTANCE_MATRIX[k][j])

#将所有点到自身的距离修改为0
for i in range(NODES_NUM):
    DISTANCE_MATRIX[i][i] = 0
print(DISTANCE_MATRIX)