- 目标函数或约束条件包含非线性函数，则为非线性规划问题
- 若线性规划的最优解存在，则其最优解只能在其可行域的边界上达到，而非线性规划的最优解则可能在其可行域的任意一点达到
- 例1.
import numpy as np
from scipy.optimize import minimize
# 目标函数
def objective(x):
return x[0] ** 2 + x[1]**2 + x[2]**2 +8
# 约束条件
def constraint1(x):
return x[0] ** 2 - x[1] + x[2]**2 # 不等约束
def constraint2(x):
return -(x[0] + x[1]**2 + x[2]**2-20) # 不等约束
def constraint3(x):
return -x[0] - x[1]**2 + 2
def constraint4(x):
return x[1] + 2*x[2]**2 -3 # 不等约束
# 边界约束
b = (0.0, None)
bnds = (b, b ,b)
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
con3 = {'type': 'eq', 'fun': constraint3}
con4 = {'type': 'eq', 'fun': constraint4}
cons = ([con1, con2, con3,con4]) # 3个约束条件
# 计算
x0=np.array([0,0,0])#定义初始值
solution = minimize(objective, x0, method='SLSQP', \
bounds=bnds, constraints=cons)
x = solution.x
print('目标值: ' + str(objective(x)))
print('答案为')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))- 例2.梯度下降法
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
def Fun(x, y): # 原函数
return x - y + 2 * x * x + 2 * x * y + y * y # 方程为z=x1^2+2x2^2-4x1-2x1x2
def PxFun(x, y): # x偏导
return 1 + 4 * x + 2 * y
def PyFun(x, y): # y偏导
return -1 + 2 * x + 2 * y
# 初始化
fig = plt.figure() # figure对象
ax = Axes3D(fig) # Axes3D对象
X, Y = np.mgrid[-2:2:40j, -2:2:40j] # 取样并作满射联合
Z = Fun(X, Y) # 取样点Z坐标打表
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap="rainbow")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# 梯度下降
step = 0.0008 # 下降系数
x = 0
y = 0 # 初始选取一个点
tag_x = [x]
tag_y = [y]
tag_z = [Fun(x, y)] # 三个坐标分别打入表中，该表用于绘制点
new_x = x
new_y = y
Over = False
while Over == False:
new_x -= step * PxFun(x, y)
new_y -= step * PyFun(x, y) # 分别作梯度下降
if Fun(x, y) - Fun(new_x, new_y) < 7e-9: # 精度
Over = True
x = new_x
y = new_y # 更新旧点
tag_x.append(x)
tag_y.append(y)
tag_z.append(Fun(x, y)) # 新点三个坐标打入表中
# 绘制点/输出坐标
ax.plot(tag_x, tag_y, tag_z, 'r.')
plt.title('(x,y)~(' + str(x) + "," + str(y) + ')')
plt.show()- 例3.粒子群法
import math
import random
import numpy as np
import matplotlib.pyplot as plt
class PSO:
def __init__(self, dimension, time, size, low, up, v_low, v_high):
# 初始化
self.dimension = dimension # 变量个数
self.time = time # 迭代的代数
self.size = size # 种群大小
self.bound = [] # 变量的约束范围
self.bound.append(low)
self.bound.append(up)
self.v_low = v_low
self.v_high = v_high
self.x = np.zeros((self.size, self.dimension)) # 所有粒子的位置
self.v = np.zeros((self.size, self.dimension)) # 所有粒子的速度
self.p_best = np.zeros((self.size, self.dimension)) # 每个粒子最优的位置
self.g_best = np.zeros((1, self.dimension))[0] # 全局最优的位置
# 初始化第0代初始全局最优解
temp = -1000000
for i in range(self.size):
for j in range(self.dimension):
self.x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j])
self.v[i][j] = random.uniform(self.v_low, self.v_high)
self.p_best[i] = self.x[i] # 储存最优的个体
fit = self.fitness(self.p_best[i])
# 做出修改
if fit > temp:
self.g_best = self.p_best[i]
temp = fit
def fitness(self, x):
"""
个体适应值计算
"""
x1 = x[0]
x2 = x[1]
y = math.floor(x1**2+2*x2*x2-4*x1-2*x1*x2)
# print(y)
return y
def update(self, size):
c1 = 2.0 # 学习因子
c2 = 2.0
w = 0.8 # 自身权重因子
for i in range(size):
# 更新速度(核心公式)
self.v[i] = w * self.v[i] + c1 * random.uniform(0, 1) * (
self.p_best[i] - self.x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.x[i])
# 速度限制
for j in range(self.dimension):
if self.v[i][j] < self.v_low:
self.v[i][j] = self.v_low
if self.v[i][j] > self.v_high:
self.v[i][j] = self.v_high
# 更新位置
self.x[i] = self.x[i] + self.v[i]
# 位置限制
for j in range(self.dimension):
if self.x[i][j] < self.bound[0][j]:
self.x[i][j] = self.bound[0][j]
if self.x[i][j] > self.bound[1][j]:
self.x[i][j] = self.bound[1][j]
# 更新p_best和g_best
if self.fitness(self.x[i]) > self.fitness(self.p_best[i]):
self.p_best[i] = self.x[i]
if self.fitness(self.x[i]) > self.fitness(self.g_best):
self.g_best = self.x[i]
def pso(self):
best = []
self.final_best = np.array([1, 2, 3, 4, 5])
for gen in range(self.time):
self.update(self.size)
if self.fitness(self.g_best) > self.fitness(self.final_best):
self.final_best = self.g_best.copy()
print('当前最佳位置：{}'.format(self.final_best))
temp = self.fitness(self.final_best)
print('当前的最佳适应度：{}'.format(temp))
best.append(temp)
t = [i for i in range(self.time)]
plt.figure()
plt.plot(t, best, color='red', marker='.', ms=15)
plt.rcParams['axes.unicode_minus'] = False
plt.margins(0)
plt.xlabel(u"迭代次数") # X轴标签
plt.ylabel(u"适应度") # Y轴标签
plt.title(u"迭代过程") # 标题
plt.show()
if __name__ == '__main__':
time = 50
size = 100
dimension = 5
v_low = -1
v_high = 1
low = [1, 1, 1, 1, 1]
up = [25, 25, 25, 25, 25]
pso = PSO(dimension, time, size, low, up, v_low, v_high)
pso.pso()- 例4.遗传算法
import math, random
class Population:
# 种群的设计
def __init__(self, size, chrom_size, cp, mp, gen_max):
# 种群信息合
self.individuals = [] # 个体集合
self.fitness = [] # 个体适应度集
self.selector_probability = [] # 个体选择概率集合
self.new_individuals = [] # 新一代个体集合
self.elitist = {'chromosome': [0, 0], 'fitness': 0, 'age': 0} # 最佳个体的信息
self.size = size # 种群所包含的个体数
self.chromosome_size = chrom_size # 个体的染色体长度
self.crossover_probability = cp # 个体之间的交叉概率
self.mutation_probability = mp # 个体之间的变异概率
self.generation_max = gen_max # 种群进化的最大世代数
self.age = 0 # 种群当前所处世代
# 随机产生初始个体集，并将新一代个体、适应度、选择概率等集合以 0 值进行初始化
v = 2 ** self.chromosome_size - 1
for i in range(self.size):
self.individuals.append([random.randint(0, v), random.randint(0, v)])
self.new_individuals.append([0, 0])
self.fitness.append(0)
self.selector_probability.append(0)
# 基于轮盘赌博机的选择
def decode(self, interval, chromosome):
'''将一个染色体 chromosome 映射为区间 interval 之内的数值'''
d = interval[1] - interval[0]
n = float(2 ** self.chromosome_size - 1)
return (interval[0] + chromosome * d / n)
def fitness_func(self, chrom1, chrom2):
'''适应度函数，可以根据个体的两个染色体计算出该个体的适应度'''
interval = [-10.0, 10.0]
(x, y) = (self.decode(interval, chrom1),
self.decode(interval, chrom2))
n = lambda x, y: math.sin(math.sqrt(x * x + y * y)) ** 2 - 0.5
d = lambda x, y: (1 + 0.001 * (x * x + y * y)) ** 2
func = lambda x, y: 0.5 - n(x, y) / d(x, y)
return func(x, y)
def evaluate(self):
'''用于评估种群中的个体集合 self.individuals 中各个个体的适应度'''
sp = self.selector_probability
for i in range(self.size):
self.fitness[i] = self.fitness_func(self.individuals[i][0], # 将计算结果保存在 self.fitness 列表中
self.individuals[i][1])
ft_sum = sum(self.fitness)
for i in range(self.size):
sp[i] = self.fitness[i] / float(ft_sum) # 得到各个个体的生存概率
for i in range(1, self.size):
sp[i] = sp[i] + sp[i - 1] # 需要将个体的生存概率进行叠加，从而计算出各个个体的选择概率
# 轮盘赌博机（选择）
def select(self):
(t, i) = (random.random(), 0)
for p in self.selector_probability:
if p > t:
break
i = i + 1
return i
# 交叉
def cross(self, chrom1, chrom2):
p = random.random() # 随机概率
n = 2 ** self.chromosome_size - 1
if chrom1 != chrom2 and p < self.crossover_probability:
t = random.randint(1, self.chromosome_size - 1) # 随机选择一点（单点交叉）
mask = n << t # << 左移运算符
(r1, r2) = (chrom1 & mask, chrom2 & mask) # & 按位与运算符：参与运算的两个值,如果两个相应位都为1,则该位的结果为1,否则为0
mask = n >> (self.chromosome_size - t)
(l1, l2) = (chrom1 & mask, chrom2 & mask)
(chrom1, chrom2) = (r1 + l2, r2 + l1)
return (chrom1, chrom2)
# 变异
def mutate(self, chrom):
p = random.random()
if p < self.mutation_probability:
t = random.randint(1, self.chromosome_size)
mask1 = 1 << (t - 1)
mask2 = chrom & mask1
if mask2 > 0:
chrom = chrom & (~mask2) # ~ 按位取反运算符：对数据的每个二进制位取反,即把1变为0,把0变为1
else:
chrom = chrom ^ mask1 # ^ 按位异或运算符：当两对应的二进位相异时，结果为1
return chrom
# 保留最佳个体
def reproduct_elitist(self):
# 与当前种群进行适应度比较，更新最佳个体
j = -1
for i in range(self.size):
if self.elitist['fitness'] < self.fitness[i]:
j = i
self.elitist['fitness'] = self.fitness[i]
if (j >= 0):
self.elitist['chromosome'][0] = self.individuals[j][0]
self.elitist['chromosome'][1] = self.individuals[j][1]
self.elitist['age'] = self.age
# 进化过程
def evolve(self):
indvs = self.individuals
new_indvs = self.new_individuals
# 计算适应度及选择概率
self.evaluate()
# 进化操作
i = 0
while True:
# 选择两个个体，进行交叉与变异，产生新的种群
idv1 = self.select()
idv2 = self.select()
# 交叉
(idv1_x, idv1_y) = (indvs[idv1][0], indvs[idv1][1])
(idv2_x, idv2_y) = (indvs[idv2][0], indvs[idv2][1])
(idv1_x, idv2_x) = self.cross(idv1_x, idv2_x)
(idv1_y, idv2_y) = self.cross(idv1_y, idv2_y)
# 变异
(idv1_x, idv1_y) = (self.mutate(idv1_x), self.mutate(idv1_y))
(idv2_x, idv2_y) = (self.mutate(idv2_x), self.mutate(idv2_y))
(new_indvs[i][0], new_indvs[i][1]) = (idv1_x, idv1_y) # 将计算结果保存于新的个体集合self.new_individuals中
(new_indvs[i + 1][0], new_indvs[i + 1][1]) = (idv2_x, idv2_y)
# 判断进化过程是否结束
i = i + 2 # 循环self.size/2次，每次从self.individuals 中选出2个
if i >= self.size:
break
# 最佳个体保留
# 如果在选择之前保留当前最佳个体，最终能收敛到全局最优解。
self.reproduct_elitist()
# 更新换代：用种群进化生成的新个体集合 self.new_individuals 替换当前个体集合
for i in range(self.size):
self.individuals[i][0] = self.new_individuals[i][0]
self.individuals[i][1] = self.new_individuals[i][1]
def run(self):
'''根据种群最大进化世代数设定了一个循环。
在循环过程中，调用 evolve 函数进行种群进化计算，并输出种群的每一代的个体适应度最大值、平均值和最小值。'''
for i in range(self.generation_max):
self.evolve()
print(i, max(self.fitness), sum(self.fitness) / self.size, min(self.fitness))
if __name__ == '__main__':
# 种群的个体数量为 50，染色体长度为 25，交叉概率为 0.8，变异概率为 0.1,进化最大世代数为 150
pop = Population(50, 24, 0.8, 0.1, 150)
pop.run()
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
