9.5.2　遗传算法在科学研究中的应用实例

1．函数的优化

例9.1　利用遗传算法求解下列函数在区间［-60，60］的极大值：

此函数的二维图像如图9.3所示，可以看出此函数在大部分区域的值为0.5，在对角线上有多个局部极大值，全局极大值为1且位置不唯一。

图9.3　函数的图像

解：下面利用遗传算法函数ga求全局极大值。

首先编写目标函数并以文件名myfun存盘。

然后在MATLAB工作窗口输入下列命令：

改变参数，再进行运算：

若要改变其他参数，则按类似方法设置。

例9.2　利用遗传算法计算下面函数的最大值：

f（x）＝xsin（10πx）＋2.0，x∈［-1，2］

解：首先编写目标函数的M文件并以文件名myfun存盘。

利用遗传算法工具箱的GUI进行计算。

在MATLAB工作窗口输入：

打开遗传算法的GUI，在Fitness function窗口输入@myfun，在Number of variables窗口输入变量数目1，其他参数选缺省值，然后单击Start按钮运行遗传算法，得到如图9.4的结果。

图9.4　遗传算法运行结果

Fitness function value（函数值）：-3.849619541712781。对于本问题为3.8496。

Optimization terminated: maximum number of generations exceeded。

final point（变量值）：1.85139。

该函数的曲线如图9.5所示。

图9.5　f（x）＝xsin（10πx）＋2.0的图像

2．优化参数及优化问题

例9.3　体重约70kg的某人在短时间内喝下2瓶啤酒后，隔一段时间测量他的血液中酒精含量（mg/100mL），得到表9.3所列的数据。

表9.3　酒精在人体血液中分解的动力学数据

根据酒精在人体血液分解的动力学规律可知，血液中酒精浓度与时间的关系可表示为

c（t）＝k（e^-qt－e^-rt）

试根据表中数据求出参数k、q、r。

解：编写目标函数并以文件名myfun存盘。

然后在MATLAB工作窗口输入下列命令：

得到结果：x_min＝72.9706　　0.0943　　3.9407

由于遗传算法是一种随机性的搜索方法，所以每次运算可得到不同的结果。为了得到最终的结果，用直接搜索工具箱中的fminsearch函数求出最佳值：

图9.6为原始数据及用优化结果绘制的曲线。

图9.6　酒精在人体血液中分解的动力学曲线

从这个例子可看出，用遗传算法求解非线性最小二乘问题时，对最终的结果要用其他方法进行验证。

例9.4　沈阳南部浑河沿岸4个排污口污水处理效率非线性规划问题。

（杨晓华，陆桂华，等．自适应基因算法在环境优化问题的应用［J］．河海大学学报．2002，30（2）：39～41）

解：首先编写目标函数文件myfun。

然后在工作窗口输入以下命令：

可以多运行几次，以求得最好的结果。如果要使运算结果重复，可使用以下方法：

再运行，就可以得到与前一次同样的结果。

3．遗传算法在变量筛选中的应用

在科学研究中，经常会遇到非线性的多变量问题。目前处理非线性问题最流行的方法是人工神经网络法，然而该方法易发生过拟合现象，即建立的模型的误差很小，但对未知样本的预报误差则较大。其他方法如偏最小二乘和主成分分析也不能得到较理想的结果。变量扩维—筛选方法是一种简单实用的处理非线性相关关系问题的方法。

变量扩维—筛选方法是采用先扩维，即引入变量的非线性项，如变量的平方项、二次交叉项等形成新的变量，作为候选变量，再筛选，从大量的候选变量中选出最优的变量子集，用这些变量子集建立含非线性因子的拟线性模型。

变量扩维—筛选方法可分为两个步骤：

①变量扩维：将含有变量x₁，x₂，…，x_n的数据矩阵X扩维，引入变量的非线性项如，，…，x₁x₂，…，x₁/x₂和其他函数形式的项，这样将X扩维到X′。

②从矩阵X′的变量筛选出一些重要的变量，或最佳变量组合形成的矩阵X″来建立模型，使得所建立的模型有较强或最好的预报能力。

