阻尼最小二乘与模拟退火结合.doc

阻尼最小二乘法与模拟退火法结合实现非线性模型参数的估计收稿日期： 接受日期：*基金项目：国家自然科学基金（30960031）和江西省自然科学基金（2009GZN0076）。刘建军， 陈明锋， 叶子飘（ 井冈山大学数理学院， 江西吉安343009）摘 要：将模拟退火法与阻尼最小二乘法相结合，得到了对非线性模型参数估计的算法。该算法最大特点是克服了LM算法无法跳出局部优解和初始值要求高的难题，同时克服了模拟退火法搜索效率逐步降低的问题，同时给出了阻尼最小二乘法编程技巧。将该算法运用于植物光合作用光响应新模型对水稻的光响应曲线的拟合，得到的拟合参数与DPS拟合值极为接近，残差平方和比DPS软件更小、相关系数更大。并给出了该算法的MATLAB 7.0程序。关键词：阻尼最小二乘法；模拟退火法；算法结合；曲线拟合；MATLAB1 前言 参数估计（Parameter Estimation）是数理统计的一个重要分支，更是测量数据处理（Surveying Data Processing）理论的重要组成部分1-2。现实中的实际模型是非线性模型（Nonlinear Model）。因此，研究非线性模型参数估计就有重大的意义。求非线性强度很强的非线性模型参数时，由于线性近似将产生大于观测的模型误差，所以一般采用迭代法求解，较常用的迭代方法有牛顿类法（包括牛顿法、信赖域法、拟牛顿法）、最速下降法、高斯牛顿法、改进的高斯牛顿法和阻尼最小二乘法3，不同的迭代算法有其各自的优点与缺点。当非线性模型平差系统缺乏基准时，阻尼最小二乘法仍能顺利迭代，而阻尼最小二乘法亦不需太多地改变传统平差程序，所以，对于非线性模型的秩亏自由网平差，采用阻尼最小二乘法最好。因为阻尼最小二乘法是局部收敛4，当初值较差时，会出现迭代发散现象。如何有效地确定参数初始值始终是难以克服的瓶颈，因此，一些实际问题可能永远无法获得满意解。非经典优化算法模拟退火法5（Simulated Annealing Algorithm，SAA）通过解的暂时的恶化，跳出局部最优的“陷阱”，可得到全局最优解。但其搜索效率会随着向最优解靠拢而逐渐的降低6，而且要得到较满意的解，则运算时间就会增加。本文在阻尼最小二乘法的基础上，发展了阻尼最小二乘法与模拟退火法结合的耦合算法用模拟退火法产生阻尼最小二乘法的初值。这将是较好的非线性模型参数估计方法。将此算法应用到植物光合作用光响应新模型对水稻光响应曲线的拟合，并将数值实验结果与DPS软件比较。2 算法分析2.1最小二乘法对于非线性方程的拟合，非线性模型的相应误差方程为： (1) 参数估计值; 测量值。于是残差平方和为： (2)所有解构成解空间，在这空间中选取的一个估计值，若满足： (3)则是的最小二乘解。 2.2.1 阻尼最小二乘法（Levenberg-Marquardt算法） (4)对在泰勒展开，进行线性化的误差方程为： (5) (6) (7)在求解时采用加入阻尼的迭代方法7： (8)式中的为大于0的任意常数。2.2.2 Levenberg-Marquardt算法的编程技巧在编程计算时用差分替代微分，简化表达（5）式得到下式： （9）设定步长为 （10）那么经过n次迭代就可以得到。为第k次迭代阻尼因子，未经迭代时即为初始阻尼因子，为任意常数。将列平方和作为对角元构成，与合并，得到。dy与n维零向量合并，得到。所以有： （11）式中的是的伪逆矩阵。 （12） 将代人（2）式得到并与比较。如果大就增大，小就减小，增大减小倍数由计算的精度决定，一般定为10。重复以上过程，就可计算出。Levenberg-Marquardt算法对于非线性模型秩亏自由网平差较适合，在参数的初始值接近时，可以很快得到的值，如果的初始值与差距较大时，就无法得到全局最优解。2.3 模拟退火算法（Simulated Annealing Algorithm）模拟退火算法是模拟统计物理学里面的固体退火过程，利用计算机随机产生可能的解模拟粒子的排列状态9，模拟淬火过程中的能量。根据Metropolis准则对应的转移概率： （13）确定是否接受从当前解i到新解j的转移。在进行足够多的转移后，缓慢减小T的值(与“徐徐”降温相对应)，如此重复，直至满足某个停止准则时算法终止。但随着向最优解靠近其搜索效率会逐渐下降。2.4算法结结合由上面的分析可得，Levenberg-Marquardt算法与Simulated Annealing 算法有其各自的优点与不足，它们可以优势互补。把它们结合在一起得到更优的拟合方法。首先采用模拟退火法计算Levenberg-Marquardt算法的初始值，当温度T下降时参数在解空间中任意变化。从而跳出局部的较优解，向全局最优解靠拢。然后用Levenberg-Marquardt算法法迭代得到全局最优解得确切值。图1 算法结合的流程图Fig.1. Flowchart of algorithm combining3 拟合植物光合作用光响应新模型曲线植物光合作用光响应新模型的数学表达式为10-11： （14）式中：P(I)为净光合速率，I为光强，为光响应曲线的初始斜率，为修正系数，为一个与光强无关的系数，Rd为暗呼吸速率。由（14）式拟合光响应数据求系数的MATLAB 7.0程序为如下：function beta= yzpfit(X,y)format long;%start Simulated annealing t=6000;%initial temperaturebeta=rand(1,4);%generated initial value Randomly model=inline(beta(1)*(1-beta(2)*X).*X./(1+beta(3)*X)-beta(4),beta,X);%new model of light-response curve of photosynthesis yfit = model(X,beta);r1=yfit*yfit;while t2 beta1=rand(1)/10 rand(1)/100 rand(1)/100 rand(1)*2; yfit1 = model(beta,X);rr=(y-yfit1)*(y-yfit1);if rrr1 r1=rr; beta=beta1; elseif randexp(-(rr-r1)/t/1000) r1=rr; beta=beta1;endt=t-1;end%end Simulated annealing and beginning Marquardt methodyfit = model(beta,X);n = numel(X);%Number of variablesp = 4;%Parameter numbersqrteps = sqrt(eps);%Initial biasJ = zeros(n,p);%The establishment of the initial matrixr = y(:) - yfit(:);sse = r*r;%Deviation sum of squares Set up convergence tolerances from options.maxiter =80;%Cyclesbetatol =eps;%Accuracyrtol = eps;fdiffstep=eps(0.5);%Stepzbeta = zeros(1,4);%The initial changes of parameter zerosp = zeros(p,1);% Set initial weight for LM algorithm.lambda = .01;%Initial amount of Marquardtiter = 0;while iter maxiter iter = iter + 1; betaold = beta; sseold = sse; % Compute a finite difference approximation to the Jacobian for k = 1:p delta = zbeta; delta(k) = fdiffstep*beta(k); yplus =model(beta+delta,X); dy = yplus(:) - yfit(:); J(:,k) = dy/delta(k);%approximation to the Jacobian end % Levenberg-Marquardt step Jplus = J; diag(sqrt(lambda*sum(J.2); %r1=J*r; rplus = r; zerosp; step = Jplus rplus; beta(:) = beta(:) + step; % Evaluate the fitted values at the new coefficients and compute the residuals and the SSE. yfit = model(beta,X); r = y(:) - yfit(:);sse = r*r; % If the LM reduced SSE, reducing lambda if sse sseold lambda = 10*lambda; Jplus = J; diag(sqrt(lambda*sum(J.2); step = Jplus rplus; beta(:) = betaold(:) + step; yfit = model(beta,X); r = y(:) - yfit(:); sse = r*r; end end % Check The Cycles and the convergence of SSE.end format longz=model(beta,X);%plot(X,y,+,X,z,-o);xgxs=corrcoef(y,z);xgxc=xgxs(1,2);%Correlationsse;beta=beta,xgxc,sse; 我们以表1给出了温度为、CO2浓度为300 mol mol-1条件下水稻的光响应数据为例12，说明结合算法与DPS（Data Processing System）数据处理软件13拟合这组数据的优缺点。表1：温度为30 C，CO2浓度为300 mol mol-1条件下水稻的光响应数据实测值12Table 1 Measured data for rice at 30 C and 300 mol mol-1.12光强I200117991599139911991000800599398200净光合速率P (I)22.522.923.122.621.720.719.21713.88.15光强I180161141121806040201净光合速率P (I)7.266.445.885.052.831.590.334-0.179-0.78表2是利用上面的算法所给MATLAB程序和DPS数据软件处理表1的数据得到的拟合结果。表2：利用上面的算法所给MATLAB程序处理表1的数据得到的拟合结果与DPS对比Table 2 Comparison of results fitted by Matlab mentioned above and DPS初始斜率修正系数系数暗呼吸速率Rd相关系数r残差平方和RMatlab程序0.061 9180.000 1220.001 435 1.417 6780.999 51.340 562DPS0.063 9250.000 1080.001 4361.416 5470.998 52.534 761图2. 水稻在温度为30 C，CO2浓度为300 mol mol-1条件下的光响应Fig.2 Light-response curve of photosynthesis for rice at 30 C and 300 mol mol-1 CO2 concentrationTl为叶面温度，Ca为气室CO2浓度，Pmax为最大净光合速率，Isat为饱和光强，Ic为光补偿点，Rd为暗呼吸速率， 为测量点， + 为拟合点 从图2可知，利用此方法可以得到很好的拟合结果，且与实测值符合程度很高。因此，利用模拟退火法与阻尼最小二乘法相结合完成可以实现非线性方程的拟合。在拟合曲线时不需要初始值。可见该方法有其优点。参考文献:1成 平, 陈希孺, 陈桂景, 吴传义. 参数估计M. 上海: 上海科学技术出版社, 19852Blaha G. Non-iterative approach to nonlinear least-squares adjustment J. Manuscripta Gepdaetica.1994, 19: 199-2123Marquardt DW. An algorithm for least squares estimation of nonlinear parameters J. SIAM Journal, 1963, 11: 431-441【2】Vol.19, No.4, 1994我不知道该怎么写【3】Volume 11, Issue 24陈宝林. 最优化理论与算法M. 北京: 清华大学出版社, 19895李士勇. 模糊控制神经网络控制和智能控制M. 哈尔滨: 哈尔滨工业大学出版社，19986康立山. 非数值并行算法（第一册）模拟退火算法M. 北京: 科学出版社, 1994:50-61.7王新洲. 非线性模型参数估计理论与应用M. 武汉: 武汉大学出版社, 2002:63.8胡茂林. 矩阵计算与应用M. 北京：科学出版社, 2008: 152-156.9帅 虹. 空间网格结构风振抑制的改进模拟退火算法研究D. 上海交通大学，200810Ye ZP. A new model for relationship between light intensity and the rate of photosynthesis in Oryza sativa J. Photosynthetica, 2007, 45(4): 637-640.11叶子飘, 赵泽海. 遮光对三叶鬼针草光合作用和叶绿素含量的影响J. 生态学杂志，2009, 28(1):19-22.12叶子飘, 王健林. 植物光合-光响应模型的比较分析J. 井冈山学院学报, 2009, 30(4): 9-13.13唐启义, 冯明光. DPS数据处理系统2实验设计、统计分析及数据挖掘M. 北京: 科学技术出版社, 2006 Estimating the non-linear model parameter using combination of damped least-squares method and simulated annealing methodLIU Jian-Jun, CHEN Ming-Feng， YE Zi-Piao(Maths & Physics College, Jinggangshan University, Ji an, Jiangxi 343009, P. R. China)Abstract: Combination the simulated annealing method with the damped minimum least-squares method, a algorithm about the function curve fitting was obtained. This algorithm not only overcomes shortcoming of the Levenberg-Marquardt algorithm, which can not jump out of local optimal solution and highly qualitative setting initi