启发式算法为克服优化过程中出现的局部最优解,因为在非凸优化中,往往会陷入局部最优。
1、传统启发式
1.1 贪心算法
贪心算法的核心是总做出当前状态下看起来最好的决策,得到最优解需要问题具有无后效性,即当前的解对之后的解没有影响,一般问题都达不到这个条件,所以需要制定贪心策略,最终得到最优解近似。以下记录几个基本问题:
1. 背包问题
# 背包问题
if __name__ == '__main__':
beg = 54 # 背包54kg
value = 0 # 已经获得的价值
choice = []
while beg > 0: # 如果背包还有空位,则递归
if beg >= 8: # 选择当前这一步的最优解,既选择B商品
beg = beg - 8
value = value + 13
choice.append("B")
elif beg >= 10: # 要是B商品选择不了,则选择第二单位价值的物品,即C物品
beg = beg - 10
value = value + 15
choice.append("C")
elif beg >= 6:
beg = beg - 6
value = value + 8
choice.append("A")
else: # 当所有的物品都选择不了,则退出
break
print("剩余的背包重量:", beg)
print("获得的总价值:", value)
print("选择的物品的类型及顺序:", choice)
剩余的背包重量: 0
获得的总价值: 86
选择的物品的类型及顺序: ['B', 'B', 'B', 'B', 'B', 'B', 'A']
里面的数据可以自行调整。
2.分糖果问题
# 分糖果
# 孩子的需求因子更小则其更容易被满足,故优先从需求因子小的孩子尝试,可以得到正确的结果
# 因为我们追求更多的孩子被满足,所以用一个糖果满足需求因子较小或较大的孩子都是一样的)
# g 需求因子 s 糖果大小数组
class Solution:
def findContentChild(self, g, s):
g = sorted(g)
s = sorted(s)
child = 0
cookie = 0
while child < len(g) and cookie < len(s):
if g[child] <= s[cookie]:
child += 1
cookie += 1
return child
if __name__ == "__main__":
g = [5, 10, 2, 9, 15, 9]
s = [6, 1, 20, 3, 8]
S = Solution()
result = S.findContentChild(g, s)
print(result)
结果:3
可以满足三个孩子。
3.摇摆排序
# 摇摆序列
# 当序列有一段连续递增(或递减)时,为形成摇摆子序列,我们只需要保留这段连续递增(或递减)的首尾元素
# 这样更有可能使得尾部的后一个元素成为摇摆子序列的下一个元素
class Solution:
def maxlength(self, nums, newnums):
if len(nums) < 2:
return len(nums)
BEGIN = 0
UP = 1
DOWN = 2
STATE = BEGIN
max_length = 1
vision = [UP, BEGIN, DOWN]
for i in range(0, len(nums)-1): # i从0开始
if STATE == 0: # 在开始状态
if nums[i] < nums[i+1]: # 不论上升还是下降,长度都会加1,不同的是当前状态
STATE = 1
max_length += 1
newnums.append(nums[i])
elif nums[i] > nums[i+1]:
STATE = 2
max_length += 1
newnums.append(nums[i])
if STATE == 1: # 上升状态
if nums[i] > nums[i+1]:
STATE = 2
max_length += 1
newnums.append(nums[i])
if STATE == 2: # 下降状态
if nums[i] < nums[i+1]:
STATE = 1
max_length += 1
newnums.append(nums[i])
if len(newnums) < max_length:
newnums.append(nums[len(nums) - 1])
return max_length
if __name__ == "__main__":
S = Solution()
g = [1, 17, 5, 10, 13, 15, 10, 5, 16, 8]
newnums = []
result = S.maxlength(g, newnums)
print(result)
print("摇摆序列为:", newnums)
7
摇摆序列为: [1, 17, 5, 15, 5, 16, 8]
4.移除K个数字
# 移除K个数字
class Solution:
def removeknums(self, nums, k): # nums类型为str
s = []
nums = list(map(int, nums)) # 将str转化为list
for i in range(len(nums)):
number = int(nums[i])
while len(s) != 0 and s[len(s) - 1] > number and k > 0:
s.pop(-1) # pop移除列表元素 默认移除最后一个元素
k -= 1
if number != 0 or len(s) != 0:
s.append(number) # append增加列表元素
while len(s) != 0 and k > 0:
s.pop(-1)
k -= 1
result = ''.join(str(i) for i in s) # join方法合并字符串
return result
if __name__ == "__main__":
S = Solution()
print(S.removeknums("7612397", 3))
结果:1237
https://blog.csdn.net/SweetSeven_/article/details/95197131
1.2 局部搜索
通常考察一个算法的性能通常用局部搜索能力和全局收敛能力这两个指标。局部搜索是指能够无穷接近最优解的能力,而全局收敛能力是指找到全局最优解所在大致位置的能力。
局部搜索算法是在一组可行解的基础上,在当前解的邻域内进行局部搜索产生新的可行解的过程。局部搜索中关键概念是邻域以及邻域动作。邻域是给定点附近其他点的集合,在距离空间中,邻域一般被定义为以给定点为圆心的一个圆;而在组合优化问题中,邻域一般定义为由给定转化规则对给定的问题域上每结点进行转化所得到的问题域上节点的集合。局部搜索算法的基本思想:在搜索过程中,始终选择当前点的邻居中与离目标最近者的方向搜索。
局部搜索有三大问题,一是容易陷入局部极值点而结束,二是步长的选择,三是起始点的选择,都对最终的优化结果影响很大。解决以上问题的方法如下,每次并不一定选择邻域内最优的点,而是依据一定的概率,从邻域内选择一个点;可变步长;随机生成一些初始点,从每个初始点出发进行搜索,找到各自的最优解,再从这些最优解中选择一个最好的结果作为最终的结果。
1.3 爬山算法
# y = sin(x^2) + 2*cos(2*x)
import numpy as np
import matplotlib.pyplot as plt
import math
xpoints = np.arange(5, 8, 0.001)
ypoints = np.sin(xpoints * xpoints) + 2*np.cos(2*xpoints)
plt.plot(xpoints, ypoints)
plt.show()
# 搜索步长
DELTA = 0.01
# 定义域x从5到8闭区间
BOUND = [5, 8]
# 随机取乱数100次
GENERATION = 100
def F(x):
return math.sin(x * x) + 2.0 * math.cos(2.0 * x)
def hillClimbing(x):
while F(x + DELTA) > F(x) and x + DELTA <= BOUND[1] and x + DELTA >= BOUND[0]:
x = x + DELTA
while F(x - DELTA) > F(x) and x - DELTA <= BOUND[1] and x - DELTA >= BOUND[0]:
x = x - DELTA
return x, F(x)
def findMax():
highest = [0, -1000]
for i in range(GENERATION):
x = np.random.rand() * (BOUND[1] - BOUND[0]) + BOUND[0] # 区间[5,8]
currentValue = hillClimbing(x)
print('current value is :', currentValue)
if currentValue[1] > highest[1]:
highest[:] = currentValue
return highest
[x, y] = findMax()
print('highest point is x :{},y:{}'.format(x, y))
python实现爬山算法 - 腾讯云开发者社区-腾讯云 (tencent.com)
2、元启发式
2.1 禁忌搜索算法
将最优解放入禁忌表中,更新禁忌条件,可以避免重复搜索已有的最优解,从而跳出局部最优解。
2.2 模拟退火算法
求解下列函数最优值
# 模拟退火
import math # 数学运算
from random import random # 随机数
import matplotlib.pyplot as plt # 画图
def func(x, y): # z值计算公式,传递入x,y,返回res
res = 4 * x ** 2 - 2.1 * x ** 4 + x ** 6 / 3 + x * y - 4 * y ** 2 + 4 * y ** 4
return res
# __init__方法实例化时会自动调用,可以有参数,通过__init__传递到类的实例化操作上
# self代表类的实例,而非类
# x为公式里的x1,y为公式里面的x2
class SA:
def __init__(self, func, iter=1000, T0=100, Tf=0.01, alpha=0.99): # 方程 迭代次数 温度上限 温度下限 学习率
self.func = func # 往类中传递进函数
self.iter = iter # 内循环迭代次数,即为L =100
self.alpha = alpha # 降温系数,alpha=0.99
self.T0 = T0 # 初始温度T0为100
self.Tf = Tf # 温度终值Tf为0.01
self.T = T0 # 当前温度为T0
self.x = [random() * 12 - 6 for i in range(iter)] # 随机生成100个x的值
self.y = [random() * 12 - 6 for i in range(iter)] # 随机生成100个y的值
self.most_best = []
self.history = {'f': [], 'T': []}
def generate_new(self, x, y): # 扰动产生新解的过程
while True:
x_new = x + self.T * (random() - random()) # self.T
y_new = y + self.T * (random() - random()) # 一个类不存在继承,继承只在类和类之间
if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
break # 重复得到新解,直到产生的新解满足约束条件
return x_new, y_new
def Metrospolis(self, f, f_new): # Metropolis准则
if f_new <= f:
return 1
else:
p = math.exp((f - f_new) / self.T) # 接收新解的概率
if random() < p:
return 1
else:
return 0
def best(self): # 获取最优目标函数值
f_list = [] # f_list数组保存每次迭代之后的值
for i in range(self.iter):
f = self.func(self.x[i], self.y[i])
f_list.append(f)
f_best = min(f_list) # 在一组迭代中选择最小值
idx = f_list.index(f_best) # 找到索引
return f_best, idx # f_best,idx分别为在该温度下,迭代L次之后目标函数的最优解和最优解的下标
def run(self): # 2.运行此函数,__init__自动运行,得到初值
count = 0 # 外循环迭代计数
# 外循环迭代直到不满足while条件
# 直到当前温度小于等于终止温度的阈值
while self.T > self.Tf:
# 内循环迭代100次(每次内循环在同一温度下)
for i in range(self.iter):
f = self.func(self.x[i], self.y[i]) # f为迭代一次后的值 x y为随机生成的100个数的数组
x_new, y_new = self.generate_new(self.x[i], self.y[i]) # 使用函数generate_new产生新解
f_new = self.func(x_new, y_new) # 产生新值
if self.Metrospolis(f, f_new): # 使用函数Metrospolis判断是否接受新值
self.x[i] = x_new # 如果接受新值,则把新值的x,y存入x数组和y数组替代原有的数
self.y[i] = y_new
# 迭代L次记录在该温度下最优解
ft, _ = self.best() # return f_best, idx 返回最优解和最优解下标
self.history['f'].append(ft) # 最优解
self.history['T'].append(self.T) # 此时的温度 绘图
# 温度按照一定的比例下降(冷却)
self.T = self.T * self.alpha # 影响产生新解(generate_new)、接受新解的概率(metrospolis)
count += 1
# 得到最优解
f_best, idx = self.best()
print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")
sa = SA(func)
sa.run() # 1.将func传递进类SA运行类中的run(self)函数
plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis() # 翻转
plt.show()
关于代码:认识到__init__,self,return的用法。
关于模拟退火算法
外层循环:退火过程,降温,温度值下降迭代
内层循环:每次内层循环时温度相同,可以控制迭代次数,如100次,首先生成的x,y数值是一组100维度的数组,产生x_new和y_new,在100次迭代过程中寻找最优值,每次迭代得到的目标函数优于之前的,将有一定概率替代原有数值,也有极小的概率不替代,这就是模拟退火的关键所在,极小的不替代概率可以帮助跳过局部最优。
ps.产生新解和接受新解均受温度的影响,随温度下降,接受使目标函数上升的概率逐渐增大(求最大值),和退火过程相似,当温度逐渐下降,分子运动的活跃度也将下降,变化的可能性就会降低。
https://blog.csdn.net/BubbleCodes/article/details/119492432
https://blog.csdn.net/weixin_48241292/article/details/109468947
2.3 遗传算法
迭代次数可以提前给定,也能根据上下两次迭代适应度的差值判断是否收敛。
# 画图
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
fig = plt.figure() #定义新的三维坐标轴
ax = plt.axes(projection='3d')
x = np.arange(-3, 3, 0.01)
y = np.arange(-3, 3, 0.01)
x, y = np.meshgrid(x, y)
z = 3*(1-x)**2*np.exp(-(x**2)-(y+1)**2)- 10*(x/5 - x**3 - y**5)*np.exp(-x**2-y**2)- 1/3**np.exp(-(x+1)**2 - y**2)
# 背景空白
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
#作图
ax.plot_surface(x,y,z,rstride = 1, cstride = 1, cmap='rainbow')
plt.savefig('plot123_2.png', dpi=300)
plt.show()
# 遗传算法
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
DNA_SIZE = 24 # DNA编码长度
POP_SIZE = 200 # 种群规模
CROSSOVER_RATE = 0.8 # 交叉率
MUTATION_RATE = 0.005 # 变异率
N_GENERATIONS = 50 # 最大迭代次数
X_BOUND = [-3, 3]
Y_BOUND = [-3, 3]
# 迭代50次,每次都在200个个体的种群中迭代
def F(x, y):
return 3*(1-x)**2*np.exp(-(x**2)-(y+1)**2)- 10*(x/5 - x**3 - y**5)*np.exp(-x**2-y**2)- 1/3**np.exp(-(x+1)**2 - y**2)
def plot_3d(ax): # 三维图
X = np.linspace(*X_BOUND, 100)
Y = np.linspace(*Y_BOUND, 100)
X,Y = np.meshgrid(X, Y)
Z = F(X, Y)
ax.plot_surface(X,Y,Z,rstride=1,cstride=1,cmap=cm.coolwarm)
ax.set_zlim(-10,10)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.pause(3)
plt.show()
def get_fitness(pop):
x,y = translateDNA(pop)
pred = F(x, y)
return (pred - np.min(pred)) + 1e-3 # 减去最小的适应度是为了防止适应度出现负数,通过这一步fitness的范围为[0, np.max(pred)-np.min(pred)],最后在加上一个很小的数防止出现为0的适应度
def translateDNA(pop): # pop表示种群矩阵,一行表示一个二进制编码表示的DNA,矩阵的行数为种群数目
x_pop = pop[:,1::2] # 奇数列表示X
y_pop = pop[:,::2] # 偶数列表示y
# pop:(POP_SIZE,DNA_SIZE)*(DNA_SIZE,1) --> (POP_SIZE,1)
x = x_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(X_BOUND[1]-X_BOUND[0])+X_BOUND[0]
y = y_pop.dot(2**np.arange(DNA_SIZE)[::-1])/float(2**DNA_SIZE-1)*(Y_BOUND[1]-Y_BOUND[0])+Y_BOUND[0]
return x,y
def crossover_and_mutation(pop, CROSSOVER_RATE = 0.8):
new_pop = []
for father in pop: # 遍历种群中的每一个个体,将该个体作为父亲
child = father # 孩子先得到父亲的全部基因(这里我把一串二进制串的那些0,1称为基因)
if np.random.rand() < CROSSOVER_RATE: # 产生子代时不是必然发生交叉,而是以一定的概率发生交叉
mother = pop[np.random.randint(POP_SIZE)] # 再种群中选择另一个个体,并将该个体作为母亲
cross_points = np.random.randint(low=0, high=DNA_SIZE*2) # 随机产生交叉的点
child[cross_points:] = mother[cross_points:] # 孩子得到位于交叉点后的母亲的基因
mutation(child) # 每个后代有一定的机率发生变异
new_pop.append(child)
return new_pop
def mutation(child, MUTATION_RATE=0.003):
if np.random.rand() < MUTATION_RATE: # 以MUTATION_RATE的概率进行变异
mutate_point = np.random.randint(0, DNA_SIZE*2) # 随机产生一个实数,代表要变异基因的位置
child[mutate_point] = child[mutate_point]^1 # 将变异点的二进制为反转
def select(pop, fitness): # nature selection wrt pop's fitness
idx = np.random.choice(np.arange(POP_SIZE), size=POP_SIZE, replace=True,
p=(fitness)/(fitness.sum()) )
return pop[idx]
# 主要是使用了choice里的最后一个参数p
# 参数p描述了从np.arange(POP_SIZE)里选择每一个元素的概率,
# 概率越高约有可能被选中,最后返回被选中的个体即可
def print_info(pop):
fitness = get_fitness(pop)
max_fitness_index = np.argmax(fitness)
print("max_fitness:", fitness[max_fitness_index])
x,y = translateDNA(pop)
print("最优的基因型:", pop[max_fitness_index])
print("(x, y):", (x[max_fitness_index], y[max_fitness_index]))
if __name__ == "__main__":
fig = plt.figure()
ax = Axes3D(fig)
plt.ion() # 将画图模式改为交互模式,程序遇到plt.show不会暂停,而是继续执行
plot_3d(ax) # def plot_3d
pop = np.random.randint(2, size=(POP_SIZE, DNA_SIZE*2))
# matrix(POP_SIZE, DNA_SIZE)
# 种群矩阵
for _ in range(N_GENERATIONS): # 迭代N代
x,y = translateDNA(pop)
if 'sca' in locals():
sca.remove()
sca = ax.scatter(x, y, F(x,y), c='black', marker='o');plt.show();plt.pause(0.1)
pop = np.array(crossover_and_mutation(pop, CROSSOVER_RATE))
# 交叉变异 返回new_pop
# F_values = F(translateDNA(pop)[0], translateDNA(pop)[1])#x, y --> Z matrix
fitness = get_fitness(pop)
# 计算适应度
pop = select(pop, fitness)
# 选择生成新的种群
print_info(pop)
plt.ioff()
plot_3d(ax)
遗传算法详解 附python代码实现_重学CS的博客-CSDN博客_python遗传算法