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

阻尼最小二乘法与模拟退火法结合实现非线性模型参数的估计收稿日期： 接受日期：
*基金项目：国家自然科学基金（30960031）和江西省自然科学基金（2009GZN0076）。

刘建军， 陈明锋， 叶子飘
（ 井冈山大学数理学院， 江西吉安343009）
摘 要：将模拟退火法与阻尼最小二乘法相结合，得到了对非线性模型参数估计的算法。该算法最大特点是克服了LM算法无法跳出局部优解和初始值要求高的难题，同时克服了模拟退火法搜索效率逐步降低的问题，同时给出了阻尼最小二乘法编程技巧。将该算法运用于植物光合作用光响应新模型对水稻的光响应曲线的拟合，得到的拟合参数与DPS拟合值极为接近，残差平方和比DPS软件更小、相关系数更大。并给出了该算法的MATLAB 7.0程序。
关键词：阻尼最小二乘法；模拟退火法；算法结合；曲线拟合；MATLAB

1 前言
参数估计（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 combining
3 拟合植物光合作用光响应新模型曲线
植物光合作用光响应新模型的数学表达式为[10-11]：
（14）
式中：P(I)为净光合速率，I为光强，α为光响应曲线的初始斜率，β为修正系数，γ为一个与光强无关的系数，Rd为暗呼吸速率。
由（14）式拟合光响应数据求系数的MATLAB 7.0程序为如下：
function beta= yzpfit(X,y)
format long;%start Simulated annealing
t=6000;%initial temperature
beta=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 t>2
beta1=[rand(1)/10 rand(1)/100 rand(1)/100 rand(1)*2];
yfit1 = model(beta,X);
rr=(y-yfit1)*(y-yfit1)';
if rr<r1
r1=rr;
beta=beta1;
elseif rand<exp(-(rr-r1)/t/1000)
r1=rr;
beta=beta1;
end
t=t-1;
end%end Simulated annealing and beginning Marquardt method
yfit = model(beta,X);
n = numel(X);%Number of variables
p = 4;%Parameter number
sqrteps = sqrt(eps);%Initial bias
J = zeros(n,p);%The establishment of the initial matrix
r = y(:) - yfit(:);
sse = r'*r;%Deviation sum of squares Set up convergence tolerances from options.
maxiter =80;%Cycles
betatol =eps;%Accuracy
rtol = eps;
fdiffstep=eps^(0.5);%Step
zbeta = zeros(1,4);%The initial changes of parameter
zerosp = zeros(p,1);% Set initial weight for LM algorithm.
lambda = .01;%Initial amount of Marquardt
iter = 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 = 0.1*lambda; % f the LM increased SSE, increased lambda t
else
while 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 long
z=model(beta,X);
%plot(X,y,'+',X,z,'-o');
xgxs=corrcoef(y,z);
xgxc=xgxs(1,2);%Correlation
sse;
beta=[beta,xgxc,sse];
我们以表1给出了温度为、CO2浓度为300 μmol mol-1条件下水稻的光响应数据为例[12]，说明结合算法与DPS（Data Processing System）数据处理软件[13]拟合这组数据的优缺点。
表1：温度为30 °C，CO2浓度为300 μmol mol-1条件下水稻的光响应数据实测值[12]
Table 1 Measured data for rice at 30 °C and 300 μmol mol-1.[12]
光强I
2001
1799
1599
1399
1199
1000
800
599
398
200
净光合速率
P (I)
22.5
22.9
23.1
22.6
21.7
20.7
19.2
17
13.8
8.15
光强I
180
161
141
121
80
60
40
20
1

净光合速率
P (I)
7.26
6.44
5.88
5.05
2.83
1.59
0.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
残差平方和R
Matlab程序
0.061 918
0.000 122
0.001 435
1.417 678
0.999 5
1.340 562
DPS
0.063 925
0.000 108
0.001 436
1.416 547
0.998 5
2.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 concentration
Tl为叶面温度，Ca为气室CO2浓度，Pmax为最大净光合速率，Isat为饱和光强，Ic为光补偿点，Rd为暗呼吸速率， ○ 为测量点， —+— 为拟合点
从图2可知，利用此方法可以得到很好的拟合结果，且与实测值符合程度很高。因此，利用模拟退火法与阻尼最小二乘法相结合完成可以实现非线性方程的拟合。在拟合曲线时不需要初始值。可见该方法有其优点。
参考文献:
[1]成 平, 陈希孺, 陈桂景, 吴传义. 参数估计[M]. 上海: 上海科学技术出版社, 1985
[2]Blaha G. Non-iterative approach to nonlinear least-squares adjustment [J]. Manuscripta Gepdaetica.1994, 19: 199-212
[3]Marquardt 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 2

[4]陈宝林. 最优化理论与算法[M]. 北京: 清华大学出版社, 1989
[5]李士勇. 模糊控制·神经网络控制和智能控制[M]. 哈尔滨: 哈尔滨工业大学出版社，1998
[6]康立山. 非数值并行算法（第一册）模拟退火算法[M]. 北京: 科学出版社, 1994:50-61.
[7]王新洲. 非线性模型参数估计理论与应用[M]. 武汉: 武汉大学出版社, 2002:63.
[8]胡茂林. 矩阵计算与应用[M]. 北京：科学出版社, 2008: 152-156.
[9]帅 虹. 空间网格结构风振抑制的改进模拟退火算法研究[D]. 上海交通大学，2008
[10]Ye 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 method

LIU 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 initial value, but also does the question of simulated annealing algorithm, which its searching efficiency gradually reduces. Also the algorithm of MATLAB programs was given in this paper. At last, the algorithm was used in the new model of light-response curve of photosynthesis to fit photosynthetic response to irradiance for rice. The obtained parameters and fitted values are very close to those obtained by DPS, whereas the residual sum of squares is smaller than that of DPS, and the correlation coefficient is greater than that of DPS. At last, the combination algorithm program of MATLAB 7.0 was given in this paper.
Key words: damped least-squares algorithm; simulated annealing algorithm; algorithm combination; curve fitting; MATLAB
9