变量扩维较为简单，关键是变量筛选。变量筛选问题，特别是当变量的数目比较大时，是十分复杂的问题。解决这个问题可以采用多种方法，其中遗传算法是其中的一种。

在处理变量筛选问题时，遗传算法的编码一般采用二进制编码。对变量数为n的问题，可用一个含有n个0或1的字符串表示一个变量组合，1和0分别表示此变量选中和未选中，1在字符串的位置表示变量的序号。如00110110，表示有8个变量，其中第3、4、6和7变量被选中。

编码结束后，再利用一般的遗传算法的基本步骤，就可以求出最佳个体，即变量数及含义。

适应度函数用PRESS值。此值的含义如下：将m样本中m－1个样本用作训练样本，剩下的一个样本作检验样本。利用m－1样本建模，用检验样本代入模型，可求得一个估计值y₁。然后换另外一个样本作为检验样本，用其余样本建模，检验样本检验，得到第二个估计值y₂。如此循环m次，每次都留下一个样本做估计，最后可求得m个估计值，并可求出m个预报残差y_i－y_i－1，再将这m个残差平方求和，即为PRESS。此值越小，表示模型的预报能力越强。

为了减少计算量，在实际中可以通过普通残差来求PRESS，即

其中，e_i为普通残差；h_ii为第i个样本点到样本点中心的广义化距离，h_ii＝（X^TX）^-1。X为数据矩阵，x_i为X中的某一行矢量。

例9.5　某钢铁公司炼钢转炉的炉龄按30天炉／天炼钢规模，大约一个月就需等炉一次进行检修。为了减少消耗，厂方希望建立炉龄的预测模型，以便适当地调节参数，以延长炉龄。通过实际测定，得到表9.4的数据，其中x₁为喷补料量，x₂为吹炉时间，x₃为炼钢时间，x₄为钢水中含锰量，x₅为渣中含铁量，x₆为作业率，目标变量y为炉龄（炼钢炉次／炉）。

表9.4　转炉炉龄数据

解：由于自变量与目标变量间可能存在非线性关系，因此采用变量扩维—筛选方法处理。变量扩维引入原变量的非线性因子，即引入自变量的平方项及二次项交叉项参加建模，原自变量加上下列的非线性因子共27个因子，其编号如下：

根据以上数据就可以通过遗传算法筛选最终的变量数，即哪些变量对炉龄的影响最大。

首先编写一个适应度函数：

然后打开遗传算法的GUI就可以计算，注意此时界面中的Options的Population type要选择Bit string，其余参数按情况设定，有些可以采用默认值。

通过运算得到其中的一次结果如下：x＝［0 0 1 1 0 0 1 0 0 1 1 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0］，即变量序号为3、4、7、10、11、12、13、14、15和22，其PRESS值为10.1293。

对求出的变量的原始数据（即不进行归一化，这样在实际中应用更方便）进行多元线性回归，可得到以下的关系式：

在实际工作中，可以通过逐步回归（Stepwise）或其他方法验证上述的结果。

4．基于遗传算法的聚类分析

例9.6　在科学研究中，聚类分析是非常重要的方法。利用遗传算法也可以进行聚类分析。

人类对二维、三维图像有很强的识别能力，如果有可能将高维空间数据分布的结构特征用二维（或三维）图像显示，利用人类对二维（或三维）图像的识别能力考察高维空间数据分布结构的特征，就可能构成一种极方便的模式识别方法。

设有高维空间数据点X_i（x_i1，x_i2，…，x_im），其二维显示的对应点是Y_i（y_i1，y_i2），则y_i1，y_i2应是x_i1，x_i2，…，x_im的某种函数。如果y值是各x的某一线性组合，则二维图像是高维图像的投影，如果y值和x值是非线性函数，则二维图像是高维图像的非线性映射（Non-Linear Mapping，NLM）。现利用NLM方法分析表9.5所列的数据。

表9.5　15个标准中国茶叶样品的化学成分

解：根据NLM方法，映射时的误差函数为

