混料D-最优近似设计的改进模拟退火算法

许秀钿,张崇岐

(广州大学 经济与统计学院,广东广州 510006)

§1 引言

混料试验设计自1958年Scheffé[1]提出以来发展十分迅速,当前已被广泛应用于工业,食品,医药,金融等领域.例如,工业领域,孙世昌等[2]在2021年利用混料试验设计优化烟用胶囊壁材配方,提高了工业生产的烟用胶囊质量;食品领域,张丽霞等[3]在2021年通过最优混料试验设计优化风味蛋白酶,碱性蛋白酶和复合蛋白酶的组合,为花生饼粕制备花生呈味基料提供一种高效复合酶;医药领域,张雪等[4]在2022年通过混料D-最优设计优选出适宜的挤出搓圆法制备丸剂的复合辅料;金融领域,燕飞和张崇岐[5]用混料模型G-最优设计来研究证券投资组合.

通过近年来混料的应用研究[6-7]可以发现,如何科学的设计试验方案,以较少的工作量获得最优的试验结果是当前混料设计的研究热点.由于传统的计算方法对于求解高阶混料模型或是约束域较为复杂的混料模型难度较大,利用优化算法求解混料最优设计[8]逐渐成为该领域关注的对象.

模拟退火算法[9](Simulated Annealing)又称SA算法,是在1953年由Metropolis提出的一种理论上的全局优化算法.在1983年Kirkpatrick[10]成功将其应用到求解组合优化问题上.当前应用较为广泛,比较具有代表性的智能优化算法[11]是遗传算法[12],禁忌搜索算法,模拟退火算法[13-14].与其他种类的优化算法相比,模拟退火算法理论较为完善.本文介绍一种改进的SA算法,该算法保留了SA算法易于求出全局最优解的优点,同时在求解精度方面相比一般的SA算法有一定程度的提升.

§2 预备知识

最基本的混料模型为

其中y是响应变量,β是未知参数组成的d维列向量,f(x)是在x ∈χ上给定的d维函数列向量,χ是一个封闭紧集,ε为随机误差,设其服从均值为0,方差为常数的正态分布.

假设设计ξ由设计点zi和对应的测度ri组成,

其中r1+r2+···+rN=1,Ξ是所有设计构成的集合.

设计ξ的信息矩阵[15]定义为

称信息矩阵的行列式最大化即max detM(ξ)的设计ξ为D-最优设计.

混料试验是研究各个分量比例x1,x2,···,xq的设计问题,q分量混料系统各种成份的占比xi(i=1,2,···,q)必须满足约束条件

即在无附加约束的混料问题上,试验的可行域是一个q-1维正规单纯形

在实际的问题中,有时出于经济与技术方面的考虑,混料试验往往还需满足其他约束条件,常见的附加约束类型[16]有

(1) 同时具有上,下界约束条件(SCCs):C′s:0≤ai ≤xi ≤bi ≤1,i=1,2,···,q;

(2) 具有线性约束条件(MCCs):C′s:0≤ai ≤≤bi ≤1,i=1,2,···,q;

(3) 具有非线性约束条件(NCCs):C′s:ai ≤φi(x1,x2,···,xq)≤bi ≤1,i=1,2,···,q;其中ai,bi都是已知的常数,φi(x1,x2,···,xq)是各变量的非线性函数.

1958年,Scheffé首次提出单纯形格子设计,奠定了混料试验设计的基础.目前许多设计方法是由该设计衍生出来的,其中混料规范多项式应用较为广泛.在寻找最优设计的过程中,常使用的q分量m阶混料规范多项式有

(1)m=1时,

(2)m=2时,

(3)m=3时,

显然当q或m增加时,模型的参数个数也会随之增加.

如果试验设计的试验点数与相应的混料规范多项式模型中未知参数的个数相等,则设计为饱和设计.易知饱和设计是试验点数最少情况下的最优设计.

§3 改进SA算法

在实际的混料最优设计的问题求解中,由于模型参数个数较多或是约束域较为复杂,模型的D-最优设计的解析解往往较难求出.此时借助计算机编程算法求解是一种很好的解决方法.

由于D-最优设计具有等测度的性质,故求解D-最优近似最优设计的过程也是求解测度为参数个数的倒数的最优支撑点的过程.

模拟退火算法是借鉴固体退火过程原理设计而成的算法.固体温度很高时,内能很大,固体内部粒子的无规则运动很快,此时整个试验域的点都有一定概率得以取到.随着温度参数慢慢降低,内能逐渐减小,固体内部粒子逐渐趋于有序,即算法中的解趋于稳定.但是,这个稳定的解可能只是局部最优解.所以在解逐渐收敛的过程中,模拟退火算法会以一定的概率跳出这个局部最优解.即固体内能增加,取到的解远离目标函数的时候,算法以一定概率的接纳,而不是一味拒绝.这也是模拟退火算法的核心思想.

3.1 改进SA算法的实现过程

基于改进的SA算法求解混料模型最优设计的基本步骤如下.

1.设定初始值.设置迭代的初始设计点,为确保迭代能遍历整个可行域,初始设计点一般随机生成并均匀分布在整个可行域上.

