function zf=Levenberg_Marquat_Model( )
clear;
% 2. LM算法
% 初始猜测s
%a0=10; b0=0.5;
%y_init = a0*exp(-b0*data_1);
disp('初始弹性模量');
data_1=[34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
34.5
]
disp('实测挠度');
obs_1=[1.293604
3.44961
5.174415
7.00975
10.901479
14.793208
16.195984
16.769108
17.198952
15.353562
13.508171
11.662781
8.838121
5.523826
3.314295
-3.755037
-6.258395
-10.013431
-14.229802
-19.368841
-24.507879
-29.646918
-35.193007
-42.587791
-49.67558
-55.84238
-62.00918
-63.55088
-62.00918
-55.84238
-49.67558
-42.587791
-35.193007
-29.646918
-24.507879
-19.368841
-14.229802
-10.013431
-6.258395
-3.755037
3.314295
5.523826
8.838121
11.662781
13.508171
15.353562
17.198952
16.769108
16.195984
14.793208
10.901479
7.00975
5.174415
3.44961
1.293604
]
% 数据个数
Ndata=length(obs_1);
% 参数维数
Nparams=length(data_1);
% 迭代最大次数
n_iters=10;
% LM算法的阻尼系数初值
lamda=1;
% step1: 变量赋值
updateJ=1;
data_est=data_1;
%it=1;
%e_lm=1;
%eps=1e-4;
% step2: 迭代
%while e_lm>eps
for it=1:n_iters
% it=it+1;
if updateJ==1
% 根据当前估计值，计算雅克比矩阵
J=zeros(Ndata,Nparams);
%调用ANSYS程序
system(' "J:\Program Files\Ansys Inc\v100\ANSYS\bin\intel\ansys100.exe" -b -p ane3fl -i "J:\ParaMSB\ZF.txt" -o "J:\ParaMSB\vm5.out" -j it');
delete('it.err');
delete('it.lock');
fid = fopen('B.txt', 'r');
J = fscanf(fid, '%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g ', [29 55]); % It has two rows now.
J = J' ;
fclose(fid);
% 根据当前参数，得到函数值
fid = fopen('Output.txt', 'r');
y_est = fscanf(fid, '%g ', [55 1]); % It has two rows now.
fclose(fid);
% 计算误差
d=obs_1-y_est;
% 计算（拟）海塞矩阵
H=J'*J;
% 若是第一次迭代，计算误差
if it==1
e=dot(d,d);
end
end
% 根据阻尼系数lamda混合得到H矩阵
H_lm=H+(lamda*eye(Nparams,Nparams));
% 计算步长dp，并根据步长计算新的可能的\参数估计值
dp=inv(H_lm)*(J'*d(:));
g = J'*d(:);
disp('更新后的弹性模量 ')
data_lm=data_est+0.1*dp %核心迭代代码
data_lm=0.95*data_lm;
zft(:,it:it)=data_lm;
fid=fopen('input.txt','wt');
fprintf(fid,'%16.5f \n',data_lm );
fclose(fid);
system(' "J:\Program Files\Ansys Inc\v100\ANSYS\bin\intel\ansys100.exe" -b -p ane3fl -i "J:\ParaMSB\ZFB.txt" -o "J:\ParaMSB\vm5.out" -j it');
delete('it.err');
delete('it.lock');
% 计算新的可能估计值对应的y和计算残差e
fid = fopen('Output.txt', 'r');
y_est_lm = fscanf(fid, '%g ', [55 1]); % It has two rows now.
disp('更新后的挠度')
y_est_lm
zf(:,it:it)=y_est_lm;
fclose(fid);
disp('迭代次数')
it
disp('误差值')
d_lm=obs_1-y_est_lm
disp('误差和值')
e_lm=dot(d_lm,d_lm)
% 根据误差，决定如何更新参数和阻尼系数
% if e_lm<0.01*e
% break;
% end
if e_lm<=e
lamda=lamda/3;
data_est=data_lm;
e=e_lm;
updateJ=1;
else
updateJ=0;
lamda=lamda*3;
end
end
%显示优化的结果
data_est
end