其中，、d_ij分别为高维数据和二维数据的欧氏距离。据此可利用遗传算法对该函数进行最小化处理，找到合适的二维数据结构，完成高维数据到二维数据的非线性映射。

首先编写目标函数并以文件名myfun存盘。

然后利用遗传算法的GUI就可计算，其中变量数设置为30（样品数×映射维数）。一次的计算结果如下：

得到图9.7所示的结果，从图中可看出各样本的聚类情况，如1、2、8三个样本在二维空间距离比较接近，可以认为是一类，以此类推。此分类情况与普通的聚类分析的结果相同。

图9.7　高维数据映射到二维空间的结果图像

例9.7　在模式识别中，其中有一类是无管理方法，即只有一大批已知样本，事先没有规定分类标准，也没有规定要分成几类，却要通过信息处理将样本分为合适的若干类。

如对15个人的三类人群进行某4项指标的测定，结果如表9.6所列。可以认为他们可分为健康、亚健康和不健康三类，但不知道具体哪一人对应的类别，试对他们进行自动归类。

表9.6　原始数据　mg/kg

解：由于有15个样品，3个种类，所以设计的编码是长度为15，数字为1、2、3组成的随机组合。如（1 2 3 3 3 2 1 3 3 3 3 3 2 2 1）表明序号为1、7和15的为第一类，第2、6、13、14序号的为第二类。适度值函数采用类间距离之和及类内距离之和。

由于本题的编码与一般问题不同，所以采用下列自编的程序进行计算。限于篇幅，其中一些函数的源程序就不再列出。其中选择算子并不是按适应度函数的大小来进行选择的，而是将函数值转化为适应性概率，并且函数值越大，对应的概率越大。计算公式如下：

其中，P为个体的适应性概率；q为几何排序常数，一般取0.08；r个体按从大到小排列序号，函数值最大的个体序号为1，概率最大；N为种群中个体数目。另外，在常规的变异算子中要进行适当的修改，使之不出现0。

类间距离则是求两类间中心的距离，类内距离则是求样品到类中心的距离。类中心是指类的各特征向量的平均值。

运行以上程序，得到其中的一次结果：

y＝［1　1　1　1　1　1　3　3　3　3　2　2　2　2　2］，即最终样品的分类情况，其函数值为41.3261。

与其他算法求得的结果相同。

5．模拟退火算法的应用

例9.8　试用模拟退火算法求解不重复历经我国31个省会城市，并且路径最短的线路，即邮递员难题（TSP问题）。其中城市坐标为

解：根据模拟退火算法的原理，编写如下函数。

运行便可得到结果，其中一次的结果，如图9.8所示。

图9.8　城市坐标及最短路径

需要说明的是，标准的TSP问题解应该是封闭的曲线，在此例中由于没有对一些条件进行约束，所以求出的并不是一条封闭曲线。

例9.9　用模拟退火算法求解下列函数的极小值：

解：利用较高版本MATLAB中优化工具箱的模拟退火函数求解。首先编写下列函数：

然后在命令窗口输入下列命令：

结果为：在x＝（-0.0904　0.7126）处有最小值fval＝-1.0316。

例9.10　利用模拟退火算法对表9.7所列的某年我国20个地区的三次产业产值数据进行聚类分析。

表9.7　某年我国20个地区的三次产业产值

解：对题中的数据作图，可以看出分成二类或者三类。与例9.8类似，编制相应的程序。限于篇幅，不再列出全部的程序，其中目标函数为

其中，X为样本向量；ω为聚类划分；X^（ω_i）为第i个聚类的中心，即类内所有样本特征值的平均值；为样本到对应聚类中心距离。

首先将前三个样本分为三类，再计算其余样本与各类中心的距离，按距离远近划分三类（也可以随机对所有样本进行分类）。然后用模拟退火算法对初始分类的结果进行修正，得到最终的分类结果。其中一次将数据分成三类的结果为

第1类：1　3　6　8　9　11　12　13　16　17

第2类：10　15　18

第3类：2　4　5　7　14　19　20