2.产生随机扰动项.由一个产生函数在初始设计点附近产生一个位于试验可行域上的新设计点.本文提出的改进的SA算法相比传统的SA算法,考虑了随机扰动项波动范围的处理,将步长设置为温度的平方根.温度越低,得到的解约趋于稳定时,步长越短,波动范围越小,能更好的提高搜索设计点的效率.

3.判断新解是否能被接受.最常用的判断依据是Metropolis准则,即当新设计点的目标函数小于当前设计点的目标函数时,接受新设计点作为新的当前设计点;否则则以一定的概率接受新设计点为新的当前设计点.

4.生成最优设计点集.当“温度”降到“冷却温度”时,迭代停止,最终保留下来的目标函数最小的设计点即为所得最优支撑点.本文提出的改进SA算法将马尔科夫链长度设置为温度的函数,使得温度越低时,马尔科夫链越长.相较于一般的SA算法,极大的提高了搜索和迭代的效率.

将最优支撑点与其对应的测度组合即可得到该问题的D-最优近似设计.此外,对于求解出的混料模型,也可应用该算法求出最优解.

3.2 改进SA算法求解D-最优近似设计的工作流程

第1步 初始化

设置初始温度T0,最终冷却温度Tf,温度衰减速率α ∈(0,1),马尔科夫链长度L,步长Step.

k=L/t为每一温度的迭代次数,显然温度降低时,迭代的次数将会增加.

第2步 产生随机扰动项,随机生成一个领域N(i)

第3步 计算目标值增量∆f=f(j)-f(i),其中j ∈N(i)

第4步 判断新解是否被接受,若∆f<0或∆f>0且exp(-)>U(0,1),则令i=j

第5步 若达到热平衡,即某一温度的迭代结束时,则继续降温,直至降为最终冷却温度,否则返回第2步

第6步 生成最优设计点集

将最终冷却温度时保留的最优支撑点集与取值为参数个数倒数的测度组合即可得到所求D-最优近似设计.

将最优支撑点集进行试验,得到相应的响应值.代入混料模型中求解出各参数估计值,即可得到经验回归方程.通过改进的SA算法,将上述步骤中目标函数改为经验回归方程,初始化时只随机抽取1个设计点,便可类似求解出试验的最优设计点.

3.3 改进SA算法求解的伪代码

§4 混料试验的应用:纱线伸长度[18]

纺纱线前需要将三种原料:聚乙烯(x1),聚苯乙烯(x2)和聚丙烯(x3)混合制成纤维.固定在作用力下纱线的伸长度是试验者感兴趣的响应.

此次试验采用三分量二阶混料规范多项式模型,其表达式为

其可行域为S3-1,分别利用改进SA算法和传统的SA算法求出其D-最优近似设计,结果如表1和图1所示,图1中改进的SA算法得出的设计点标记为三角形,传统的SA算法得出的设计点标记为正方形,理论上的设计点标记为圆形.

图1 两种算法得出的解与理论解的三元图

表1 改进的SA算法与传统的SA算法得出的D-最优近似设计解析解比较

易知此时参数个数为6,当给定初始温度T0=100,最终冷却温度Tf=1,温度衰减速率α=0.99,马尔科夫链长度L=1000000时,2种算法的三分量二阶混料规范多项式模型D-最优近似饱和设计的结果如表1所示.对于三分量二阶混料规范多项式,D-最优设计的理论解[19-20]是单纯形上3个顶点与3个边界中点的组合,即

因此从图1中的结果可知,对于D-近似最优设计而言,用改进的SA算法比传统的SA算法求出的D-最优近似设计数值解跟理论解更为接近.

从图2中可以看到,随着迭代次数的增加,温度逐渐降低.改进的模拟退火算法在目标函数值降为-4.7162×10-9之后达到平衡,此时代码运行总时间为3233.792秒;传统的模拟退火算法在目标函数值降为-3.9007×10-11之后达到平衡，此时代码运行总时间为13465.094秒.可以发现改进SA算法的目标函数值更低,即由改进的SA算法得到的设计较传统SA算法更优.同时,根据这两种算法的代码运行时间可知,改进的SA算法相比传统的SA 算法求解速度更快.

图2 改进(左图)和传统(右图)的SA算法的D-最优近似设计目标函数值与迭代情况

本次试验的部分试验数据和设计点的响应平均值分别如表2和图3所示.

图3 设计点的响应平均值

表2 纱线伸长度部分试验数据

通过最小二乘法即可算出本例所拟合的经验回归方程为

该模型的响应等值线图如图4所示.

图4 纱线伸长度的等值线图

从算法运行结果可以看出,通过改进的SA算法得到最优设计点为(0.3014,0.0082,0.6904),其对应的响应为17.2905,如图5所示.

图5 响应目标函数值与迭代情况

§5 结束语

在给定混料模型和可行域的情况下,通过改进的SA算法来求解混料D-最优近似设计不失为一个好方法.本文给出了一般情况下的混料规范多项式求解D-最优近似设计解析解的改进SA算法,通过比较可知,该算法相比传统的SA算法,求解的速度更快,精度更高.此外,本文提出的改进的SA算法也可用于求解其他最优设计的设计点,由于篇幅限制,本文仅以D-最优设计为例进行展开.