代码收藏家 利用scikit中的遗传算法求解(整数01)约束规划实例详解教程+利用scipy.optimize求解约束规划问题
注意标准形式
下面两个方法约束规划的一般标准形式为：
利用scikit-opt的遗传算法求解约束规划问题
先放上链接：scikit-opt网址
主要四个步骤：
下面依照此题多约束为例
可知该题有5个不等式约束，且决策变量为01整数，后面将具体讲解如何将目标函数的约束条件加入GA模型中
一：import scikit-opt库
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sko.GA import GA ###########此行为调出scikit-opt中的遗传算法
###########
#若安装
#则： pip install scikit-opt
#####################################
二：定义目标函数
看线面的代码我们发现里面有个参数p，这个参数我理解就是GA中的某一个种群，真正调用这个函数的时候我们不需要写参数和括号
def schaffer(p):
'''
目标函数: min(xi)求和
决策变量维度：n
决策变量类型：0 1 整数
:param p:
:return:
'''
x=np.array(p)###转化为数组
return sum(x)###返回将目标函数的值
三：定义约束
注意这里遗传算法中的约束应该默认为≤0的形式，而且不管有多少个约束，最终将约束的值及你给过一个max函数累加后返回，同样这里也需要决策变量，这里也是参数P，我们真正调用时，不需要写括号和参数
def constraint_yueshu(p):
####注意这里的p就是种群，对应本约束问题，就是x，但是p是个一维的数据，若想利用矩阵计算，可将p reshape，即进行尺寸的转化
x=np.array(p)###转化为一维数组
#注意需要将原约束转化为标准形式，<=0的形式。即求和1000-求和x*a<=0,
he = 0
for j in range(5):###由于这里面涉及到5个约束，且形式相似，因此利用for循环计算，也可单独计算，而决策变量的约束系数aij是一个二维数组
tem=0
for i in range(n):###这里的n就是决策变量的维度（个数）
tem+=x[i]*a[i][j]
##这一层for循环结束后计算出了一个约束条件的值
he+=max(0,1000-tem)
#这里用1000-tem就相当于转化为了标准形式
###这一层的for循环结束后就计算出了所有的约束条件的值（不是严格意义的数值，而是大于或小于0）
#这里的max（0，tem）是关键，因为第二层for循环结束后就是计算出了一个约束条件的值，如果该值＜0证明满足约束，max（0，tem）将返回0，对结果没有影响，若＞0则返回该值，最终的约束条件值得和->he必然＞0，即只要有一个约束不成立，那么最终返回值就是＞0得，也即这个种群对应的决策变量值不满足约束
print(he)###这个只是便于检测约束和，看是否满足约束
return he
constraint= [constraint_yueshu]
###注意这里的调用计算约束函数的函数并没有写括号，传参数
四：代入模板函数，设定相关参数
n=10##设定决策变量x的个数
ga = GA(func=schaffer, n_dim=n, size_pop=50, max_iter=800, prob_mut=0.001, lb=[0]*n, ub=[1]*n, precision=[1]*n，constraint_ueq=constraint)
###注意这里的目标函数schaffer并没有写括号，传参数
best_x, best_f = ga.run()
###这里的best_x就会返回最终的最优决策变量值，所有决策变量的取值
###这里的best_f就会返回最终的最优目标函数值
print('best_x:', best_x, '\n', 'best_y:', best_f)
参数说明：（官方文档有）
五：结果可视化
Y_history = pd.DataFrame(ga.all_history_Y)
fig, ax = plt.subplots(2, 1)
ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
Y_history.min(axis=1).cummin().plot(kind='line')
plt.show()
####这里是分别画出了每代Y_history.index中的各种群目标函数值Y_history.values和最优值Y_history.min
六：完整代码：
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sko.GA import GA ###########此行为调出scikit-opt中的遗传算法
###########
#若安装
#则： pip install scikit-opt
#####################################
def schaffer(p):
'''
目标函数: min(xi)求和
决策变量维度：n
决策变量类型：0 1 整数
:param p:
:return:
'''
x=np.array(p)###转化为数组
return sum(x)###返回将目标函数的值
def constraint_yueshu(p):
####注意这里的p就是种群，对应本约束问题，就是x，
####但是p是个一维的数据，若想利用矩阵计算，可将p reshape，即进行尺寸的转化
x=np.array(p)###转化为一维数组
#注意需要将原约束转化为标准形式，<=0的形式。即求和1000-求和x*a<=0,
he = 0
for j in range(5):###由于这里面涉及到5个约束，且形式相似，因此利用for循环计算，也可单独计算，而决策变量的约束系数aij是一个二维数组
tem=0
for i in range(n):###这里的n就是决策变量的维度（个数）
tem+=x[i]*a[i][j]
##这一层for循环结束后计算出了一个约束条件的值
he+=max(0,1000-tem)
#这里用1000-tem就相当于转化为了标准形式
###这一层的for循环结束后就计算出了所有的约束条件的值（不是严格意义的数值，而是大于或小于0）
#这里的max（0，tem）是关键，因为第二层for循环结束后就是计算出了一个约束条件的值，
#如果该值＜0证明满足约束，max（0，tem）将返回0，对结果没有影响，
#若＞0则返回该值，最终的约束条件值得和->he必然＞0，
#即只要有一个约束不成立，那么最终返回值就是＞0得，也即这个种群对应的决策变量值不满足约束
print(he)###这个只是便于检测约束和，看是否满足约束
return he
constraint= [constraint_yueshu]##这一块是真正代入模板的
n=10##设定决策变量x的个数
ga = GA(func=schaffer, n_dim=n, size_pop=50, max_iter=800, prob_mut=0.001, lb=[0]*n, ub=[1]*n, precision=[1]*n，constraint_ueq=constraint)
###注意这里的目标函数schaffer并没有写括号，传参数
best_x, best_f = ga.run()
###这里的best_x就会返回最终的最优决策变量值，所有决策变量的取值
###这里的best_f就会返回最终的最优目标函数值
print('best_x:', best_x, '\n', 'best_y:', best_f)
Y_history = pd.DataFrame(ga.all_history_Y)
fig, ax = plt.subplots(2, 1)
ax[0].plot(Y_history.index, Y_history.values, '.', color='red')
Y_history.min(axis=1).cummin().plot(kind='line')
plt.show()
###注意这里并没有给出二维矩阵a的值，只是举例子
利用scipy.optimize求解
详细见代码
参考自：代码来源
'''参考自：https://blog.csdn.net/qq_39241986/article/details/104832513'''
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
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])
##此处的多约束也可仿照上面GA算法中的思想计算约束值之和，最终相当于只输出一个约束
##而等式约束可以考虑为tem+=man(0,绝对值)
##决策变量的符号约束
b=(0,None)
bnds=(b,b,b)
###求解
x=np.array([0,0,0])#定义初始值
solution=minimize(objective,x,method='SLSQP',bounds=bnds,constraints=cons)
##打印结果
x=solution.x
print('目标值: ' + str(objective(x)))
print('最优解为')
print('x1 = ' + str(round(x[0],2)))
print('x2 = ' + str(round(x[1],2)))
print('x3 = ' + str(round(x[2],2)))
print(solution)
利用geatpy
注意，其实仔细看geatpy官方文档是比较好改的。需要改的也只有自定义的目标函数了，虽然看着后面很多代码很复杂，但是其实都不用动。
可以看官方文档中的《数据结构》和《快速入门》来改函数。
看数据结构方便自己自定义目标函数的约束条件。
文档网址
import numpy
import pip
import pandas as pd
import numpy as np
import geatpy as ea
#########################
weidu=402
#########################
#############读取数据
path='E:\桌面\gyx.xlsx'
dataframe=pd.read_excel(path,'供')
data=dataframe.iloc[0:402,1:dataframe.shape[1]].values
########################################################################
class MyProblem(ea.Problem): # 继承Problem父类
def __init__(self):
name = 'MyProblem' # 初始化name（函数名称，可以随意设置）
M = 1 # 初始化M（目标维数）
maxormins = [1] # 初始化目标最小最大化标记列表，1：min；-1：max
Dim = weidu # 初始化Dim（决策变量维数）
varTypes = [1] * Dim # 初始化决策变量类型，0：连续；1：离散
#这是一个list,初始化varTypes(决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的)
lb = [0]*Dim # 决策变量下界
ub = [1]*Dim # 决策变量上界
lbin = [1]*Dim # 决策变量下边界(0表示不包含该变量的下边界，1表示包含)
ubin = [1]*Dim # 决策变量上边界(0表示不包含该变量的下边界，1表示包含)
# 调用父类构造方法完成实例化
ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb,ub, lbin, ubin)
def aimFunc(self, pop): # 目标函数，pop为传入的种群对象
Vars = pop.Phen # 得到决策变量矩阵,行表示一完整个体，列表示一个决策变量（参见geatpy数据结构文档）
############################################################################################
##############改动之处##################
##############改动之处##################
##############改动之处##################
##############改动之处##################
# x1 = Vars[:, [0]] # 取出第一列得到所有个体的x1组成的列向量
# x2 = Vars[:, [1]] # 取出第二列得到所有个体的x2组成的列向量
# x3 = Vars[:, [2]] # 取出第三列得到所有个体的x3组成的列向量
#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处
#取出第i列得到所有个体的x[i]组成的列向量
#######################################
##############改动之处##################
##############改动之处##################
##############改动之处##################
##############改动之处##################
x=[[] for i in range(weidu)]
for i in range(weidu):
x[i]=Vars[:,[i]]
##############
###########################################################################################
##############改动之处##################
##############改动之处##################
##############改动之处##################
##############改动之处##################
#计算目标函数值，赋值给pop种群对象的ObjV属性
# pop.ObjV = 4 * x1 + 2 * x2 + x3
# print('Pop.objv:',pop.Objv)
# # 采用可行性法则处理约束，生成种群个体违反约束程度矩阵
# pop.CV = np.hstack([2 * x1 + x2 - 1, # 第一个约束
# x1 + 2 * x3 - 2, # 第二个约束
# np.abs(x1 + x2 + x3 - 1)]) # 第三个约束
#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处#改动之处
#######################################
##############改动之处##################
##############改动之处##################
##############改动之处##################
##############改动之处##################
# # 计算目标函数值，赋值给pop种群对象的ObjV属性
##如何改动可参看geatpy数据结构ObjV这样就知道了
pop.ObjV=x[0]
for i in range(1,weidu):
pop.ObjV +=x[i]
# 采用可行性法则处理约束，生成种群个体违反约束程度矩阵
list = []
for i in range(24):
tem = 0
for j in range(weidu):
if data[j][0] == 'A':
tem += data[j][i + 1] * x[j] / 0.6
elif data[j][0] == 'B':
tem += data[j][i + 1] * x[j] / 0.66
else:
tem += data[j][i + 1] * x[j] / 0.72
list.append(2*1e4-tem)
pop.CV = np.hstack(list)
###########################################################################################
"""============================实例化问题对象========================"""
problem = MyProblem() # 实例化问题对象
"""==============================种群设置==========================="""
Encoding = 'RI' # 编码方式
NIND = 100 # 种群规模
Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges,
problem.borders) # 创建区域描述器
population = ea.Population(Encoding, Field, NIND) #
#实例化种群对象（此时种群还没被真正初始化，仅仅是生成一个种群对象）
"""===========================算法参数设置=========================="""
myAlgorithm = ea.soea_DE_best_1_L_templet(problem, population) #
#实例化一个算法模板对象
myAlgorithm.MAXGEN = 1500 # 最大进化代数
myAlgorithm.mutOper.F = 0.6 # 差分进化中的参数F
myAlgorithm.recOper.XOVR = 0.8 # 设置交叉概率
myAlgorithm.logTras = 100 #
#设置每隔多少代记录日志，若设置成0则表示不记录日志
myAlgorithm.verbose = True # 设置是否打印输出日志信息
myAlgorithm.drawing = 1 #
#设置绘图方式（0：不绘图；1：绘制结果图；2：绘制目标空间过程动画；3：绘制决策空间过程动画）
"""==========================调用算法模板进行种群进化==============="""
[BestIndi, population] = myAlgorithm.run() #
#执行算法模板，得到最优个体以及最后一代种群
BestIndi.save() # 把最优个体的信息保存到文件中
"""=================================输出结果======================="""
print('评价次数：%s' % myAlgorithm.evalsNum)
print('时间已过 %s 秒' % myAlgorithm.passTime)
if BestIndi.sizes != 0:
print('最优的目标函数值为：%s' % BestIndi.ObjV[0][0])
print('最优的控制变量值为：')
for i in range(BestIndi.Phen.shape[1]):
print(BestIndi.Phen[0, i])
else:
print('没找到可行解。')
注：此题的例子是2021数学建模国赛的C题066论文约束规划问题的简化
