290-机器学习算法全家桶(python代码),包含智能优化算法,单/双重分解算法,残差修正,降噪,统计模型,机器学习算法

#!/usr/bin/env python# coding: utf-8# In[1]:import time  # 导入 time 模块，用于时间相关操作#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制import copy  # 导入 copy 模块，用于复制对象import numpy  # 导入 numpy 模块，用于数值计算import random  # 导入 random 模块，用于生成随机数import datetime  # 导入 datetime 模块，用于日期和时间相关操作import math  # 导入 math 模块，用于数学计算import pandas as pd  # 导入 pandas 模块，用于数据处理import numpy as np  # 导入 numpy 模块，用于数值计算import matplotlib.pyplot as plt  # 导入 matplotlib.pyplot 模块，用于绘图from PyEMD import EMD, EEMD, CEEMDAN  # 从 PyEMD 模块中导入 EMD、EEMD、CEEMDAN 类，用于经验模态分解from sampen import sampen2  # 从 sampen 模块中导入 sampen2 函数，用于计算样本熵from vmdpy import VMD  # 从 vmdpy 模块中导入 VMD 类，用于变分模态分解import tensorflow as tf  # 导入 tensorflow 模块，用于深度学习from sklearn.cluster import KMeans  # 从 sklearn.cluster 模块中导入 KMeans 类，用于聚类分析from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, mean_absolute_percentage_error  # 从 sklearn.metrics 模块中导入评估指标函数from sklearn.preprocessing import MinMaxScaler  # 从 sklearn.preprocessing 模块中导入 MinMaxScaler 类，用于数据归一化from tensorflow.keras.models import Sequential  # 从 tensorflow.keras.models 模块中导入 Sequential 类，用于构建模型from tensorflow.keras.layers import Dense, Activation, Dropout, LSTM, GRU  # 从 tensorflow.keras.layers 模块中导入 Dense、Activation、Dropout、LSTM、GRU 类，用于构建神经网络层from keras.layers import Conv1D, Add, Dropout, Flatten, MaxPooling1D  # 从 keras.layers 模块中导入 Conv1D、Add、Dropout、Flatten、MaxPooling1D 类，用于构建卷积神经网络层from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping  # 从 tensorflow.keras.callbacks 模块中导入 ReduceLROnPlateau、EarlyStopping 类，用于回调函数import warnings  # 导入 warnings 模块，用于忽略警告信息from scipy.fftpack import hilbert, fft, ifft  # 从 scipy.fftpack 模块中导入 hilbert、fft、ifft 函数，用于傅里叶变换from math import log  # 导入 math 模块中的 log 函数，用于数学计算中的对数运算from typing import List  # 导入 typing 模块中的 List 类型，用于函数参数注解#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制warnings.filterwarnings("ignore")  # 忽略警告信息# In[2]:# 调用GPU加速gpus = tf.config.experimental.list_physical_devices(device_type='GPU')for gpu in gpus:    tf.config.experimental.set_memory_growth(gpu, True)# In[3]:plt.rcParams['font.sans-serif'] = ['Times New Roman']  # 设置图形的字体为'Times New Roman'，用于正常显示中文plt.rcParams['axes.unicode_minus'] = False# 用来解决图形中负号显示异常的问题，设置为False表示不使用unicode编码来显示负号，以确保负号显示正常# In[4]:df_raw_data = pd.read_csv('焦作.csv', usecols=[0,1])  # 从名为'焦作.csv'的CSV文件中读取数据，仅选择第一列和第二列进行读取，并存储在DataFrame对象df_raw_data中series_close = pd.Series(df_raw_data['AQI'].values, index=df_raw_data['time'])  # 将df_raw_data中列名为'AQI'的数据转换为Series对象，其中'AQI'列的值将作为Series的值，而列名为'time'的数据将作为Series的索引test = df_raw_data['AQI'].values[int(len(df_raw_data['AQI'].values)*0.8):]  # 从'AQI'列的值中根据切片操作获取80%之后的部分，并将其赋值给test变量# In[5]:timestep = 30# 定义时间步数，用于创建训练集和测试集tau = 0.# 设置VMD分解的参数tau，用于控制模态函数的带宽DC = 0# 设置VMD分解的参数DC，用于控制是否提取直流分量init = 1# 设置VMD分解的参数init，用于设置初始模态函数数量tol = 1e-7# 设置VMD分解的参数tol，用于控制停止迭代的阈值# In[6]:classProblemModel:def__init__(self, bounds=None):        self.bounds = bounds  # 初始化类的实例时，将参数bounds赋值给实例变量self.boundsdefgetIndependentVar(self):if self.bounds isnotNone:  # 如果self.bounds不为None            independentVar = []for bound in self.bounds:# 在每个bound范围内生成一个随机数，并添加到independentVar列表中                independentVar.append(bound[0] + random.random() * (bound[1] - bound[0]))return independentVar  # 返回生成的独立变量列表else:pass# 如果self.bounds为None，则不执行任何操作#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制defgetNewVar(self, var_1, var_2):if self.bounds isnotNone:  # 如果self.bounds不为None            newVar = []#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制            step_random = random.random()  # 生成一个0到1之间的随机数for v_1, v_2 in zip(var_1, var_2):# 根据var_1和var_2的值生成一个新的变量，并添加到newVar列表中                newVar.append(v_1 + step_random * (v_2 - v_1))return newVar  # 返回生成的新变量列表else:pass# 如果self.bounds为None，则不执行任何操作defgetValue(self, variable):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制if len(variable) == 2:  # 如果variable列表长度为2            K = int(variable[0])  # 将variable的第一个元素转换为整数赋值给变量K            alpha = variable[1]  # 将variable的第二个元素赋值给变量alphaif K < self.bounds[0][0]:  # 如果K小于bounds的第一个元素的第一个值                K = self.bounds[0][0]  # 将K赋值为bounds的第一个元素的第一个值if K > self.bounds[0][1]:  # 如果K大于bounds的第一个元素的第二个值                K = self.bounds[0][1]  # 将K赋值为bounds的第一个元素的第二个值if alpha < self.bounds[1][0]:  # 如果alpha小于bounds的第二个元素的第一个值                alpha = self.bounds[1][0]  # 将alpha赋值为bounds的第二个元素的第一个值if alpha > self.bounds[1][1]:  # 如果alpha大于bounds的第二个元素的第二个值                alpha = self.bounds[1][1]  # 将alpha赋值为bounds的第二个元素的第二个值            u, u_hat, omega = VMD(df_raw_data['AQI'], alpha, tau, K, DC, init, tol)# 调用VMD函数，传入相应参数，并将返回的结果分别赋值给变量u、u_hat和omega            EP = []for i in range(K):  # 对于范围为K的循环迭代                H = np.abs(hilbert(u[i, :]))  # 对u[i, :]进行希尔伯特变换并取绝对值，赋值给变量H                e1 = []for j in range(len(H)):  # 对于范围为H长度的循环迭代                    p = H[j] / np.sum(H)  # 计算H[j]占H所有元素之和的比例，赋值给变量p                    e = -p * log(p, 2)  # 计算信息熵并赋值给变量e                    e1.append(e)  # 将e添加到e1列表中                E = np.sum(e1)  # 计算e1列表中所有元素的和，并赋值给变量E                EP.append(E)  # 将E添加到EP列表中            s = np.sum(EP) / K  # 计算EP列表中所有元素的和除以K，并赋值给变量sreturn s  # 返回变量s作为结果else:return1# 如果variable列表长度不为2，则返回1作为结果# In[7]:classNectarSource:    problem_src = None# 类变量problem_src，用于存储问题源def__init__(self, position):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制        self.position = position  # 实例变量position，存储蜜源的位置        self.value = self.problem_src.getValue(position)  # 实例变量value，根据问题源的getValue方法计算蜜源的值if self.value >= 0:  # 如果蜜源的值大于等于0            self.fitness = 1 / (1 + self.value)  # 计算蜜源的适应度else:            self.fitness = 1 + math.fabs(self.value)  # 计算蜜源的适应度（值小于0时）        self.trail = 0# 初始化实例变量trail为0，用于存储蜜源的信息素浓度# In[8]:classABCAlgor:    LIMIT = 10# 类变量LIMIT，用于指定蜜源的最大尝试次数def__init__(self, problem, employedNum, onlookerNum, maxIteration):        NectarSource.problem_src = problem  # 将问题源赋值给类变量problem_src        self.problem = problem  # 实例变量problem，存储问题源        self.employedNum = employedNum  # 实例变量employedNum，指定雇佣蜜蜂数量        self.onlookerNum = onlookerNum  # 实例变量onlookerNum，指定观察蜜蜂数量        self.maxIteration = maxIteration  # 实例变量maxIteration，指定最大迭代次数        self.nectarSrc = []  # 实例变量nectarSrc，存储蜜源列表        self.bestNectar = NectarSource(self.problem.getIndependentVar())  # 实例变量bestNectar，存储最优蜜源        self.resultRecord = []  # 实例变量resultRecord，存储迭代过程中的结果记录for i in range(self.employedNum):            self.nectarSrc.append(NectarSource(self.problem.getIndependentVar()))  # 初始化蜜源列表#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制defupdateNectarSrc(self, index):        src = self.nectarSrc[index]  # 获取指定索引的蜜源        src_another = random.choice(self.nectarSrc)  # 从蜜源列表中随机选择另一个蜜源while src_another is src:            src_another = random.choice(self.nectarSrc)  # 如果随机选择的蜜源与当前蜜源相同，则重新选择        src_new = NectarSource(self.problem.getNewVar(src.position, src_another.position))  # 基于当前蜜源和另一个蜜源生成新的蜜源if src_new.fitness > src.fitness:  # 如果新蜜源的适应度大于当前蜜源的适应度            self.nectarSrc[index] = src_new  # 更新蜜源列表中的蜜源为新蜜源else:            self.nectarSrc[index].trail += 1# 否则，当前蜜源的信息素浓度加1defemployedProcedure(self):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制        length = len(self.nectarSrc)  # 获取蜜源列表的长度for i in range(length):            self.updateNectarSrc(i)  # 更新蜜源列表中的每个蜜源defonlookerProcedure(self):        sum_fitness = 0# 初始化总适应度为0for src in self.nectarSrc:            sum_fitness += src.fitness  # 计算总适应度        length = len(self.nectarSrc)  # 获取蜜源列表的长度for i in range(length):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制            probability_fit = self.nectarSrc[i].fitness / sum_fitness  # 计算每个蜜源的适应度概率for onlookerBee in range(self.onlookerNum):if random.random() < probability_fit:  # 根据适应度概率决定是否更新蜜源                    self.updateNectarSrc(i)  # 更新蜜源列表中的蜜源defupdateBestNectar(self):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制for src in self.nectarSrc:if src.fitness > self.bestNectar.fitness:  # 如果蜜源的适应度大于最优蜜源的适应度                self.bestNectar = src  # 更新最优蜜源为当前蜜源defscoutProcedure(self):        length = len(self.nectarSrc)  # 获取蜜源列表的长度for i in range(length):if self.nectarSrc[i].trail >= self.LIMIT:  # 如果蜜源的信息素浓度超过限制                self.nectarSrc[i] = NectarSource(self.problem.getIndependentVar())  # 重新生成蜜源defrun(self):for i in range(self.maxIteration):  # 在最大迭代次数范围内进行迭代            print('迭代次数：', i)  # 输出当前迭代次数            self.employedProcedure()  # 执行雇佣蜜蜂阶段            self.onlookerProcedure()  # 执行观察蜜蜂阶段            self.updateBestNectar()  # 更新最优蜜源#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制            self.scoutProcedure()  # 执行侦查蜜蜂阶段            self.updateBestNectar()  # 再次更新最优蜜源            self.resultRecord.append(self.bestNectar.value)  # 将当前最优蜜源的值记录下来return self.bestNectar.position, self.bestNectar.value  # 返回最优蜜源的位置和值        print('best solution:', int(self.bestNectar.position[0]), self.bestNectar.position[1])  # 输出最优蜜源的位置        print('best value:', self.bestNectar.value)  # 输出最优蜜源的值defshowResult(self):for result in self.resultRecord:            print(result)  # 输出每次迭代的结果记录        plt.plot(self.resultRecord)  # 绘制结果曲线        plt.title('result curve')  # 设置标题        plt.tight_layout()  # 调整布局        plt.show()  # 显示结果曲线#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制# In[9]:beesNum = 1# 蜜蜂数量为1employedNum = int(beesNum/2)  # 雇佣蜜蜂数量为蜜蜂数量的一半（向下取整）onlookerNum = int(beesNum/2)  # 观察蜜蜂数量为蜜蜂数量的一半（向下取整）maxIteration = 1# 最大迭代次数为1problem = ProblemModel(bounds=([3, 25], [100, 3500]))  # 创建问题模型对象，设置变量范围abcSolution = ABCAlgor(problem, employedNum, onlookerNum, maxIteration)  # 创建ABC算法对象，传入问题模型和参数abcSolution.run()  # 运行ABC算法abcSolution.showResult()  # 显示结果# In[10]:defceemdan_decompose(series=None, trials=10, num_clusters=3):    decom = CEEMDAN()  # 创建CEEMDAN对象    decom.trials = trials  # 设置分解的试验次数    df_ceemdan = pd.DataFrame(decom(series.values).T)  # 对数据进行CEEMDAN分解并转换为数据框    df_ceemdan.columns = ['imf'+str(i) for i in range(len(df_ceemdan.columns))]  # 为每一列设置列名为'imf' + 对应的索引号return df_ceemdan  # 返回分解后的数据框# In[11]:defsample_entropy(df_ceemdan=None, mm=1, r=0.1):    np_sampen = []  # 存储样本熵的列表for i in range(len(df_ceemdan.columns)):        sample_entropy = sampen2(list(df_ceemdan['imf'+str(i)].values), mm=mm, r=r, normalize=True)  # 计算样本熵        np_sampen.append(sample_entropy[1][1])  # 将样本熵的值添加到列表中    df_sampen = pd.DataFrame(np_sampen, index=['imf'+str(i) for i in range(len(df_ceemdan.columns))])  # 创建样本熵的数据框，设置行索引为'imf' + 对应的索引号return df_sampen  # 返回样本熵的数据框# In[12]:defkmeans_cluster(df_sampen=None, num_clusters=3):    np_integrate_form = KMeans(n_clusters=num_clusters, random_state=9).fit_predict(df_sampen)  # 使用K均值聚类进行聚类操作    df_integrate_form = pd.DataFrame(np_integrate_form, index=['imf'+str(i) for i in range(len(df_sampen.index))], columns=['AQI'])  # 创建聚类结果的数据框，设置行索引为'imf' + 对应的索引号，列名为'AQI'return df_integrate_form  # 返回聚类结果的数据框# In[13]:defintegrate_imfs(df_integrate_form=None, df_ceemdan=None):    df_tmp = pd.DataFrame()  # 创建一个空的数据框用于存储临时结果for i in range(df_integrate_form.values.max()+1):        df_tmp['imf'+str(i)] = df_ceemdan[df_integrate_form[(df_integrate_form['AQI']==i)].index].sum(axis=1)  # 对每个聚类簇内的IMF分量进行求和，得到综合的IMF分量    df_integrate_result = df_tmp.T  # 对临时结果进行转置    df_integrate_result['sampen'] = sample_entropy(df_tmp).values  # 计算综合的IMF分量的样本熵，并将其作为新的列添加到结果数据框中    df_integrate_result.sort_values(by=['sampen'], ascending=False, inplace=True)  # 根据样本熵降序排列综合的IMF分量    df_integrate_result.index = ['co-imf'+str(i) for i in range(df_integrate_form.values.max()+1)]  # 为综合的IMF分量设置新的行索引，命名规则为'co-imf' + 对应的索引号    df_integrate_result = df_integrate_result.drop('sampen', axis=1, inplace=False)  # 移除样本熵这一列return df_integrate_result.T  # 返回结果数据框的转置# In[14]:defvmd_decompose(series=None, alpha=math.floor(abcSolution.bestNectar.position[1]), tau=tau,                   K=math.floor(abcSolution.bestNectar.position[0]), DC=DC, init=init, tol=tol, draw=True): #更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制    imfs_vmd, imfs_hat, omega = VMD(series, alpha, tau, K, DC, init, tol)  # 使用VMD方法对时间序列进行分解    df_vmd = pd.DataFrame(imfs_vmd.T)  # 将VMD分解得到的IMFs转换为数据框    df_vmd.columns = ['imf'+str(i) for i in range(K)]  # 为每个IMF分量设置列名，命名规则为'imf' + 对应的索引号return df_vmd  # 返回VMD分解得到的IMFs数据框# In[15]:defcreate_train_test_set(data=None, timestep=timestep, co_imf_predict_for_fitting=None):if isinstance(data, pd.DataFrame):  # 如果数据是数据框类型        dataY = data['sum'].values.reshape(-1, 1)  # 将数据框中名为'sum'的列作为目标变量，并将其转换为二维数组形式        dataX = data.drop('sum', axis=1, inplace=False)  # 去除数据框中的目标变量列，得到输入特征矩阵else:  # 如果数据是数组类型        dataY = data.values.reshape(-1, 1)  # 将数据转换为二维数组形式作为目标变量        dataX = dataY  # 输入特征矩阵与目标变量相同    scalarX = MinMaxScaler(feature_range=(0,1))  # 创建MinMaxScaler对象，将输入特征矩阵归一化到0~1之间    dataX = scalarX.fit_transform(dataX)  # 对输入特征矩阵进行归一化if co_imf_predict_for_fitting isnotNone:  # 如果提供了用于拟合的共同IMF预测结果        co_imf_predict_for_fitting = scalarX.transform(co_imf_predict_for_fitting)  # 对共同IMF预测结果进行归一化    scalarY = MinMaxScaler(feature_range=(0,1))  # 创建MinMaxScaler对象，将目标变量归一化到0~1之间    dataY = scalarY.fit_transform(dataY)  # 对目标变量进行归一化    trainX, trainY = [], []  # 初始化训练集的输入特征和目标变量列表for i in range(len(dataY)-timestep):        trainX.append(np.array(dataX[i:(i+timestep)]))  # 将前timestep个时间步作为输入特征添加到训练集的输入特征列表        trainY.append(np.array(dataY[i+timestep]))  # 将第timestep个时间步的目标变量添加到训练集的目标变量列表if co_imf_predict_for_fitting isnotNone:  # 如果提供了用于拟合的共同IMF预测结果if i < (len(dataY)-timestep-len(co_imf_predict_for_fitting)):                trainX[i] = np.insert(trainX[i], timestep, dataX[i+timestep], 0)  # 将真实的输入特征添加到训练集的输入特征中else:                trainX[i] = np.insert(trainX[i], timestep, co_imf_predict_for_fitting[i-(len(dataY)-timestep-len(co_imf_predict_for_fitting))], 0)  # 将共同IMF预测结果添加到训练集的输入特征中return np.array(trainX), np.array(trainY), scalarY  # 返回训练集的输入特征、目标变量和目标变量的归一化器对象# In[16]:''' 种群初始化函数 '''definitial(pop, dim, ub, lb):    X = np.zeros([pop, dim])  # 创建一个形状为 (pop, dim) 的全零数组 X，用于存储种群的初始位置for i in range(pop):for j in range(dim):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制            X[i, j] = random.random() * (ub[j] - lb[j]) + lb[j]  # 生成一个位于 lb[j] 和 ub[j] 之间的随机数，赋值给 X[i, j]return X, lb, ub'''边界检查函数'''defBorderCheck(X, ub, lb, pop, dim):for i in range(pop):for j in range(dim):#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制if X[i, j] > ub[j]:  # 如果 X[i, j] 大于上界 ub[j]，将其设置为上界                X[i, j] = ub[j]elif X[i, j] < lb[j]:  # 如果 X[i, j] 小于下界 lb[j]，将其设置为下界                X[i, j] = lb[j]return X'''计算适应度函数'''defCaculateFitness(X, fun):    pop = X.shape[0]  # 种群大小    fitness = np.zeros([pop, 1])  # 创建一个形状为 (pop, 1) 的全零数组 fitness，用于存储适应度值for i in range(pop):        fitness[i] = fun(X[i, :])  # 计算第 i 个个体的适应度值，将结果赋值给 fitness[i]return fitness'''适应度排序'''defSortFitness(Fit):    fitness = np.sort(Fit, axis=0)  # 按列排序适应度值数组 Fit，将结果赋值给 fitness    index = np.argsort(Fit, axis=0)  # 返回按列排序后的索引数组，将结果赋值给 indexreturn fitness, index'''根据适应度对位置进行排序'''defSortPosition(X, index):    Xnew = np.zeros(X.shape)  # 创建一个与位置数组 X 相同形状的全零数组 Xnewfor i in range(X.shape[0]):        Xnew[i, :] = X[index[i], :]  # 根据索引数组 index 对位置数组 X 进行排序，并将结果赋值给 Xnewreturn Xnew# In[17]:defDBO(pop, dim, lb, ub, MaxIter, fun):# 参数设置    PballRolling = 0.2# 滚球蜣螂比例    PbroodBall = 0.4#产卵蜣螂比例    PSmall = 0.2# 小蜣螂比例    Pthief = 0.2# 偷窃蜣螂比例    BallRollingNum = int(pop*PballRolling) #滚球蜣螂数量    BroodBallNum = int(pop*PbroodBall) #产卵蜣螂数量    SmallNum = int(pop*PSmall) #小蜣螂数量    ThiefNum = int(pop*Pthief) #偷窃蜣螂数量    X, lb, ub = initial(pop, dim, ub, lb)  # 初始化种群    fitness = CaculateFitness(X, fun)  # 计算适应度值# 记录全局最优    minIndex = np.argmin(fitness)  # 找到适应度最小值的索引    GbestScore = copy.copy(fitness[minIndex])  # 复制最小适应度值作为全局最优分数    GbestPositon = np.zeros([1, dim])  # 创建一个全0矩阵来存储全局最优位置    GbestPositon[0, :] = copy.copy(X[minIndex, :])  # 复制最小适应度对应的位置为全局最优位置    Curve = np.zeros([MaxIter, 1])  # 创建一个全0数组用于记录迭代过程中的最优分数    Xl = copy.deepcopy(X)  # 用于记录上一代的种群位置# 记录当前代种群    cX = copy.deepcopy(X)  # 复制当前种群位置    cFit = copy.deepcopy(fitness)  # 复制当前种群适应度值for t in range(MaxIter):   # 迭代次数循环，从0到MaxIter-1        print("第" + str(t) + "次迭代")# 蜣螂滚动 文献中式（1），（2）更新位置# 获取种群最差值        maxIndex = np.argmax(fitness)  # 找到适应度最大值的索引        Wort = copy.copy(X[maxIndex, :])  # 复制最大适应度对应的位置为Wort        r2 = np.random.random()  # 生成一个随机数r2for i in range(0,BallRollingNum): # 循环迭代变量 i 在从 0 到 BallRollingNum（不包括）的范围内if r2<0.9: # 如果 r2 的值小于 0.9if np.random.random()>0.5:  # 如果 np.random.random() 生成的随机数大于 0.5                    alpha=1# 设置alpha为1else:  # 如果 np.random.random() 生成的随机数不大于 0.5                    alpha=-1# 设置alpha为-1                b = 0.3# 设置b为0.3                k = 0.1# 设置k为0.1#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制                X[i,:]=cX[i,:]+b*np.abs(cX[i,:]-Wort)+alpha*k*Xl[i,:]# 对 X 数组中的第 i 行进行赋值计算，计算结果由 cX[i, :], b, np.abs(cX[i, :] - Wort), alpha 和 k * Xl[i, :] 组合得到else:  # 如果 r2 的值不小于 0.9                theta = np.random.randint(180)# 生成一个0到179之间的随机整数if theta==0or theta == 90or theta == 180: # 如果theta的值为0、90或180度                    X[i,:]=copy.copy(cX[i,:])# 将当前位置复制给新位置else:#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制                    theta = theta*np.pi/180# 将theta转换为弧度制                    X[i,:]=cX[i,:]+np.tan(theta)*np.abs(cX[i,:]-Xl[i,:])# 对 X 数组中的第 i 行进行赋值计算，计算结果由 cX[i, :], np.tan(theta), np.abs(cX[i, :] - Xl[i, :]) 组合得到for j in range(dim): # 循环迭代变量 j 在从 0 到 dim 的范围内if X[i,j]>ub[j]:  # 如果 X 数组中的第 i 行、第 j 列的值大于上界 ub[j]                    X[i,j]=ub[j]  # 将 X 数组中的第 i 行、第 j 列的值设置为上界 ub[j]if X[i,j]<lb[j]:  # 如果 X 数组中的第 i 行、第 j 列的值小于下界 lb[j]                    X[i,j]=lb[j]  # 将 X 数组中的第 i 行、第 j 列的值设置为下界 lb[j]            fitness[i]=fun(X[i,:]) # 计算第 i 行的适应度值，计算结果由 fun 函数根据 X[i, :] 得到if fitness[i]<GbestScore:  # 如果第 i 行的适应度值小于全局最优适应度值 GbestScore                GbestScore=copy.copy(fitness[i]) # 将全局最优适应度值 GbestScore 更新为第 i 行的适应度值的副本                GbestPositon[0,:]=copy.copy(X[i,:])  # 将全局最优位置 GbestPositon 的第一行设置为 X 数组中第 i 行的副本# 当前迭代最优#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制        minIndex=np.argmin(fitness)  # 使用 np.argmin 函数找到适应度值数组 fitness 中的最小值的索引        GbestB = copy.copy(X[minIndex,:])  # 将 GbestB 设置为 X 数组中索引为 minIndex 的行的副本# 蜣螂产卵 ，文献中式（3）        R=1-t/MaxIter # 根据当前迭代次数 t 和最大迭代次数 MaxIter 计算 R 的值        X1=GbestB*(1-R)  # 根据全局最优解 GbestB 和 R 计算 X1 的值        X2=GbestB*(1+R)  # 根据全局最优解 GbestB 和 R 计算 X2 的值        Lb = np.zeros(dim)  # 创建长度为 dim 的零数组 Lb 和 Ub        Ub = np.zeros(dim)#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制for j in range(dim):  # 循环迭代变量 j 在从 0 到 dim 的范围内            Lb[j]=max(X1[j],lb[j])  # Lb[j] 的值为 X1[j] 和 lb[j] 中的较大值            Ub[j]=min(X2[j],ub[j])  # Ub[j] 的值为 X2[j] 和 ub[j] 中的较小值for i in range(BallRollingNum,BallRollingNum+BroodBallNum):# 循环迭代变量 i 在从 BallRollingNum 到 BallRollingNum + BroodBallNum 的范围内            b1=np.random.random()  # 生成一个随机数 b1            b2=np.random.random()  # 生成一个随机数 b2            X[i,:]=GbestB+b1*(cX[i,:]-Lb)+b2*(cX[i,:]-Ub)# 根据公式更新 X 数组中的第 i 行的值for j in range(dim): # 循环迭代变量 j 在从 0 到 dim 的范围内if X[i,j]>ub[j]:  # 如果 X 数组中的第 i 行、第 j 列的值大于上界 ub[j]                    X[i,j]=ub[j]  # 将 X 数组中的第 i 行、第 j 列的值设置为上界 ub[j]if X[i,j]<lb[j]:  # 如果 X 数组中的第 i 行、第 j 列的值小于下界 lb[j]                    X[i,j]=lb[j]  # 将 X 数组中的第 i 行、第 j 列的值设置为下界 lb[j]            fitness[i]=fun(X[i,:])  # 计算第 i 行的适应度值，计算结果由 fun 函数根据 X[i, :] 得到if fitness[i]<GbestScore:  # 如果第 i 行的适应度值小于全局最优适应度值 GbestScore                GbestScore=copy.copy(fitness[i]) # 将全局最优适应度值 GbestScore 更新为第 i 行的适应度值的副本                GbestPositon[0,:]=copy.copy(X[i,:]) # 将全局最优位置 GbestPositon 的第一行设置为 X 数组中第 i 行的副本# 小蜣螂更新#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制#文献中(5),(6)        R=1-t/MaxIter # 根据当前迭代次数 t 和最大迭代次数 MaxIter 计算 R 的值        X1=GbestPositon[0,:]*(1-R) # 根据全局最优位置 GbestPositon 和 R 计算 X1 的值        X2=GbestPositon[0,:]*(1+R) # 根据全局最优位置 GbestPositon 和 R 计算 X2 的值        Lb = np.zeros(dim) # 创建长度为 dim 的零数组 Lb 和 Ub        Ub = np.zeros(dim)#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制for j in range(dim):  # 循环迭代变量 j 在从 0 到 dim 的范围内            Lb[j]=max(X1[j],lb[j])  # Lb[j] 的值为 X1[j] 和 lb[j] 中的较大值            Ub[j]=min(X2[j],ub[j])  # Ub[j] 的值为 X2[j] 和 ub[j] 中的较小值for i in range(BallRollingNum+BroodBallNum,BallRollingNum+BroodBallNum+SmallNum):  # 循环迭代变量 i 在从 BallRollingNum + BroodBallNum 到 BallRollingNum + BroodBallNum + SmallNum 的范围内            C1 = np.random.random([1,dim]) # 生成一个随机数组 C1，形状为 [1, dim]            C2 = np.random.random([1,dim]) # 生成一个随机数组 C2，形状为 [1, dim]            X[i,:]=GbestPositon[0,:]+C1*(cX[i,:]-Lb)+C2*(cX[i,:]-Ub) # 根据公式更新 X 数组中的第 i 行的值for j in range(dim): # 循环迭代变量 j 在从 0 到 dim 的范围内if X[i,j]>ub[j]: # 如果 X 数组中的第 i 行、第 j 列的值大于上界 ub[j]                    X[i,j]=ub[j] # 将 X 数组中的第 i 行、第 j 列的值设置为上界 ub[j]if X[i,j]<lb[j]: # 如果 X 数组中的第 i 行、第 j 列的值小于下界 lb[j]                    X[i,j]=lb[j] # 将 X 数组中的第 i 行、第 j 列的值设置为下界 lb[j]#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制            fitness[i]=fun(X[i,:]) # 计算第 i 行的适应度值，计算结果由 fun 函数根据 X[i, :] 得到if fitness[i]<GbestScore: # 如果第 i 行的适应度值小于全局最优适应度值 GbestScore                GbestScore=copy.copy(fitness[i]) # 将全局最优适应度值 GbestScore 更新为第 i 行的适应度值的副本                GbestPositon[0,:]=copy.copy(X[i,:]) # 将全局最优位置 GbestPositon 的第一行设置为 X 数组中第 i 行的副本# 当前迭代最优        minIndex=np.argmin(fitness) # 使用 np.argmin 函数找到适应度值数组 fitness 中的最小值的索引        GbestB = copy.copy(X[minIndex,:]) # 将 GbestB 设置为 X 数组中索引为 minIndex 的行的副本# 偷窃蜣螂更新 # 文献中式（7）for i in range(pop-ThiefNum,pop): # 循环迭代变量 i 在从 pop - ThiefNum 到 pop 的范围内            g=np.random.randn() # 生成一个符合标准正态分布的随机数 g            S=0.5# 设置 S 的值为 0.5            X[i,:]=GbestPositon[0,:]+g*S*(np.abs(cX[i,:]-GbestB)+np.abs(cX[i,:]-GbestPositon[0,:])) # 根据公式更新 X 数组中的第 i 行的值for j in range(dim): # 循环迭代变量 j 在从 0 到 dim 的范围内if X[i,j]>ub[j]: # 如果 X 数组中的第 i 行、第 j 列的值大于上界 ub[j]                    X[i,j]=ub[j] # 将 X 数组中的第 i 行、第 j 列的值设置为上界 ub[j]if X[i,j]<lb[j]: # 如果 X 数组中的第 i 行、第 j 列的值小于下界 lb[j]                    X[i,j]=lb[j] # 将 X 数组中的第 i 行、第 j 列的值设置为下界 lb[j]            fitness[i]=fun(X[i,:]) # 计算第 i 行的适应度值，计算结果由 fun 函数根据 X[i, :] 得到if fitness[i]<GbestScore: # 如果第 i 行的适应度值小于全局最优适应度值 GbestScore                GbestScore=copy.copy(fitness[i]) # 将全局最优适应度值 GbestScore 更新为第 i 行的适应度值的副本                GbestPositon[0,:]=copy.copy(X[i,:])  # 将全局最优位置 GbestPositon 的第一行设置为 X 数组中第 i 行的副本# 记录t代种群        Xl= copy.deepcopy(cX) # 将 Xl 设置为 cX 的深拷贝，即创建一个与 cX 一样的副本#更新当前代种群#更多模型搜索“清朝简单的饮料”或者“安仁坊天蝎座果蔬”，支持模型定制for i in range(pop):  # 循环迭代变量 i 在从 0 到 pop 的范围内if fitness[i]<cFit[i]:  # 如果第 i 个个体的适应度值小于当前代个体的适应度值 cFit[i]                cFit[i]=copy.copy(fitness[i]) # 将当前代个体的适应度值 cFit[i] 更新为第 i 个个体的适应度值的副本                cX[i,:]=copy.copy(X[i,:])  # 将当前代个体 cX 的第 i 行设置为 X 数组中第 i 行的副本        Curve[t] = GbestScore # 将 Curve 数组的第 t 个元素设置为全局最优适应度值 GbestScorereturn GbestScore, GbestPositon, Curve # 返回全局最优适应度值 GbestScore、全局最优位置 GbestPositon 和适应度曲线 Curve